Global peatland area and carbon dynamics from the Last Glacial Maximum to the present – a process-based model investigation

Peatlands are an essential part of the terrestrial carbon cycle and the climate system. Understanding their history is key to understanding future and past land–atmosphere carbon fluxes. We performed transient simulations over the last 22 000 years with a dynamic global peat and vegetation model forced by Earth system model climate output, thereby complementing data-based reconstructions for peatlands. Our novel results demonstrate a highly dynamic evolution with concomitant gains and losses of active peatland areas. Modeled 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 they are partly lost in midlatitude regions. In the tropics, peatlands are partly lost due to the flooding of continental shelves and are regained through nonlinear responses to the combined changes in temperature, precipitation, and CO₂. Large north–south shifts of tropical peatlands are driven by shifts in the position of the intertropical convergence zone associated with the abrupt climate events of the glacial termination. Time slice simulations for the Last Glacial Maximum (LGM) demonstrate large uncertainties in modeled peatland extent (global range from 1.5 to 3.4 Mkm², million square kilometers) stemming from uncertainties in climate forcing. The net uptake of atmospheric CO₂ by peatlands, modeled at 351 GtC since the LGM, considers decay from former peatlands. Carbon uptake would be misestimated, in particular during periods of rapid climate change and subsequent shifts in peatland distribution, when considering only changes in the area of currently active peatlands. Our study highlights the dynamic nature of peatland distribution and calls for an improved understanding of former peatlands to better constrain peat carbon sources and sinks.

1 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 multiple forms from minerotrophic fens to ombrotrophic bogs and forested 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 they store a share of the total global soil organic carbon that is up to an order of magnitude higher (Page et al., 2011; Yu, 2012). At the same time, they are a significant carbon sink (e.g., Gorham et al., 2012; Lähteenoja et al., 2012; Leifeld et al., 2019) and a large natural source of methane (e.g., Frolking and Roulet, 2007; LAI, 2009; Korhola et al., 2010; Yu et al., 2013; Packalen et al., 2014); thus, they are an integral part of the terrestrial carbon cycle (Gorham, 1991; Yu, 2011; Page et al., 2011). Most of today's peatlands, formed over the past 12 000 years as a result of deglacial climate change and ice sheet retreat (e.g., Halsey et al., 2000; Gajewski et al., 2001; MacDonald et al., 2006; Gorham et al., 2007; Yu et al., 2010; Ruppel et al., 2013; Morris et al., 2018; Treat et al., 2019). Since then, northern peatlands alone have sequestered about 500 GtC (Yu et al., 2010; Yu, 2012), resulting in a net cooling effect on the climate (Frolking and Roulet, 2007). The drainage and conversion of existing peatlands to plantations or other forms of land use leads to release of peat carbon into the atmosphere, adding to the ongoing global warming trend (Dommain et al., 2018; Leifeld et al., 2019). Additionally, global warming will likely diminish the net carbon sink of remaining global peatlands (Spahni et al., 2013; Gallego-Sala et al., 2018; Wang et al., 2018; Leifeld et al., 2019; Ferretto et al., 2019), despite a possible increase in the sink of some northern peatlands (Swindles et al., 2015; Chaudhary et al., 2017a).

Despite their global importance, peat research has long focused almost exclusively on northern high-latitude peatlands with about 80 % of dated peat cores taken in Europe and North America, which only covers roughly 40 % of global peat area (Xu et al., 2018; Treat et al., 2019). Although research on tropical peatlands has increased in recent years (e.g., Page et al., 2011; Dommain et al., 2011, 2014; Lawson et al., 2015; Silvestri et al., 2019; Gumbricht et al., 2017; Cobb and Harvey, 2019; Leng et al., 2019; Illés et al., 2019), our understanding about tropical peatlands, including their dynamics and life cycles, is still limited. This also entails ongoing new discoveries of previously unknown peatland complexes such as in the Congo Basin (Dargie et al., 2017). The tendency to search for the deepest core within a peatland (Loisel et al., 2017) and the acute lack of information about the fate of old and buried peat (Treat et al., 2019) represent additional sampling biases that contribute to our limited understanding of peatland evolution and its drivers. These gaps in our understanding are also reflected in the large ranges of estimates of today's peatland area (e.g., Yu et al., 2010; Page et al., 2011; Loisel et al., 2017; Xu et al., 2018) and peatland carbon (e.g., Tarnocai et al., 2009; Yu et al., 2010; Yu, 2012; Page et al., 2011; Gumbricht et al., 2017) in the literature. Only recently, a highly contested study proposed a doubling of conventional northern high-latitude peat carbon stock estimates (Nichols and Peteet, 2019; Yu, 2019). Refining our understanding and estimates of peatland carbon dynamics is timely, as the potential past and future effects of peatlands on the global carbon cycle are substantial, and knowledge of the amount, timing, and speed of carbon removal and release is crucial to constrain them.

Results from process-based models can offer an independent perspective on the transient evolution of global peatlands and peat carbon stocks, complementing data-based reconstructions of global peatland expansion and carbon accumulation (e.g., MacDonald et al., 2006; Gorham et al., 2007; Yu et al., 2010; Ruppel et al., 2013; Dommain et al., 2014; Loisel et al., 2017; Treat et al., 2019). Efforts to model peatlands and processes within them exist at the site level (e.g., Frolking et al., 2010; Morris et al., 2011; Baird et al., 2012; Morris et al., 2012; Kurnianto et al., 2015; Cresto Aleina et al., 2015; Chaudhary et al., 2017b; Cobb and Harvey, 2019) as well as on regional to global scales (e.g., Wania et al., 2009a, b; Kleinen et al., 2012; Spahni et al., 2013; Gallego-Sala et al., 2016; Alexandrov et al., 2016; Chaudhary et al., 2017a; Stocker et al., 2017; Largeron et al., 2018; Qiu et al., 2019; Swinnen et al., 2019). Although still small, the number of dynamic global vegetation models (DGVMs) with integrated peatland modules and dynamic peatland area is increasing (Kleinen et al., 2012; Stocker et al., 2014; Largeron et al., 2018; Qiu et al., 2018) enabling, for the first time, a hindcast of past and a prediction of future peatlands on large spatial and temporal scales. Representations of peatlands have also been developed for the inclusion in the land modules of complex Earth system models (Lawrence and Slater, 2008; Schuldt et al., 2013). However, peatlands are generally still prominently missing from the newest generation of Earth system models (ESMs) taking part in the sixth phase of the Climate Model Intercomparison Project (CMIP6), which is the main source for future climate and carbon cycle projections used for the determination of international climate mitigation targets (Eyring et al., 2016). Thus, rigorous testing and improvement of the existing peat modules not only has the potential to yield further insights into peatland dynamics but can also pave the way for the integration of peat into the next generation of ESMs for improved climate projections.

Peatlands and their carbon stocks evolve dynamically through time and over glacial cycles. Peatlands may disintegrate or be buried by mineral sediments when climatic conditions become locally unfavorable for peat growth or local hydrologic conditions change (e.g., Talbot et al., 2010; Tchilinguirian et al., 2014; Campos et al., 2016; Lähteenoja et al., 2012; Tipping, 1995). Peatlands on exposed coastal shelves may be flooded during periods of rising sea levels (Kreuzburg et al., 2018), and new peatlands may form in areas previously covered by continental ice sheets or in areas that were previously too cold or too dry for peat establishment. Therefore, net changes in peat extent are the difference of concomitant gains and losses in peatland area. Similarly, the net flux of CO₂ from the atmosphere to peat carbon is the sum of complex changes. Peat carbon accumulates on active and expanding peatlands. Dying peatlands may lose some of the accumulated carbon to the atmosphere through degradation whereas another part might be buried and, thus, conserved on long timescales. Estimating peat carbon stocks for today's active peatlands is an important but insufficient step to fully constrain the influence of peat carbon changes on the atmospheric carbon balance. At the same time, peatlands are, in the absence of abrupt anthropogenic disturbance, slowly reacting systems with process timescales ranging from years to millennia. Thus, the present distribution of peatland and peat carbon and their future fate depend on past peatland dynamics and legacy effects from the last glacial–interglacial climate transition as well as the current interglacial (Holocene). However, model studies that thoroughly investigate the establishment and the disintegration of global peatlands constraining the total carbon balance transiently over the deglaciation are still lacking.

Here, our goal is to present a rigorous model investigation of peatland area and carbon dynamics since the Last Glacial Maximum (LGM), 21 000 years before present (BP). We use a DGVM to simulate the LGM peatland distribution and assess uncertainties stemming from the climate forcing. Transient model and factorial simulations from the LGM to the present are analyzed to learn about past peatland dynamics, underlying drivers, and the net peatland carbon balance. Model results are compared to available data for the present and the LGM as well as to reconstructions of modern-day peatland initiation and development.

2 Methods

2.1 Model description

The simulations presented here were performed with the Land surface Processes and eXchanges (LPX-Bern) dynamic global vegetation model (DGVM) version 1.4 (Lienert and Joos, 2018). It includes an interactive carbon, water, and nitrogen cycle and simulates dynamic vegetation composition with plant functional types (PFTs), which compete for water, light, and nutrients (Sitch et al., 2003; Xu-Ri et al., 2012; Spahni et al., 2013). The implementation of permafrost and peatlands as long-term carbon stores are based on the Lund–Potsdam–Jena Wetland Hydrology and Methane (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 “peat”, “mineral soils”, and “old peat” (formerly active peat now treated as mineral soils) land classes. In this study, anthropogenic land use and land use change and corresponding land classes are not considered. Estimates of peat carbon loss through land use only exist for the industrial period. Leifeld et al. (2019) estimate that about 22 ± 5 GtC of peat carbon was lost globally between 1850 and 2015. Houghton and Nassikas (2017) used results from Randerson et al. (2015) and Hooijer et al. (2010) to estimate carbon loss from the draining and burning of peatlands for oil palm plantations in Southeast Asia. Losses are negligible before 1980, and they amount to about 6 GtC for the period from 1980 to 2015. These estimates of the loss of peat carbon through land use, although substantial, are still small compared with the total pool sizes. Peatland vegetation is represented by five peat PFTs: Sphagnum and flood tolerant graminoids, as mostly indicative of high-latitude peatlands, and flood tolerant tropical evergreen and deciduous tree PFTs as well as a flood tolerant version of the C4 grass PFT, as mostly indicative of 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 a fluctuating water table position (WTP; Spahni et al., 2013). For the determination of the carbon flux between acrotelm and catotelm a fixed average acrotelm carbon (C) density (18.7 kg C m⁻³) is used, along with a fixed acrotelm depth of 0.3 m. The difference between this target acrotelm density and the actual average acrotelm density, determined by carbon influx from the litter pools and heterotrophic respiration within the acrotelm, is used to determine the size and sign of the daily flux between acrotelm and catotelm (Spahni et al., 2013). Decay rates are modulated by temperature in the catotelm and by temperature and the WTP in the acrotelm (Wania et al., 2009a).

The area fraction covered by peat (f_peat) in a given grid cell is determined dynamically with the DYPTOP (Dynamical Peatland Model Based on TOPMODEL) module (Stocker et al., 2014). The TOPMODEL approach (Beven and Kirkby, 1979) is used to predict the monthly inundated area fraction given sub-grid-scale topographic information and mean grid cell WTP. Here, the WTP calculation of mineral soils has changed slightly with respect to Stocker et al. (2014), with drainage runoff excluded from the calculation. The area potentially available for peatlands (f_pot) is then determined by inundation persistency. Peatlands expand or shrink towards a changing f_pot at a rate of 1 % of the current f_peat per year. The grid cell fraction lost during peatland retreat is treated as a separate land use class named “old peat”; it inherits the carbon stocks of the dying peat and is subsequently treated in the same way as the mineral soils regarding vegetation, hydrology, and carbon cycling. Growing active peatlands first expand on eventual old peat, inheriting the remaining carbon there.

As vegetation growth and carbon cycling continues normally on the old peat fraction, the carbon inherited by the former peatland, which would form distinct organic soil layers in the real world, can not be distinguished from new carbon accumulated by new non-peatland vegetation in the model. The same is true for grid cells that get flooded by rising sea levels. Given the evidence of coastal peat carbon deposits (Kreuzburg et al., 2018; Treat et al., 2019), we assume that most of the carbon is buried within sediments rather than released to the atmosphere during flooding. In case of flooding in the model, carbon from all land use classes in the respective cell is combined into a single “flooded” land use class, where it slowly decays with a mean lifetime of 15 kyr. Despite this mixing of carbon from different sources, we can track peat carbon in post-processing using the transient model output for peatland area changes, decay rates of slow pool carbon, and carbon input into the catotelm of the active peatlands. Area changes are used to transfer carbon between active, old, and flooded peatlands. Transient decay rates are used to decay the carbon in the respective pools. Thus, carbon is tracked from the entry into the catotelm of an active peatland until decay in either an active peatland, an old peatland, or peatland flooded by ocean. This approach can not take the acrotelm carbon into account; however acrotelm carbon constitutes only a small part of total peat carbon (5 % in the preindustrial (PI) period, and we can assume that this carbon at the peat surface is quickly respired after peatland transformation. The “old peat carbon” calculated this way represents the remaining peatland carbon after peatland death and is used in the calculation of the peatland carbon balance (see Sect. 3.3.5).

Peatland existence, beyond a small peatland “seed” ($f_{\mathrm{peat}}={\mathrm{10}}^{-\mathrm{5}}$) in every grid cell, is further limited by criteria on its carbon (C) and water balance. In this study, the evapotranspiration for peatland tree PFTs is now calculated analogously to non-peatland tree PFTs using demand and supply functions (Sitch et al., 2003). The determination of the criterion of a positive water balance (precipitation over evapotranspiration ratio in peat >1), however, was kept functionally unchanged to Stocker et al. (2014). The C criteria were slightly improved in this study. The peat establishment and persistence criterion on the C balance during the spin-up is a positive net ecosystem production (NEP) and an acrotelm to catotelm flux higher than 10 $\mathrm{g}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$, or C stocks of the peat seed exceeding 50 kg m⁻² as in Stocker et al. (2014). During the transient run this criterion is changed so that peat establishment depends on the acrotelm to catotelm flux alone. For peat persistence, the sharp C stock threshold is softened. For a peat C stock from 50 kg m⁻² to about 45 kg m⁻², f_pot is reduced to an actual potential peatland fraction (f_apot) according to a sigmoid function:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{f}_{\mathrm{apot}}& =f_{\mathrm{pot}}\times {\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\mathrm{20}{e}^{-\mathrm{2.4}\left({\mathrm{C}}_{\mathrm{peat}}-\mathrm{46}\right)}}}|\mathrm{if}{\mathrm{C}}_{\mathrm{peat}}\\ & <\mathrm{50}{\mathrm{kg}\mathrm{m}}^{-\mathrm{2}},\end{array}\end{array}$$

where C_peat represents the peatland soil carbon pool in units of kilograms per square meter (kg m⁻²). This avoids peatland collapse due to a sharp threshold. Peatlands can now endure short periods of carbon loss, even with C pools falling below the threshold of 50 kg m⁻², but they have to suffer area losses as a consequence, as f_peat now approaches f_apot.

The representation of peatlands in the LPX described above is a simplification in many respects. The absence of local processes and information like lateral water flow, the influence of sea level variations on the water balance, local soil features, or the influence of animals by grazing and river damming can limit the ability of the TOPMODEL approach to predict peatlands on a regional to local scale. Further, direct anthropogenic influences such as land use, drainage, or peat mining are not considered. The lack of a distinction and transition between different peatland types, like fens, bogs, blanket bogs, or marshes, neglects possible differences in the constraints on their formation and evolution. The treatment of acrotelm and catotelm as single carbon pools and the absence of strong disturbances such as peat fires constitute limits on the comparability of the model results to peat core carbon profiles. This simplified representation, however, has been shown to reproduce peatland area and carbon accumulation well within the observational constraints (Wania et al., 2009a; Spahni et al., 2013; Stocker et al., 2014, 2017) while using a minimal set of free parameters. Our efficient representation allows for long transient paleo-simulations and sensitivity studies such as those we present here.

2.2 Simulation setup

The transient LPX simulations from the Last Glacial Maximum (LGM), 22 kyr before present (BP), until present were run with a model resolution of 2.5^∘ latitude × 3.75^∘ longitude and were forced with CO₂ (Joos and Spahni, 2008), temperature, and precipitation fields as well as transient evolving orbital parameters influencing available photosynthetic active radiation. Temperature and precipitation anomalies are taken from the transient CCSM3 run TraCE21k (Liu et al., 2009). The TraCE21k experiment constitutes a unique climate forcing, not only because it is currently the only published transient simulation over the deglaciation using a fully coupled general circulation model (GCM) but also because the meltwater forcing in TraCE21k was chosen, using sensitivity experiments, to best reproduce the abrupt climate events such as the Bølling–Allerød (BA) and the Younger Dryas (YD) (He, 2011). The TraCE21k anomalies are imposed on the CRU TS 3.1 (Mitchell and Jones, 2005) base climate from 1960 to 1990. Thus, interannual variability is adopted from TraCE21k. The land–sea–ice mask is changes 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.

In addition to the standard LPX transient simulation, five transient factorial simulations were performed using the same setup but keeping one of the five transient forcings (land–sea–ice mask, orbital, CO₂, precipitation, and temperature) constant at LGM levels for each simulation. These were used to identify the dominant drivers and driver contributions through time by comparing the factorial and standard runs (see Sect. 3.3).

To investigate the uncertainty stemming from the choice of climate forcing, seven additional LGM time slice simulations were performed. Mean LGM climate anomalies from six different PMIP3 models (including CCSM4, COSMOS-ASO, IPSL-CM5A-LR, MIROC-ESM, MPI-ESM-P, and MRI-CGCM3; Braconnot et al., 2012) and the mean LGM anomaly of the TraCE21k simulation were imposed on the CRU 3.1 climatology from 1901 to 1931. Thus, interannual variability is adopted from CRU. CO₂, ice-sea-land mask, and orbital parameters are set to LGM levels (in this case 21 kyr BP). Two of the eight available LGM simulations from phase 3 of the Paleoclimate Modelling Intercomparison Project (PMIP3) were not used (FGOALS-g2 and CNRM-CM5) due to their 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 for an additional 2 kyr under unchanged conditions. In the analysis, the temporal mean over the last 2 kyr is used.

2.3 Validation data

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 the tropical forests or the Arctic tundra. Even estimates for the relatively well-studied northern high-latitude peatlands have a range from 2.4 to 4.0 M km² (see Loisel et al., 2017 for a review). The total area of tropical peatlands is even less well defined and estimates range from 0.37 to 1.7 M km² (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 modeling, satellite imaging, and topographic data and, thus, tries to also account for currently undiscovered peatlands (Gumbricht et al., 2017). The most extensive and comprehensive compilation of known peatlands to date is the recent PEATMAP by Xu et al. (2018). PEATMAP shows a distribution shifted more towards the tropics than previous estimates from the literature. For example, Yu et al. (2010) estimates the area of northern peatlands to be 4 M km² and the area of tropical peatlands to be 0.37 M km², whereas PEATMAP gives 3.18 and 0.99 M km², respectively. In Fig. 1, Table 1, and Sect. 3.1, the LPX present-day peatland extent 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 $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ grid cell (MGK13G), as a proxy for peatland initiation. The MGK13G dataset is used in this study for comparison with 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 grid cells 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 utilizing 954 basal dates from 138 sites, with at least three dated cores per site. Their expansion reconstruction (KOR10) shows delayed expansion compared with Loisel et al. (2017) and the fastest expansion between 3 and 5 kyr BP. Both MGK13S and KOR10 are compared to the expansion simulated by LPX for currently existing northern peatlands (see Sect. 3.3.4), thereby not including area changes of peatlands that had previously existed but have already disappeared.

.

Figure 1 Global present-day peatland distribution according to PEATMAP (Xu et al., 2018) in a 0.5^∘ × 0.5^∘ gridded version (a) and simulated by LPX-Bern after the transient “standard” setup simulation from 22 kyr BP to present (b). The colored rectangles show three of the regions listed in Table 1: northern Asia (red), the western Siberian lowland (orange) region, and Southeast Asia (green)

3 Results and discussion

3.1 Distribution and carbon inventories of present-day peatlands

3.1.1 Peatland area

The modern peatland distribution simulated by LPX-Bern (standard run) compares well to the distribution given by PEATMAP (see Fig. 1 and Table 1). LPX and PEATMAP yield a very similar global peatland area, with respective values of 4.37 and 4.23 M km². The same is true for the latitudinal distribution. LPX simulates 3.2 M km² in the high latitudes (>30^∘ N) and 1.15 M km² in the tropics (30^∘ S–30^∘ N), whereas PEATMAP gives 3.18 and 0.99 M km², respectively. This broad-scale agreement between LPX and PEATMAP notably emerges without any tuning of LPX against PEATMAP. These results are similar to previous results using the LPX (Stocker et al., 2014) with slightly larger tropical peatland coverage in the current study. Differences are due to a new model version after data assimilation, LPX v1.4, (Lienert and Joos, 2018) and the additional model changes described in Sect. 2.1.

Minor and major differences in peatland area between LPX and PEATMAP are seen on the local to regional scale (Fig. 1 and Table 1). In the tropics, LPX simulates more peat in South America and Southeast Asia than PEATMAP indicates. Compared to the estimate of Gumbricht et al. (2017), the LPX peatland extent is similar for South America and a factor of 2 smaller for Southeast Asia. The vast peatland complex in the Congo Basin is almost absent in LPX. In the northern mid to high latitudes, LPX seems to underestimate European peatland area by a factor of 2 and slightly overestimates peatland area in northern Asia, mostly west and east of the western Siberian lowland (WSL) peat complex. In North America, LPX simulates more peat in Alaska and Quebec and less in Western Canada than PEATMAP.

Other modeling studies present results from prognostic simulations of Northern Hemisphere peatlands. Kleinen et al. (2012) simulated peatland dynamics and carbon accumulation over the past 8000 years using the coupled climate carbon cycle model CLIMBER2-LPJ. However, no quantitative results in terms of peatland area were reported. Qiu et al. (2019) used the ORCHIDEE-PEAT DGVM (Qiu et al., 2018) to simulate northern (>30^∘ N) peat expansion over the Holocene. Their simulated northern present-day peatland area is, at 3.9 M km², slightly larger than in LPX. They find similar regional discrepancies between simulated and observation-based peat area in North America, northern Europe, and Asia, as described above for LPX. Peat area dynamics in ORCHIDEE-PEAT are also using the TOPMODEL approach following Stocker et al. (2014), with some different expansion criteria. This might indicate that these discrepancies could have their source in the TOPMODEL approach and its limitations. Another major source of uncertainties is in the climate data used to force LPX (see also Sect. 3.2.2). In particular, precipitation data show large discrepancies between available observational products (Sun et al., 2018).

Table 1 Peatland area for different regions and latitudinal bands as given by PEATMAP (Xu et al., 2018) for today and peatland area and their carbon stocks as simulated by LPX-Bern for the preindustrial (PI) period and the Last Glacial Maximum (LGM) in the transient “standard” setup simulation from 22 kyr BP to present. The extent of the northern Asia, western Siberian (WS) lowland, and Southeast Asia regions are shown in Fig. 1.

3.1.2 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. Therefore, the ranges of both carbon estimates and area estimates are large (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 estimates in the tropics and northern Asia. Estimates for northern peatlands range from 270 to 604 GtC and were obtained using various methods and area estimates (see Yu, 2012 and Yu et al., 2014 for a review). The modern carbon inventory of northern peatlands simulated by LPX at the end of the transient standard run from the LGM until present is, at 361 GtC, well within this observational range. In the tropics, LPX simulates a peat carbon inventory of 136 GtC which is substantially larger than classical estimates from the literature that range from 50 to 87 GtC (Yu et al., 2010; Page et al., 2011). These, however, also assume a substantially smaller tropical peatland area than LPX or PEATMAP suggest (see Sect. 3.1.1). Including estimates for the newly discovered peat in the Congo Basin, Dargie et al. (2017) estimate a tropical peat inventory of 69.6–129.8 GtC, which is closer to the LPX results. 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 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 spin-up scheme, whereas the pools in this study are filled over a transient run. Spahni et al. (2013) also report northern peatland carbon stocks after a transient LPX run from the LGM, although with prescribed not prognostic peatland 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., 2019) in northern peatlands.

3.2 Peatlands during the Last Glacial Maximum

3.2.1 Peatland distribution and carbon storage

Under LGM conditions, global simulated peatland area and carbon inventories are reduced compared with the preindustrial period (Table 1). Globally, simulated peatland area and the peat C inventory are 38 % and 45 % smaller at the LGM than during the PI period, respectively. This reduction is dominated by the northern extra-tropics, where peat extent and the C inventory are almost 60 % smaller at the LGM than during the PI period. In contrast, the peat C inventory in the tropics is only about 3 % smaller, and the tropical peat area is even 7 % larger at the LGM than during the PI period. 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 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.

Simulated peatland coverage in northern mid and high latitudes is smaller and is shifted southwards at the LGM compared with the PI period. Ice sheets covered large parts of Europe and North America during the LGM preventing vegetation and peat from growing. However, peat is also mostly absent in northern Asia and the WSL due to the substantially colder and dryer conditions compared with today. On the other hand, large peatland complexes are simulated along the southern ice sheet margins in North America and in Europe (Fig. 2), in regions where modeled peatlands are mostly absent under current conditions. This even leads to a simulated net increase of peatland area in Europe (+43 %) compared with the present. Verifying the existence of these extensive LGM peatlands that do not exist under present conditions (compare Fig. 1) is difficult, as existing compilations of peat core dates focus almost exclusively on today's existing peatlands (MacDonald et al., 2006; Gorham et al., 2007; Yu et al., 2010; Loisel et al., 2017). In a recent study, Treat et al. (2019) presented a compilation of dated buried peat deposits and simulated peatland area and carbon stocks. Their simulation also suggests large midlatitude peatlands in North America in agreement with our results. Their peat deposits data for the LGM (Fig. 2; dots), together with pollen analyses suggesting the presence of at least some Sphagnum in eastern North America (Halsey et al., 2000; Gajewski et al., 2001), provide plausible evidence for the existence of midlatitude LGM peatlands in North America and Europe. Their extent, however, is probably overestimated in our simulation. Comparisons between the North American LGM hydroclimate in TraCE21k and proxy reconstructions have resulted in a poor skill score, especially in eastern North America (Lora and Ibarra, 2019). Bad performance in this region is shared with all PMIP3 models (Lora and Lora, 2018).

Figure 2 Peatland distribution at the Last Glacial Maximum (LGM) as simulated by LPX-Bern in the standard setup (a), agreement (as number of models simulating peat in given grid cell) between LGM time slice simulations run with LPX and forced with different climates anomalies from six PMIP3 models as well as the TraCE21k anomaly (b), and the squared correlation coefficient for a linear regression between physical properties of the different time slice simulations (precipitation minus evapotranspiration, P−E, and growing degree days above 0 ^∘C) and peat fraction in the respective cells. Only cells with significant correlation (p>0.05) are plotted. The colored shading in panel (b) indicates how many time slice simulations show a peat fraction of >0.05 in the respective cell. The color code in panel (c) denotes the dominant predictor in the respective cell. Dots in panels (a) and (b) show buried and still active peat deposits that indicate active peat accumulation during the LGM (24.5–17.5 kyr BP). Peat core data are from Treat et al. (2019).

3.2.2 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 parameterizations, 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) to constrain model parameters, is an approach to improve model performance in the light of uncertain key parameters. Another source of uncertainty stems from uncertainties in the prescribed forcings. The orbital parameters, atmospheric CO₂ mixing ratio, and land–sea–ice mask for the LGM and their deglacial evolution are all well constrained for the purpose of peat modeling, in contrast to the climate anomalies. Although there are paleoclimate reconstructions for the LGM (Bartlein et al., 2011; Schmittner et al., 2011; Annan and Hargreaves, 2013), 6 k (Bartlein et al., 2011) and the last millennium (Hakim et al., 2016; Tardif et al., 2019), they lack the temporal resolution and/or spatial coverage needed for a global transient simulation from the LGM to present. Climate models can fill these gaps; however, climate anomalies are model dependent, and model performance differs between variables, regions, and the simulated time period (Harrison et al., 2014). These differences in climate models have been shown to propagate large uncertainties into carbon cycle projections (Stocker et al., 2013; 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 M km²), peat area (1.5–3.4 M km²), and peat carbon (147–347 GtC). Interestingly, simulated wetland and peatland area and the peat C inventory for the 21 kyr BP period are also substantially different between the standard transient simulation using temporally evolving climate anomalies from the TraCE21k simulation compared with the time slice simulation forced with TraCE21k anomalies (Fig. 3, gray star versus gray dot). This highlights both the influence of different methods of input preparation, with slightly different treatment of anomalies and an interannual variability taken from TraCE21k in the transient simulation and from CRU 3.1 for the time slice, as well as the importance of memory effects for a slowly reacting system such as peatlands.

Agreement on the simulated peat extent among the seven simulations differs among regions (Fig. 2b). It tends to be higher in the tropics and East Asia and lower in North America, Europe, and western Siberia. Differences in temperature and precipitation anomalies propagate into differences in the water balance and productivity, partly limited by growing season length, and thus 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. 2c). Temperature, precipitation, precipitation minus evapotranspiration (P−E), and growing degree days over 0 ^∘C (GDD₀) are considered as climatic predictor variables for the peat fraction within a grid cell. We correlated, for each grid cell, the seven climatological mean values of a selected predictor with the modeled peat fraction from each of the seven time slice simulations. P−E and GDD₀ show significant correlations (p<0.05) in more grid cells than precipitation and temperature, respectively. Both moisture balance and GDD₀ have been shown to be among the most important predictors of northern peat initiation and carbon accumulation in the past (Morris et al., 2018; 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 differences in local precipitation anomalies driving P−E. The same 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 largest tropical peatlands. However, in parts of central South America, temperature is the dominant predictor signaling a fragile carbon balance, where peat presence in some of the model climates 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 differences. The time slice forced with TraCE21k climate shows the peatland distribution in North America more shifted to the east compared with most other PIMP3 forcings alongside warmer and wetter conditions (see also Lora and Lora, 2018). The 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 time slice, temperature anomalies with respect to the preindustrial period are lowest and peat is subsequently shifted northwards compared with other time slices. Similar is true for northern and East Asia where lower temperature anomalies allow for more peatlands. Large areas in central Europe, East Asia, and South America show differences in peatland extent induced by differences in climate forcings, but no significant correlations between peat fraction and predictor variables are found. This might be the result of nonlinear interactions and threshold behaviors not captured by our linear regression approach. Taken together, these findings demonstrate a strong sensitivity of simulated peat extent and the C inventory to the prescribed climate fields and a strong dependence of the results on the choice of climate model output used to force LPX. In other words, caution is warranted when interpreting model results for times and regions in which proxy records or observations are sparse and have limited power to constrain the actual climate conditions. This holds not only in the context of this study but also for global peat and carbon cycle model studies in general.

Figure 3 Spread of the global annual mean inundated area (meanIA), peatland area, and peatland carbon during the LGM (21 kyr BP) in time slices forced with climate anomalies from seven different climate models (see Sect. 2.2). The box, the solid line, and the dashed line indicate the interquartile range, the median, and the mean from the seven simulations, respectively. The star indicates the values of the transient simulation for the same time period using climate anomalies from TraCE21k.

3.3 Transient peat evolution

Figure 4a shows the peatland evolution in the transient model run. The model simulates the establishment and expansion of peatlands under favorable conditions as well as 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. 4b). 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 (Figs. 4c, 6).

Figure 4 Simulated peatland area over time (a), gross positive and negative peatland area changes in 0.5 kyr bins, the evolution of today's simulated peatland area and old peat area (b), and driver contributions to the same changes (c), calculated using factorial simulations (see Sect. 2.2). Contributions by regions (b) and by drivers (c) are plotted cumulatively. Vertical bars indicate the Last Glacial Maximum (LGM) period, the Heinrich Stadial 1 (HS1) Northern Hemisphere cold phase, the Bølling–Allerød (BA) Northern Hemisphere warm phase, and the Younger Dryas (YD) Northern Hemisphere cold phase.

Figure 5 Maps showing the changes in the peatland fraction during three periods of past substantial climate change: (a) the Heinrich Stadial 1 (HS1, 17.43–14.63 kyr BP), (b) the Bølling–Allerød Northern Hemisphere warm period (BA, 14.63–12.85 kyr BP), and (c) the Younger Dryas Northern Hemisphere cold period (YD, 12.85–11.65 kyr BP).

3.3.1 22–17.43 kyr BP

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 grid cells still approaching equilibrium after the spin-up (Fig. A1). North America already sees an accelerating carbon accumulation with unchanging area before the HS1, driven mostly by increasing temperature. Carbon and area also increase in Europe with large temperature-driven fluctuations

Figure 6 Drivers of the change in peatland area from LGM to present. Colors indicate the most important driver, and colored shading shows the contribution of the respective driver on a scale from 0 (no contribution) to 1 (only contributor).

3.3.2 17.43–11.65 kyr BP

Three main features characterize the peat area evolution over the last glacial termination: (i) a northward shift in the distribution of northern extra-tropical peatlands, including peat expansion in northern Asia; (ii) dipole-like north–south shifts in tropical South America, associated with north–south shifts of the rain belts of the intertropical convergence zone (ITCZ); and (iii) flooding of peatlands on continental shelves, mostly in Southeast Asia, due to rising sea levels.

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₂ (Monnin et al., 2001; Shakun and Carlson, 2010; Ritz et al., 2013). The termination is divided into the Heinrich Stadial 1 (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 reorganization of ocean circulation, which is 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 and with it the high precipitation zones in the tropics (McGee et al., 2014; Shi and Yan, 2019; Cao et al., 2019). These climate dynamics are well captured by the transient TraCE21k simulation (Liu et al., 2009).

The responses of peatlands in LPX to these climatic changes are drastic. Large shifts in peatland area start to set in at the 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. 4b). In the northern mid and high latitudes, peatlands shift north and eastward (see Figs. 5, A2). Peatlands disappear in midlatitude North America and Europe, and new peatlands emerge at higher latitudes and in cold continental regions of Asia. These new peatlands include the large peat complex in the western Siberian lowland (WSL) region. Some of the peatlands established in northern Europe during HS1 vanish again during the BA. The described changes are driven by the TraCE21k climate which shows a substantial warming and wetting of the Northern Hemisphere beginning during the HS1. Temperature is the dominant driver for peat loss and expansion in Europe and North America. The loss of old peat is especially abrupt in North America (see Fig. A1b). Here, precipitation decreases and temperature increases abruptly over southwestern North America at 13 870 BP. Both changes decrease the water balance given by P−E which leads to a decrease in potential peat area and, thus, loss of the previous extensive peat complexes. As this abrupt climate change occurs at a discontinuity of the trace boundary conditions (changes in ice sheet configuration and freshwater forcing), the speed of this change is probably drastically overestimated. In northern Asia both temperature and precipitation drive the peatland expansion. This expansion sees a pronounced halt during the YD where the Northern Hemisphere climate briefly returns to more glacial conditions (see Fig. A2e).

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. 5). 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 Africa half the peatland area is lost during the BA mostly driven by drying. In Southeast 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 with the LGM and at the onset of the HS1. However, much less carbon is stored in active peatlands at the beginning of the Holocene 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.

3.3.3 11.65–0 kyr BP

Modeled peatlands in the Holocene show a continuous net expansion in the northern extra-tropics, with newly forming peatlands more than balancing the loss of peatlands elsewhere. The Holocene experienced relative stability in climate and CO₂ levels compared with the termination. The early to mid Holocene was likely characterized by warmer summer temperatures than the preindustrial period 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 at a rate of up to 12 mm yr⁻¹ (Henton et al., 2006).

For northern peatlands, positive area changes are consistently larger than negative changes throughout most of the Holocene (see Fig. 4a). This leads to a large continuous area expansion. Old peatland area is also simulated to increase continuously during the Holocene (Fig. 4b), 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 slowdown and 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 (Figs. A1 and A2). Here, both negative and positive dynamics 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 midlatitude 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. A1b).

The tropical peatland area is simulated to stay relatively stable throughout the Holocene, with positive changes balancing negative changes. Southeast 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 nonlinear effects. Peat area in South America increases slightly over the Holocene with mostly precipitation-driven fluctuations in between. On the other hand, fluctuations in region-integrated 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₂ concentrations. The new area gradually degrades again until 3 kyr BP with another peak at 0.5 kyr BP.

3.3.4 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; 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 ceased 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. Thus, we limit most of the model–data comparison of the transient behavior to today's existing peatlands. Figures 7 and 8c show modeled initiation date frequency, area, and carbon dynamics of northern peatlands that are still active at present.

Figure 7a compares LPX results to a gridded “oldest age” dataset compiled by Loisel et al. (2017), and Fig. 7b and c compare LPX results to two different reconstructions for lateral expansion (Loisel et al., 2017; Korhola et al., 2010) based on similar methods but different underlying peat core datasets (see Sect. 2.3). The two reconstructions for peat expansion agree on a limited pre-Holocene expansion, but they disagree substantially on the timing of fastest expansion during the Holocene (Fig. 7). Both simulated initiation and peat expansion have peaks about 4 kyr earlier than the reconstructions. The model simulates early initiation of today's northern peatlands, already beginning in HS1, and a large expansion during BA. The reconstructions, on the other hand, suggest lateral peat expansion picking up only with the transition into the Holocene. Agreement between model and reconstructions becomes good in the mid to late Holocene.

The early expansion in the model also propagates to the carbon balance for presently active peatlands. The model simulates 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. 7c). The summed simulated carbon increase from the LGM to the PI period in today's northern peatlands amounts to 313 GtC (Fig. 8c).

The early expansion of northern peatlands in the simulation is mostly dominated by peat establishment in western Siberian lowland (WSL) region and northern Asia in general (see Sect. 3.3.2 and Fig. A1). The dominant drivers of this expansion are temperature and precipitation, which, according to TraCE21k, both increase substantially over northern Asia during the HS1 and BA. A similar simulated early expansion into the WSL was reported by Treat et al. (2019), with the coupled CLIMBER2-LPJ setup. Morris et al. (2018) investigated possible climatic drivers for peat initiation in a modeling study using the HadCM3 model. They suggest that the WSL responded to an increase in effective precipitation at about 11.5 kyr BP, instead of the early warming. One source of the model data mismatch could lie in the uncertainties in climate anomalies discussed in Sect. 3.2.2. Especially at high latitudes, climate anomalies can vary greatly between climate models, and model performance at one point in time does not always correspond to performance at another point in time (Harrison et al., 2014). Although the freshwater change in TraCE21k was designed to capture the rapid climate events during the glaciation, the magnitude and timing of regional or even hemispheric changes can still have large biases. To date, TraCE21k is the only available transient GCM simulation, but new simulations under the umbrella of PMIP4 might shed more light on the model dependence of the warming pattern in question (Ivanovic et al., 2016). Another source of the mismatch could lie in the simple representation of peatlands in the model, which might be unsuitable to reproduce specific initiation and expansion pathways, like terrestrialization and fen–bog transition that might have been important controlling factors in that time and region (Kremenetski et al., 2003). One example could be the relative weakness of the initiation criteria on the moisture balance (precipitation over evapotranspiration >1), which is almost always weaker than the indirectly mediated condition on inundation persistence. This might pose a problem, especially in the WSL where moisture balance might have been the driving factor for peat initiation (Morris et al., 2018). Lastly, although the WSL is relatively densely sampled and reconstructions of peat initiation are robust, other areas of northern Asia are vastly under sampled and reconstructions are less reliable.

Throughout the tropics dated buried and active peat cores already show peatland presence during and preceding the LGM. However, peatland extent or evolution towards the presence are not well constrained and are subject to large uncertainties. The 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 the largest expansion rates between 8 and 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 that inferred by basal ages alone, with 90 % of today's peat establishing after 7 kyr BP and 60 % after 3 kyr BP. In this study, the dominant control on peatland area was found to be local sea level. Rising sea level during the termination and the early Holocene triggered the establishment of inland peatlands through alterations in moisture availability and the hydrological gradient, and the stabilization of sea level and subsequent sea level regression after 4 kyr BP prompted the establishment of coastal peatlands. In contrast to these reconstructions, the transient simulation shows that 60 % of today's tropical peatland area is already present in the LGM and only small expansion during the recent millennia. The sparsity of the data warrants caution when comparing reconstructions to model results. However, in Southeast Asia, this discrepancy could indicate the importance of the feedback of sea level on local hydrology, which is missing in LPX.

Figure 7 Simulated and reconstructed dynamics of today's existing northern peatlands: (a) peatland initiation frequency (two overlapping histograms), (b) peatland area expansion and (c) expansion rate, and (d) the net carbon balance (NCB) normalized by respective estimates of today's carbon pool. The reconstruction datasets are described in Sect. 2.3. The background colors indicate different time periods, as in Fig. 4

3.3.5 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. Thus, we quantify the overall contribution of peatland soils and peat carbon to the changes in the global 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 buried 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.

Figure 8 Carbon on global land (soil + litter + vegetation) (a); on active peatland areas, old peatland areas, and on flooded continental shelves (b); and in today's active peatlands (c). The brown line in panel (b) represents global carbon that originated from peatlands, even if it is not part of an active peatland anymore. The background colors indicate different time periods, as in Fig. 4

Figure 8b shows the temporal evolution of carbon stored in soils of active peatlands, old peat soil carbon remaining on former peatlands, and old peat carbon stored on flooded continental shelves. Here, the old peat pools presented exclusively include the carbon from the organic-rich layers of formerly active peatlands that remain after accounting for decomposition over time (see Sect. 2.1). In the PI period, 499, 139, and 22 GtC of peat carbon are stored in their respective pools. The total simulated increase in peat carbon from the LGM to the PI period within these three pools is 351 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 a value of 224 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 (Fig. 8c), and we would overestimate the net peat accumulation since the LGM. While the latter difference in net peat 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 as 12 GtC when including old peat, versus 102 GtC when only looking at the carbon in today's active peat for the period from 14.6 to 10 kyr BP.

Peatlands contribute about 40 % to the total land biosphere carbon increase of 893 GtC. The result for the total land carbon increase between the LGM and the PI period is in good accordance with a recent estimate, integrating multiple proxy constraints (median: 850 GtC; 450–1250 GtC ± 1 standard deviation) (Jeltsch-Thömmes et al., 2019). The model also simulates the total change of the land biosphere carbon inventory between the beginning of the Holocene and the preindustrial period in reasonable agreement with the reconstruction by Elsig et al. (2009). The simulated temporal evolution, however, is different, with a rapid uptake in the early Holocene in the reconstruction compared with a delayed uptake in the mid to late Holocene in the simulation.

4 Conclusions

We used the LPX-Bern dynamic global vegetation model to produce an in-depth model analysis of the transient area and carbon dynamics of global peatlands from the Last Glacial Maximum (LGM) to the present. For the LGM, peatland area, reduced to the tropics and northern midlatitudes, is predicted at 2.687 Mkm² in the transient run, storing 275.6 GtC of carbon. Under LGM climatic conditions, LPX-Bern predicts peatlands in areas with low or no peat cover at present or on currently 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 Paleoclimate Modelling Intercomparison Project (PMIP3). This results in a peat area range of 1.5–3.4 Mkm² with a carbon storage of 147–347 GtC. This large range illustrates the dependence 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 on their carbon storage make it difficult to further constrain this range. At the same time, there are currently 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 from the LGM to present.

A driver attribution of the simulated transient evolution of peatlands using factorial simulations showed regional and temporal differences. Modeled changes in the tropics were dominated by shifts in the position of the intertropical convergence zone and associated precipitation changes during the last glacial termination as well as by rising sea level. Changes in the northern high latitudes are mostly driven by temperature and precipitation increases. The 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 of area and carbon, with the fastest positive and negative changes happening during the termination (Heinrich Stadial 1 and Bølling–Allerød). This reveals a different perspective from 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. Thus, when assessing the net carbon balance of global peatlands over time, accounting for paleo-peatlands and their remains becomes essential. In our transient simulation the LGM to the PI period net peat carbon balance is predicted at 351 GtC, with 499, 139, and 22 GtC stored in the PI period 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 estimates from the literature, whereas predictions for the tropics are larger than most estimates. However, data constraints in the tropics are significantly weaker, as peatland science has long focused on the northern high latitudes and has only began accelerating its effort in the tropics over the last few decades. Even fewer data are 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 the uncertainties on a global scale. It contributes to a foundation for a better understanding of past peat dynamics and emphasizes the importance of treating and understanding peatlands as dynamic and evolving systems. In a next step, the results presented here can serve as a starting point for projections of future peat dynamics under different scenarios.

A growing database of buried peat and knowledge emerging from the growing literature on anthropogenically drained peatlands might shed more light on the fate of old peat carbon and inform future modeling studies. New time slice and transient climate model simulations under PMIP4 (Ivanovic et al., 2016; Kageyama et al., 2017) as well as an increased effort from the peat community to fill gaps in sample coverage both for today's peatlands and buried peat layers, especially in North America, northern Asia, and the tropics, might help to constrain past peat dynamics further and to test the robustness of the results presented here. At the same time, there is potential for improvements to the LPX-Bern that could decrease model data mismatches, especially on the regional scale. Future improvements could include refining the moisture balance criteria on peat initiation, improving hydrology and boundary conditions on continental shelves, and finding key processes that might benefit from a more complex representation, such as a multilayer peat profile and distinctions between different peatland types.

Appendix A

Figure A1 Simulated global and regional peatland area and carbon dynamics over time, relative to PI levels. PI levels are given in millions of square kilometers (M km²) for peat area and in gigatons of carbon (GtC) for peat carbon. The extent of the northern Asia, western Siberian (WS) lowlands, and Southeast Asia regions are shown in Fig. 1. The background colors indicate different time periods, as in Fig. 4

Figure A2 Simulated global and regional gross positive and negative changes in peatland area in 0.5 kyr bins, relative to PI levels. PI levels are given in millions of square kilometers (M km²) for peat area and in gigatons of carbon (GtC) for peat carbon. Colors indicate driver contributions to changes attributed using factorial simulations. The extent of the northern Asia, western Siberian (WS) lowland, and Southeast Asia regions are shown in Fig. 1. The background colors indicate different time periods, as in Fig. 4