Peatland area and carbon over the past 21,000 years – a global process based model investigation

Peatlands are an essential part of the terrestrial carbon cycle and the climate system. Understanding their history is key to understand future and past land-atmosphere carbon fluxes. We performed transient simulations over the past 22,000 years with a dynamic global peat and vegetation model forced by Earth System Model climate output, thereby complementing databased reconstructions for peatlands. Our novel results demonstrate a highly dynamic evolution with concomitant gains and 5 losses of active peatland areas. Modelled gross area changes exceed net changes several fold, while net peat area increases by 60 % over the deglaciation. Peatlands expand to higher northern latitudes in response to warmer and wetter conditions and retreating ice sheets and are partly lost in mid-latitude regions. In the tropics peatlands are partly lost due to flooding of continental shelves and regained by non-linear interactions between temperature, precipitation and CO2. Large north-south shifts of tropical peatlands are driven by shifts in the position of the Inter Tropical Convergence Zone associated with the abrupt 10 climate events of the glacial termination. Time slice simulations for the Last Glacial Maximum (LGM) demonstrate large uncertainties in modelled peatland extent (global range: 1.5 to 3.4 Mkm) stemming from uncertainties in climate forcing. Net uptake of atmospheric CO2 through peatlands, modelled at 350 GtC since the LGM, includes decay from former peatlands. Carbon uptake would be misestimated, in particular during periods of rapid climate change and subsequent peatland area shifts, when considering only changes in the area of currently active peatlands. Our study highlights the dynamic nature of peatland 15 distribution and calls for an improved understanding of former peatlands to better constrain peat carbon sources and sinks.


Introduction
Peatlands are a wetland landscape type that is characterized by permanently waterlogged conditions, resulting in accumulation of dead plant material as peat (Gorham, 1957;Moore, 1989;Blodau, 2002). Peatlands are globally distributed and can take 20 multiple forms from minerotrophic fens to ombrotrophic bogs and forrested tropical peat swamps (Rydin and Jeglum, 2013;Page and Baird, 2016;Lindsay, 2018). Peatlands cover less than 3% of the global land area (Xu et al., 2018), but store a share (Sitch et al., 2003;Xu-Ri et al., 2012;Spahni et al., 2013). Implementation of the long term terrestrial carbon stores, permafrost and peatlands, are based on the LPJ-WHyMe model (Wania et al., 2009a, b) and a module to simulate peat area dynamically (Stocker et al., 2014). Peatlands are represented as a separate land class within a grid cell. The area of each grid cell is split into a fraction covered by the land classes "peat", "mineral soils" and "old peat" (former active peat now treated as mineral soils).
In this study, anthropogenic land use and land use change and corresponding land classes are not considered. Peat vegetation 5 is represented by five peat PFTs: Sphagnum and flood tolerant graminoids as indicative mostly for high latitude peatlands, and flood tolerant tropical evergreen and decidious tree PFTs as well as a flood tolerant version of the C4 grass PFT, as indicative mostly for tropical peatlands (Stocker et al., 2014). Carbon cycling in peat soils is based on the distinction between a lower, fully water saturated slow overturning pool (catotelm) and an upper fast overturning pool (acrotelm) with fluctuating water table position (WTP) . Acro-to catotelm flux is determined by an average acrotelm bulk density and the 10 acrotelm mass balance, determined by organic influx from the litter pools and heterotrophic respiration within the acrotelm.
Decay rates are modulated by temperature in the catotelm and by temperature and WTP in the acrotelm.
The area fraction covered by peat (f peat ) in a given grid cell is determined dynamically with the DYPTOP module (Dynamical Peatland Model Based on TOPMODEL) (Stocker et al., 2014). The TOPMODEL approach (Beven and Kirkby, 1979) is sea-ice mask is changing every 1 kyr according to Peltier (2004) and is interpolated in between. The model is spun up under LGM conditions for 2.5 kyr before starting the transient simulations.
Additional to the standard LPX transient simulation, five transient factorial simulations were performed using the same setup but for each keeping one of the five transient forcings (land-sea-ice mask, orbital, CO 2 , precipitation, and temperature) constant at LGM levels. These were used to identify the dominant drivers and driver contributions through time, by comparing 5 the factorial and standard runs. (see Sect. 3.3) To investigate the uncertainty stemming from the choice of climate forcing, seven additional LGM timeslice simulations were performed. Mean LGM climate anomalies from six different PMIP3 models (Braconnot et al., 2012) and the mean LGM anomaly of the TraCE21k simulation were imposed on the CRU 3.1 climatology from 1901-1931. Inter annual variability thus is adopted from CRU. CO 2 , ice-sea-land mask and orbital parameters are set to LGM levels (in this case 21 kyr BP). Two of the 10 eight available PMIP3 LGM simulations were dismissed (FGOALS-g2 and CNRM-CM5), because of the poor performance compared to observational data, especially in the variables of temperature and precipitation (Harrison et al., 2014). Simulations are spun up for 2.5 kyr and run an additional 2 kyr under unchanged conditions. In the analysis the temporal mean over the last 2 kyr is used.

Validation data 15
Even estimates for the current global peatland area are still subject to large uncertainties as peatlands often lie in remote, inaccessible or understudied regions, such as tropical forests or the Arctic tundra. Even estimates for the relatively well studied northern high-latitude peatlands have a range of 2.4 -4.0 Mkm 2 (see Loisel et al. (2017) for a review). Total area of tropical peatlands is even less well defined and estimates range from 0.37-1.7 Mkm 2 (Yu et al., 2010;Page et al., 2011;Gumbricht et al., 2017). The upper end of this range is given by an estimate that uses an expert system method, combining hydrological 20 modelling, satellite imaging, and topographic data and thus tries to also account for still undiscovered peatlands (Gumbricht et al., 2017). The to date most extensive and comprehensive compilation of known peatlands is the recent PEATMAP by Xu et al. (2018). PEATMAP shows a distribution shifted more towards the tropics, than previous literature estimates. For example Yu et al. (2010) estimates the area of northern peatlands to 4 Mkm 2 and of tropical peatlands to 0.37 Mkm 2 , whereas PEATMAP gives 3.18 Mkm 2 and 0.99 Mkm 2 , respectively. In Fig. 1 25 and global distribution are compared against a 0.5°× 0.5°gridded version of PEATMAP.
Measured peat core basal dates have long been used to estimate northern peat initiation and lateral expansion through time. Yu et al. (2013) compiled a dataset containing 2808 basal dates combining published datasets from MacDonald et al. (2006), Gorham et al. (2007) and Korhola et al. (2010). Loisel et al. (2017) used this dataset (MGK13) to produce a version with only the oldest date per 1°x1°gridcell (MGK13G), as a proxy for peatland initiation. The MGK13G dataset is used in this 30 study to compare to simulated northern peat initiation (see Sect. 3.3.4). Multiple local basal dates are needed to disentangle lateral expansion from initiation. Loisel et al. (2017) compiled a reconstruction based on the gridded MGK13 dataset, but only gridcells with three or more peat cores were considered (MGK13S). Expansion curves were built regionally and then stacked to compensate for regional sampling bias. Korhola et al. (2010) used a similar approach using 954 basal dates from 138 sites,  Table 1: Northern Asia (red), Western Siberian lowlands (orange), and Southeast Asia (green) .

Peatland carbon
Total peat carbon estimates are closely linked to the estimates for area and thus inherit their uncertainties. Additional assumptions on bulk density and peat depth introduce additional uncertainties. The range of carbon estimates therefore is similarly large as for area (Gorham, 1991;Turunen et al., 2002;Yu et al., 2010). The research bias allows for more constrained estimates in well studied regions such as Europe and North America and less constrained in the tropics and Northern Asia. Estimates 5 for Northern peatlands range from 270-604 GtC, obtained with various methods and area estimates (see Yu (2012) andYu et al. (2014) for review). The modern carbon inventory of northern peatlands simulated by LPX at the end of the transient standard run from the LGM till present is with 361 GtC well within this observational range. In the tropics, LPX simulates a peat carbon inventory of 135 GtC which is substantially larger than classical literature estimates that range from 50-87 GtC (Yu et al., 2010;Page et al., 2011). These, however, also assume a substantially smaller tropical peatland area than LPX or 10 PEATMAP suggest (see Sect. 3.1.1). Gumbricht et al. (2017) calculate an even larger area than LPX and combining their area estimate with the peat properties assumed by Page et al. (2011) would result in a tropical peat inventory of 350 GtC. Previous studies with LPX-Bern reported somewhat different carbon inventories than given here. Stocker et al. (2014) reported 460 GtC and 88 GtC for northern and tropical peatlands respectively. Differences stem from an updated model version, also resulting in different areas as mentioned in the previous section. Additionally, their carbon stocks were the results of an accelerated spinup scheme, whereas in this study the pools are filled over a transient run. Spahni et al. (2013) also reports northern peatland carbon stocks after a transient LPX run from the LGM, however with prescribed not prognostic peatland 5 area. Their simulation resulted in 365 GtC stored in northern peatlands.
Other model studies with dynamic peatland area reported 317 GtC after an 8 kyr Holocene run (Kleinen et al., 2012) and 463 GtC after a 12 kyr Holocene run (Qiu et al., 2018b) in northern peatlands.

Area distribution and carbon inventories of peatlands during the Last Glacial Maximum
Under LGM conditions global simlulated peatland area and carbon inventories are reduced compared to preindustrial (Table 1).

10
Globally, simulated peatland area and peat C inventory are 38% and 45% smaller at LGM than PI respectively. This reduction is dominated by the northern extratropics, where peat extent and C inventory are by almost 60% smaller at LGM than PI. In contrast, peat C inventory in the tropics is only about 3% smaller and tropical peat area is even 7% larger at LGM than PI.
This difference in the tropics is mostly linked to large peatlands simulated on flat exposed continental shelves in Southeast Asia at the LGM, which were subsequently flooded during the deglaciation. Another modeling study by Kaplan (2002) also 15 suggests extensive wetlands on the flat Sunda Shelf, but reconstructions of Indonesian peatlands suggest that vast peat presence in Indonesia during the LGM is unlikely (Dommain et al., 2014). Establishment of now existing inland peatlands seems to be connected to rising sea level (Dommain et al., 2011). In sediment cores from the now submerged Sunda Shelf, there is little evidence of peatlands during the LGM (Hanebuth et al., 2011). Dommain et al. (2014) suggest that the shelf, although with a small topographic gradient, had an effective drainage system with deeply incised river valleys, preventing the formation of large wetlands. Both, the hydrological feedback of rising sea level and deep river systems, are not represented in LPX and thus might limit the models ability to reproduce peat and wetland dynamics in this region correctly.

5
Simulated peatland coverage in northern mid and high-latitudes is smaller and shifted southwards at LGM compared to PI.
Ice-shields covered large parts of Europe and North America during the LGM preventing vegetation and peat to grow. But also in Northern Asia and the WSL, peat is mostly absent due to the substantially colder and dryer conditions compared to today. On the other hand, large peatland complexes are simulated along the southern ice-shield margins in North America and in Europe ( shared with all PMIP3 models (Lora and Lora, 2018).

Uncertainties from climate forcing
The peat distribution as simulated by LPX-Bern for the LGM and the past 20,000 years is subject to many uncertainties. Uncertainties arise from model parameterisations, not only in the peat module, but through all components of the model, and are often hard to quantify. Data assimilation, as done recently for the LPX in Lienert and Joos (2018)  uncertainties into carbon cycle projections Ahlström et al., 2017).
We assess the uncertainty in peatland area and peat carbon stemming from climate forcing uncertainties. Climate anomalies from seven different models are used to force the LPX into seven different LGM states (see Sect. 2.2). This yields a very wide range for global mean inundated area (2.6-3.6 Mkm 2 ), peat area (1. 5-3.4 Mkm 2 ) and peat carbon (144-343 GtC). Interest-5 ingly, simulated wetland and peatland area and peat C inventory for the 21 kyrBP period are also substantially different between the standard transient simulation using temporally evolving climate anomalies from the Trace21k simulation compared to the timeslice simulation forced with Trace21k anomalies ( Fig. 3; red star versus magenta dot). This highlights both the influence of different methods of input preparation, with slightly different treatment of anomalies and an inter annual variability taken form TraCE21k in the transient simulation and from CRU 3.1 for the timeslice, as well as the importance of memory effects 10 for a slowly reacting system such as peatlands.
Agreement on simulated peat extent among the seven simulations differs among regions ( Fig. 2 (b)). It tends to be higher in the tropics and East Asia, and lower in North America, Europe and West Siberia. Differences in temperature and precipitation anomalies propagate into differences in the water balance and productivity, partly limited by growing season length, and thus 15 into differences in peat abundance and extent.
A statistical analysis of the differences in climatic drivers and simulated peat area reveals regionally different mechanisms ( Fig. 2 (c)). Temperature, precipitation, precipitation minus evapotranspiration (P-E), and growing degree days over 0 • C (GDD 0 ) are considered as climatic predictor variables for the peat fraction within a gridcell. We correlated, for each grid cell, the seven climatological mean values of a selected predictor with the modelled peat fraction from each of the seven time 20 slice simulations. P-E and GDD 0 show significant correlations (p<0.05) in more gridcells than precipitation and temperature, respectively. Both moisture balance and GDD 0 have been shown to be among the most important predictors of northern peat initiation and carbon accumulation in the past (Morris et al., 2018a;Charman et al., 2013). In LPX the water balance, influenced by P-E, and the carbon balance, influenced by temperature and growing season length, define thresholds on peatland existence and size. In eastern Europe, differences in peat extent between the seven LGM time slice simulations are mostly driven by 25 differences in local precipitation anomalies driving P-E. Similar is true in the tropics, with MRI-CGCM3 and IPSL-CM5A-LR being the driest models with the least tropical peatlands and TraCE21k and COSMOS-ASO being the wettest with the most tropical peatlands. In parts of central South America however, temperature is the dominant predictor signaling a fragile carbon balance, where peat presence in some models is possible because of cooler conditions and thus reduced respiration. In the south of North America, moisture balance, with contributions of both P and E, is the dominant determinant of the inter model 30 differences. The timeslice forced with TraCE21k climate shows the peatland distribution in North America more shifted to the east compared to most other PIMP3 forcings alongside warmer and wetter conditions (see also Lora and Lora (2018)).
Peatland extent in the north of North America is sensitive to temperature differences with longer growing season allowing for increased productivity and therefore peat formation. In the MRI-CGCM3 timeslice temperature anomalies with respect to preindustrial are lowest and peat is subsequently shifted northwards compared to other time slices. Similar is true for Northern  peatlands under favorable conditions, but also the decay and disappearance of peatlands under unfavorable conditions. Both processes can happen simultaneously on a global as well as a regional scale (Fig. A1). To treat carbon storage in a consistent manner, we distinguish between the active peatlands, which are treated as peatlands in the LPX, and old peatlands, which are treated as mineral soils. Old peatlands inherit the carbon stocks of the peatlands that are shrinking or vanishing. Similarly, growing active peatlands first expand onto the area of old peatlands inheriting the remaining carbon stored there (see also Sect. 3.3.5). In the analysis, we decompose the net changes of peatland area into gross positive and negative changes. This allows for a deeper insight into the underlying temporal dynamics (Fig. 4 (b)). Transient factorial runs, performed over the same time period as the standard setup (see Sect. 2.2), allow us to attribute driver contributions to the simulated changes ( Fig. 4 (c) and 5). Global changes in peatland area and carbon before the onset of the Heinrich Stadial 1 (HS1) are small, due to the relatively small changes in the main drivers. There is initial carbon loss in some regions of the tropics, due to some gridcells still approaching equilibrium after the spinup (Fig A1). North America sees an accelerating carbon accumulation with unchanging area already before the HS1, driven mostly by increasing temperature. Carbon and area also increase in Europe with large temperature  of peatlands on continental shelves, mostly in South-East Asia, due to the beginning of rising sea levels. 5 The last termination represents the transition of the climate system from the last glacial to the current interglacial, accompanied by large warming, ocean circulation changes, and an increase in atmospheric CO 2 (Monnin et al., 2001;Shakun and Carlson, 2010;Ritz et al., 2013). The termination is divided into the HS1 (17.43 -14.63 kyr BP) northern hemisphere (NH) cold period, the Bølling-Allerød (BA, 14.63 -12.85 kyr BP) NH warm period, and the Younger Dryas (YD, 12.85 -11.65 kyr BP) NH cold period (Rasmussen et al., 2014), These NH cold-warm swings are associated with a large-scale reorganiza-10 tion of ocean circulation, thought to have been provoked by freshwater release from ice sheet melting leading to changes in the ocean heat transport (Stocker and Johnsen, 2003). With changing low-to high-latitude temperature gradients the ITCZ shifted The responses of peatlands in LPX to these climatic changes are drastic. Large shifts in peatland area start to set in at the 15 onset of the HS1 and increase into the BA. During the BA, peatlands show the fastest gross positive and negative area changes throughout the simulation (see Fig. 4 (b)). In the northern mid and high latitudes, peatlands shift north and eastward (see Fig.  and in cold continental regions of Asia. These new peatlands include the large peat complex in the Western Siberian Lowlands (WSL). Some of the peatlands established in northern Europe during HS1 vanish again during the BA. These changes are driven by the Trace21k climate which shows a substantial warming and wettening of the Northern Hemisphere already starting during the HS1. Temperature is the dominant driver for peat loss and expansion in Europe and North America. In Northern Asia both temperature and precipitation drive the peatland expansion. This expansion sees a pronounced halt during the YD 5 where northern hemisphere climate is briefly returning to more glacial conditions (see Fig. A2 (e)).

Transient peat evolution
In the tropics, the area and carbon changes are mostly driven by precipitation changes. Largest changes are simulated in South America, where precipitation patterns respond to changes in ITCZ position (see Fig. 6). During the HS1 and the YD where the Atlantic meridional overturning circulation (AMOC) is in a reduced state, peatland area shifts to the south following the southward shift of the ITCZ. During the BA the AMOC is strong and precipitation and peatland area shift back north. In 10 Africa half the peatland area is lost during the BA mostly driven by drying. In South East Asia peatlands are lost over the whole termination due to precipitation changes and the onset of sea level rise, which starts to flood the large continental shelves at about 16 kyr BP.

The shifts in peatland distribution result in a similar global peatland area at the beginning of the Holocene compared to the
LGM and at the onset of the HS1. However, much less carbon is stored in active peatlands at the beginning of the Holocene 15 than during the LGM and at the onset of the HS1. Thus, the carbon density per unit area is much lower for the newly established peatlands than for the lost LGM peatlands.

11.65 kyr BP -0 kyr BP
Modelled peatlands in the Holocene show a continuous net expansion in the northern extratropics, with new forming peatlands more than balancing the loss of peatlands elsewhere. The Holocene experienced relative stability in climate and CO 2 levels 20 compared to the termination. The early to mid Holocene was likely characterized by warmer summer temperatures than preindustrial with a larger seasonality in the northern hemisphere (Marcott et al., 2013;Liu et al., 2014;Samartin et al., 2017). Ice sheet retreat and sea level rise lagged behind the deglacial temperature increase and was mostly completed at about 7 kyr BP (Peltier, 2004). Locally new land keeps emerging to this day due to isostatic rebound. This effect is especially pronounced in the Hudson Bay Lowlands, where new land emerges with a rate of up to 12 mm yr −1 (Henton et al., 2006). 25 For northern peatlands, positive area changes are consistently larger than negative changes throughout most of the Holocene (see Fig. 4 (a)). This leads to a large continuous area expansion. Old peatland area is also simulated to increase continuously during the Holocene (Fig. 4 (b)), showing that the parallel positive and negative changes are more than mere fluctuations of existing peat but that there is actual continuous peatland loss and growth. Net area increase picks up at about 9 kyr BP, decreases in the late Holocene and turns into a net area reduction in the last millennium. This late Holocene slow down and 30 reversal of net peatland area growth is most pronounced in Northern Asia where increasing negative changes start to balance and eventually offset the still large positive changes ( Fig. A1 and A2). Both negative and positive dynamics here are driven by temperature and precipitation. The early fast expansion in Northern Asia is offset by a temperature driven net area loss in Europe which is recovered partly towards the late Holocene. Net area increase in North America is delayed by continued loss of mid latitude peat and the slow retreat of the Laurentide Ice Sheet, which limits the establishment of new peatlands. Today's peatlands in North America start to establish after about 9 kyr BP with most of today's peatlands forming between 7 and 2 kyr BP. Carbon stocks follow these regional trends but with larger relative increases especially towards the late Holocene. As the timescale for building up carbon pools is generally much longer than the timescale of potential area changes, fluctuations in area, mostly by young peatlands, are smoothed in the carbon stocks (see e.g. Fig. A1 (b)). 5 The tropical peatland area is simulated to stay relatively stable throughout the Holocene, with positive changes balancing negative changes. South East Asia sees a reduction in area in the early Holocene due to continued sea level rise and a subsequent gradual recovery of integrated peat area driven by precipitation and non linear effects. Peat area in South America increases slightly over the Holocene with mostly precipitation driven fluctuations in between. On the other hand, fluctuations in regionintegrated peat carbon stocks are largely absent in South America, as carbon, with changing area, is shifted between peat and old peat pools. Peat carbon stocks in South America show a large relative increase following a near linear path. Africa sees an increase in area at about 10 kyr BP driven by precipitation and enabled by high CO 2 concentrations. The new area gradually 5 degrades again until 3 kyr BP with another peak at 0.5 kyr BP.

Model versus reconstructions
The study of peatland initiation, life cycle, dynamics and responses to external forcing has been focused on today's existing and active peatlands. This work includes large compilations of peat core basal dates (MacDonald et al., 2006;Gorham et al., 2007;Yu et al., 2010) which are used to reconstruct initiation dates and lateral expansion (Yu et al., 2010;Korhola et al., 2010;10 Dommain et al., 2014;Loisel et al., 2017) of the sampled active peatlands. This approach, however, does not include earlier peatlands that dried out, were buried or flooded or otherwise seized to be active accumulating peatlands. Treat et al. (2019) presented a first compilation of dated buried peat layers, but the small sample size make quantitative reconstructions difficult.
We thus limit most of the model-data comparison of the transient behaviour to the today's existing peatlands. Figures 7 and 8 (c) show modeled initiation date frequency, area and carbon dynamics of northern peatlands that are still active at present. The early expansion in the model also propagates to the carbon balance for presently active peatlands. The model simulates 25 earlier accumulation extending into the HS1 and slower accumulation during the early Holocene than suggested by net carbon balance (NCB) reconstructions by Yu (2011) (Fig. 7 (c)). The summed simulated carbon increase from LGM to PI in today's northern peatlands amounts to 343 GtC (Fig. 8 (c)). Throughout the tropics dated buried and active peat cores show peatland presence already during and preceding the LGM.
But peatland extent or evolution towards the presence are not well constrained and are subject to large uncertainties. The 10 small number of available dated tropical peat cores impedes a statistical approach. Applied nevertheless, it indicates a more or less continuous growth of today's peatlands since about 19 kyr BP with largest expansion rates between 8 -4 kyr BP (Yu et al., 2010). Dommain et al. (2014) reconstructed the evolution of Indonesian peatlands using a combination of dated cores and a transfer function between depth and age. They argue for a peat expansion much later than infered by basal ages alone, with 90% of todays peat establishing after 7 kyr BP and 60% after 3 kyr BP. In this study the dominant control on 15 peatland area was found to be local sea level. Rising sea level during the termination and the early Holocene triggering the establishment of inland peatlands through alterations in moisture availability and the hydrological gradient, and stabilization of sea level and subsequent sea level regression after 4 kyr BP prompting the establishment of coastal peatlands. In contrast to these reconstructions, the transient simulation shows 60% of today's tropical peatland area already present in the LGM and only small expansion during the last millennia. The sparsity of the data warrant's caution when comparing to model results.

20
However, in South East Asia, this discrepancy could indicate the importance of the feedback of sea level on local hydrology, missing in LPX.

Transient carbon balance of peatland soils and the land biosphere as seen by the atmosphere
In this section, we address how carbon stored in soils of active peatlands and carbon stored in the remains of former peat soils changed over time. Thereby, we quantify the overall contribution of peatland soils and peat carbon to the changes in the global 25 carbon inventory of the land biosphere.
When trying to quantify the net effect of peatlands on the atmosphere, looking only at carbon stored in today's active peatlands can be misleading. Former active peatlands have transformed into other landscapes. Organic rich peat layers may now be burried under mineral soils on land or in coastal ocean sediments (Treat et al., 2019;Kreuzburg et al., 2018). When analyzing the transient carbon balance of global peatlands such "old peat carbon" pools have to be considered. simulated increase in peat carbon from LGM to PI within these three pools is 350 GtC. This represents the simulated net carbon accumulation of global peat and thus the net amount of carbon sequestered from the atmosphere by peat. When only considering carbon stored in active peatlands we would underestimate the deglacial peat carbon change with 223 GtC. On the other hand if we only consider the carbon stored in today's active peatlands, the inferred deglacial change amounts to 365 GtC (Figure 8 (c)), and we would overestimate the net peat accumulation since the LGM. While the latter difference in net peat 5 carbon accumulation between the complete and incomplete accounting scheme appears small for the total deglacial change, the difference can be substantial and relevant for other periods. For example, a particular large difference is identified for the phase of high peatland expansion and loss rates as simulated from the Bølling/Allerød to the Preboreal in our model. Here the carbon balance is given by the complete accounting, including old peat, as 12 GtC versus 102 GtC when only looking at the carbon in today's active peat for the period from 14.6 to 10 kyr BP.  LGM climatic conditions, LPX-Bern predicts peatlands in areas with low or no peat cover at present, or on now submerged continental shelves. Uncertainty from the climate forcing was assessed by using, in addition to the TraCE21k, climate anomalies from six different time slice simulations for the LGM from phase 3 of the Paleomodel Intercomparison Project PMIP3. This results in a peat area range of 1. 5-3.4 Mkm 2 with a carbon storage of 144-343 GtC. This large range illustrates the dependence 10 of results on, uncertain, LGM climate conditions and the sensitivity of simulated peatlands to these differences. Sparse data on paleo peatlands, on their extent, and their carbon storage make it difficult to further constrain this range. At the same time there 21 https://doi.org/10.5194/bg-2020-110 Preprint. Discussion started: 26 April 2020 c Author(s) 2020. CC BY 4.0 License. are, to date, only a few coupled climate simulations for the LGM and only one transient simulation with an atmosphere-ocean general circulation model available for the period LGM to present.
A driver attribution of the simulated transient evolution of peatlands using factorial simulations showed regional and temporal differences. Modelled changes in the tropics were dominated by shifts in the position of the Inter Tropical Convergence Zone and associated precipitation changes during the last glacial termination as well as by rising sea level. Changes in the 5 northern high latitudes are mostly driven by temperature and precipitation increases. Largest model mismatches to available area reconstructions can be seen in the onset and timing of the earliest expansion of today's northern peatlands. A strong warming in the climate forcing during Heinrich Stadial 1 and the Bølling/Allerød triggers a first expansion into Northern Asia, which according to reconstructions only starts during the Preboreal, about 4 kyr later.
The simulated transient evolution of peatlands is characterized by continuous and simultaneous increases and decreases 10 of area and carbon, with fastest positive and negative changes happening during the termination (Heinrich Stadial 1 and Bølling/Allerød). This reveals a different perspective than the commonly assumed linear and continuous growth of global peatlands. Instead peatlands become a dynamic, growing, dying and shifting landscape. Carbon in soils of formerly active peatlands can be trapped in mineral soils or ocean sediments. When assessing the net carbon balance of global peatlands over time, accounting for paleo peatlands and their remains thus becomes essential. In our transient simulation the LGM to PI net 15 peat carbon balance is predicted at 350 GtC, with 499 GtC, 139 GtC, and 22 GtC stored at pre-industrial in soils of still active peatlands, in the remains of former peat soils on land, and in the remains on submerged shelves, respectively. For today's active northern peatlands, simulated peat area and carbon is in good accordance with the range of literature estimates, whereas predictions for the tropics are larger than most estimates. However, data constraints in the tropics are significantly weaker as peat sciences has long focused on the northern high latitudes and only in the last decades is accelerating its effort in the tropics.

20
Even fewer data is available to constrain old peat carbon that remains outside of today's active peatlands.
Taken together our study provides an in depth model analysis of peatland development, the associated drivers, and uncertainties on a global scale. It contributes to a foundation for a better understanding of past peat dynamics and emphasises the importance of treating and understanding peatlands as dynamic and evolving systems.   Figure A1. Simulated global and regional peatland area and carbon dynamics over time, relative to PI levels. PI levels are given in Mkm 2 for peat area and GtC for peat carbon. The extent of the regions Northern Asia, Western Siberian (WS) Lowlands, and Southeast Asia are shown in Fig. 1.