Transient simulations of the carbon and nitrogen dynamics in northern peatlands: from the Last Glacial Maximum to the 21st century

. The development of northern high-latitude peatlands played an important role in the carbon (C) balance of the land biosphere since the Last Glacial Maximum (LGM). At present, carbon storage in northern peatlands is substantial and estimated to be 500 ± 100 Pg C (1 Pg C = 10 15 g C). Here, we develop and apply a peatland module embedded in a dynamic global vegetation and land surface process model (LPX-Bern 1.0). The peatland module features a dynamic nitrogen cycle, a dynamic C transfer between peatland acrotelm (upper oxic layer) and catotelm (deep anoxic layer), hydrology- and temperature-dependent respiration rates, and peatland speciﬁc plant functional types. Nitrogen limitation down-regulates average modern net primary productivity over peatlands by about half. Decadal acrotelm-to-catotelm C ﬂuxes vary between − 20 and + 50 g C m − 2 yr − 1 over the Holocene. Key model parameters are calibrated with reconstructed peat accumulation rates from peat-core data. The model reproduces the major features of the peat core data and of the observation-based modern circumpolar soil carbon distribution. Results from a set of simulations for possible evolutions of northern peat development and areal extent show that soil C stocks in modern peatlands increased by 365–550 Pg C since the LGM, of which 175–272 Pg C accumulated between 11 and 5 kyr BP. Furthermore, our simulations suggest a persistent C sequestration rate of 35–50


Introduction
Northern high-latitude peatlands represent a substantial carbon (C) pool of the terrestrial biosphere. Peat mainly consists of partially decomposed, plant-derived organic matter. Peatlands often form under wet conditions, where the water table typically is close to the surface limiting the supply of oxygen to the soil. The anoxic conditions in the waterlogged soil significantly reduce the decomposition rates of the plant material, and thus allow for peat accumulation over long timescales. Due to its extremely slow decay, accumulated C in present-day peatland dates back to the end of the Last Glacial Maximum (LGM), about 16 500 years (yr) before present (BP) (MacDonald et al., 2006). Over this long period, peatlands can form organic soil layers of several meters in depth. This amounts to an exceptionally high soil C density (C storage per unit area) and a total of 500 ± 100 Pg C (1 Pg C = 10 15 g C) stored in the northern high-latitude peatlands in spite of their limited extent (Yu, 2012). Peat soils (histosols) can also occur in areas of permafrost (histels) and contribute substantially to the total circumpolar soil C inventory (Tarnocai et al., 2009).
The fate of this significant C pool under future climate change is unclear (Davidson and Janssens, 2006 hand, projected warmer temperatures prolong the growing season for peatland vegetation, and thus increase peatland net primary productivity (NPP) and C accumulation (Loisel et al., 2012). The projected increase in precipitation over peatlands may further favor peat accumulation. On the other hand, warmer temperatures enhance water loss through evaporation and transpiration, which in turn affects water table position and eventually NPP. This process seems to be an important factor for peatland modeling (Clark et al., 2010). In addition, warmer temperatures increase soil C decomposition and emissions of the greenhouse gases (GHG), including carbon dioxide (CO 2 ) and methane (CH 4 ). Many other processes and characteristics, e.g. peatland hydrology, nutrient availability, plant species composition, and microbial community, can have an impact on the net C balance in peatlands. Peatland modeling offers a way to analyze the dynamic interaction of these processes in the past, and allows to project changes in peatland C stocks in the future.
Peat growth modeling has a long history and most recent peatland models are related to the simple conceptual model by Clymo (1984). For example, recent models representing peat evolution over the Holocene (Bauer, 2004;Heinemeyer et al., 2010;Frolking et al., 2010) include feedbacks and influences of hydrology, plant communities, and peat properties on peat accumulation. These models are applied at individual peatland sites and simulation results were then compared with peat-core data. However, in these models the ecohydrological feedback is missing, that is, effects of hydrology on peatland vegetation, litter quality and peat decomposition that in turn affect water table position. Global dynamic vegetation models like the Lund-Potsdam-Jena (LPJ) model have been extended for the simulation of a peatland C cycle (Wania et al., 2009a;Wania et al., 2009b;Kleinen et al., 2012). These LPJ models assign specific soil C pools to the acrotelm, a near surface soil layer with varying water table, and the catotelm, a deeper soil layer with permanent water saturation. For northern peatlands an important aspect is also the simulation of freezing and thawing in peat soils, especially in the permafrost area (Wania et al., 2009a). In the case of LPJ-WHy used in Wania et al. (2009b), the peatland vegetation includes two specific plant functional types (PFTs) that are adapted for water saturated or inundated environments. However, none of these models included a dynamic N cycle that could limit plant production in these often nutrient-poor environments . Ecophysiological theory and remote-sensing data of evapotranspiration and plant production suggest that global terrestrial plant productivity is reduced on average of 16-28 % owing to N limitations (Fisher et al., 2012).
Here we present a new peatland module with the implementation of a dynamic N cycle and a dynamic acrotelm-tocatotelm C transfer process into the framework of the Land surface Processes and eXchanges (LPX) model (Sect. 2.2). LPX peatland parameters are then calibrated against reconstructed peat accumulation rates since the LGM from peat-cores in northern high-latitudes regions (Yu et al., 2009, Sect. 4). We apply LPX in a transient simulation of global peatland development since the LGM in the circum-Arctic region that illustrates the evolution of the persistent C sink in the terrestrial biosphere (Sect. 5). From these simulations we evaluate present-day peatland N fluxes and pools (Sect. 5.1) and peatland C pools (Sect. 5.3) against observations. Finally, we apply LPX to future scenarios until 2100 AD, driven by simulated climate based on two representative concentration pathways (RCP 2.6 and RCP 8.5;Sect. 5.5). Results are discussed in Sect. 6 and concluded in Sect. 7.
The peatland module in LPX is based on LPJ-WHy developed by Wania et al. (2009a) and Wania et al. (2009b), and further improved as described in Sect. 2.2. Figure 1 shows a scheme of C fluxes and pools in this module. C is assimilated mainly by two PFTs: peat mosses (Sphagnum spp.) and flood tolerant C3 graminoids (e.g., Carex spp., Eriophorum spp., Juncus spp., Typha spp.). Other PFTs can grow, but are limited in gross primary productivity (GPP), due to stress of standing water in peatlands (Wania et al., 2009b), which affects all major C fluxes of these PFTs. After subtracting plant growth and maintenance respiration, a fraction of C from NPP is transformed to exudates in the rhizosphere (exufrac = 17.5 % of NPP; Wania et al., 2009b). The reduced NPP (82.5 %) is then allocated to plant C pools as leaves and stems for Sphagnum mosses, leaves and roots for graminoids, and leaves, roots, sapwood and heartwood for trees (Sitch et al., 2003;Wania et al., 2009b).  The water table depth  (d wt ) can vary in the range (light blue area) between the acrotelm depth (d acro = 0.3 m) and the maximum of standing water above the surface (0.1 m). The catotelm is considered to be always water saturated (dark blue area). For the litter pools "ag" and "bg" denote aboveground and belowground, respectively. While in LPX leaf and wood turnover C enters the aboveground litter pool, root turnover C is added to the belowground litter pool. Thus, the fractional plant cover and related NPP of Sphagnum mosses versus graminoids directly modulate the proportion of C input to these two litter pools (Fig. 1). Both litter pools undergo oxic and anoxic organic matter decomposition, where C is partially respired to the atmosphere (atmfrac; Sitch et al., 2003) and partially moved to the soil C pools (soilfrac = 100 % − atmfrac). In mineral soils, the majority of the soil organic matter (SOM) input from litter decomposition enters the intermediate C pool (fastfrac = 98.5 %) with a relatively fast turnover rate, while the remainder enters the C pool with a relatively slow turnover rate (slowfrac = 1.5 %; Sitch et al., 2003). For peatlands we follow the approach of Wania et al. (2009b), and identify the acrotelm as the intermediate soil C pool, while the anoxic catotelm is represented by the slow soil C pool (Fig. 1). In this approach acrotelm and catotelm C layers are not treated as peat cohorts, but simulated as uniform C pools. Aerobic and anaerobic decomposition of exudates, litter and soil C is parameterized by individual rates (k) that are modified by soil temperature (R T ) and water-filled pore space (R water ): (1) Table 1 lists an overview of the maximum C decomposition rates at 10 • C (k 10 ) for the individual C pools and of different land units used in LPX. Note that in the litter and soil pools, we do not further distinguish C components by their original PFT. All plant type remnants have the same decomposition rate. Also, we do not consider the priming effect of fresh litter C on decomposition rate of peat entering the catotelm. With the addition of a dynamical N cycle to LPX (Stocker et al., 2013), the fraction of litter C respired to the atmosphere has been reduced from 70 to 60 % and more C enters the soil pools (soilfrac = 40 %). In order to match observed patterns of C content in soils (Tarnocai et al., 2009;Batjes, 2008), the decomposition rates were thus slightly tuned again for natural and agricultural mineral soils (Table 1). Decomposition rates for peatlands are calibrated in this study (see Sect. 4). The temperature response function R T is given according to Lloyd and Taylor (1994) and Sitch et al. (2003) for all C pools. For peatlands, we modify R water in the acrotelm according to the water table depth and soil moisture (see Sect. 2.2.1), and keep it constant for exudates, litter and catotelm C (R water = R moist = 0.35; Wania et al., 2009b). In addition to model changes for acrotelm decomposition, we further developed the peatland module in LPX including several new processes and mechanisms. A simple N cycle in peatlands is described in Sect. 2.2.2. Section 2.2.3 describes the process that acrotelm-to-catotelm C transfer rates dynamically depend on acrotelm C pool size.

Dependency of acrotelm decomposition rates on hydrology
The parameterization of acrotelm decomposition includes a temperature modifier, based on soil temperature at 0.25 m depth, and a constant soil moisture modifier of 0.35 as in LPJ-WHy (Wania et al., 2009b To include the effect of varying water table depth on acrotelm decomposition in peatlands, we modify the acrotelm decomposition rate as a weighted average of moisture dominated (R moist , aerobic) and water saturated (R sat , anaerobic) regimes. Parameterizations follow the approach by Frolking et al. (2010): where d l and w l are individual acrotelm layer depth intervals (above the water table depth) and soil moisture contents, respectively. w opt (= 0.45) is the soil moisture content with optimal conditions for aerobic heterotrophic respiration (Frolking et al., 2010). c moist (= 2.15) is a scaling parameter that results in R moist = 0.35 for a completely saturated acrotelm (Wania et al., 2009a) and in R moist = 0.56 for w l = 0.
where R ca (= 0.001) is the ratio of catotelm to acrotelm decomposition rates, and c sat (= 1.15 m) is a parameter of the scaling length, here representing the distance from surface to mid-depth of the catotelm in the model. If the water table drops to the lowest level in the acrotelm (d wt = d acro ), R sat approaches its maximum of 0.35, which is equal to the ratio of anaerobic to aerobic respiration (Wania et al., 2009b). Combining both respiration modifiers results in the weighted mean acrotelm respiration modifier by water-filled pore space (R water ): The parameterizations in the above equations shows that R water will vary between ∼ 0.30 (water table above surface) and 1.0 (water table at bottom of acrotelm and optimal soil moisture content).

Nitrogen cycle in peatlands
Peatlands can be classified as bogs and fens according to their nutrient sources. Bogs mainly depend on the atmosphere as a nutrient source (ombrotrophic), whereas fens are fed mainly by surface water or groundwater (minerotrophic). In our large-scale modeling approach presented here (with grid-cell area of ∼ 54 × 10 3 km 2 ) we do not simulate explicit lateral water transport within or between grid cells, and thus do not explicitly distinguish between fens and bogs. However, we include both atmospheric N deposition and an implicit N source as N inputs to peatlands. The latter can be interpreted as N 2 fixation and other N input processes that allow LPX to simulate conditions as in minerotrophic ecosystems . Thus, the model represents average large-scale landscape C and N cycling. Furthermore, it dynamically calculates the competition among PFTs, also depending on N cycling and subsequent N availability. The simulation of N cycling in peatlands follows the approach implemented for mineral soils (Xu-Ri and Prentice, 2008;Xu-Ri et al., 2012;Stocker et al., 2013). N uptake is governed by NPP and PFT-specific C : N ratios of new biomass production. It is calculated as the minimum of N uptake capacity and N availability, where the former has a soil temperature dependency (Xu-Ri and Prentice, 2008). Although PFTs can be N limited, they do not directly compete for N. Thus, all PFTs are equally limited and no PFT has advantage to access mineral N. This limitation has then a feedback on the C balance of each PFT, in addition to the limitation by water and light. Specifically for peatlands we adapted the C : N mass ratios of new production to 30.0 and 54.29 for Sphagnum and flood-tolerant C3 graminoids, respectively, and the C : N ratios of SOM to 55.0 when derived from Sphagnum mosses and 50.0 from flood-tolerant C3 graminoids (McGuire et al., 1992;Melillo et al., 1993;Heijmans et al., 2001;Limpens et al., 2006). Note that the ultimate C : N ratio of SOM is an average, weighted by the relative inputs from these PFTs. N allocation to leaf, roots and sapwood follows Xu-Ri and Prentice (2008). N is then cycled through plant, litter and SOM pools in parallel with C, but soil C : N ratios are held constant over time. This implies that an additional N source to the peatland ecosystem exists as N 2 fixation by bacteria or algae in peatlands. N is mineralized in parallel to litter and soil decomposition ( Fig. 1) and enters the ammonium (NH + 4 ) pool, where it is subject to volatilization, leaching, nitrification to nitrates (NO − 3 ) and subsequent denitrification (Li et al., 1992;Xu-Ri and Prentice, 2008).

Dynamic acrotelm-to-catotelm C transfer rates
In LPJ-WHy the average acrotelm-to-catotelm C transfer rate is held constant at 12 g C m −2 yr −1 (Wania et al., 2009b). Within the LPX peatland module presented here we now define a dynamic C transfer depending on actual accumulation and decomposition rates in the acrotelm. While the water table varies in the acrotelm (top 0.3 m), the deeper catotelm is considered to be water saturated at all times in the model. Thus, the net C balance of the acrotelm leads to a C flux between the acrotelm and catotelm C pools (F AC ): where C acro is the acrotelm C, F LA is the litter input to the acrotelm ( Fig. 1) and F AR the acrotelm respiration as in Eqs. (1)-(4): If accumulated C acro exceeds average acrotelm C pool size (C AAS ) over the course of a year, F AC is transferred from the acrotelm to the catotelm at a daily rate the following year. C AAS is calculated from a fixed acrotelm depth (d acro = 0.3 m) and an average carbon density (ρ acro = 18.7 kg C m −3 derived from Belyea and Clymo, 2001): If acrotelm C respires due to dry climatic conditions, F AC is negative and represents the deficit C (difference to C AAS ) that is transferred from the catotelm to the acrotelm. So this implies that a peatland can shrink in its pool size and could even disappear if conditions are too dry. We are aware that peat bulk density and decomposition rates of C returning from the catotelm are different from the overlying acrotelm layers. For example, peat compacts when it is submerged under the water table, resulting in greater bulk density, and SOM becomes more recalcitrant over time. This is clearly a model simplification that plays a role for the rate at which peatlands can disappear (Leifeld et al., 2012).

Climate input data for simulations from the LGM to present
In order to simulate peatland development in the northern high-latitudes, we perform simulations with LPX forced by climate anomalies relative to present-day (surface temperature, total precipitation, and cloud cover) from snapshot simulations since the LGM with the coupled oceanatmosphere Hadley Centre climate model (HadCM3; Singarayer and Valdes, 2010). Climate anomaly data are available for every 1000 yr. The simulated climate anomalies are linearly interpolated to individual years between every 1000 yr time slices, and added present-day observed gridded climate (CRU; Mitchell and Jones, 2005), annual time series for these climate variables are generated. The climate simulation does not include millennial-scale variability like the Bølling-Allerød or the Younger Dryas. We keep the number of wet days per month unchanged, as in the CRU data. The LPX simulations were performed with prescribed atmospheric CO 2 (Joos and Spahni, 2008) and variable orbital forcing (Berger, 1978) since the LGM. N deposition plays a crucial role for the peatland N availability and its productivity. Since peatlands are often N limited, atmospheric N deposition is an important N source. At present, the N deposition in peatlands is in the range of 0.2-1.2 g N m −2 yr −1 in North America (Bubier et al., 2007), 0.3-0.8 g N m −2 yr −1 in eastern Canada and Sweden, but higher for central (3-4 g N m −2 yr −1 ) and western Europe (up to 8 g N m −2 yr −1 ; Turunen et al., 2004). For the simulations over the LGM and the Holocene we assume N deposition was considerably lower, and thus we prescribe a constant atmospheric deposition of NH x and NO y according to the preindustrial (1850 AD) values of Lamarque et al. (2011), which are estimated to be on average 0.10 g N m −2 yr −1 and maximum values of 0.85 g N m −2 yr −1 for northern peatlands.
Simulations are started with a spinup of 1500 yr, and an analytical solution to ensure that all C pools have established equilibrium conditions at the beginning of the LGM (21 kyr BP; 1 kyr = 1000 yr; BP = before present, and present = 1950 AD). In our main scenarios, no peatlands and peat carbon exist at the LGM and peat carbon build-up is simulated transiently. The influence of the initial conditions simulated at the end of the spin-up on results for non-peatland carbon is decreasing relatively rapidly as slowest decomposition time scales for carbon on mineral soils are typically less than a few centuries. LPX then runs forward on a daily time step until present (2000 AD), which requires the interpolation from monthly input data to quasi daily values. LPX has a resolution of 2.5 • latitude and 3.75 • longitude. For a global simulation from the LGM to present, the code is parallelized on 40 CPUs, and it takes about 4 days of computation time.

Setup for global simulations
For global simulations, we prescribe changes in polar ice sheets and land area, affected by global sea level rise, including effects of land isostatic uplift in response to ice melting (Peltier, 2004). It is assumed that at the start of the simulation, at the LGM, no peatlands existed in present-day peatland areas. These areas are either covered by the ice sheets, or climatic conditions were not conductive for peat formation. Since LPX has no bioclimatic restriction of peat initiation, peatlands could grow even under LGM conditions in some regions. We thus prescribe the time for the oldest peat initiation possible for major peatland regions as 17 kyr BP for North America and 16 kyr BP for Europe and Siberia on the basis of basal age data in MacDonald et al. (2006). During the simulation peatlands naturally develop on grid cell land area that becomes available after the retreat of northern ice sheets. As peatlands are simulated on a sub grid cell fraction, carbon and water are transferred between different pools using a scheme developed for anthropogenic land-use change (Strassmann et al., 2008;Stocker et al., 2011). With this scheme it is possible to prescribe the expanding peatland area and allow for C accumulation over non-peatland vegetation in a framework that conserves the C mass. We assume a strictly monotonic expansion of peatland area fraction over time (f (t), t) using a sigmoid function: where f present is the peatland fraction per grid cell at present and t S the time with maximum frequency of peatland initiation. We choose t S = 10 kyr BP, the time when peatland initiation peaks (MacDonald et al., 2006;Yu et al., 2010). During the simulation peatland expansion deviates from the analytical function f (t) due to changes in the land mask. The latter is updated every 1000 yr by linearly interpolating the area changes (1 ‰ grid cell area per year). For present-day northern high-latitude peatland fractions we used histel and histosol fractions from the Northern Circumpolar Soil Carbon Database (NCSCD; Tarnocai et al., 2007) with a total area of 2.71 × 10 6 km 2 (1.42 × 10 6 km 2 and 1.29 × 10 6 km 2 in North America and Eurasia, respectively). These presentday and scaled peatland maps of the past are shown in Fig. 2.
It is important to note that according to the ICE-5G ice-sheet reconstruction (Peltier, 2004) large areas of present peatlands were ice covered in North America at 10 kyr BP. Total simulated peatland area thus increases from 0 at the LGM to 2.71 × 10 6 km 2 at the present. Our calculations of peatland area in the model roughly match the cumulative curve of peat initiation dates (Fig. 3). However, there are some noticeable differences. Before the Holocene, total peatland area prescribed in LPX (Eq. 8) increased slower than the rate of increase in the number of individual peatlands, especially compared to Alaska (Fig. 3). That could be simply because many peatlands in their initial states were very small in size. Also, in the early Holocene until about 6 kyr BP, LPX peatland area increases faster than the rate of peatland initiation (MacDonald et al., 2006). Thus, the prescribed increase in peatland area in the early Holocene implies that this increase is mostly controlled by the lateral expansion of existing peatlands. This is plausible if one assumes that each initial peatland expands itself over time and climate conditions during the Holocene thermal maximum ( Fig. 4c; also see Jones and Yu, 2010;Yu et al., 2010) may favor peat area expansion. In summary, during the glacial-interglacial transition only peatlands in Alaska ( Fig. 3; Jones and Yu, 2010) and Siberia (Smith et al., 2004) contribute to circumpolar peatland area. With the retreat of the Laurentide and the Fennoscandian ice sheets in the early Holocene, the peatlands develop in these areas and C accumulation starts (MacDonald et al., 2006).

Setup for site simulations
For simulations at individual sites we do not change the peatland area over time as the analysis is based on peat-core data per unit area (m 2 ). However, we do adjust peat initiation dates of individual grid cells to match basal peat data of the corresponding peat-core site. Other than that, the setup for these simulations is identical to that of the global simulations.

Setup for future simulations
Simulations with LPX for future scenarios have been conducted by Stocker et al. (2013). A total of 31 simulations were run until 2100 AD with prescribed climate from CMIP5 model outputs (CMIP5, 2009) for scenarios RCP 2.6 and RCP 8.5 (numbers indicate radiative forcing target in W m −2 at 2100 AD; RCP, 2009;Meinshausen et al., 2011). The climate changes averaged over peatland area (2100-2010 AD) range from a warming of 0.5 • C (RCP 2.6) to 9.2 • C (RCP 8.5). Average precipitation changes in peatland area for the 31 simulations range from −2 (RCP 2.6) to +202 mm yr −1 (RCP 8.5). These LPX runs have been initialised from the end state of a transient run since the LGM (Sect. 3.2.1) at pre-industrial, and then forced with the climate output from the range of Earth System Models mentioned above until 2100 AD (Stocker et al., 2013;CMIP5, 2009).

Calibration of model parameters to site data
In this study, we use the term "apparent peat accumulation rate" for the accumulation rate of peat layers, as observed at the present, that underwent decomposition since their formation. This is in contrast to the "peat accumulation rate" or net ecosystem production (NEP) that represents the C balance of a peatland in the actual year of peat formation. In addition For the LPX maps the light blue areas denote the ICE-5G northern ice sheets, and green areas the ice-free land mass, expanding into present-day oceans due to the lower sea level in the past (Peltier, 2004). we use the term "gross peat accumulation rate" for peatland accumulation, or C balance, in the year of peat formation, excluding peatland decomposition in the catotelm since or in the year of peat formation.

Peatland apparent accumulation rate data
As a tuning target, we use reconstructed Holocene apparent peat accumulation rates from a compilation of 33 peatland sites in the northern high-latitudes (Yu et al., 2009). The data include fens and bogs as well as permafrost peatlands from the following major peatland regions: Alaska, Canada, Scotland, Finland and Siberia. Each region itself represents an average of 3 or more records with basal ages between 5.5 and 14.2 kyr BP. Measured apparent C accumulation rates range between 8.4 and 38.0 g C m −2 yr −1 for these regions over the Holocene, with an overall time-weighted average rate of 18.6 g C m −2 yr −1 for all 33 sites (Yu et al., 2009). The locations of the 33 reconstructed apparent peat accumulation rate sites fall into 24 grid cells of the LPX model grid. For the comparison with model output, accumulation rate data are binned at 1000 yr intervals over the entire Holocene. We use these data as targets to find the best combinations of model parameters.

Parameters regulating catotelm accumulation rates
The simulation of past peatland accumulation rates very much depends on the balance of C transfer from the acrotelm to the catotelm (F AC ) and loss by heterotrophic respiration from the catotelm to the atmosphere (F CR ). We do not consider any lateral transport or other processes for C change in the catotelm. On the landscape scale also erosion may become important for peatland ecosystems. Here, erosion processes are not considered. balance as described in Eq. (5). The change of the catotelm C pool (C cato ) thus can be written as where F LC is the C flux from the litter directly to the catotelm from roots (Fig. 1). In this equation F AC is the dominating term, defining the long-term peat accumulation rate and its variability over time. F LC represents in general a very small flux, and F CR scales almost linearly with C cato , analogous to Eq. (6). From Eqs. (5) and (9) one can conclude that the simulated peat accumulation rates depend on the input (F LA , and F LC ) and respiration (F AR , and F CR ) C fluxes in the acrotelm and catotelm. Parameters used to configure these fluxes are thus subject to tuning. In the following, we limit our tuning to three parameters: an additional globally uniform N input rate (N in ), and decomposition rate constants for both the acrotelm and the catotelm (k 10,acro , and k 10,cato ; see Table 1). The input C fluxes directly relates to NPP, which is explicitly calculated each day in LPX as a function of temperature, water availability, incoming solar shortwave radiation, atmospheric CO 2 concentration and the N availability (Xu-Ri and Prentice, 2008). Of these influences the N availability and in general the N cycle in peatlands is a new feature compared to a previous model version (Wania et al., 2009b). As mentioned earlier, N input is crucial for N availability and finally peatland NPP. For the past 21 kyr BP the natural N input, in particular N fixation, is highly uncertain. With the setup described so far, we might underestimate total N input to the peatland ecosystem in our simulations. Therefore, we explore the impact of an additional globally uniform N input rate (N in ) in peatlands as a tuning parameter with values of 0.1, 0.2, 0.3, or 0.4 g N m −2 yr −1 . Important tuning parameters for the acrotelm and catotelm C balance are the decomposition rate constants (k 10,acro and k 10,cato ). Beside the direct effects on the C balance, the decomposition also has implications for the N cycling as it dynamically affects the N mineralization rate from litter, acrotelm and catotelm peat, and thus the N availability ( Fig. 1; Xu-Ri and Prentice, 2008). Previous modeling studies using simple peat schemes (no simulation of plant cohorts) used decomposition rate constants for the acrotelm at 0.03 yr −1 (k 10,acro ; Wania et al., 2009b) or 0.067 yr −1 (k acro ; Kleinen et al., 2012), and for the catotelm at 0.001 yr −1 (k 10,cato ; Wania et al., 2009b) or 0.0000335 yr −1 (k cato ; Kleinen et al., 2012). We slightly varied decomposition rate constants to 0.02, 0.03, 0.04, or 0.05 yr −1 for the acrotelm and to 0.0004, 0.0007, or 0.001 yr −1 for the catotelm in our tuning exercises.

Metric for apparent accumulation rate history
We used the simulated C balance of the catotelm (Eq. 9) for model calibration with peat-core data. While the model simulates actual C accumulation for individual years in the past, the peat-core data measure C increase in peat profiles that underwent decomposition since it was first accumulated (Yu, 2011). We thus use simulated changes in catotelm gross C accumulation (mainly acrotelm-to-catotelm C transfer, F AC ) summed over 1000 yr intervals, and let this accumulated C decompose until present with the time dependent catotelm decomposition rate as calculated by LPX, averaged over the same intervals. In this way we derive an apparent accumulation rate history, which can be compared to the peatcore data in Yu et al. (2009). As the metric for evaluating the comparison we use the root-mean-square deviation (RMSD) of simulated versus reconstructed apparent C accumulation rates. The RMSD is calculated for each 1000 yr interval at each site, and then averaged over all 33 sites. The tuning of the model parameters (N in , k 10,acro and k 10,cato ) is thus optimized (with the smallest RMSD) for all 33 records and for all time periods together, and not for the present day peatland C storage (kg C m −2 ). Intervals with negative rates in LPX have been excluded, as this can not be observed by peat core data.

Results of parameter sensitivity simulations
In a total of 48 simulations that combine all parameter values of N in , k 10,acro and k 10,cato , the simulated apparent peat accumulation rates have an average RMSD to the peat-core data of 15.1 g C m −2 yr −1 (ranging from 12.7 to 19.7 g C m −2 yr −1 . This RMSD value is quite large, as it is of the same order of magnitude as the C accumulation rate itself. Uncertainties in reconstructed peat-core data, a scaling effect from the grid cell to the site levels, or local differences in modelled and real past climate contribute to large deviations. Fig. 4. Simulated gross (blue) and apparent (red) C accumulation rates for peatland sites in (a) Scotland and (b) Canada compared with reconstructed C accumulation rates (green; Yu et al., 2009;Anderson, 2002;Gorham et al., 2003). Top left values denote RMSD of simulated (red) versus reconstructed accumulation rates (green). See Appendix Fig. A1 for comparisons of simulated and reconstructed C accumulation rates from all other grid cells and sites. Panel (c) shows the annual average temperature (solid lines) and changes in annual precipitation relative to present (dashed lines) for the corresponding grid cell in Scotland (orange) and Canada (purple) in the HadCM3 model (Singarayer and Valdes, 2010).
The best agreement, i.e., with the lowest RMSD, is obtained for N in = 0.2 g N m −2 yr −1 , k 10,acro = 0.04 yr −1 and k 10,cato = 0.0004 yr −1 (RMSD = 12.7 g C m −2 yr −1 ). However, the chosen parameters are not unambiguously the best. There are 9 other combinations that have a RMSD of 14.0 g C m −2 yr −1 or less. Simulations with N in = 0.1 and 0.2 g N m −2 yr −1 have nearly an identical low RMSD, that is always lower than in simulations with N in = 0.3 and 0.4 g N m −2 yr −1 , irrespective of the other two parameters. In a similar way, the RMSD for k 10,acro = 0.05 yr −1 is always higher and for k 10,cato = 0.0004 yr −1 always lower than the RMSD for the other decomposition rate constants.
In Fig. 4 we highlight a comparison of simulated versus reconstructed apparent C accumulation rates for two sites in Scotland and one site in Canada (as compiled in Yu et al., 2009;Anderson, 2002;Gorham et al., 2003). The apparent accumulation rates at the Scottish sites, both located in the same model grid cell, show an early Holocene increase and a linear decrease towards the present. LPX captures the behavior of the Glen Carron bog, but overestimates present day accumulation rates. The Canadian site in Fig. 4b shows a local minimum in accumulation rate in the mid-Holocene, where LPX simulates no peat accumulation at all during the mid-Holocene. At this particular point in the simulation, precipitation decreases (Fig. 4c), the water table drops and acrotelm decomposition is increased. This leads to a gap in C accumulation of the peatland. Overall, LPX simulates lower rates at this particular site. Such deviations are to be expected for a simulation where parameters are optimized for all 33 sites.
Generally, simulated apparent C accumulation rates in peatlands clearly show less variability over the Holocene than reconstructed peatland data (see Fig. A1 in the Appendix for comparisons at each of these sites/grids). Partly this can be explained by an underestimation of centennial to millennial variability in the climate input data as they were interpolated from 1000 yr snapshot climate simulations. For two fens in Northern Canada the model simulates very little NPP because annual temperature is too low (below −12 • C). Plant growth is limited either directly by cold climate or indirectly by simulated N limitations due to very slow decomposition and mineralization under extremely cold climate. As a consequence no peatland C is accumulated (Fig. A1 in the Appendix). These two sites have been excluded for the RMSD calculation.
The impact of the N cycle was tested separately by running the model without the dynamic N cycle, that is, atmfrac = 70 %, N in = 0.0 g N m −2 yr −1 , but all other parameters www.clim-past.net/9/1287/2013/ Clim. Past, 9, 1287-1308, 2013 are identical. The optimal decomposition rate constants in this case are k 10,acro = 0.1 yr −1 and k 10,cato = 0.001 yr −1 as given in Table 1. For this tuning peatland NPP is larger without N limitation, thus higher decomposition rates make up for a larger C input to the acrotelm and catotelm. The resulting correlation of simulated versus reconstructed peat C accumulation rates averaged by region and over 1000 yr time intervals is shown in Fig. 5. Average apparent accumulation rates for Scotland and Finland are considerably closer to reconstructions, when simulations include the dynamic N cycle. However, it should be noted that the two model versions (LPX with and without dynamic N cycle) are tuned separately and peat accumulation behaves differently. An interpretation of N cycle dynamics from Fig. 5 is thus not straightforward. The average accumulation rates in Siberia are underestimated by LPX in both cases with or without N dynamics. The bias likely originates from a peat core near Tomsk, Russia (Yu et al., 2009, Cell 3 in Fig. A1 in the Appendix), which has very high accumulation rates that are not reproduced in LPX with the optimized parameters.  . LPX values denote the averages weighted by peatland area and the simulated range of grid cell values in brackets. For calculations, organic matter (OM) contains 50 % carbon (C).

Peatland
Prescribed or Observations parameter simulated in LPX  N input fluxes

Global simulations of peatland development
Global simulations of dynamical vegetation and peatland development from the LGM to present are performed with the climate forcing described in Sect. 3 and the parameters optimized as described in Sect. 4. Compared to simulations with previous versions of the model (e.g., Joos et al., 2004) the inclusion of the dynamic N cycle is considered as a major model change, strongly affecting the land C cycle. We thus check simulated C and N fluxes and pools against observations, before we present results from transient simulations since the LGM.

Validation of N cycle in peatlands
In Table 2 we summarize simulated N fluxes, N pools and C : N ratios in peatlands that were averaged over the period 1950-2000 AD for all PFTs and the northern peatland area. In general LPX simulates values that fall in the range of Fig. 6. Comparison of (a) gross primary production, (b) net primary production, (c) heterotrophic respiration and (d) net ecosystem production in northern high-latitude peatlands for LPX simulations with and without dynamic N cycle (DyN). Each point represents a 10 yr average. Solid lines are linear fits through points (R 2 , coefficient of determination, is given in the headers) and dashed lines represent the 1 : 1 line.
observations . Average N input and output in the peatland ecosystem are smaller in the simulations than in observations. The difference can be explained by a low N deposition in LPX, which uses pre-industrial (1850 AD) N deposition value, whereas observations are based on present-day fluxes with much higher N deposition rates .
The weighted averages of C : N ratios for global peatland area and fractional plant cover in LPX are simulated as 44.7, 41.8, 51.0 and 34.0 for vegetation, aboveground litter, belowground litter and soil, respectively (corresponding N content is given in Table 2). Since mosses have no roots, below ground litter in the model is calculated entirely from root turnover of graminoids. The higher C : N ratio of below ground compared to above ground litter is inherited from the higher C : N ratio of plant production for graminoids compared to mosses. Simulated peat (both acrotelm and catotelm) C : N ratios are in the range of 20-42 (11.9-20.0 mg N g −1 dry biomass), which is at the lower end of observed C : N ratios of 30-55 .
How do the N cycle processes affect plant growth and ecosystem C balance in peatlands? In the LPX version with a dynamic N cycle photosynthesis and GPP are not very different to the version without a dynamic N cycle (Fig. 6), but NPP is reduced considerably if not enough nitrogen is available to support plant growth. This additional limitation affects the overall growth of plants and thus all subsequent C fluxes including heterotrophic respiration (HR) and NEP. It is hypothesised that in northern high-latitudes N mineralization is slow and N availability is low in spite of high SOM content (Chapin et al., 1995;Nadelhoffer et al., 1991). This effect can be seen in Fig. 6; some sites have no growth (in either GPP or NPP) in LPX simulations with dynamic N cycle, and thus also no HR or net ecosystem production (NEP). NEP reflects the balance between NPP and HR (NEP = NPP − HR). Thus, NEP depends on how much carbon is transferred to the catotelm where it is subject to very low respiration rates. An important difference between LPJ-WHy and the current LPX version is that the acrotelm-tocatotelm flux is set constant in LPJ-WHy, whereas it is computed from the acrotelm C balance in LPX. A lower NPP in LPX than in LPJ-WHy does therefore not necessarily translate into a lower NEP. Comparisons of NEP with different driving variables are shown in Fig. A2 in the Appendix.
The reduction of maximum simulated NPP from 800 to 450 g C m −2 yr −1 due to the inclusion of a dynamic N cycle is an important aspect for peat growth and crucial for the C balance in the acrotelm and peat accumulation rates. In simulations without dynamic N cycle the resulting NPP is clearly overestimated as discussed in detail for LPJ-WHy by Wania et al. (2009b). They find that in LPJ-WHy the excessive NPP values lead to an overestimation of the seasonal amplitude of NEP compared to data. In their analysis Wania et al. (2009b) hypothesise that this is due to the lack of nitrogen or phosphorus limitation of NPP. However, their annual NEP values had no offset to measured data. An explanation for this is given by the actual state of peatland soil C pools. In simulations without dynamic N cycle peatland soil C pools are much closer to equilibrium than in simulations with a dynamic N cycle, due to the higher NPP over many thousands of years. According to Eq. (6) for the acrotelm and analogously for other C pools, the larger peatland soil C pools lead to higher HR. This means that NEP in both simulations can be in the same range, despite NPP being considerably lower in simulations with dynamic N cycle. Furthermore, because peatlands in these two simulations have a different accumulation rate history since the LGM, their peatland NEP at a chosen location are not necessarily correlated. This is shown by the weak correlation in Fig. 6d.

Variability of acrotelm-to-catotelm carbon transfer
The global area-weighted average acrotelm-to-catotelm C transfer (F AC ), including positive and negative rates in LPX, is 22.3 g C m −2 yr −1 at present. This average value is substantially larger than the constant flux of 12 g C m −2 yr −1 assumed in LPJ-WHy (Wania et al., 2009b). The acrotelmto-catotelm C transfer shows a large spatio-temporal variability. A frequency distribution of F AC (10 yr averages of individual grid cells for the peatland land unit) over the last 10 kyr shows a wide range, even spreading out to negative rates (Fig. 7).
The way we implemented F AC (Eq. 5) very often leads to slightly negative values whenever acrotelm respiration is larger than C input to acrotelm on an annual basis. So, for balanced or decaying peatlands, F AC oscillates around 0 g C m −2 yr −1 with a higher probability for negative rates. Thus, there are lots of grid cells and time periods that fall within the same bin (−1 < 0 g C m −2 yr −1 ) that shows up as a frequency peak in Fig. 7 (a total of 72 187 occurrences). The positive acrotelm-to-catotelm rates can be approximated with a normal distribution centered at 20 g C m −2 yr −1 and a standard deviation of 7 g C m −2 yr −1 . However, the tails of the simulated distribution deviate significantly from the normal distribution. Some grid cells show transfer rates higher than 50 g C m −2 yr −1 . This is quite a remarkable high transfer rate, as it represents a grid cell average. It has been shown that peatland C accumulation rates, and thus also acrotelmto-catotelm transfer rates, for individual sites indeed can exceed 50 g C m −2 yr −1 (Yu et al., 2009).
Although acrotelm decomposition (F AR , Eq. 6) is now dependent on water table depth (3.7 cm on average) and plant available soil moisture content (0.41 on average), the inferred (Eqs. 2-4) average respiration modifier (R water = 0.47) is only slightly higher than in LPJ-WHy (R water = 0.35; Wania et al., 2009b). However, the variability in R water increases the variability in F AR , which has also a direct impact on the distribution of F AC . Fig. 7. Frequency histogram of acrotelm-to-catotelm transfer rates (10 yr averages) from simulated grid cells in LPX over the Holocene (last 10 kyr and for bins with a width of 1 g C m −2 yr −1 ). Negative values on the x-axis denote a transfer of carbon from the catotelm to acrotelm, implying C loss. Rates below zero (−1 < 0 g C m −2 yr −1 ) occur most frequent and are outside the frequency scale as indicated. A normal distribution (µ, σ ) is shown as an approximative distribution of positive rates (red line).

Comparison of modeled results with observed soil C data
At present, the peatland area in LPX is prescribed to 2.71 × 10 6 km 2 according to the NCSCD (see Sect. 3.2.1).
Other studies use much larger present-day peatland areas of ∼ 4 × 10 6 km 2 in total , whereof mainly Eurasian peatland area is larger than in our study. Therefore, in an initial test, we compare simulated soil C densities with the NCSCD data (Tarnocai et al., 2007) for North America. The NCSCD has detailed C inventories for North America at soil depths of 30 cm, 100 cm, and total soil C (down to 300 cm). Although LPX only simulates soil hydrology and thermal dynamics of the top 200 cm, there is no explicit limit of simulated C in soil pools. We thus compare LPX soil C with total soil C in the NCSCD. For this, the high-resolution soil C data were regridded to the LPX grid using the fractional peatland area; deviation in total soil C stocks between the regridded and original data set are small (Fig. 8).
The simulated total C inventories in peat soils, in permafrost soils, and all soil classes are well within the ±10 % uncertainty range of the observation-based estimates (Fig. 8). Simulated soil C in North America is 332 Pg C for all soils (observation based: 344 Pg C), 177 Pg C for peatlands (167 Pg C), 155 Pg C for mineral soils (177 Pg C), whereof 105 Pg C are associated with permafrost soils (115 Pg C; Fig. 8. Present-day distribution of peatland and all land (peatlands, permafrost and mineral soils) soil C densities in North America: (a) peatland C density from NCSCD data (Tarnocai et al., 2007), (b) peatland C density from NCSCD data averaged on LPX grids, (c) peatland C density simulated by LPX, (d) all soils C density from NCSCD data, (e) all soils C density from NCSCD data averaged on LPX grids, (f) all soils C density simulated by LPX. Numbers denote the area-weighted total C pools. Tarnocai et al., 2009). The simulated spatial distribution of peatland C and soil C densities is also in broad agreement with observations (Fig. 8). The observed band of high C densities stretching northwest from the US lake region to the Arctic Ocean is well represented in the model. However, compared to the NCSCD, the LPX peatland and total soil C stocks are overestimated in Alaska, and LPX total soil C is underestimated in northern Canada. A possible reason for the latter may be a misrepresentation of temperature effects on the dynamic N cycle in high northern latitudes. Simulations without dynamic N cycle show slightly higher permafrost soil C in northern Canada (+10 kg C m −2 ). The LPX peatland soil C overestimation in Alaska may possibly be related to uncertainties in the prescribed climate evolution or in peat initiation. An alternative explanation could be that effects of litter quality on decomposition, not included in LPX, are important. Also, the flooding disturbance might have been important in the mid-and late Holocene at some regions in Alaska that was not explicitly simulated in our model. A regression of non-zero values between peatland soil C data and modelled soil C densities on peatland yields a coefficient of determination of R 2 = 0.54.
Total present-day mineral and peat soil C and N pools simulated in LPX are 1154 Pg C and 63 Pg N in the Northern Hemisphere (30-90 • ), and 1714 Pg C and 108 Pg N globally. Thereof 365 Pg C and 11 Pg N are located in peatlands in the Northern Hemisphere using the peatland distribution by Tarnocai et al. (2009). Of the global soil C pool, 608 Pg C and 29 Pg N are located in northern circumpolar permafrost area with peatlands therein contributing 275 Pg C and 8 Pg N. Albeit total permafrost C is slightly lower than in Tarnocai et al. (2009), it is important to note that peatlands within permafrost area (histels) are a significant part of the soil C pool, and may behave differently under climate change than mineral soil permafrost. Further, the model does not simulate soil C stored in Yedoma or deltaic deposits, which are estimated to be 407 and 241 Pg C, respectively (Tarnocai et al., 2009), but it does take into account soil C that is preserved in continental shelves due to rising sea levels during the glacialinterglacial transition. However, the latter does not contribute significantly at present as we assume that this C pool mostly respired over the Holocene.  Figure 9 shows simulated changes in average water table depth, permafrost thaw depth, average NEP (net ecosystem production), and in peatland soil C (both the acrotelm and catotelm) since the LGM. The water table rises in parallel to increased precipitation of +240 mm yr −1 over peatland areas and since the LGM. The rise in water table position probably lead to a long-term decrease in acrotelm respiration and thus an increase in acrotelm-to-catotelm C transfer rate. Thaw depth increases due to an average change in annual mean peatland temperature from −16 to −4 • C from the LGM to present. The highest rate of increase in thaw depth at 12-7 kyr BP may be related to the Holocene thermal maximum (HTM) in high northern latitudes (Kaufman et al., 2004), where abundant peatlands exist. The slight decrease in thaw depth over the last 7 kyr, especially in North America, may be caused by the late Holocene climate cooling, after the HTM. Both changes in water table and thaw depth increase the saturation of liquid water in the soil and the amount of water available to plants; this tends to increase NPP and NEP.

Results for peatland NEP and C pools
A major observation from peat-core data is high apparent peat accumulation rates at the beginning of the Holocene and a decreasing trend towards the present . Yu et al. suggest that in the early Holocene high summer insolation and strong summer-winter climate seasonality lead to a longer and more intensive growing season, and thus higher annual NPP and finally peat C accumulation. This feature is not present in all peat cores (Fig. A1 in the Apppendix), and the robustness of this conclusion may be affected by limited data availability as only data from 33 sites are currently available.
For comparison, we examined the evolution of the area weighted northern high-latitude average peatland NEP (Fig. 9b). NEP in g C m −2 yr −1 represents the average peat C accumulation rate per unit area and should not be confused with total peat accumulation, which is the product of average NEP and peat area. We recall that the Bølling-Allerød and Younger Dryas climate swing are not present in the climate data of the Hadley Centre model used to force LPX in this study (see Fig. 4c). This likely affects the simulated NEP during the transition and the early Holocene in Europe and Siberia. In Europe and Siberia, NEP shows a broad maximum between 13 and 10 kyr BP, followed by a millennial-scale decrease during the Holocene. In North America, on average NEP increases over the transition and the early Holocene to peak around 8 kyr BP, after the major retreat of the Laurentide Ice sheet. In brief, our results are compatible with the suggestion of an early Holocene peak in peat NEP, given the uncertainties in our climate forcing data.
The prescribed peatland area (T09; Tarnocai et al., 2009) is increasing over time (Figs. 2 and 3) and this translates into higher overall C uptake by peatlands. Total peatland C pools start to increase significantly at about 12 kyr BP in Europe and Siberia, and at about 10 kyr BP in North America (Fig. 9c). Rates of C pool increase are linear, but show different slopes in the two regions (Fig. 9d). These results suggest that northern peatlands have been a persistent land C sink since the early Holocene.

Alternative peatland development scenarios
Total northern peatland C as simulated by LPX for the peatland area T09 (Tarnocai et al., 2009) is 365 Pg C at the present, which is within the range of previously published estimates (Tarnocai et al., 2009;Yu et al., 2010). However, since the total amount changes with the area and time since peatland initiation, we simulated peatland development for two additional scenarios. First, assuming that some peatlands are much older than what peat basal dates have suggested, we carried out a test simulation with no restriction to peat initiation (T09 LGM ). Under these conditions, LPX simulates peatland C accumulation already in the LGM. As a consequence, catotelm C pools are much larger and closer to the steady state at present. Second, we simulated peatland development using a total northern (30-90 • N) peatland area of 4 × 10 6 km 2 as in Yu et al. (2010, Y10). The distribution of northern high-latitude peatlands affects the area weighted mean of peatland soil C, NEP, and peatland C uptake rate (Table 3). While present-day peatland NEP averages to similar values for all scenarios, peatland C amount and uptake rates at present and in the past are considerably different from each other (Fig. 9). Below we analyse the simulation results for the three scenarios T09, T09 LGM and Y10.
In scenario T09 LGM present-day peatland soil C is nearly double the value of T09 (Table 3), despite the fact that the increase from the LGM to present is only marginally larger for the northern total and is almost identical for North America. In LPX, additional peatland soil C in T09 LGM is located in the West Siberian Lowlands (WSL) during the LGM. For this specific region, the T09 LGM scenario simulates 400 Pg C at present and thus clearly overestimates GIS-based observations of 70 Pg C . Scenarios T09 and Y10 with peat initiation after the LGM (earliest at 16 kyr BP) have much lower present-day peat C inventories of 97 and 90 Pg C for West Siberia, respectively. From the T09 LGM simulation it becomes clear that even though present day peatland NEP is not that different from that in T09, the peatland initiation and history are very relevant for total C storage in these ecosystems.
The Y10 scenario suggests northern total peatland soil C of 550 Pg C, which approximately proportionally related to the larger area compared to T09 (Table 3). Surprisingly, this value is very close to the mean value (547 Pg C) independently estimated by Yu et al. (2010), despite the fact that LPX results have not been tuned to the total C stock of northern peatlands. It is also noted that the tuning of model parameters to accumulation rates for only a few sites has a significant impact on the total amount of simulated peatland soil C. Also, the inclusion of the dynamic N limitation has a large effect on the amount of peatland soil C, as simulations without dynamic N limitation yield values of more than 720 Pg C for the Y10 scenario. Irrespective of the scenarios T09 and Y10, peatlands in North America and Europe/Siberia each contribute about 50 % to the northern peatland total. In the Y10 scenario, North American peatland soil C is about ∼ 100 Pg C larger than for T09, which is not supported by the NCSCD data that show even lower soil C than in T09 (Fig. 8). However, the peatland area of scenario Y10 used in LPX simulations is also proportionally larger in North America. One can thus conclude that differences in the geographical distribution of peatland area between scenarios T09 and Y10 are probably less important than the difference in total peatland area for simulated peatland soil C.  Tarnocai et al. (2009) and Y10 by Yu et al. (2010). T09 LGM has no restriction on peat initiation and peatland areas during the LGM.

Transient simulation results for future scenarios
The model was calibrated on the millennial timescales resolved by the peat accumulation data. It appears justified to apply the model also on the decadal-to-centennial time scales that are typical for anthropogenic climate and environmental changes as relevant processes such as NEP, water table depth variations, or acro-to-catotelm carbon transfer or decomposition rates are allowed to vary rapidly within the model. Further, projected changes in climate and atmospheric composition over this century are comparable in magnitude to the changes that peatland experienced over the past 20 kyr. LPX simulations were performed for future scenarios (Sect. 3.2.3). Results for peatland NEP and C pool changes from 2010 to 2100 AD are shown in Fig. 10. Overall, changes in peatland C stocks over the 21st century are very small and less than 6 Pg C, irrespective of the forcing scenario applied. Peatland NEP is still positive today and C pools are linearly increasing until mid 21st century in all runs. In RCP 2.6, peatland C stocks steadily increase towards 2100 AD. After 2050 AD, the median of simulations under RCP 8.5 show a decrease in average peatland NEP, becoming negative by 2100 AD. Negative peatland NEP represents Fig. 10. Results of 31 LPX simulations forced by CMIP5 model outputs using scenarios RCP 2.6 (blue) and RCP 8.5 (red) for the period 2010-2100 AD as described in Stocker et al. (2013). Plots show that (a) peatland water table depth and thaw depth increase over time, (b) peatland NEP and (c) soil C start to decrease around 2050 AD in some simulations using RCP 8.5. Plots show the 5 yr running average of model output.
an actual loss of peat C as shown by the decrease in peatland C stock of up to ∼ 1 Pg C relative to 2050 AD. There is a large spread in results among the CMIP5 models and thus in LPX when these different forcings are applied. The HadGEM2 model (CMIP5, 2009) simulates large climatic changes in the northern high-latitudes and its forcing reduces peatland C pools after 2050 AD in LPX, while the CCSM4 model (CMIP5, 2009) yields comparably smaller climatic changes translating into a positive C balance also under RCP 8.5. The average peatland water table simulated in LPX drops by a few centimeters. In all LPX simulations the average peatland thaw depth is increasing consistently, between 7 and 74 cm until 2100 AD. The increase in thaw depth also represents an increase in active layer depth, which accelerates microbial induced SOM decomposition in thawed peat.
Despite the large warming and the increase in thaw depth, no simulations show a net C loss in peatland C pools in 2100 AD compared to 2010 AD (Fig. 10b). LPX simulation results thus suggest that northern high-latitude peatlands as a whole are quite resistant to C loss under the prescribed climate change alone until 2100 AD, which is in contradiction to earlier estimates (Davidson and Janssens, 2006;Ise et al., 2008;Dorrepaal et al., 2009). All RCP 8.5 simulations with decreasing peatland soil C after 2050 AD yet show increased peatland NPP at 2100 AD that partly compensates the increased soil decomposition rates.

Discussion
Calibrating model parameters with reconstructed peat accumulation rates is an important step towards robust simulations of peatland C dynamics. As common to large-scale environmental modeling, the upscaling from site data to large regions and to the global scale is problematic. It is thus not surprising that LPX reproduces C fluxes for regional averages better (Fig. 5) than at individual sites (Fig. A1). Also, the inclusion of the dynamic N cycle process resulted in a better agreement of simulated NPP and apparent peatland C accumulation rates from the observations. On the other hand, the simulated C and N cycles also strongly depend on the accuracy of GCM simulated climate used as the input data for peatland simulations. Any peat accumulation rate parameterization is tied directly to input climate and must eventually be recalibrated for deviations on centennial to millennial scales. For the simulations in this study, we used average millennial climate output from simulations with the coupled ocean-atmosphere Hadley-Centre climate model (HadCM3; Singarayer and Valdes, 2010) that does not fully reflect millennial-scale climatic features such as the the Bølling-Allerød (BA) warming period or the Younger Dryas (YD) cold event (Severinghaus et al., 1998). It is expected that the BA warm period likely favored peatland development, while the cold temperatures (decrease of ∼ 10 • C) during the YD most probably delayed and slowed down peat growth. While the BA and YD events stand out as northern hemispheric synchronious climate events there was probably also more disconnected regional millennial scale variability during the Holocene. Such variability probably did not affect the catotelm C pool directly as the effect on decomposition of deep peat is small. Significant C pool changes are possible through the mechanism of acrotelm-to-catotelm C transfer as discussed in Sect. 2.2.3. However, generally increased peatland NPP under warmer temperatures and associated higher litter inputs compensate effects of accelerated acrotelm decomposition under warm conditions.
From the LGM to the pre-industrial period, LPX simulates a net increase in total northern peatland and global mineral soil C of 440 Pg C and global vegetation C of 100 Pg C. This is within the uncertainty of recent estimates (Ciais et al., 2012) and much less than in simulations with a previous version of the model (Joos et al., 2004). The attribution of the increase considerably changed since then as the model has undergone significant changes, such as the addition of permafrost, peatland and dynamic N cycle modules.
The amount and timing of early peatland C uptake is important for the explanation of atmospheric CO 2 evolution and the global C budget during the glacial-interglacial transition (Yu, 2011;Elsig et al., 2009;Menviel and Joos, 2012). Our results suggest a C uptake by currently existing peatlands (from T09 and Y10 areas) of 175-272 Pg C between 11 and 5 kyr BP, and an uptake of 175-253 Pg C during the past 5 kyr (Fig. 9c). These carbon fluxes are substantial. For comparison, Elsig et al. infer from their ice core CO 2 and carbon isotope data a total terrestrial uptake of 290 ± 36 Pg C over the period from 11 to 5 kyr BP and a small release of 36 ± 37 Pg C thereafter. Taken at face value, the differences imply an additional uptake of roughly 20-120 Pg C in the early Holocene and a release of roughly 210-290 Pg C after 5 kyr BP due to other terrestrial processes such as changes in the Monsoon system and the greening of the Sahara (Indermühle et al., 1999;Schurgers et al., 2006), land use emissions (Stocker et al., 2011), early Holocene forest regrowth, C release linked to continental shelve flooding due to sea-level rise, or the loss of carbon from former peatlands.
However, our understanding of land C stock changes on peatlands since the LGM and over the Holocene is still incomplete. Peatland data and our model results imply that most of present-day peatlands did not exist or had small areal extents during the LGM. Our simulations do not include other peatlands that may have existed during the LGM or during the course of the Holocene, perhaps at lower latitudes as indicated by paleo records of Sphagnum (Halsey et al., 2000), but then disappeared over time due to changes in climates. By omitting "lost peatlands" we likely overestimate the northern net peatland C sink, as do all estimates based on peat-core data. This omitted C loss to the atmosphere has implications for the terrestrial C balance over the course of the last glacial-interglacial transition. One could argue that peatland disappearance at one site and appearance at another site had a balancing effect on C accumulation on peatlands. In fact, some peatlands have basal ages of 2-3 kyr BP, whereas other peat-core records are missing peat for the last 4 kyr (e.g., Peteet et al., 1998), suggesting a shift in peatland regions over time. Also, the LPX simulations suggest that at some sites, peatland C accumulation rate were negative during the Holocene. This cannot be confirmed directly, but is possible within the uncertainty of apparent C accumulation rate measurements.
Assuming that simulated peatland C accumulation is well represented for today, LPX simulates present-day accumulation rates of 35 to 50 Pg C per 1000 yr (Table 3), which is in the range of previous estimates (Yu, 2011). This represents a long-term persistent C sink in the land biosphere, previously under-represented in global carbon cycle models. Simulations of LPX for future climate change scenarios suggest that peatlands remain a global C sink at least until 2050 AD (Fig. 10). In these LPX simulations, anthropogenic N deposition slightly amplifies peatland NPP, but it is not taken into account that possibly heterotrophic respiration is also directly increased by N deposition (Bragazza et al., 2006). In addition, the effects of frequent droughts on microbial growth and C decomposition in peatlands (Fenner and Freeman, 2011) are not included in the LPX simulations. This might lead to enhanced peatland C loss by microbial decomposition until 2100 AD. The difference of simulated soil C in 2100 AD between the two RCPs is only ∼ 1.5 Pg C on average, and thus negligible for the global C budget. However, the important point here is that peatland ecosystems could completely change sign from a net C sink to a C source. A previous study on peatland model intercomparisons suggests that the bioclimatic space associated with peat could decline considerably under a changing climate (Clark et al., 2010). The importance of high-latitude processes is further underlined by Stocker et al. (2013), who show that annual CH 4 emissions from peatlands increase considerably (+120-200 %) until 2100 AD in the case of RCP 8.5. Therefore, present-day peatlands and related emissions of GHG will undergo substantial changes, if anthropogenic C emissions are not reduced.

Conclusions
Present-day peatlands represent a substantial land C pool that developed since the last glacial-interglacial transition. Simulations of peatland C accumulation with the LPX model have been calibrated to reconstructed apparent C accumulation rates (Yu et al., 2009) and compared to present-day North America soil C inventories (Tarnocai et al., 2009). Beside the representation of C fluxes, LPX simulations include a fully dynamic N cycle, which considerably improved the agreement with observations and reconstructions. Results of simulations from the LGM to present show that ∼ 50 % of present-day peatland C established before 5 kyr BP. Present-day northern total peatland C is simulated as 365-550 Pg C depending on the peatland distribution and total area (Tarnocai et al., 2009;Yu et al., 2010), partly located in permafrost-affected regions. Simulations for future projections show that northern high-latitude peatlands remain a net C sink at least until 2050 AD, then becomes a small net C source depending on the prescribed RCP and CMIP5 model. From our analysis we conclude that improving the spatial coverage of peatland C accumulation rate and inventory data can greatly narrow down the uncertainty in modeling the evolution of this important land C pool in the past, and improve projections of its response to future climate change. In addition, estimates of the area and of peat C stocks in peatlands that disappeared in the past are missing. The neglection of these peatlands limits the assessment of the role of peatland dynamics in the context of C cycle changes since the LGM. Figure A1 shows the comparison of simulated versus reconstructed C accumulation rates for all peat core sites (Yu et al., 2009). Figure A2 shows NEP sensitivity against chosen peatland variables in LPX simulations.  A1. LPX simulations of gross (blue) and apparent (red) peat accumulation rates for grid cells (1-24), where peat-core data of apparent C accumulation rates (green) are available. Multiple sites may be located within the same model grid cell. RMSD for the 33 individual sites (D and site no. as in Table 1 of Yu et al., 2009) are given in g C m 2 yr −1 .

Fig. A2.
Characteristics of simulated peatland NEP as a function of (a) acrotelm-to-catotelm transfer (AtoC), (b) NPP, (c) ratio of precipitation to equilibrium evapotranspiration (P/EET), and (d) water table position (WTP). Colored points are values for the period 1950-2000 AD, and black dots from the 50 yr period centered at 10 kyr BP. The strong correlation of NEP (peat C accumulation) with NPP indicates that C accumulation is driven by plant production as suggested by Charman et al. (2013). That study also highlights that peatland ecosystems have a precipitation to equilibrium evapotranspiration ratio larger than one as shown in (c). According to LPX simulations, NEP seems to have a relative maximum at a water table of ∼ 10 cm below peat surface.