Modelling past, present and future peatland carbon accumulation across the pan-Arctic region

Correspondence to: N. Chaudhary (nitin.chj@gmail.com) 10 Abstract. Most northern peatlands developed during the Holocene, sequestering large amounts of carbon in terrestrial ecosystems. However, recent syntheses have highlighted the gaps in our understanding of peatland carbon accumulation. Assessments of the long-term carbon accumulation rate and possible warming driven changes in these accumulation rates can therefore benefit from processbased modelling studies. We employed an individualand patchbased dynamic global ecosystem 15 model with dynamic peatland and permafrost functionalities and vegetation dynamics to quantify longterm carbon accumulation rates and to assess the effects of historical and projected climate change on peatland carbon balances across the pan-Arctic region. Our results are broadly consistent with published regional and global carbon accumulation estimates. A majority of modelled peatland sites in Scandinavia, Europe, Russia and central and eastern Canada change from carbon sinks through the 20 Holocene to potential carbon sources in the coming century. In contrast, the carbon sink capacity of modelled sites in Siberia, Far East Russia, Alaska and western and northern Canada was predicted to increase in the coming century. The greatest changes were evident in eastern Siberia, northwest Canada and in Alaska, where peat production, from being hampered by permafrost and low productivity due the cold climate in these regions in the past, was simulated to increase greatly due to warming, wetter 25 climate and higher CO2 levels by the year 2100. In contrast, our model predicts that sites that are expected to experience reduced precipitation rates and are currently permafrost free will lose more carbon in the future.


Introduction
The majority of the northern peatlands developed during the Holocene ca.8-12 thousand years (kyr) ago after the deglaciation of the circum-Arctic region (MacDonald et al., 2006).The availability of new land surfaces owing to ice retreat (Dyke et al., 2004;Gorham et al., 2007), climate warming following deglaciation (Kaufman et al., 2004), increased summer insolation (Berger and Loutr, 2003), more pronounced seasonality (Yu et al., 2009), greenhouse gas emissions (MacDonald et al., 2006) and elevated moisture conditions (Wolfe et al., 2000) are some of the factors that promoted the rapid expansion of the northern peatlands.Moderate plant productivity together with depressed decomposition due to saturated conditions led to a surplus of carbon (C) input relative to output, resulting in the accumulation of peat (Clymo, 1991).Peatlands of the Northern Hemisphere are estimated to have sequestered approximately 350-500 PgC during the Holocene (Gorham, 1991;Yu, 2012).
Peatlands share many characteristics with upland mineral soils and non-peat wetland ecosystems.However, they constitute a unique ecosystem type with many special characteristics, such as a shallow water table depth, C-rich soils, a unique vegetation cover dominated by bryophytes (hereinafter referred to as "mosses"), spatial heterogeneity, anaerobic biogeochemistry and permafrost in many regions.Due to their high C density and the sensitivity of their C exchange with the atmosphere to temperature changes, these systems are an important component in the global C cycle and the coupled Earth system (MacDonald et al., 2006).Lately, considerable effort has been made to incorporate peatland accumulation processes into models with the purpose of understanding the role of peatlands in sequestering C, thereby lowering the radiative forcing of past climates (Frolking and Published by Copernicus Publications on behalf of the European Geosciences Union.Roulet, 2007;Wania et al., 2009a;Frolking et al., 2010;Kleinen et al., 2012;Tang et al., 2015) and how they might affect future climate warming and C cycling (Ise et al., 2008;Swindles et al., 2015).Clymo (1984) developed a simple one-dimensional peat accumulation model and described the main processes and mechanisms involved in peat growth and its development.This model became the starting point for later work in many peat growth modelling studies.Hilbert et al. (2000) developed a theoretical peat growth model with an annual step, modelling the interaction between peat accumulation and water table depth using two coupled non-linear differential equations.Using a similar approach, Frolking et al. (2010) developed a complex Holocene peat model by combining the dynamic peat accumulation model of Hilbert et al. (2000) with a peat decomposition model (Frolking et al., 2001).They showed that the model performed fairly well in simulating the long-term peat accumulation, vegetation and hydrological dynamics of a temperate ombrotrophic bog in Ontario, Canada.Though the models mentioned above are detailed enough to capture the peat accumulation and decomposition processes quite robustly, they lack soil freezingthawing processes, and this limits their application over regions where such processes occur.Wania et al. (2009a) were the first to account for peat dynamics in a model for largearea application by incorporating peatland functionality in the LPJ-DGVM model, which was designed for regional and global simulation of ecosystem responses to climate change (Sitch et al., 2003).Their approach included a number of novel features, such as a detailed soil freezing-thawing scheme, peatland-specific plant functional types (PFTs) and a vegetation inundation stress scheme, but it employed a twolayer representation of the peat profile, which is not as detailed as the process-based dynamic multilayer approaches taken by Bauer et al. (2004), Heinemeyer et al. (2010) and Frolking et al. (2010).Other model representations have also included peatland processes in their frameworks (Morris et al., 2012;Alexandrov et al., 2016;Wu et al., 2016) and been shown to perform reasonably at different sites.In addition, some of these models have been applied over large areas (Kleinen et al., 2012;Schuldt et al., 2013;Stocker et al., 2014;Alexandrov et al., 2016) to simulate regional peatland dynamics.
Though much information is available about the past and present rates of C accumulation in the literature, recent syntheses have highlighted the existing spatial gaps in data availability across the pan-Arctic (45-75 • N) region (Yu et al., 2009;Loisel et al., 2014).The extent and remoteness of many locations present challenges for the reliable estimation of total C, basal ages and accumulation rates of peat carbon.This demands the use of process-based modelling for upscaling and interpolation across the pan-Arctic distribution area.We employed LPJ-GUESS, an individual-and patch-based dynamic ecosystem model (Smith et al., 2001) extended to represent the characteristic vegetation, biogeochemical and hydrological dynamics of high-latitude peatlands to simulate C accumulation of peatlands across the pan-Arctic region under past, present and future climates (Chaudhary et al., 2017a).The model accounts for the close intercoupling between peatland and permafrost dynamics, which is critical for the evolution of these ecosystems and their carbon dynamics in the warming regional climate.We assess the potential effects of historical and projected climate and atmospheric CO 2 on peatland C balances and permafrost distribution at the regional scale across the pan-Arctic region.

Model description
LPJ-GUESS (Lund-Potsdam-Jena General Ecosystem Simulator) is a process-based model of vegetation dynamics, plant physiology and the biogeochemistry of terrestrial ecosystems.It simulates vegetation structure, composition and dynamics in response to changing climate and soil conditions based on an individual-and patch-based representation of the vegetation and ecosystems of each simulated grid cell and is optimized for regional and global applications (Smith et al., 2001(Smith et al., , 2003;;Miller and Smith, 2012).The model has been evaluated in comparison to independent datasets and other models in numerous studies; see e.g.McGuire et al. (2012), Piao et al. (2013), Smith et al. (2014) and Ekici et al. (2015).
We employed a customized Arctic version of the model (Miller and Smith, 2012) that has been developed to include dynamic, multilayer peat accumulation functionality and permafrost dynamics.The model represents the major physical and biogeochemical processes in upland and wetland arctic ecosystems, including an expanded set of plant functional types (PFTs) specific to these areas (McGuire et al., 2012;Miller and Smith, 2012).The revised model is described in outline below, while a full description can be found in Chaudhary et al. (2017a).In our approach, vegetation and peatland C dynamics are simulated on multiple connected patches to account for the functional and spatial heterogeneity in peatlands.The simulated PFTs have varied structural and functional characteristics and can establish in each connected patch and compete for soil resources, space and light.The composition in terms of relative PFT abundance and the physical structure of the plant community are emergent outcomes of this competition.The model is initialized with a random surface comprised of 10 patches of uneven height.Heterogeneity in the height of adjacent patches is a precondition for hydrological redistribution between them, which mediates vegetation succession and affects the peat accumulation rate, as described below.The soil-peat column is represented by four different vertically resolved layers.A dynamic single snow layer overlays the peat column, represented by a dynamic litter-peat layer consisting of a number of sublayers, updated yearly, that depends on thickness.Un-derneath the peat column is a fixed 2 m deep mineral soil column consisting of 0.1 m thick sublayers, which is underlain by a 48 m deep "padding" column consisting of relatively thicker sublayers.The soil temperature is updated daily for each sublayer at different depths, enabling the simulation of a dynamic soil thermal profile as a basis for the representation of permafrost (Wania et al., 2009a).The fractions of ice and water as well as the mineral and peat fractions in each layer govern the heat capacities and thermal conductivities and affect the freezing and thawing processes of soil water in peat and mineral soil layers (Wania et al., 2009a).The fractions of water and ice in the sublayers are updated each day depending upon variation in soil temperature and fractional mineral content, following Hillel (1998).A detailed description of the permafrost and soil temperature scheme is available in Chaudhary et al. (2017a), Miller and Smith (2012) and the references therein.
A water bucket scheme was used to simulate peatland hydrology in which the assumption is made that precipitation (rain) and snowmelt are the main input of water.Evapotranspiration, drainage, surface and base runoff are the major water balance processes in the peat layers (Gerten et al., 2004).The model also includes lateral flow of water between patches, an important governing process of vegetation and C dynamics of peatlands that is lacking in most peatland models (Chaudhary et al., 2017b).A simple lateral flow scheme connects higher elevated patches (hummocks) to lower depressions (hollows).The water table position (WTP) of individual patches is reset to the mean landscape WTP on each daily time step, effecting lateral flow from patches with a higher WTP following the current day's rainfall, snowmelt and evapotranspiration fluxes to those with lower WTP.This in turn affects the plant productivity and decomposition rates in each patch and results in dynamic surface conditions over time.
Five PFTs are used to represent the main functional elements of peatland vegetation: graminoids (Gr), mosses (M), high summergreen shrubs (HSS), low summergreen shrubs (LSS) and low evergreen shrubs (LSE).PFTs differ in the physiological, morphological and life history characteristics that govern their interactions and responses to climate and an evolving system state.Key PFT parameters in the present study include C allocation, phenology, rooting depth, tolerance for waterlogging and decomposability of PFT-derived litter (Miller and Smith, 2012).Prescribed bioclimatic limits (Miller and Smith, 2012) and favoured annual average WTP (aWTP) ranges determine PFT presence or absence (see Table A1 in the Appendix) and reflect their bioclimatic distribution.Shrubs are favoured in dry conditions (Malmer et al., 2005) where aWTP is below −25 cm (we use a sign convention in which a negative value of WTP signifies a water table below the peat surface).Conversely, mosses and graminoids are more vulnerable to dry conditions.Graminoids favour saturated conditions and establish when aWTP is above −10 cm, while mosses establish when the aWTP is between +5 and −50 cm.The establishment function is implemented annually and dependent on aWTP.
Peat accumulation arises from the balance between the annual addition of new litter layers on top of the mineral soil column and the daily decomposition rate.C originating from different PFTs accrues as litter in the peat layers at variable rates depending on differences in PFT mortality, productivity and leaf turnover.The accumulated peat decomposes on a daily time step based on the plant litter types in each layer of a patch with decomposition rates that are controlled by soil physical and hydrological properties in each layer.Differences in peat decomposition rates among PFTs arise from their intrinsic properties and structure, parameterized using an initial decomposition rate, k o (see Table A1; Aerts et al., 1999;Frolking et al., 2001;Chaudhary et al., 2017a), which is assumed to decline over time (Clymo et al., 1998).
The way plants access water from the mineral soil and dynamic peat layers in each patch, which is dependent on the combined depth of dynamic peat layers and the mineral soil layers, necessitated a readjustment of the soil layer representation relative to the standard version of LPJ-GUESS.In the modified water uptake scheme, there are two static underlying mineral soil layers: an upper mineral soil (UMS) layer and a lower mineral soil (LMS) layer at 0.5 and 1.5 m of depth, respectively.The fraction of roots in these two layers in the absence of peat is prescribed for each PFT and determines the daily plant uptake of water from the mineral soil (Table A1; Chaudhary et al., 2017a).We assigned rooting depth fractions of 0.7 and 0.3 to the shrub PFTs UMS and LMS, respectively, while graminoids were assumed to have relatively shallow rooting depths with fractions of 0.9 and 0.1 in the UMS and LMS, respectively (Bernard and Fiala, 1986;Malmer et al., 2005;Wania et al., 2009b).During the initial stages of peat accumulation, plant roots are still present in both in UMS and LMS, but as peat builds up part of the root fraction is transferred to the growing peat layers, allowing plants to access water from the peat soil.Mosses are assumed to take up water from the top 50 cm of peat (Shaw et al., 2003;Wania et al., 2009b) once peat height exceeds 50 cm.Before this, mosses take water only from the mineral soil.All other PFTs can take up water from both mineral soil layers and peat layers until peat height reaches 2 m, after which they can only access water from the peat soil layers.

Hindcast experiments
To initialize the model with vegetation in equilibrium with early Holocene climate, the model was run from bare ground surface conditions for the first 500 years by repeatedly recycling the first 30 years of the Holocene climate dataset (see below).The mineral and peat layers were forced to remain saturated for the entire initialization period.The peat decomposition, soil temperature and water balance calculawww.biogeosciences.net/14/4023/2017/Biogeosciences, 14, 4023-4044, 2017 Location of 180 randomly selected simulation sites spread across 10 geographical zones between 45 and 75 • N.
tions began when the peat column reached a minimum thickness of 0.5 m.We adopted this model initialization strategy to avoid a sudden collapse of the peat column in very dry conditions: continuous dry periods tend to increase temperaturedependent decomposition, particularly for shallow peat layers, reducing the accumulation rate.
To adequately represent the peatland history and dynamics across the major bioclimatic domains of the pan-Arctic region, the model was applied at 180 grid points (referred to as "sites" below) distributed among 10 geographical zones spanning the circum-Arctic from 45 to 75 • N (Fig. 1); each zone is represented by 10-20 randomly selected points (see Fig. 1 and Table 1).While peatland initiation started at ca. 12-13 kyr BP in high-latitude areas, the majority of peatlands formed after 10 kyr BP (MacDonald et al., 2006).Therefore, each simulation was run for 10 100 years and was comprised of three distinct climate forcing periods.The first phase, the Holocene, lasted from 10 kyr before present (BP) until 0 BP.During this period, the model was forced with daily climate fields (temperature, precipitation and cloudiness) constructed by interpolating between monthly values from 10 000 calendar years before present (cal BP) until 1900.The monthly Holocene climate forcing data were prepared by the delta-change method by applying relative monthly anomalies in temperature and precipitation for the nearest GCM grid cell (see Sect. 2.3.2) to the site location to their average monthly values from the CRU TS 3.0 global gridded climate dataset (Mitchell and Jones, 2005) from the period 1901 to 1930.We then linearly interpolated the values between the millennium time slices to get values for each year of the simulation.This method conserves interannual variability for temperature and precipitation from the baseline historical climate  throughout the simulation.Finally, the monthly Holocene temperature values were interpolated to daily values, while total monthly precipitation was distributed randomly among the number (minimum 10) of rainy days per month.For cloudiness, the monthly CRU values from the years 1901-1930 were re- peated for the entire simulation period.The second historical phase ran from 1901 until 2000.During this period, we forced the model with the CRU data.Finally, the future scenario phase (see Sect. 2.3.2) ran from 2001 until 2100, applying anomalies extracted for the RCP8.5-forcedGCM climate fields (Sect.2.3.2) for each location.Annual CO 2 concentration values to force our model from 10 kyr BP to 1850 AD were interpolated from the millennial values used as a boundary condition in the Hadley Centre Unified Model (UM; Miller et al., 2008) time slice experiments that were run for each millennium from 10 kyr BP to 1850 AD.From the year 1850 to 2000, we used CO 2 values from atmospheric or ice core measurements.
Accurate prediction of total C accumulation at any particular location depends on selecting the right inception period, the C content and lability of the peat material, its bulk density over time and depth and local hydro-climatic conditions (Clymo, 1992;Clymo et al., 1998).Bulk density and C fraction values vary widely among different peatlands, and reliable estimates are often lacking (Clymo et al., 1998).Basal ages, which are proxies for peatland initiation history, are often hard to determine and are not available for many key peatland types.For example, eastern Siberia and European Russia are regions that have not been well studied in this regard (Loisel et al., 2014;Yu et al., 2014a).We therefore started simulations at the same time (10 kyr BP) for all 180 sites and fixed initial bulk densities to 40 kg C m −3 .T8.5 RCP8.5 temperature only 3.
FTPC8.5 RCP8.5 including all treatments The carbon accumulation rate (CAR) of a peatland is the balance between biological inputs (litter addition) and outputs (decomposition and leaching); both input and output fluxes are quite sensitive to climate variability (Clymo, 1991).The long-term (apparent) rate of C accumulation (LARCA) expresses the rate of C accumulated in a peatland since its inception (Clymo et al., 1998) and is a useful metric of the sequestration capacity of peatlands because the current C uptake rate (ARCA; here specified as the recent 30 years) is a snapshot in time that is not expected to reflect the C balance dynamics through the history of the peatland (Lafleur et al., 2001;Roulet et al., 2007).We calculated the rate of C accumulation as LARCA and as the actual (net) rate of C accumulation (ARCA; see Fig. 2).We also calculated the near future rate of C accumulation (NFRCA) from 2001 to 2100 for the 10 studied zones (see below).

Climate change experiments
To investigate the sensitivity of CAR to climate change, future experiments were performed (see Table 2) by extending the base experiment (BAS) covering the Holocene and recent past climate (to year 2000) for an additional century to the year 2100 (Table 2).Climate output from the Coupled Model Intercomparison Project Phase 5 (CMIP5) RCP8.5 (Moss et al., 2010) runs performed with the Hadley Global Environment Model 2 (HadGEM2-ES; Collins et al., 2011) was used to provide anomalies for future climate forcing.HadGEM2-ES is an updated version of the same model chosen for the Holocene anomaly fields.It is in the middle of the range of models contributing to the CMIP5 ensemble in terms of simulated temperature change across the Arctic region (Andrews et al., 2012;Klein et al., 2014).Atmospheric CO 2 concentrations for model input were taken from the RCP8.5 scenario extracted from the International Institute for Applied Systems Analysis website (IIASA; http: //tntcat.iiasa.ac.at/RcpDb/; page visited 14 June 2017).In the first three experiments, the single-factor effect of temperature (T8.5), precipitation (P8.5) and CO 2 (C8.5) was examined, followed by a combined experiment (FTPC8.5) in which change in all three drivers was used to force the model.The model output variables examined here include total CAR, net primary productivity (NPP), net ecosystem C exchange www.biogeosciences.net/14/4023/2017/Biogeosciences, 14, 4023-4044, 2017 ) for each zone (Fig. 1) for each 1000-year period for the last 10 000 years (NEE), permafrost distribution, active layer depth (ALD) and regional soil C balance.

Model evaluation
To evaluate the model, we compared simulated CAR with regional Holocene C accumulation records synthesized across the pan-Arctic region, hereinafter referred to as the "literature range".We also compared the model results for millennial time slices with Holocene LARCA values based on the 127 sites analysed by Loisel et al. (2014) and the 33 sites analysed by Yu et al. (2009).The Loisel et al. (2014) dataset is more comprehensive and contains more basal points compared to Yu et al. (2009).In Yu et al. (2009), many key regions, such as the Hudson Bay Lowlands, western Europe, and western and eastern Siberia, are not present, while the Loisel et al. (2014) dataset omits some regions, such as eastern Siberia and European Russia.Furthermore, the points in these two datasets were limited to areas south of 69 • N (< 69).
3 Results In Sect.3.1, we discuss the simulated temporal and spatial patterns of peatland C accumulation across the pan-Arctic re- gion.Drivers and response mechanisms underlying the simulated patterns are discussed in Sect.3.2.

Hindcast experiment
The mean modelled CAR among all 180 sites was 35.9 g C m −2 yr −1 , after which it followed a similar temporal pattern to observed CAR values (Fig. 3a; Yu et al., 2009;Loisel et al., 2014).The observed rate calculated by Yu et al. (2009) shows a dip after 5 kyr BP, but the modelled result exhibited no such deviation (Fig. 3a).The observed rate reported by Loisel et al. (2014) is a little higher than the simulated rate before 4 kyr BP and for the present climate.Modelled CAR was higher at the beginning of the simulation except in Zone J (Fig. 3b).Zones A and B covering the Scandinavian and European regions had high CAR in the beginning of the Holocene, which then declined through the Holocene, while Zone E covering eastern Siberia displays a peak suggesting an accelerated rate of C accumulation by the year 1900.Almost all regions exhibited similar CAR for 7-8 kyr BP and followed different trajectories thereafter.Scandinavia (Zone A), Europe (B), south-western Siberia (D), central Canada (H) and northern Canada (J) exhibit lower LARCA values compared to the pan-Arctic average (Fig. 4 I and bars; positive bar value means C source) with northern Canadian (J) and European (B) sites accumulating the lowest amounts of carbon (14.5 ± 14.8 and 14.2 ± 3.7 g C m −2 yr −1 , respectively) through the Holocene.The other five zones (C, E, F, G and I) showed relatively higher mean LARCA values, and the peatlands in eastern Siberia (E), Alaska (F) and western Canada (G) had the highest mean LARCA values (26.8 ± 13.8, 26.4 ± 16.3 and 26.6 ± 14.7 g C m −2 yr −1 , respectively).The global mean LARCA (black dashed line) for the 10 zones was 20.8 ± 12.3 g C m −2 yr −1 (Fig. 4I and Table 1).
Comparing mean ARCA for each zone with the respective LARCA values indicates that the majority of sites accumulated relatively more C in the last 30 years except Scandinavia (A), while in Europe (B) the changes were almost negligible (Fig. 4II).The global mean ARCA for the last 30 years was 29.4 ± 27.8 g C m −2 yr −1 , suggesting an upward trend in CAR since the beginning of the Holocene (Fig. 4II and Table 1).
Interpolated values of permafrost (characterized in this study by ice fraction in the peat soil), ALD, CAR and accumulated litter are presented for the recent past and future climate in Figs. 5, 6 and A1. Figure 5a shows that permafrost was widely distributed from Siberia to Canada and in parts of northern Scandinavia around the end of the 20th century according to our model.The majority of these permafrost areas were associated with shallow active layers (ALD < 100 cm), while in the southern parts of Siberia and Canada the active layers are relatively deeper (Fig. 5d).The presence of permafrost shows no simple relationship to peatland CAR (Fig. 6a), ranging from moderate to high litter accumulation in different permafrost areas (Fig. A1a).Large parts of western Canada, Alaska and Siberia accumulated relatively high amounts of C by the year 2000 (Fig. A1a) according to our model.

Climate change experiment
In the FTPC8.5 experiment, in which all the drivers were combined, the global mean FLARCA (20.78 g C m −2 yr −1 ) was largely unchanged from the mean LARCA (20.8 g C m −2 yr −1 ; see Figs. 2, 4I and 7I).However, the change in CAR was quite evident in certain geographic zones (Fig. 7I and bars; positive bar value means C source).Some regions showed an increase in C accumulation, while others become C neutral or sources of C. While Scandinavian (A), European (B) and central and eastern Canadian (H, I) sites are projected to become C sources (Fig. 7I and bars), the remaining zones are projected to become stronger sinks in this scenario.For example, the uptake capacity of northern Canadian (J) sites is projected to increase fourfold, to 52.3 ± 19.2 g C m −2 yr −1 from (its LARCA value of) 14.5 ± 14.8 g C m −2 yr −1 (Table 1 and Fig. 7I).All zones showed a decline in CAR in the T8.5 experiment relative to the recent historical climate (Fig. 7II); the positive effects of temperature on soil organic matter decomposition rates explain this change.An exception to this general pattern is seen for northern Canada (Zone J) where warming has a positive effect on CAR (Fig. 7II and bars): higher temperatures create a more suitable environment for plant growth in this region where cold weather and permafrost limit plant (and therefore litter) production under present climate conditions (see Fig. A2j).The mean modelled global NFRCA in the T8.5 experiment from 2000 to 2100 was 1.52 g C m −2 yr −1 (Fig. 7II; black dashed line).This was a significant drop when compared to modelled LARCA and ARCA.In this experiment, the ESM-derived (Collins et al., 2011) surface air temperature anomalies used to force our model increase by approximately 5 • C by 2100 relative to 2000.Higher temperature is associated with elevated decomposition rates, leading to more C loss and higher heterotrophic respiration.Projected precipitation increases in the P8.5 experiment resulted in higher CAR in all zones (Fig. 7 III and bars).Regionally, Siberian and far eastern Russian (C, D, E), Alaskan (F) and Canadian (G, H, I, J) sites showed the largest changes, while very little change was seen for Scandinavia (A) and Europe (B).Elevated atmospheric CO 2 enhanced photosynthesis, which led to higher CAR in the C8.5 experiment in all zones (Fig. 7IV and bars).
Our simulations suggest that the significant temperature increase implied by the RCP8.5 future scenario will lead to the disappearance or fragmentation of permafrost from the peat soil and deeper active layers (Fig. 5b and e).Additional soil water changes resulting from the effects of higher temperatures on evapotranspiration rates could then either suppress or accelerate the decomposition rate at many peatland locations (Fig. 7II).Effects of precipitation changes and rising CO 2 concentrations on plant productivity can offset decomposition changes in terms of effects on the peat accumulation rate.In the Siberian (C, D and E) and Alaskan (F) zones, the projected higher decomposition rates are compensated for by higher plant productivity due to increases in soil moisture and CO 2 fertilization (Fig. 7III and IV; bars), leading to a net increase in CAR by 2100 in this scenario.
From Fig. 5b, it is evident that permafrost area declines, remaining limited to central and eastern parts of Siberia and the northern Canadian region under the future experiment in our model (FTPC8.5).Permafrost disappears from large parts of western Siberia and southern parts of Canada with very little remaining presence in Scandinavia (Fig. 5b).This degradation (Fig. 5f) leads to wetter conditions initially in large areas of peatlands currently underlain by permafrost.Wetter conditions together with CO 2 fertilization lead to high CAR in these areas with high C build-up.In contrast, non-permafrost peatlands showed a decline in CAR and in total litter accumulation due to higher decomposition rates (Figs.6b, c and  A1b, c) as a result of increases in evapotranspiration, which draw down WTP.

Discussion
Recent CAR tends to be higher compared to LARCA because older peat would have experienced more decay losses, leaching and erosion (Lafleur et al., 2001).This is clearly reflected in our result (Table 1) where LARCA < ARCA in most cases, even though in our study only decay losses were considered.The variability in LARCA among sites within a region with relatively similar climate highlights the influence of local factors (Borren et al., 2004).If climate was the major driving factor behind observed variations in LARCA, then all the peatland types within one climate zone would be ex-pected to have similar LARCA values.LARCA is highly influenced by local hydrology, topography, climate conditions, permafrost, fire events, substrate, microtopography and vegetation succession (Clymo, 1984;Robinson and Moore, 2000;Beilman, 2001;Turunen et al., 2002;Turetsky et al., 2007).Some studies attribute differences in LARCA values to the overrepresentation of terrestrialized peatlands and an underrepresentation of paludified or shallow peatlands (Botch et al., 1995;Tolonen and Turunen, 1996;Clymo et al., 1998) in estimations of this metric.Our model initialization allowed vegetation to reach an equilibrium with the climate of 10 kyr ago, but the model ignores the presence of ice over some parts of the study area at this time, thus overestimating the vegetation cover at the beginning of the simulation and leading to higher CAR than observed (Fig. 3a, b).In addition, the underlying topography is a major factor for peat initiation and lateral expansion of any peatland complex, but no such data are available for regional simulations.Therefore, we assumed a moist and on average uneven horizontal soil surface upon which peatland could potentially form at each of our 180 simulation points, ignoring the role of underlying topography and its effects on water movement within a basin (Tang et al., 2015).However, the lateral exchange between higher and lower patches within an overall horizontal landscape was included in our model (see Sect. 2).
The mean modelled LARCA across the pan-Arctic study area was 20.8 ± 12.3 g C m −2 yr −1 , which is a value that falls within the reported range for northern peatlands, namely 18.6-22.9g C m −2 yr −1 (Yu et al., 2009;Loisel et al., 2014).However, the Loisel et al. (2014) dataset is not completely representative of the pan-Arctic region, and data from some key regions are missing, such as eastern Siberia and European Russia (Yu et al., 2014a).The Loisel et al. (2014) dataset includes points that are mainly from deep or central parts and shallow peat basins are underrepresented (Mac-Donald et al., 2006;Gorham et al., 2007;Korhola et al., 2010).Furthermore, the dataset is limited to areas south of 69 • N. Inclusion of shallow peatland complexes and more subarctic and arctic sites in the syntheses might conceivably bring down the mean observed pan-Arctic LARCA value.Nevertheless, the overall trend of the modelled pan-Arctic averaged CAR (n = 180) for the last 10 kyr is quite similar to these published syntheses (Fig. 3a and b and Table 1).
Suitable climate and optimal local hydrological conditions influenced by favourable underlying topographical settings accelerated CAR, which led to the formation of large peatland complexes in the pan-Arctic region (Yu et al., 2009).High CAR is associated with high plant productivity and a moist climate, leading to shorter residence time in acrotelm layers with generation of recalcitrant peat or a combination of any of these factors (Yu, 2006).In many regions, CAR is also influenced by the presence of permafrost.Under stable or continuous permafrost conditions, CAR slows down or ceases (Zoltai, 1995;Blyakharchuk and Sulerzhitsky, 1999) due to low plant productivity.CAR may also become negative due to wind abrasion and thermokarst erosion, but these factors are not considered in our simulations.In contrast, areas underlain by sporadic and discontinuous permafrost sequester relatively more C (Kuhry and Turunen, 2006).
Significant increases in temperature are expected at high latitudes in the coming century, even under the most optimistic emissions reduction scenarios.Under these conditions, some peatlands could sequester more C (Charman et al., 2013), while others could turn into C sources and degrade (Ise et al., 2008;Fan et al., 2013).Permafrost peatlands are sensitive ecosystems and respond quite rapidly to temperature change and other aspects of climate (Christensen et al., 2004).The formation of thermokarst lakes, degradation of palsa, flooding and subsidence of the land surface are key features that might indicate and result from rapid warming and permafrost decay.Soil subsidence-driven pond formation has been observed to lead to a total shift from a recalcitrant moss-dominated vegetation community to dominance by non-peat-forming taxa, such as Carex spp.(Malmer et al., 2005).However, the complex physical dynamics inducing such changes are not included in our model.
In our scenario simulations (Table 2), we find that higher temperature leads to thawing of permafrost that in turn increases the moisture availability, at least initially.The rise in temperature also results in early spring snowmelt and a longer growing season (Euskirchen et al., 2006), while in the same time frame atmospheric CO 2 concentration will also increase.These factors lead to increases in plant productivity, leading to higher CAR (Klein et al., 2013;Chaudhary et al., 2017a), even in cases where moisture-and temperaturedriven peat decomposition also speeds up.
High temperature and limited moisture conditions with limited or no permafrost have been generally found to accelerate peat decomposition (Franzén, 2006;Ise et al., 2008;Bragazza et al., 2016).This will also result in the drawdown of water position and dominance of woody shrubs.The latter trend, namely an expansion of shrubs across the Arctic and   (Sturm et al., 2005;Loranty and Goetz, 2012).
Conversely, warmer and wetter future climate conditions, in combination with CO 2 fertilization, could lead to increased CAR in areas projected to have a higher precipitation rate, compensating for the temperature enhancement of decomposition.
We now go on to discuss the simulated responses of peatland to the differential climate conditions of the studied regions in relation to available literature.
The modelled LARCA (14.2 ± 3.7 g C m −2 yr −1 ) for central and eastern Europe (Zone B) is relatively low.However, while some sites in this region are reported as being quite productive (21.3 ± 3.7 g C m −2 yr −1 ; Anderson, 2002), long-term CAR estimates are available for relatively few sites (Charman, 1995;Anderson, 1998), making a comparison difficult.The points that fall in the British Isles showed lower modelled LARCA (12-14 g C m −2 yr −1 ) values than the observed literature range, indicating shortcomings in the simulation of local hydrological conditions or a possible bias in the climate forcing of our model.A decline in CAR in Scandinavia and Europe over recent decades is apparent in our simulations.Some observational studies also point to a reduced rate of C accumulation in recent years for this region (Clymo et al., 1998;Klarqvist et al., 2001b;Gorham et al., 2003).This slowing has been attributed to an increase in decay rates due to climate and hydrological changes, the development of a stable structure (Malmer and Wallen, 1999), divergence in the rate of nutrient supply or a combination of these factors (Franzén, 2006).Our model predicts that the C balance of Scandinavian peatlands will decrease after 2050 and become C neutral, with peatland in the European region becoming a C source in the same time frame (Fig. 7I zones A  and B).The simulated future C losses are associated with an increase in the decomposition rate due to higher temperatures and a lower soil water table, the latter resulting from the combination of marginal or no increase in precipitation and soil water loss due to higher evapotranspiration.

Siberia (zones C and D) and far eastern Russia
(zone E) Large peatland complexes were formed in western Siberia during the Holocene and around 40 % of the world's peat deposits are found in this region, covering more than 300 million ha (Turunen et al., 2001;Bleuten et al., 2006).LARCA for western Siberia has been estimated at 5.4 to 38.1 g C m −2 yr −1 (Beilman et al., 2009).The modelled LARCA for the north-west and south-west region is 24.6 ± 14.6 and 16.7 ± 8.6, respectively (Fig. 4I Zones C and D and Table 3).The combined average modelled LARCA for the northern and south-western Siberian (C + D) zones is 20.6 g C m −2 yr −1 .Turunen et al. (2001) report average LARCA from 11 sites in north-western Siberia at 17.3 g C m −2 yr −1 (range from 12.1 to 23.7 g C m −2 yr −1 ).Botch (1995) estimated relatively higher LARCA (31.4-38.1 g C m −2 yr −1 ) for the raised string bogs in western Siberia.These observations are in line with our modelled range of 24.6 ± 14.6 g C m −2 yr −1 for the north-western sites.Borren et al. (2004) found LARCA values between 19 and 69 g C m −2 yr −1 for the southern taiga zones of southwestern Siberia.The modelled LARCA value for the southwestern zone (D) is 16.7 ± 8.6 g C m −2 yr −1 .The apparent underestimation by our model could be explained by the relatively larger area encompassed by our simulations, extending into warmer southerly areas with limited peat accumulation compared to the aforementioned study (Fig. 5 Zones C and D and Table 3).Borren and Bleuten (2006) modelled a LARCA range of 10-85 g C m −2 yr −1 (mean 16 g C m −2 yr −1 ) for a large mire complex in south-western Siberia, and our value falls within this range.
The mean observed LARCA was 10.6 ± 5.5 g C m −2 yr −1 for a permafrost polygon peatland of far eastern Russia (Gao and Couwenberg, 2015).Botch et al. (1995) cite CAR values of 44.8 g C m −2 yr −1 for both the Kamchatka and Sakhalin regions and 33.6 g C m −2 yr −1 for far eastern regions.Our modelled estimate of 26.8 ± 13.8 g C m −2 yr −1 is broadly comparable to the range of these observations.
Our model predicted that the sink capacity (22.7 g C m −2 yr −1 ) of the entire Russian region (C, D and E) was higher than the pan-Arctic average (Fig. 4 and Table 3).In the future, higher temperature and precipitation, together with increases in snow depth, result in permafrost degradation that will lead to a deeper active layer in the western part of Siberia (Fig. 5b, e).Plants experience improved hydrological conditions due to a deeper ALD.Thawing of the permafrost in the peat and mineral soils coupled with a longer growing season and CO 2 fertilization leads to higher plant productivity, offsetting the higher decomposition rate and leading to an increase in CAR (Fig. 6b, c).Hence, this region is projected to act as a C sink in the future (Fig. 7I).It is notable in our simulations that temperature increases in the T8.5 experiment have a very limited overall effect on decomposition rates in Russia (Zones C, D and E), while precipitation and CO 2 fertilization have a positive effect on C build-up (Fig. 7II, III and IV).
4.3 Canada (Zones G to J) and Alaska (Zone F) Canada's Mackenzie River basin and the Hudson Bay Lowlands are two of the largest peatland basins in the world (Beilman et al., 2008).The individual observed C accumulation rates vary considerably across Canada, and the LARCA for the entire Canadian region ranges from 0.2 to 45 g C m −2 yr −1 (see Table 3).The modelled mean LARCA value averaged among Zones G to J (the entire Canadian region) is 21.2 g C m −2 yr −1 .Most observational studies have been carried out in the western and central regions of Canada (Halsey et al., 1998;Vitt et al., 2000;Beilman, 2001;Yu et al., 2003;Sannel and Kuhry, 2009).However, in recent years, studies have been conducted in the Hudson Bay Lowlands and the James Bay Lowlands of eastern Canada (Loisel and Garneau, 2010;van Bellen et al., 2011;Bunbury et al., 2012;Lamarre et al., 2012;Garneau et al., 2014;Holmquist and MacDonald, 2014;Packalen and Finkelstein, 2014).Observed LARCA in Zone I is relatively low, as peatlands initiated later in this region due to a late Holocene thermal maximum (5.0-3.0 kyr; Yu et al., 2009) and the presence of the remnants of the Laurentide ice sheet (Gorham et al., 2007).In our model simulations, all peatlands were initiated at the same time and we have not considered the influence of ice sheet cover, which explains the higher modelled CARs (25.3 ± 11.8 g C m −2 yr −1 ) in the eastern region.The observed LARCA of the three main eastern regions in Canada is as follows: Quebec (26.1 g C m −2 yr −1 ; Garneau et al., 2014), Hudson Bay Lowlands (18.5 g C m −2 yr −1 ; Packalen and Finkelstein, 2014) and James Bay Lowlands (23.9 g C m −2 yr −1 ; Holmquist and MacDonald, 2014).Other studies in the area have similar values (see Table 3).Our simulations suggest that permafrost will disappear from large areas of southern Canada under the RCP8.5 climate change scenario, leading to deeper ALD (Fig. 5b, e).While western and northern Canadian regions sequester C at higher rates from 2001 to 2100 in our simulations, central and eastern parts turn into a C source over the same time period (Fig. 6c).Decomposition rates will increase due to higher temperatures, overriding the positive gains due to precipitation and C fertilization in central and eastern regions (Fig. 7 Zones H and I).
The majority of simulated points in northern Canada (Zone J) are in the continuous or discontinuous permafrost region (Sannel and Kuhry, 2009).Observed LARCA values in this zone vary from 0.2 to 16.5 g C m −2 yr −1 (see Ta-ble 3).Similarly, the modelled CAR of the northern Canadian sites was lowest (14.5 ± 14.8 g C m −2 yr −1 ) as a result of cold climate conditions (Table 4).The mean temperature in this zone is around −15 • C with a short growing season and low precipitation, the majority of which falls as snow.In some sites, negligible CARs were noticed due to extremely cold climate conditions that limited plant productivity.In other subarctic regions, similar effects of cold climate and permafrost conditions have been observed.For instance, LARCA ranges from 12.5 to 16.5 g C m −2 yr −1 for the central polygon peatlands in western Canada (Vardy et al., 2000) and 11 g C m −2 yr −1 in the northern Yukon (Ovenden, 1990).Similarly, polygon peat plateaus in eastern Siberia have sequestered C at low rates (10.2 g C m −2 yr −1 ; Gao and Couwenberg, 2015).Lately, owing to recent climate warming and permafrost thaw, bioclimatic conditions have changed in these peatlands and many of them have seen twofold to threefold increases in CAR (Ali et al., 2008;Loisel and Garneau, 2010), indicating a recent shift toward an increased C sink capacity.A fourfold increase in CAR associated with permafrost thaw and increased primary productivity was simulated under future warming by our model (Table 1 and Fig. 7 Zone J).
Alaska hosts around 40 million ha of peatland area (Kivinen and Pakarinen, 1981).Studies show that LARCA in this region ranges from 5 to 20 g C m −2 yr −1 (see Table 3).Our modelling results (26.4 ± 16.3 g C m −2 yr −1 ) may be overestimations (Table 1 and Fig. 4 Zone F).The higher CAR values in our simulations are caused by high plant productivity, moist climate conditions, the generation of recalcitrant peat or a combination of these factors.This overestimation of CAR in Alaska casts doubt on the simulated large future sink capacity of the study area (55.5 ± 16.3 g C m −2 yr −1 ) under the RCP8.5 scenario.

Future climate impacts on peatlands
Our simulations under the RCP8.5 future forcing indicate a sharp reduction in the area underlain by permafrost, for example in western Siberia and western Canada, leading to an initial increase in moisture conditions or wet surfaces there.The increase in moisture conditions can dampen the amplifying effects of temperature on decomposition rates, leading to net increase in CAR (Figs. 5,6 and 7).By 2100, our model indicates that permafrost areas will be limited to eastern Siberia, northern and western Canada and parts of Alaska (Fig. 5b).
In the future, areas currently devoid of permafrost, mainly Europe and Scandinavia, eastern parts of Canada and European Russia, could lose a substantial amount of C due to the drying of peat in conjunction with a deeper WTP (Figs. 6  and 7).In a modelling study, Ise et al. (2008) Borren and Bleuten (2006), using a three-dimensional dynamic model with imposed artificial drainage to simulate the Bakchar bog in western Siberia, indicated that LARCA will drop from 16.2 to 5.2 g C m −2 yr −1 during the 21st century due to higher decomposition linked to reduced peat moisture content.Our simulations are based on climate forcing derived from the RCP8.5 scenario output from one Earth system model (HadGEM2-ES).We expect that simulated changes in permafrost and C accumulation would be more moderate and slower if the model were forced with more moderate levels of climate change.
Overall, we found that Scandinavia, Europe, Russia and central and eastern Canadian sites could turn into C sources, while C uptake could be enhanced at other sites (Figs. 6  and 7).The greatest changes were evident in eastern Siberia, north-western Canada and in Alaska.Peat production was initially hampered by permafrost and low productivity due to the cold climate in these regions, but initial warming coupled with a moisture-rich environment and greater CO 2 levels could lead to rapid increases in CAR by 2100 in this scenario.In contrast, sites that experience reduced precipitation rates and that are currently without permafrost could lose more C in the future.

Conclusion and outlook
Our model, which among large-scale models of high-latitude peatlands uniquely accounts for feedbacks between hydrology, peat properties, permafrost and dynamics of vegetation across a heterogeneous peatland landscape, is able to reproduce broad, observed patterns of peatland C and permafrost dynamics across the pan-Arctic region.Under a business-as-usual future climate scenario, we showed that non-permafrost peatlands may be expected to become a C source due to soil moisture limitations, while permafrost peatlands gain C due to an initial increase in soil moisture, which suppresses decomposition while enhancing plant production.We also demonstrate that the extant permafrost area will be reduced and limited to central and eastern parts of Siberia and the northern Canadian region by the late 21st century, disappearing from large parts of western Siberia and southern parts of Canada with very little presence in Scandinavia.Our modelling approach contributes to an understanding of long-term peatland dynamics at a regional and global scale.As such it complements empirical research in this field but also synthesizes the implications of current empirical knowledge and understanding, on the basis of which our model was constructed and evaluated.We plan to incorporate methane biogeochemistry and nutrient dynamics in the next model update.In the future, the model will be coupled to the atmospheric component of a regional Earth system model to examine the role of peatland-mediated biogeochemical and biophysical feedbacks to climate change in the Arctic and globally.
Code and data availability.Model code can be inspected by contacting the corresponding lead author, Nitin Chaudhary, or Paul Miller (paul.miller@nateko.lu.se).Readers who would like to use our code in their own research can contact Paul Miller directly for information on conditions of use.

Figure 2 .
Figure 2. Commonly used measures of peat accumulation rate: long-term (apparent) rate of C accumulation (LARCA), recent rate of C accumulation (RERCA), actual (true) rate of C accumulation (ARCA), simulated future long-term (apparent) rate of C accumulation (FLARCA) and near future rate of C accumulation (NFRCA; adapted from Rydin and Jeglum, 2013).

Figure 3 .
Figure 3. (a) Simulated and observed mean C accumulation rate (g C m −2 yr −1 ) for each 1000-year period for the last 10 000 years.Red: simulated mean (and standard error of the mean) CAR based on 180 random sites.Blue and black points are observed C accumulation rates (g C m −2 yr −1 ) based on 127 (Loisel et al., 2014; blue points) and 33 sites (Yu et al., 2009; black points) across the northern peatlands with error bars showing the standard errors of the means.(b) Mean C accumulation rate (g C m −2 yr −1 ) for each zone (Fig. 1) for each 1000-year period for the last 10 000 years

Figure 4 .
Figure 4. Simulated Holocene peat accumulation rates across the 10 zones considered in this study (blue dots) and for the pan-Arctic region as a whole (dashed black line).The x axis shows the number of sites partitioned into 10 zones.The black dashed line is the pan-Arctic average with standard deviation (black line outside the y axes) and the red dashed line is the average among zones with the standard deviation as a light red patch.(I) Simulated long-term (apparent) rate of C accumulation (LARCA); (II) simulated actual (true) rate of C accumulation (ARCA) for the last 30 years.Blue bars show the difference between ARCA and LARCA mean values for the respective zone (II − I).

Figure 5 .
Figure 5. Modelled September ice fraction (0-1) in the peat soil (as a proxy for permafrost distribution) interpolated among simulation points averaged over (a) 1990-2000 and (b) 2090-2100.(c) Continuous and discontinuous permafrost zones and the modelled mean September active layer depth (ALD in cm) interpolated among simulation points for (d) 1990-2000 and (e) 2090-2100.(f) Net change in total ALD (e-d).

Figure 7 .
Figure 7. Simulated C accumulation rate (blue lines) for each zone (refer to Figs. 1 and 4) and across the pan-Arctic region.The black dashed line is the pan-Arctic average with standard deviation (black line outside); the red dashed line is the average for the respective zone with the standard deviation as a light red patch.(I) Average simulated near future rate of C accumulation (NFRCA) for 2001-2100 in the FTPC8.5 experiment; (II) simulated NFRCA in the T8.5 experiment, (III) simulated NFRCA in the P8.5 experiment and (IV) simulated NFRCA in the C8.5 experiment.Blue bars show the difference between the FLARCA and LARCA values for each zone.

Table 1 .
Mean modelled C accumulation rates at different timescales in 10 geographical zones.

Table 2 .
Summary of hindcast and global change experiments.

Table 3 .
Observed regional long-term rate of peatland C accumulation across northern latitude areas.

Table 3 .
Continued.Chaudhary et al.:Modelling past, present and future peatland carbon accumulation beyond in the next half of the 21st century, is in keeping with other studies a FSU is the former Soviet Union.b CAR over the past 4000 years.www.biogeosciences.net/14/4023/2017/Biogeosciences, 14, 4023-4044, 2017 4034 N.
used a coupled physical-biogeochemical soil model at a site in northern Manitoba, Canada and found that peatlands could respond quickly to warming, losing labile soil organic carbon dur-

Table A1 .
Plant functional types (PFTs) simulated in this study, showing representative taxa, phenology, bio-climatic limits, water table position (WTP) threshold for establishment, prescribed root fractions in mineral soil layers and initial decomposition rate for different litter fractions.