Modelling northern peatland area and carbon dynamics since the Holocene with the ORCHIDEE-PEAT land surface model (SVN r5488)

Abstract. The importance of northern peatlands in the global carbon cycle has been
recognized, especially for long-term changes. Yet, the complex interactions
between climate and peatland hydrology, carbon storage, and area dynamics
make it challenging to represent these systems in land surface models. This
study describes how peatlands are included as an independent sub-grid
hydrological soil unit (HSU) in the ORCHIDEE-MICT land surface model. The
peatland soil column in this tile is characterized by multilayered vertical
water and carbon transport and peat-specific hydrological properties. The
cost-efficient version of TOPMODEL and the scheme of peatland initiation and
development from the DYPTOP model are implemented and adjusted to simulate
spatial and temporal dynamics of peatland. The model is tested across a
range of northern peatland sites and for gridded simulations over the
Northern Hemisphere (>30∘ N). Simulated northern
peatland area (3.9 million km2), peat carbon stock (463 Pg C), and peat
depth are generally consistent with observed estimates of peatland area (3.4–4.0 million km2), peat carbon (270–540 Pg C), and data compilations
of peat core depths. Our results show that both net primary production (NPP)
and heterotrophic respiration (HR) of northern peatlands increased over the
past century in response to CO2 and climate change. NPP increased more
rapidly than HR, and thus net ecosystem production (NEP) exhibited a
positive trend, contributing a cumulative carbon storage of 11.13 Pg C since
1901, most of it being realized after the 1950s.



Introduction
Northern peatland carbon (C) stock is estimated between 270 and 540 Pg C across an area of 3.4-4 million km 2 (Gorham, 1991;Turunen et al., 2002;Yu et al., 2010), amounting to approximately one-fourth of the global soil C pool (2000-2700 Pg C) and one-half of the current atmospheric C pool (828 Pg C) (Ciais et al., 2013;Jackson et al., 2017). Due to waterlogged, acidic, and low-temperature conditions, plant litter production exceeds decomposition in northern peatlands. More than half of northern peat carbon was accumulated before 7000 years ago during the Holocene (Yu, 2012). While being one of the most effective ecosystems at sequestering CO 2 from the atmosphere over the long term, northern peatlands are one of the largest natural sources of methane (CH 4 ), playing a pivotal role in the global greenhouse gas Published by Copernicus Publications on behalf of the European Geosciences Union.
The carbon balance of peatlands is sensitive to climate variability and climate change (Chu et al., 2015;Lund et al., 2012;Yu et al., 2003a). Projected climate warming and precipitation changes press us to understand the mechanisms of peat growth and stability and further to assess the fate of the substantial amount of carbon stored in peatlands and its potential feedbacks on the climate. Several land surface models (LSMs) have included representations of the biogeochemical and physical processes of peatlands to simulate the observed past extent and carbon balance of peatlands and predict their responses to future climate change (Chaudhary et al., 2017a, b;Frolking et al., 2010;Kleinen et al., 2012;Spahni et al., 2013;Stocker et al., 2014;Wania et al., 2009a, b;Wu et al., 2016). Water table is one of the most important factors controlling the accumulation of peat because it limits oxygen supply to the saturated zone and reduces decomposition rates of buried organic matter (Kleinen et al., 2012;Spahni et al., 2013). It is highlighted by observed and experimental findings that variations in ecosystem respiration (ER) depend on water table depth (Aurela et al., 2007;Flanagan and Syed, 2011). However, some studies showed that changes in soil water content could be very small while the water table was lowering; the drawdown of the water table caused only small changes in soil air-filled porosity and hence exerted no significant effect on ER (Lafleur et al., 2005;Parmentier et al., 2009;Sulman et al., 2009). Therefore, while studying the interactions between peatland water and carbon balances, the dynamics of soil moisture deserve special attention.
The two-layered (acrotelm-catotelm) conceptual framework was chosen by many Earth system models (ESMs) to describe peatland structures. The peat profile was divided into an upper layer with a fluctuating water table (acrotelm) and a lower, permanently saturated layer (catotelm) -using depth in relation to a drought water table or a constant value (a widely used depth is 0.3 m below the soil surface) as the discrete boundary of these two layers (Kleinen et al., 2012;Spahni et al., 2013;Wania et al., 2009a). This diplotelmic model assumes that all threshold changes in peatland soil ecological, hydrological, and biogeochemical processes occur at the same depth, causing the lack of generality and flexibility in the model, and thus possibly hindering the representation of the horizontal and vertical heterogeneity of peatlands (Fan et al., 2014;Morris et al., 2011).
To our knowledge, only two models attempted to simulate peatland area dynamics for large-scale gridded applications (Kleinen et al., 2012;Stocker et al., 2014). Kleinen et al. (2012) modelled wetland extent and peat accumulation in boreal and arctic peatlands over the past 8000 years using the LPJ model. In their study, simulated summer mean, maximum and minimum wetland extent by TOPMODEL are used as surrogates for peatland area, from the assumption that peatland will only initiate and grow in frequently inundated areas. Stocker et al. (2014) extended the scope of Kleinen et al. (2012) in the DYPTOP model. In their model, soil water storage and retention were enhanced and runoff was reduced by accounting for peatland-specific hydraulic properties. A positive feedback on the local water balance and on peatland expansion was therefore exerted by the peatland water table and peatland area fraction within a grid cell. Areas that are suitable for peatland development were distinguished from wetland extent according to temporal persistency of inundation, water balance, and peatland C balance. While both studies made pioneering progress in the modelling of peatland ecosystems, they adopted a simple bucket approach to model peatland hydrology and peatland C accumulation, and neither of them resolved the diel cycle of surface energy budget.
To tackle these above-mentioned discrepancies and estimate the C dynamic as well as the peat area, we used the ORCHIDEE-MICT land surface model incorporating peatland as a sub-grid hydrological soil unit (HSU). The vertical water fluxes and dynamic carbon profiles in peatlands are simulated with a multilayer scheme instead of a bucket model or a diplotelmic model (Sect. 2.1). Peatland extent is modelled following the approach of DYPTOP (Stocker et al., 2014) but with some adaptions and improvements (Sect. 2.2). The aim of this study is to model the spatial extent of northern peatlands since the Holocene and to reproduce peat carbon accumulation over the Holocene.

Model description
ORCHIDEE-MICT is an updated version of the ORCHIDEE land surface model with an improved and evaluated representation of high-latitude processes. Soil water freezing and melting, and subsequent changes in thermal and hydrological properties, as well as latent heat release and consumption involved in the freeze-thaw processes are all simulated by this model (Guimberteau et al., 2018). The model simulates a more rapid thermal signal propagation and a reduction in soil water infiltration and movement in frozen soil . The model calculates the active layer thickness (ALT) from simulated soil temperatures and adjusts root distribution and soil carbon inputs relative to the ALT to represent impacts of permafrost physics on plant water availability and soil carbon profiles. It is worth mentioning that the model resolves one energy budget for all soil tiles in one grid cell, but soil thermal properties of a specific grid cell are defined as the weighted average of mineral soil and pure organic soil in that grid, with C content of the grid cell derived from the soil organic C map from NCSCD (Hugelius et al., 2013) for permafrost regions and from HWSD (FAO, 2012) for nonpermafrost regions (Guimberteau et al., 2018). This makes it possible to include the impacts of peat carbon on the grid cell soil thermodynamics.
Based on ORCHIDEE-MICT, ORCHIDEE-PEAT is specifically developed to dynamically simulate northern Geosci. Model Dev., 12, 2961-2982, 2019 www.geosci-model-dev.net/12/2961/2019/ peatland extent and peat accumulation. ORCHIDEE-PEAT version 1 was evaluated and calibrated against eddycovariance measurements of CO 2 and energy fluxes, water table depth, and soil temperature from 30 northern peatland sites . Parameterizations of peatland vegetation and water dynamics are unchanged from ORCHIDEE-PEAT version 1: vegetation growing in peatlands is represented by one C 3 grass plant functional type (PFT) with shallow roots (see dedicated Sect. 2.2.1 of Qiu et al., 2018, for additional discussion on peatland PFT); surface runoff of nonpeatland areas in the grid cell is routed into peatland; vertical water fluxes in peatland HSU are modelled with peatspecific hydraulics (Text S1 in the Supplement). The large porosity (0.9 m 3 m −3 ) and the large saturated water conductivity (2120 mm d −1 ) of the peatland HSU, as well as the addition of an above-surface water reservoir, reduce runoff and increase soil water storage and retention . Therefore, the occurrence and expansion of peatland increase the grid cell mean water table and enhance inundation.
In ORCHIDEE-PEAT, the hydrology of peatland is resolved by an 11-layer physically based diffusion scheme . Compared to the two-layer bucket approach, this multilayer diffusion scheme allowed a more realistic representation of surface water fluxes and showed better performance at simulating soil water storage and soil water storage variations (Guimberteau et al., 2014;De Rosnay et al., 2002). Here, we improve peatland C dynamics by replacing the diplotelmic peatland C model in ORCHIDEE-PEAT version 1 with a multilayered one. The 32-layered thermal and C models in the standard ORCHIDEE-MICT are used to simulate peatland C accumulation and decomposition (Sect. 2.1). With fine resolution in the soil surface (10 layers for the top 1 m), this 32-layer model better represents the effects of soil temperature, soil freezing, and soil moisture on carbon decomposition continuously within the peat profile than a diplotelmic model. Furthermore, the approach proposed by Stocker et al. (2014) is incorporated into the model to simulate dynamics of peatland area (Sect. 2.2). This model simulating the dynamics of peatland extent and the vertical buildup of peat is hereinafter referred to as ORCHIDEE-PEAT v2.0.

Modelling peat accumulation and decomposition
The model has two litter C pools (metabolic and structural) and three soil C pools (active, slow, and passive); all pools are vertically discretized into 32 layers, with exponentially coarser vertical resolution as depth increases and a total depth of 38 m. Decomposition of the C in each pool and the C fluxes between the pools are calculated at each layer, with each pool having a distinct residence time. A detailed description of the litter and soil C pools and carbon flows between them can be found in the Supplement Text S2.

Peat carbon decomposition
Decomposition of peat soil C is calculated at each layer, controlled by base decomposition rates of different pools modified by soil temperature, moisture, and depth: where k i,l is the decomposition rate of pool i at layer l, k 0,i is the base decomposition rate of pool i, f T ,l is the temperature modifier at layer l, f M,l is the moisture modifier, and f Z,l is a depth modifier that further reduces decomposition at depth. For unfrozen soils, the temperature modifier is an exponential function of soil temperature, while below 0 • when liquid water enabling decomposition disappears, respiration linearly drops to zero at −1 • (Koven et al., 2011). The soil moisture modifier is prescribed from the meta-analysis of soil volumetric water content (m 3 m −3 ) and respiration relationship for organic soils conducted by Moyano et al. (2012). See Supplement Text S3 for a more detailed description of the temperature and moisture modifier. Following Koven et al. (2013), we implement a depth modifier (f Z,l ) to represent unresolved depth controls (i.e. priming effects, sorption of organic molecules to mineral surfaces) on C decomposition. This depth modifier decreases exponentially with depth: where z l (m) is the depth of the layer l, and z 0 (m) is the e-folding depth of base decomposition rate.

Vertical buildup of peat
Waterlogging and cold temperatures in northern peatland regions prevent complete decomposition of dead plant material, causing an imbalance between litter production and decay (Parish et al., 2008). The un-decomposed plant residue accumulates as peat, and consequently the peat surface shows an upward growth. Instead of modelling this upward accumulation of peat, we simulate a downward movement of C by adapting the method that Jafarov and Schaefer (2016) used to build up a dynamic surface organic layer. We first calculate the empirical carbon content at each model layer (C obs,l ) according to measured data from 102 peat cores from 73 sites (Lewis et al., 2012;Loisel et al., 2014;McCarter and Price, 2013;Price et al., 2005;Tfaily et al., 2014;Turunen et al., 2001;Zaccone et al., 2011). C obs,l is calculated as where BD l (kg m −3 ) is the soil bulk density at model layer l, which is the median observed bulk density after compiling all bulk density measurements into model depth bins (Fig. S1a). α c,l is the mass fraction of carbon in the soil (% weight) for the layer, derived from a regression of measured carbon fraction on measured bulk density from 39 cores from 29 sites (Fig. S1b). Z l (m) is the thickness of the layer. We then model the vertical downward movement of C between soil layers to mimic the aggradation of carbon in the peat as follows. If carbon in layer l (C l ) exceeds a threshold amount (C th,l ), a prescribed fraction (f ) of the carbon is moved to the layer below (l + 1). Here, the carbon flux from layer l to the layer below (l + 1) is calculated as where C l (kg m −2 ) is the carbon content of layer l. The threshold amount of carbon of layer l (C th,l ) is a prescribed fraction (f th ) of the empirically determined C obs,l : The values of model parameters f and f th do not change with soil depth. Finally, the total peat depth is defined as the depth to which carbon can be transferred: where k is the deepest soil layer where carbon content is greater than 0, C k (kg m −2 ) is the carbon content of layer k, C obs,k (kg m −2 ) is the empirical amount of carbon that layer k can hold, and Z k (m) is the thickness of layer k.

Simulating dynamic peatland area extent
In grid-based simulations, each grid cell is characterized by fractional coverages of PFTs. The dynamic coverage of each non-peatland PFT is determined by the dynamic global vegetation model (DGVM) equations as functions of bioclimatic limitations, sapling establishment, light competition, and natural plant mortality (Krinner et al., 2005;Zhu et al., 2015).
Here, a cost-efficient TOPMODEL from the DYPTOP model (Stocker et al., 2014) is incorporated and calibrated for each grid cell by present-day wetland areas that are regularly inundated or subject to shallow water tables to simulate wetland extent (Sect. 2.2.1). Then, the criteria for peatland expansion are adapted from DYPTOP to distinguish peatland from wetland (Sect. 2.2.2).

The cost-efficient TOPMODEL
Concepts of TOPMODEL (Beven and Kirkby, 1979) have been proven to be effective at outlining wetland areas in current state-of-the-art LSMs (Kleinen et al., 2012;Ringeval et al., 2012;Stocker et al., 2014;Zhang et al., 2016). Based on TOPMODEL, sub-grid-scale topography information and soil properties of a given watershed/grid cell are used to redistribute the mean water table depth to delineate the extent of sub-grid area at maximum soil water content. The empirical relationship between the flooded fraction of a grid cell and the grid cell mean water table position (WT) can be established (Fig. 1a) and approximated by an asymmetric sigmoid function, which is more computationally efficient than determining water table depth for each sub-grid pixel (Stocker et al., 2014). Here, we adopted the costefficient TOPMODEL from Stocker et al. (2014) and calibrated TOPMODEL parameters for each grid cell to match the spatial distribution of northern wetlands (see more details in Text S4). Tootchi et al. (2019) reconciled multiple current wetland datasets and generated several high-resolution composite wetland (CW) maps. The one used here (CW-WTD) was derived by combining regularly flooded wetlands (RFWs), which is obtained by overlapping three openwater and inundation datasets (ESA-CCI; Herold et al. 2015;GIEMS-D15, Fluet-Chouinard et al., 2015;and JRC, Fluet-Chouinard et al., 2015), with areas that have shallow (WT ≤ 20 cm) water tables from groundwater modelling of Fan et al. (2013). CW-WTD wetlands are static and aim at representing the climatological maximum extent of active wetlands and inundation. We therefore compare simulated maximum monthly mean wetland extent over 1980-2015 with CW-WTD to calibrate TOPMODEL parameters. Note that lakes from the HydroLAKES database have been excluded from the CW-WTD map because of their distinct hydrology and ecology compared with wetlands (Tootchi et al., 2019).

Peatland development criteria
The criteria used to constrain peatland area development are greatly inspired by DYPTOP (Stocker et al., 2014), but with some adaptions. The initiation of peatland only depends on moisture conditions of the grid cell (Fig. 1b,. First, only the sub-grid-cell area fraction that is frequently inundated has the potential to become peatland (f pot ). Stocker et al. (2014) introduced a "flooding persistency" parameter (N in Eq. 13 in Stocker et al., 2014) for the DYPTOP model to represent the temporal frequency of inundation. N is a globally uniform parameter in DYPTOP, set to 18 months during the preceding 31 years. However, the formation of peat is a function of local climate, and thus suitable formation conditions for peatland vary between geographic regions. To be specific, the accumulation of peat in arctic and northern latitudes is due to both high water table and low temperature, while it is mainly a result of waterlogging conditions in subtropical and tropical latitudes (Parish et al., 2008). Therefore, it is essential to apply different values for the flooding persistency parameter for different regions, according to local climate conditions. We redefined the requirement of persistent flooding for peatland formation as follows: the area fraction that has the potential to become peatland needs to be flooded at least Num months during the preceding 30 years, with Num being the total number of growing season months (monthly air temperature > 5 • C) in 30 years (Fig. 1b,no. 5). In this case, with the help of relatively low air temperature making shorter growing seasons, arctic and boreal latitudes need shorter inundation periods than subtropical and tropical regions to form peatland. Furthermore, as Sphagnum-dominated peatlands are sensitive to summer moisture conditions (Alexandrov et al., 2016;Gignac et al., 2000), the summer water balance of the grid cell needs to pass a specific threshold (SWB) to form peat and to achieve the potential peatland area (Fig. 1b,no. 7). The summer water balance is calculated as the difference between total precipitation (P ) and total potential evapotranspiration (PET) of May-September. We consider SWB as a tunable parameter in the model and run simulations with SWB = −6, 0, 3, and 6 cm. SWB = 6 cm is selected so that the model captures the southern frontier of peatland in Eurasia and western North America (Text S5). Note that the definitions of summer (May-September) and SWB are not applicable for tropical regions and the Southern Hemisphere.
After the initiation, the development of peatland area is controlled by both moisture conditions of the grid cell and the long-term carbon balance of the peatland HSU (Fig. 1c, nos. 9-17). If the climate becomes drier and the calculated potential peatland area is smaller than the current peatland area, the peatland HSU area will contract to the new potential peatland area fraction (Fig. 1c,no. 12). Otherwise (Fig. 1c, no. 13), the peatland has the possibility to expand when the summer water balance threshold is passed. If these criteria are satisfied, the final decision depends on the carbon density of the peatland (C peat ): the peatland can expand only when long-term input exceeds decay and a certain amount of C (C lim ) has accumulated (Fig. 1c,no. 17). C lim is defined here as a long-term peatland C balance condition, it is a product of a mean measured peat depth (1.07 m) from 40 peat cores (with peat age greater than 1.8 ka but smaller than 2.2 ka) from North American peatland (Gorham et al., 2007(Gorham et al., , 2012 and from the West Siberian lowlands (Kremenetski et al., 2003), with a dry bulk density assumption of 100.0 kg m −3 and a mean C fraction of 47 % in total peat (Loisel et al., 2014). Our estimation for C lim is 50.3 kg C m −2 , which matches well with the C density criterion (50 kg C m −2 ) chosen by Stocker et al. (2014) to represent typical peatland soil.
The moisture conditions are evaluated every month throughout the simulation, while C peat is checked only in the first month after the SubC in Spin-up1 and is checked every month in Spin-up2 and the transient simulation (see Sect. 3.2). The peatland area fraction (f peat ) is updated every month. During the simulation, the contracted area and C are allocated to an "old peat" pool and are kept track of by the model. It should be noted that drainage (drought) may cause a decrease in porosity and saturated moisture content of peat soils (Oleszczuk and Truba, 2013) and changes in peatland vegetation compositions (Benavides, 2014). But the current model structure does not allow us to take these potential changes in peatland into consideration. Therefore, paramewww.geosci-model-dev.net/12/2961/2019/ Geosci. Model Dev., 12, 2961-2982, 2019 terizations of the old peat pool are identical to mineral soils, following the study of Stocker et al. (2014). When peatland expansion happens, the peatland will first expand into this old peat area and inherit its stored C (Stocker et al., 2014).
The difference between our model and the DYPTOP model in simulating peatland area dynamics can be summarized as follows.
(1) For the TOPMODEL calibration, TOP-MODEL parameters are globally uniform in the DYPTOP model, but grid-cell-specific in ORCHIDEE-PEAT v2.0.
(2) For criteria for peatland expansion, in DYPTOP, the flooding persistency parameter is globally uniform, being 18 months in the preceding 31 years. And the ecosystem water balance is expressed as annual precipitation over actual evapotranspiration (POAET). In ORCHIDEE-PEAT v2.0, the flooding persistency parameter is grid-cell-specific, being the total number of growing season months in the preceding 30 years. And peatland expansion is limited only by summer water balance. The relative areal change in peatland is limited to 1 % per year in DYPTOP, but not limited in our model.
(3) For peatland initiation, DYPTOP prescribes a very small peatland area fraction (0.001 %) in each grid cell to simulate peatland C balance condition. Peatland can expand from this "seed" once water and carbon balance criteria are met. In ORCHIDEE-PEAT v2.0, no seed is needed because only the flooding persistency and summer water balance criteria need to be met for the first initiation of peatland (Fig. 1b); carbon balance is only checked after initiation (Fig. 1c).

Critical model parameters
The base decomposition rates of active, slow, and passive peat soil carbon pools in the model are 1.0, 0.027, and 0.0006 a −1 at a reference temperature of 30 • C, respectively (Table 1, Sect. 5: Choice of model parameters). The e-folding depth of the depth modifier (z 0 , Eq. 2) determines the general shape of increases in soil C turnover time with depth; the prescribed threshold to allow downward C transfer between soil layers (f th , Eq. 5) and the prescribed fraction of C to be transferred (f , Eq. 4) determine movement and subsequent distribution of soil C along the soil profile. We compare simulated C vertical profiles with observed C profiles at 15 northern peatland sites (Table S1) (Loisel et al., 2014) using different combinations of parameters (z 0 = (0.5, 1.0, 1.5, 2.0), f th = (0.5, 0.7, 0.9), and f = (0.1, 0.2, 0.3)) and eventually selected z 0 = 1.5 m, f th = 0.7, and f = 0.1 based on visual examinations to match the observed C content. Model sensitivity to the selection will be discussed in Sect. 5.

Simulation protocol
We conduct both site-level and regional simulations with ORCHIDEE-PEAT v2.0 at 1 • × 1 • spatial resolution. Regional simulations are performed for the Northern Hemi-sphere (> 30 • N), while site-level simulations are performed for 60 grid cells containing at least one peat core (Table S1, Fig. S2). Peat cores used in site-level simulations are from the Holocene Perspective on Peatland Biogeochemistry (HPPB) database (Loisel et al., 2014). Both site-level and regional simulations are forced by the 6-hourly meteorological forcing from the CRUNCEP v8 dataset, which is a combination of the CRU TS monthly climate dataset and NCEP reanalysis (ftp://nacp.ornl.gov/synthesis/2009/frescati/temp/ land_use_change/original/readme.htm, last access: 10 July 2019).
All simulations start with a two-step spin-up followed by a transient simulation after the pre-industrial period (Fig. S3). The first spin-up (Spin-up1) includes N cycles of a peat carbon accumulation acceleration procedure consisting of (1) 30 years with the full ORCHIDEE-PEAT (FullO) run on a 30 min time step followed by (2) a stand-alone soil carbon sub-model (SubC) run to simulate the soil carbon dynamics in a cost-effective way on monthly steps (fixed monthly litter input, soil water, and soil thermal conditions from the preceding FullO simulation). Repeated 1961-1990 climate forcing is used in Spin-up1 to approximate the higher Holocene temperatures relative to the pre-industrial period (Marcott et al., 2013). The atmospheric CO 2 concentration is fixed at the pre-industrial level (286 ppm). Each time we run the SubC for 2000 years in the first N − 1 sets of acceleration procedures, and the value of N and the time length of the last set of acceleration procedures (X) are defined according to the age of the peat core in site-level simulations, and are defined according to the reconstructed glacial retreat in regional simulations (Figs. S4, S5). The reconstructed glacial retreat used in this study is from Dyke (2004) for North America and from Hughes et al. (2016) for Eurasia (Text S6).
In the second spin-up step (Spin-up2), the full ORCHIDEE-PEAT model was run for 100 years, forced by looped 1901-1920 climate forcing and pre-industrial atmospheric CO 2 concentration so that physical and carbon fluxes can approach the pre-industrial equilibrium. After the two spin-ups, a transient simulation is run, forced by historical climate forcing from CRUNCEP and rising atmospheric CO 2 concentration. For site-level simulations, the transient period starts from 1860 and ends at the year of coring (Table S1). For regional simulations, the transient period starts from 1860 and ends at 2009.

Evaluation datasets for site-level simulations
All peatland sites used in this study are from the HPPB database (Loisel et al., 2014). All the peat cores measured peat ages and depths (60 sites, Table S1), and are hence used to evaluate simulated peat depth, with sites being grouped into different peatland types, climate zones, and ages. For peat cores where peat ages, depths, fraction of C, and bulk  Table S1), we construct vertical C profiles with this measured information to compare with our simulated C profiles.

Soil organic carbon stocks
Simulated peatland soil organic carbon (SOC) is evaluated against (1) the WISE database (Batjes, 2016) and (2) the IMCG GPD (Joosten, 2010). All the above-mentioned datasets used to evaluate ORCHIDEE-PEAT v2.0 at a regional scale are described in the Supplement Text S7. Gorham et al. (2007Gorham et al. ( , 2012 and Kremenetski et al. (2003) collected depth and age of 1685 and 130 peat cores, respectively, from literature data on peatlands in North America (NA) and in the West Siberian lowlands (WSL). These compilations make it possible for us to validate peat depths simulated by ORCHIDEE-PEAT v2.0 at regional scales, in addition to the detailed site runs in Sect. 3.3.1. Compared to the HPPB database, these datasets lack detailed peat properties (i.e. C content, peatland type), but contain more samples and cover larger areas. Note that as this study aims to reproduce development of northern peatlands since the Holocene, peat cores that are older than 12 ka are removed from the model evaluation. Finally, 1521 out of 1685 observed peat cores in NA and 127 out of 130 observed peat cores in WSL are used in model evaluation (Sect. 4.2: Peat depth).

Site simulation
We first evaluate the performance of ORCHIDEE-PEAT v2.0 in reproducing peat depths and vertical C profiles at the 60 sites from HPPB (Table S1). Out of the 60 grid cells (each grid cell corresponding to one peat core), ORCHIDEE-PEAT v2.0 produces peatlands in 57 of them. The establishment of  (Table S1). Shapes of markers indicate peatland types (bogs, fens, others); colours of markers imply climatic zones (temperate, boreal, arctic) of site locations. peatlands at Zoige, Altay, and IN-BG-1 (Table S1) is prevented in the model by the summer water balance criterion of these grid cells. Peat depths are underestimated for most sites (Fig. 2). Simulated depth of these 57 sites ranges from 0.37 to 6.64 m and shows a median depth of 2.18 m, while measured peat depth ranges from 0.96 to 10.95 m, with the measured median depth being 3.10 m ( Table 2). The rootmean-square error (RMSE) between observations and simulations is 2.45 m.
The measured and simulated median peat depths for the 14 fen sites are 3.78 and 2.16 m, compared to 3.30 and 2.18 m, respectively, for the 33 bog sites ( Table 2). The model shows slightly higher accuracy for fens than for bogs, with the RMSE for fens being 2.08 and 2.59 m for bogs. RMSEs for peat depths of sites that are older than 8 ka are greater than those of younger sites, but are smaller than the measured mean depth (3.5 m) of all peat cores. Simulated median depth of the six arctic sites is larger than observations, but that of the 47 boreal sites and the four temperate sites is smaller than observations ( Table 2). The RMSE for temperate sites is larger than that for arctic or boreal sites.
The simulated and observed vertical profiles of soil C for the 15 sites are shown in Fig. 3, simulated C concentrations are generally within the range of measurements at most of the sites, but are underestimated at Sidney bog, Usnsk Mire 1, Lake 785, and Lake 396. In the model, the buildup of peat is parameterized by downward movement of C between soil layers, with the empirical amount of C that each layer can hold being calculated from median observed bulk density and C fraction of peat core samples (Sect. 2.1.2). High C concentration of cores that have significantly larger bulk density and/or C fraction than the median of the measurements thus cannot be reproduced. This is the case of Lake 785 and Lake 396 (Table S1), where C concentrations are underestimated and depths are overestimated (Fig. 2), while simulated total C content is close to observations (for Lake 785, measured and simulated C content is 86.14 and 96.13 kg C m −2 , respectively, while values for Lake 396 are 57.2 and 70.2 kg C m −2 ).
As shown in Fig. 4, there is considerable variability in depth and C concentration profiles among peat cores within a grid cell, even though these cores have a similar age. We rerun the model at the five grid cells where more than one peat core has been sampled, with time length of the simulation being defined as the mean age of cores in the same one grid cell. The simulated peat depth and C concentration profiles at G2, G4, and G5 are generally within the range of peat core measurements (Fig. 4). Observed C fraction at grid cell G1 and G3 is much greater than the median value of all peat core samples (Sect. 2.1.2); thus simulated C concentration along the peat profile is smaller than observations, but peat depth is still overestimated by the model. This is also the case with Lake 785 and Lake 396.

Northern peatland area and C stock
Simulated maximum inundated area of the Northern Hemisphere is 9.1 million km 2 , smaller than the wetland areas in CW-WTD (∼ 13.2 million km 2 after excluding lakes). TOP-MODEL gives an area fraction at maximum soil water content while CW-WTD includes both areas seasonally to permanently flooded and areas that are persistently saturated or near-saturated (the maximum water table shallower than 20 cm) at soil surface. Therefore, an exact match between CW-WTD and the model prediction is not expected. The model generally captures the spatial pattern of wetland areas represented by CW-WTD (Fig. 5). The multi-sensor satellitebased GIEMS dataset (Prigent et al., 2007, which provides observed monthly inundation extent over the period of 1993-2007, is used to evaluate simulated seasonality of inundation. Figure 6 shows that the seasonality of inundation is generally well captured by the model, although simulated seasonal maximum of inundation extent occurs earlier than observations (except in WSL) and simulated duration of inundation is longer than observations. While our model predicts the natural extent of peatlands under suitable climate conditions, soil formation processes and soil erosion are not included in the model. We mask grid cells that are dominated by Leptosols, which are shallow or stony soils over hard rock, or highly calcareous material (Nachtergaele, 2010) (Figs. S6, S7). Peatlands have been extensively used for agriculture after drainage and/or partial extraction worldwide (Carlson et al., 2016;Joosten, 2010;Leifeld and Menichetti, 2018;Parish et al., 2008). Intensive cultivation practices might cause rapid loss of peat C and ensuing disappearance of peatland. Additionally, agri-Geosci. Model Dev., 12, 2961-2982, 2019 www.geosci-model-dev.net/12/2961/2019/   (Table S1). The black circles indicate depths of measurements; the red circles indicate the depth of each soil layer in the model.  cultural peatlands are often classified as cropland, not as organic soils (Joosten, 2010). Therefore, we masked agricultural peatland from the results by assuming that crops occupy peatland in proportion to the grid cell peatland area (Carlson et al., 2016). The distribution and area of cropland used here is from the MIRCA2000 dataset (Portmann et al., 2010), which provides monthly crop areas for 26 crop classes around the year 2000 and includes multicropping explicitly (Fig. S8). After masking Leptosols and agricultural peatlands from the simulated peatland areas and peatland C stocks, the simulated total northern peatland area is 3.9 million km 2 (f noLEP-CR , Fig. 7d), holding 463 Pg C (C noLEP-CR , Fig. 8b). These estimates fall well within estimated ranges of northern peatland area (3.4-4 million km 2 ) and carbon stock (270-540 Pg C) (Gorham, 1991;Turunen et al., 2002;Yu et al., 2010). Simulated peatland area matches relatively well with PEATMAP data in Asian Russia but overestimates peat area in European Russia (Table 3). The simulated total peatland area of Canada is in relatively good agreement with the three evaluation datasets, though the world's second largest peatland complex at the Hudson Bay lowlands (HBL) is underestimated and a small part of the northwest Canada peatlands is missing. Packalen et al. (2014) stressed that initiation and development of HBL peatlands are driven by both climate and glacial isostatic adjustment (GIA), with initiation and expansion of HBL peatlands tightly coupled with land emergence from the Tyrrell Sea, following the deglaciation of the Laurentide ice sheet and under suitable climatic conditions. The pattern of peatlands at southern HBL was believed to be driven by the differential rates of GIA rather than climate (Glaser et al., 2004a, b). More specifically, Glaser et al. (2004a, b) suggested that the faster isostatic uplift rates on the lower reaches of the drainage basin reduce regional slope, impede drainage and shift river channels. Our model, however, cannot simulate the tectonic and hydrogeologic controls on peatland development. In addition, the development of permafrost at depth as peat grows in thickness over time acts to expand peat volume and uplift peat when liquid waterfilled pores at the bottom of the peat become ice-filled pores (Seppälä, 2006). This process is not accounted for in the model and may explain why the HBL does not show up as a large flooded area today whereas peat developed in this region during the early development stages of the HBL complex. The simulated distribution of peatland area in Alaska agrees well with Yu et al. (2010) and WISE. There is a large overestimation of peatland area in the southeastern US (Table 3, Fig. 7d). The simulated peat C stock in Russia (both the Asian and the European part) and in the US is overestimated compared to IMCG GPD and WISE, but that of Canada is underestimated (Table 4, Fig. 8b). Figure 9 shows measured and simulated peat depth in NA and WSL. Some peat cores are sampled from the Canadian Arctic Archipelago, southwestern US, and the northern tip of Quebec, where there is no peatland in peat inventories or the soil database. These sites support the notion that the formation and development of peatland are strongly dependent on local conditions, i.e. retreat of glaciers, topography, drainage, vegetation succession (Carrara et al., 1991;Madole, 1976   As shown in Fig. 4, within a grid cell, sampled peat cores can have very different depths and/or ages. We calculate the mean depth of cores in each of the grid cells and compare it against the simulated mean depth. The mean age of cores in each of the grid cells is used to determine which output of the model should be examined. For instance, the mean age of the four cores in grid cell (40.5 • N, 74.5 • W) is 2.5 ka, and accordingly we pick out the simulated depth of this grid cell right after the first run of SubC (Fig. S3) to compare with the mean depth of these cores. We acknowledge that   this is still a crude comparison since the simulation protocol implies that we can only make the comparison at 2000-year intervals. Nonetheless, it is a compromise between running the model for 1815 peat cores independently and comparing the mean depth of measured points with grid-based simulated depth. As shown in Fig. 10, for each age interval (of both the West Siberian lowlands and North America), the variation in simulated depth is smaller than that in the measurement. The two deepest simulated peat measurements in WSL belong to the fourth age group (6 < Age ≤ 8 ka) and are the result of a shallow active layer; while C is moving downward to deeper and deeper layers, the decomposition is greatly limited by cold conditions at depth. At both WSL and NA, simulated median peat depths (2.07-2.36 m at WSL, 1.02-2.15 m at NA) are in relatively good agreement with measurements (1.8-2.31 m at WSL, 0.8-2.46 m at NA) for cores younger than 8 ka (Fig. 10). For the two oldest groups (peat age > 8 ka), the simulated median depths are about 0.70 m shallower than measurements at NA and about 1.04 m shallower at WSL.

Undisturbed northern peatland carbon balance in the past century
Simulated mean annual (averaged over 1901-2009) net ecosystem production (NEP) of northern peatlands varies from −63 to 46 g C m −2 a −1 (Fig. 11). The West Siberian lowlands, the Hudson Bay lowlands, Alaska, and the China-Russia border are significant hotspots of peatland C uptake. Simulated mean annual NEP of all northern peatlands over 1901-2009 is 0.1 Pg C a −1 , consistent with the previous estimate of 0.076 Pg C a −1 by Gorham (1991) N. (c, d) The scatter plot of measured and simulated peat depth for the West Siberian lowlands (c) and North America (d). For a grid cell that has multiple measured peat cores, the median depth of all measurements is plotted against the simulated depth in the scatter plot.
heterotrophic respiration (HR) show an increasing trend, but NPP rises faster than HR during the second half of the cen- Figure 11. Simulated annual net ecosystem production (NEP), averaged over 1901-2009. Obtained by multiplying peatland NEP (g C m −2 peatland a −1 ) with peatland fraction for each grid cell. tury (Fig. 12a). The increase in NPP is caused by atmospheric CO 2 concentration and increasing air temperature (Figs. 12, S9). As air (soil) temperature increases, HR also increases but lags behind NPP (Figs. 12, S9). Simulated annual NEP ranges from −0.03 to 0.23 Pg C a −1 , with a significant positive trend over the second half of the century (Fig. 12b). NEP shows a significant positive relationship with air (soil) temperature and with atmospheric CO 2 concentration (Fig. S9). CH 4 and dissolved organic carbon (DOC) are not yet included in the model; both of them are significant losses of C from peatland (Roulet et al., 2007).

Peat depth
We found a general underestimation of peat depth (Figs. 2, 10), possibly due to the following reasons. Firstly, there is a lack of specific local climatic and topographic conditions. The surfaces of peatlands are mosaics of microforms, with accumulation of peat occurring at each individual microsite of hummocks, lawns, and hollows. Differences in vegetation communities, thickness of the unsaturated zone, local peat hydraulic conductivity, and transmissivity between microforms result in considerable variation in peat formation rate and total C mass (Belyea and Clymo, 2001;Belyea and Malmer, 2004;Borren et al., 2004;Packalen et al., 2016). Cresto Aleina et al. (2015) found that the inclusion of microtopography in the hummock-hollow model delayed the simulated runoff and maintained wetter peat soil for a longer time at a peatland of northwest Russia, thus contributing to enhanced anoxic conditions. Secondly, site-specific parameters are not included in gridded simulations. Parameters describing peat soil properties, i.e. soil bulk density and soil carbon fraction, determine the amount of C that can be stored across the vertical soil profile. Hydrological parameters, i.e. the hydraulic conductivity and diffusivity, and the saturated and residual water content regulate vertical fluxes of water in the peatland soil and expansion-contraction of the peatland area, and hence influence the decomposition and accumulation of C at the sites considered. Plant trait parameters, i.e. the maximal rate of carboxylation (V cmax ), and the light saturation rate of electron transport (J max ) determine the carbon budgets of the sites . The depth modifier, which parameterizes depth dependence of decomposition, controls C decomposition at depth and is an important control on simulated total C and the vertical C profile. A third reason is sample selection bias. Ecologists and geochemists tend to take samples from the deepest part of a peatland complex to obtain the longest possible records (Gorham, 1991;Kuhry and Turunen, 2006). In contrast, the model is designed to model an average age and C stock of peatlands in a grid location, and thus preferably the simulated C concentrations of a grid cell should only be validated against grids represented by a number of observed cores. We do try to compare the model output with multiple peat cores (Figs. 4, 10), but we need to note that shallow peat is not sufficiently represented in field measurements. A fourth source of error is that simulated initiation time of peat development at some sites is too late compared to ages of measured cores. The model multiple-spin-up strategy accounts for coarse-scale ice sheet distribution at discrete Holocene intervals (Sect. 3.2, Fig. S3), and if the modelled occurrence of peatland is too late, the accumulated soil C may be underestimated. For example, at the Patuanak site, where the core age is 9017 years, the model was run with 4 times' SubC (Table S1). However, there was no peatland before the first SubC, meaning that simulated peatland at this grid cell was 2000 years younger than the observation and that our simulation missed C accumulation during the first 2000 years at this site. This may be another source of bias associated with the model resolution, namely that local site conditions fulfilled the initiation of peatland at specific locations, but the average topographic and climatic conditions of the coarse model grid cell were not suitable for peatland initiation. Also, one has to keep in mind that a single (a few) sample(s) from a large peat complex may not be enough to capture the lateral spread of peat area, which may be an important control on accumulation of C (Charmen, 1992;Gallego-Sala et al., 2016;Parish et al., 2008). The underestimation of peat depth can also come from biased climate input data: spin-ups of the model are forced with repeated 1961-1990 climate, assuming that Holocene climate is equal to recent climate. While peatland carbon sequestration rates are sensitive to climatic fluctuations, centennial-to millennial-scale climate variability, i.e. cooling during the Younger Dryas period and the Little Ice Age period and warming during the Bølling-Allerød period, is not included in the climate forcing data (Yu et al., 2003a, b). An early Holocene carbon accumulation peak was found during the Holocene Thermal Maximum when the climate was warmer than present (Loisel et al., 2014;Yu et al., 2009). Finally, effects of landscape morphology on drainage as well as drainage of glacial lakes are not incorporated and can represent a source of uncertainty.

Vertical profiles of peatland soil organic carbon
We note that caution is needed in interpreting the comparison between simulated peat C profile and measured C profile from peat cores (Figs. 3, 4). In reality, peat grows both vertically and laterally since inception, with peat deposits tending to be deeper and basal age tending to be older at the original nucleation sites/center of the peatland complex (Bauer et al., 2003;Mathijssen et al., 2017). As mentioned earlier, field measurements tend to take samples from the deeper part of a peatland complex and shallow peat is underrepresented. The model, however, only simulates peat growth in the vertical dimension and lacks an explicit representation of the lateral development of a peatland in grid-based simulations; thus simulated peat C (per unit peatland area) is diluted when the simulated peatland area fraction in the grid cell increases.
In addition, we cannot compare the simulated peat C profile Geosci. Model Dev., 12, 2961-2982, 2019 www.geosci-model-dev.net/12/2961/2019/ against the observed profile from dated peat cores because the model does not track age bins explicitly. The above-noted discrepancies between the simulation and the observation highlight both the need for more peat core data collected with more rigorous sampling methodologies and the need to improve the model. In parallel with this study, 14 C dynamics in the soil have been incorporated into the ORCHIDEE-SOM model , which may give us an opportunity to compare simulated 14 C age-depth profiles with dated peat C profiles in the future after being merged with our model.

Simulated peatland area development
The initiation and development of peatlands in NA followed the retreat of the ice sheets, as a result of the continuing emergence of new land with the potential to become suitable for peatland formation (Gorham et al., 2007;. To take glacial extent into account for simulating the Holocene development of peatlands, we use ice sheet reconstructions in NA and Eurasia (Figs. S4,S5). Not surprisingly, when ice cover is considered, the area of peatlands that developed before 8 ka is significantly decreased, while the area that developed after 6 ka is increased (Fig. 13). We use observed frequency distribution of peat basal age from MacDonald et al. (2006) as a proxy of peatland area change over time, following the assumption proposed by Yu (2011) that peatland area increases linearly with the rate of peat initiation. We grouped the data of MacDonald et al. (2006) into 2000-year bins to compare with simulated peatland area dynamics (Fig. 13). The inclusion of dynamic ice sheet coverage triggering peat inception clearly improved the model performance in replicating peatland area development during the Holocene, though the peatland area before 8 ka is still overestimated by the model in comparison with the observed frequency distribution of basal ages (Fig. 13). In spite of the difference in peatland area expansion dynamics between the simulation that considered dynamic ice sheets and the one that did not, the model estimates of present-day total peatland area and carbon stock are generally similar (Fig. S10). Without dynamic ice sheets, the model would predict only 0.1 million km 2 more peatland area and 24 Pg more peat C over the Northern Hemisphere (> 30 • N). We are aware of two studies that attempted to account for the presence of ice sheets during the Holocene (Kleinen et al., 2012) and the Last Glacial Maximum  while simulating peatland C dynamics. Kleinen et al. (2012) modelled C accumulation over the past 8000 years in the peatland areas north of 40 • N using the coupled climate-carbon cycle model CLIMBER2-LPJ. A decrease of 10 Pg C was found when ice sheet extent at 8 ka BP (from the ICE-5G model) was accounted for. Another peatland modelling study conducted by Spahni et al. (2013) with the Land surface Processes and eXchanges (LPX) model also prescribed ice sheets and land area from the ICE-5G ice sheet reconstruction (Peltier, The peatland carbon density criterion for peatland expansion (C lim ) is an important factor impacting the simulated Holocene trajectory of peatland development. Without the limitation of C lim , a larger expansion of northern peatlands would occur before 10 ka (Fig. S11). Such a premature, "explosive" increase in peatland area would result in the overestimation of C accumulated in the early Holocene in the model. In the meantime, peatland area in regions that only have small C input, i.e. Baffin Island and northeast Russia, would be overestimated (Fig. S12).

Choice of model parameters
For the active, slow, and passive peat soil carbon pool, the base decomposition rates are 1.0, 0.027, and 0.0006 a −1 at a reference temperature of 30 • C, respectively, meaning that the residence times at 10 • C (no moisture and depth limitation) of these three pools are 4, 148, and 6470 years. In equilibrium/near-equilibrium state, simulated C in the active pool takes up only a small fraction of the total peat C, while generally 40 %-80 % of simulated peat C is in the slow C pool and about 20 %-60 % is in the passive C pool. Assuming that in a peatland, the active, slow, and passive pools account for 3 %, 60 %, and 37 % (median values from the model output of the year 2009) of the total peat C, we can get a mean peat C residence time of 2500 years. If depth modifier is considered, the C residence time will vary from 2500 years at the soil surface to 13 200 years at the 2.5 m depth. For the record, in previous published large-scale diplotelmic peatland models, at 10 • C, C residence time for the acrotelm (depth = 0.3 m) ranged from 10 to 33 years and ranged from 1000 to 30 000 years for the catotelm (Kleinen et al., 2012;Spahni et al., 2013;Wania et al., 2009b). We performed sensitivity tests to show the sensitivity of the modelled peat C to model parameters at the 15 northern peatland sites where observed vertical C profiles can be constructed (Table S1). Tested parameters are the e-folding decreasing depth of the depth modifier (z 0 , Eq. 2), the prescribed thresholds to start C transfer between soil layers (f th , Eq. 5), and the prescribed fraction of C transferred vertically (f , Eq. 4). We found that changing f th or f leads to only small effects on the vertical soil C profile (see, e.g. Burnt Village peat site in Fig. S13). The parameter z 0 , by contrast, exerts a relatively strong control over C profiles. It is noteworthy that while our model resolves water diffusion between soil layers according to the Fokker-Planck equation , simulated soil moisture does not necessarily increase with depth (Fig. S14). z 0 is therefore an important parameter to constrain peat decomposition rates at depth. With smaller z 0 , decomposition of C decreases rapidly with depth, resulting in a deeper C profile (Fig. S14). Regional-scale tests verified these behaviours of the model: when f th = 0.9 is used (instead of f th = 0.7), changes in peatland area and peat C stock are negligible (Fig. S15). Without z 0 , simulated northern peatland area will not change (3.9 million km 2 ), but northern peatland C stock will be underestimated (only 300 Pg C). If z 0 = 0.5 m is applied (instead of z 0 = 1.5 m), the simulated total peat C would triple while the total peatland area would only increase by 0.2 million km 2 (Fig. S16).

Uncertainties in peatland area and soil C estimations
There are large uncertainties in estimates of peatland distribution and C storage. Some studies prescribe peatlands from wetlands. However, in spite of the fact that there are extensive disagreements between wetland maps, it is a challenge to distinguish peatlands from non-peat-forming wetlands (Gumbricht et al., 2017;Kleinen et al., 2012;Melton et al., 2013;Xu et al., 2018). Estimates based on peatland inventories are impeded by poor availability of data, nonuniform definitions of peatlands among regions, and coarse resolution (Joosten, 2010;Yu et al., 2010). In addition, as peatlands are normally defined as waterlogged ecosystems with a minimum peat depth of 30 or 40 cm, shallow peat is underrepresented. Another approach to estimate peatland area and peat C is to use a soil organic matter map to outline organic-rich areas, such as Histosols and Histels (Köchy et al., 2015;Spahni et al., 2013). This approach overlooks local hydrological conditions and vegetation composition (Wu et al., 2017). Our model estimates of peatland area and C stock generally fall well within the range of published estimates, except in the southeastern US, where there is only 0.05-0.10 million km 2 of peatland in observations but 0.37 million km 2 in the model prediction (Fig. 7d, Table 3). From the early 1600s to 2009, ∼ 50 % of the original wetlands in the lower 48 states of the US have been lost to agricultural, urban development, and other development (Dahl, 2011;Tiner, 1984). Although wetlands are not necessarily peatlands, the reported losses of wetlands in the US indicate that a potentially large area of peatlands in the US may have been lost to land use change. However, historical losses of peatlands due to land use change and the impact of agricultural drainage of peatlands have not been taken into account by our model. Another factor that might have contributed to the overestimation is a limitation of TOPMODEL, namely that the "floodability" of a pixel in the model is determined by its compound topographic index (CTI) value regardless of the pixel's location along the stream, and thus the floodability of an upstream pixel with a large CTI might be affected by downstream pixels that have a small CTI. The model's inability to resolve small-scale streamflows might be another cause of the overestimation. The simulated mean annual NPP, HR, and NEP of northern peatlands increase from about 1950 onwards. We find positive relationships between NPP and temperature, NPP and atmospheric CO 2 concentration, and HR and temperature over the past century (Fig. S9). From a future perspective, it is unclear whether the increasing trend of NEP can be maintained. While photosynthetic sensitivity to CO 2 decreases with increasing atmospheric CO 2 concentration and photosynthesis may finally reach a saturation point in the future, decomposition is not limited by CO 2 concentration and may continue to increase with increasing temperature (Kirschbaum, 1994;Wania et al., 2009b).
Our model applies a multilayer approach to simulate process-based vertical water fluxes and dynamic C profiles of northern peatlands and highlights the vertical heterogeneities in the peat profile in comparison to previous diplotelmic models (Kleinen et al., 2012;Spahni et al., 2013;Stocker et al., 2014;Wania et al., 2009b). While simulating peatland dynamics, large-scale models used a static peatland distribution map obtained from peat inventories and soil classification maps Wania et al., 2009b, a), prescribed the trajectory of peatland area development over time , or used wetland area dynamics as a proxy (Kleinen et al., 2012). We adapt the scheme of DYP-TOP to simulate spatial and temporal dynamics of northern peatland area by combing simulated inundation and a set of peatland expansion criteria (Stocker et al., 2014).
As a large-scale LSM which is designed for large-scale gridded applications, ORCHIDEE-PEAT v2.0 cannot explicitly model the lateral development of a peatland. The model therefore aims to simulate large-scale average peat depth and C profile rather than capturing local peat inception time and age-depth profiles at the location of specific peat cores. Tracers like 14 C are not included in the model, making some site-to-site evaluation in particular regarding peat inception time and age-depth profiles of peat cores difficult. For tropical peatlands, the model needs to be improved to represent its tree dominance, oxidation of deeper peat due to pneumatophore (breather roots) of tropical trees, and the greater water table fluctuations as a result of the higher hydraulic conductivity of wood peats and tropical climates (Lawson et al., 2014). In addition, tropical peat is formed in riparian seasonally flooded wetlands with water coming from upstream river networks, whereas the TOPMODEL equations used here implicitly assume a peatland is formed in a grid cell only from rainfall water falling into that grid cell. Further work to improve this simulation framework is needed in areas such as an accurate representation of the Holocene climate, higher spatial resolution, and distinguishing bogs from fens to better parameterize water inflows into peatland. Including CH 4 emissions and leaching of DOC will be helpful to get a more complete picture of peatland C budget.

Conclusions
Multilayer schemes have been proven to be superior to simple box configurations in ESMs at realistic modelling of energy, water, and carbon fluxes over multilayer ecosystems (De Rosnay et al., 2000;Jenkinson and Coleman, 2008;Best et al., 2011;Wu et al., 2016). We apply multilayer approaches to model vertical profiles of water fluxes and vertical C profiles of northern peatlands. In addition to representations of peatland hydrology, peat C decomposition, and accumulation, a dynamic model of peatland extent is also included. The model shows good performance at simulating average peat depth and vertical C profile in grid-based simulations. Modern total northern peatland area and C stock is simulated as 3.9 million km 2 and 463 Pg C (Leptosols and agricultural peatlands have been masked), respectively. While this study investigated the capability of ORCHIDEE-PEAT v2.0 to hindcast the past, in ongoing work, the model is being used to explore how peatland area and C cycling may change under future climate scenarios.
Author contributions. CQ implemented peatland water and carbon processes into ORCHIDEE-MICT, introduced the dynamic peatland area module, and performed the simulation. DZ contributed to ensuring consistency between the peatland modules and various other processes and modules in the model. PC conceived the project. PC, BG, GK, DZ, and CQ contributed to improving the research and interpreting results. SP assisted in implementing the cost-efficient TOPMODEL. AT and AD provided the dataset of wetland areas. SP, AT, AD, and AH contributed to the calibration of the TOPMODEL. All authors contributed to the paper.