Large stocks of peatland carbon and nitrogen are vulnerable to permafrost thaw

Based on the peatland map products compiled in this study (25), we estimate northern peatland extent to be 3.7 ± 0.5 million km² (mean ± root-mean-square error [RMSE]), including 1.7 ± 0.5 million km² in permafrost (Fig. 1 A and B, Table 1, and SI Appendix, Table S1; peatlands defined as >40-cm surface organic soil material). This extent is similar to inventory-based estimates (1, 19, 26, 27) but suggests that both the global soil map WISE30sec (23) and the global PEATMAP dataset (28) underestimate northern peatland extent by ∼1 million km². Our map is relatively consistent with the global maps in areas of very high peatland cover (e.g., West Siberian Lowlands and Hudson Bay Lowlands), but we identify more small peatland complexes outside of the core peatland regions (e.g., in central and eastern Siberia and the European Arctic west of the Ural Mountains). Our map of peatland extent was derived from the mean of two independent soil maps, which are the highest-resolution maps for delineating northern peatlands we know of, that is, harmonized national soil inventory maps (29, 30) and the global digital soil map SoilGrids250m (31). These maps separate peatlands with and without permafrost, but we include a maximum threshold for permafrost occurrence at mean annual air temperature (MAAT) of greater than or equal to +1 °C. See SI Appendix, supplement section S1 for detailed information on ground-truthing the peat extent maps using local observations.

We determine peatland C stocks using a newly compiled dataset of peat cores with observations of peat depths (n = 7,111) of which some cores include peat organic C and N storage (n = 782 and 105, respectively) distributed across the Northern Hemisphere (Fig. 1C and Dataset S6). Based on the peat core data, the average peat C storage is 106 ± 66 kg C⋅m⁻², ranging from 0.4 to 593 kg C⋅m⁻² with substantial variability over short distances (Fig. 1D). From these, peat cores maps of peat C and N stocks (Fig. 1 E and F) are made by combining a machine-learning model of peat depth, the maps of peatland extent, and linear models for predicting peat C and N storage from peat depth (SI Appendix, Fig. S3). The machine-learning model of peat depth combines peat core data with spatial environmental data (summarized in SI Appendix, Table S2) and shows that underlying mineral soil texture, the extent of peatlands in the area, photosynthetically active radiation, and mean summer temperature were the most important variables for predicting peat depth (SI Appendix, Fig. S4). The maps show a mean peat depth across the circumpolar north (mean ± SD) is 249 ± 97 cm and that total peat C and N stocks are 415 ± 150 and 10 ± 7 Pg (mean ± RMSE), respectively (Table 1). We estimate that permafrost peatlands store 185 ± 70 Pg C and 7 ± 4 Pg N, a substantial part of the total stocks (Table 1). These estimates are largely consistent with most of previous estimates based on data aggregated from tables (1, 19, 26). However, the spatially explicit maps reveal patterns in peat C and peat N stock distribution that have been missed previously. We map larger extent of permafrost peatlands, but we also find these peatlands to be shallower with lower C stocks per area unit than previously assumed. Permafrost peatlands are on average ∼80 cm shallower than permafrost-free peatlands. This confirms that earlier local/regional findings of limited vertical peat accumulation in permafrost peatlands (34) are applicable across the permafrost region.

Combining our peatland C stock estimates with existing estimates for peatlands in the tropics [105 Pg C (35, 36) and the extratropical Southern Hemisphere (15 Pg C (2)], we estimate that peatlands store 530 ± 160 Pg C globally, with northern peatlands accounting for ∼80% of the total. The maps of northern peatlands reveal a very pronounced latitudinal pattern in peatland extent, with nearly half the global peatland C stored between latitudes of 60 and 70° N (SI Appendix, Fig. S7).

Our estimated peatland C stocks are difficult to reconcile with a recent estimate of >1,000 Pg C in northern peatlands (20). That study uses a conservative estimate of peatland areal extent (2.9 million km²) and argues that early onset of peatland expansion after deglaciation caused very high, sustained, C accumulation. Our approaches differ notably, e.g., in their inclusion of data from outside the northern peatland region, their lack of direct bulk density data, and their lack of any observational constraints to peat depth. To accumulate >1,000 Pg C in 2.9 million km², mean peatland depths of ∼5 to 6 m are needed—twice as deep as suggested by our >7,000 data points. There is also a difficulty in reconciling a >1,000 Pg peatland C stock within the global carbon budget constraints offered by marine and ice-core paleo records of atmospheric CO₂ concentration and isotopic composition (37), whereas our study can be reconciled with those top-down estimates. We also note that this high C stock estimate is currently being questioned elsewhere (37, 38).

Combining the peatland maps with syntheses of peatland annual flux and C accumulation observations, we calculate present-day peatland GHG fluxes as a sink of atmospheric CO₂ at 0.10 ± 0.02 Pg C⋅y⁻¹, a source of CH₄ at 0.026 ± 0.002 Pg C⋅y⁻¹, and a source of nitrous oxide (N₂O) at 0.022 ± 0.005 Tg N⋅y⁻¹. These are empirically based spatial estimates of northern peatland GHG balances, and the results are similar to previous estimates (3, 6). We estimate losses into aquatic systems (dissolved and particulate organic matter) to be 0.022 ± 0.02 Pg C⋅y⁻¹ and 0.7 ± 0.5 Tg N⋅y⁻¹. We note that the observational GHG data, for CO₂ and especially N₂O, remain very limited from northern peatlands.

Our estimated net C sink varies with MAAT. We find a significant relationship between MAAT and C accumulation rates over the past 2,000 y (logistic growth model, n = 129, P < 0.05, R² = 0.3; SI Appendix, Fig. S8). By extrapolating this relationship spatially, we found an average accumulation rate of 34 g C⋅m⁻²⋅y⁻¹ across the northern peatland region. We note that such long-term C accumulation estimates also implicitly include fire dynamics occurring naturally over time at these sites. Our meta-analyses reveal no clear climatic controls on CH₄ and N₂O fluxes; but these GHGs vary depending on peatland types (Dataset S1). The variability from different peatland types is accounted for by mapping distributions of permafrost peatlands, and estimating bog and fen cover from biome distributions (SI Appendix, Table S5).

As sinks of CO₂, but net sources of CH₄, peatlands cool the climate over long timescales. Using a snapshot of annual peatland flux, peatlands warm the climate over decadal time horizons, but cool it over longer time periods (39). Using a radiative forcing model (39), we estimate that the isolated radiative forcing from one year of present-day peatland GHG exchange peaks after 30 to 40 y at +0.075 W⋅m⁻², mainly caused by net CH₄ emissions (SI Appendix, Fig. S9). Over centuries, there is a net cooling, caused by CO₂ sequestration through photosynthesis, which is reached after ∼350 y. These time constraints on radiative forcing are sensitive to the ratio of CH₄ emission to CO₂ uptake, which changes under peatland disturbances such as permafrost thaw.

Based on an equilibrium model, we estimate that the preindustrial extent of permafrost in peatlands was ∼2 million km², with a present-day coverage of 1.7 million km². This area is projected to decrease to 1 million km² at a 2 °C global warming stabilization above the preindustrial (Fig. 2B). At 6 °C global warming, we project that almost no peatland permafrost would remain. To model these permafrost losses, we used the present-day relationship between peatland permafrost extent and MAAT, extracted from maps (SI Appendix, Fig. S10), and projected it into future scenarios of global average warming stabilization above preindustrial (from +0.5° to +6 °C) using ensemble ESMs (Fig. 2A). This approach essentially adapts the method of ref. 40 specifically to permafrost peatlands. Because our approach was based on assumptions of equilibrium rather than transient responses processes, we cannot project how long it will take permafrost to thaw, but rather what the net, long-term effect will be.

To assess the effect of permafrost thaw on peatland C and N budgets, we distinguish four main potential stages in the long-term transition from stable permafrost to nonpermafrost peatlands (Fig. 2 C and D, Dataset S1, and SI Appendix, Fig. S11A). 1) Intact permafrost peatlands are sinks of CO₂ and have near-neutral CH₄ and N₂O balances (41–44). 2) Gradual active layer warming and deepening cause releases of CO₂ and N₂O from the active layer and from newly thawed peat while CH₄ remains near neutral (45–47). If thaw progresses into ice-rich permafrost, thermokarst may occur (18). 3) Young thermokarst stage fens and bogs are CO₂ sinks and CH₄ sources (44, 48, 49). There is no evidence of strong GHG losses from thawed peat (12, 50, 51), but chronosequence studies suggest large net losses of previously frozen peat (13, 18, 52), which we suggest may occur via dissolved or particulate organic C (DOC or POC) fluxes into aquatic ecosystems or through shorter transport and reposition. Young thermokarst lakes are sources of both CO₂ and CH₄ (53, 54). 4) Stabilization of postthaw peatland stages over centuries leads to weaker CH₄ sources with time, and lakes that were CO₂ sources transition into CO₂ sinks (18, 54–56). We note that all four stages may not occur everywhere, and in many cases different stages may occur simultaneously across a peatland complex. In some sites, a permafrost peatland may experience, under drier conditions, a more extended stage 2 of thaw (active layer deepening) and then progress to a postthaw stage. The spatial model approximates this variability with probability distributions for different stages (SI Appendix, Fig. S11A). Our model framework allows us to explore postthaw C balances with simplified process representation, including the extent to which changes in peatland C stocks are attributable to active layer deepening or thermokarst expansion (18), but does not include other disturbances. Long-term data on C accumulation rates inherently include peat losses to fire, but our framework does not account for C and N losses from rapidly increasing peatland fire frequencies (8) or droughts (9). Our model does not account for landscape-scale hydrological impacts of thaw, such as increased hydrological connectivity (57) or increased evapotranspiration (58).

Our modeling projects that permafrost thaw will cause a transient period of positive radiative forcing from peatlands, which will last one to three centuries (Fig. 3). This added radiative forcing is calculated from the difference between baseline peatland GHG balances at present climate and thaw scenarios. Under warming scenarios across the full range of +1.5 to +6 °C, transient losses of 2 to 6 Pg C as CO₂ and CH₄, but no significant losses of N₂O are projected (Dataset S2 and SI Appendix, Fig. S11B). Two centuries after initial thaw, the combined radiative forcing of this transitional GHG release reaches 0.05 and 0.13 W⋅m⁻², respectively, for +2 and +4 °C global warming scenarios (Fig. 3 A and B). The radiative effect is mainly caused by CH₄ release. For CO₂, initial release from peat decomposition during active-layer deepening is compensated by a net sink effect in the thermokarst stage. We project smaller CO₂ losses, but similar fluxes of methane CH₄, compared to a previous study of abrupt thaw (18). We estimate a minimal radiative forcing contribution from N₂O, allaying concerns of potential added N₂O forcing from peatland thaw (46). After long-term stabilization of the thaw-pulse (>200 y), the peatlands will be an annual net sink of C, a source of CH₄, and near neutral for N₂O (SI Appendix, Fig. S11C). Both the strength of the CO₂ sink and the CH₄ source increase with MAAT. This is because a warmer climate increases the peatland CO₂ sink capacity (due to higher plant productivity), but also increases the fraction of peatlands that change from dry permafrost peatlands to wetter permafrost-free peatlands, thus increasing CH₄ emissions. Warmer and wetter climates in the future may also lead to the formation of new peatlands, and C accumulation, in high-latitude regions (4, 59), but this is not addressed in our model calculations.

To assess the societal relevance of the projected GHG fluxes, we compare estimated net added GHG fluxes from the permafrost thaw to full anthropogenic radiative forcing scenarios consistent with United Nations Framework Convention on Climate Change Conference of the Parties negotiations (60). Across the full range of +1.5 to 6 °C global warming, northern peatland emissions amount to adding +0.4% to +0.8% onto the forcing from projected anthropogenic emissions by the year 2050 (Fig. 3C). For the low-emission scenarios, the contributions from permafrost thaw in northern peatlands peak near +2% of human emissions, but not until the next century.

In addition to GHG losses to the atmosphere, we project additional net cumulative lateral losses to aquatic systems of 10 to 30 Pg C and 0.4 to 1.1 Pg N from permafrost thaw under +1.5 to +6 °C global warming (Fig. 3 D–F). These projections are based on spatial modeling of synthesized peatland permafrost-thaw chronosequences, which have shown rapid (decadal) postthaw losses of deep permafrost peat (13, 18, 61) (SI Appendix, Table S7). These estimates are based on only five chronosequence studies and are more speculative than our other projections. Similar data have been used to infer C losses as CO₂ (18), but because GHG flux observations do not seem to support thaw-induced losses of this magnitude (Dataset S1 and refs. 12, 50, 51), we suggest these would most likely occur laterally. This putative lateral loss may occur as DOC or POC fluxes into aquatic ecosystems or through shorter transport and reposition in adjacent ecosystems. If correct, such losses of peat C could have significant implications for aquatic biogeochemistry and ecosystems (62). Averaged over space and time, our modeled losses add up to 11 to 18 kg C⋅m⁻² per unit of thawed peatland over the full ∼100-y postthaw period. Although limited, there are observations of lateral flux from thawing permafrost peatlands against which these numbers can be compared (SI Appendix, section S4), and we conclude that fluvial losses of these magnitude can be supported by data from West Siberia (63, 64) but not from a Boreal peatland dominated catchment in western Canada (65). Our projected total net gaseous and lateral C losses are 5 to 10 Pg C higher than those estimated in ref. 18, a difference mainly attributable to a higher projected area of permafrost thaw and the inclusion of active layer deepening losses in this present study. We note that large fractions (half or more) of C transported as POC/DOC may be degassed as CO₂ directly from inland water surfaces (63, 66), potentially increasing the total atmospheric burden caused by permafrost thaw in peatlands.

A total of 7,111 geolocated peat cores with peat depth data was compiled (Dataset S6). Only sites where basal peat was reached are included. A subset of 782 cores have data on peat organic carbon content (OC% by weight) and dry bulk density. A subset of 105 cores has additional data on peat total N content (weight % N). The sources of data were from refs. 26, 71–75, and previously unpublished data.

This study based estimates of peatland spatial extent on soil classification maps. The study region is limited to the extratropical northern hemisphere (defined as north of 23° latitude). Three different map products were used or evaluated for their capacity to accurately map peatland extent: the global WISE30sec dataset (23), the global SoilGrids250m dataset (31), and harmonized national and regional soil maps (29, 30). We refer to these references for details about how the maps were made. The WISE30sec dataset was not used for peatland mapping as it had too low resolution. The SoilGrids250m and the national/regional soil maps were combined and harmonized for this study; see SI Appendix for more details. All datasets were projected using equal area projections and were resampled to 5-km grids using bilinear interpolation.

The spatial scaling of peat depths was carried out using random forest machine learning (RFML). Random forest is a tree-based machine-learning method that uses bootstrapped samples (here peat cores) to grow a large number of decision trees (n_tree) with randomized environmental predictors at each tree node (m_try). These trees are then averaged to predict new data (76, 77). A RFML model, with 1,000 trees (n_tree), was trained using the observational data. In total, 6,038 peat cores had sufficiently precise geolocation and matched point-to-pixel overlays for all of the environmental training variables (n = 12; SI Appendix, Table S2). We used a 10-fold cross-validation with five repetitions providing m_try as a tunable parameter for model training using the caret package in R (77, 78). We applied bias correction to the predicted peat depths using best angle residual rotation of the peat depth map (SI Appendix, Fig. S6, and ref. 79).

By combining the RFML model of potential peat depth with the map of peat coverage, we calculated area-weighted peat depths and peat volumes. To calculate stocks, the modeled peat depths were used to estimate peat organic C and N storage (kilograms of C or N per square meter) using linear relationships formulated based on the peat core data (SI Appendix, Fig. S3). The estimated C and N storage was then used to calculate total C and N mass per pixel. Uncertainties are reported as RMSE based on 5th/95th percentiles of residuals between modeled and observed values of peat depth.

The baseline C and N balances of peatlands, including GHGs (CO₂, CH₄, and N₂O), were estimated based on paleo-reconstructions of C balances as well as syntheses of flux measurements from permafrost- and permafrost-free peatlands (Dataset S1). Paleo-observations were used for long-term net C budgets and syntheses of GHG flux measurements for shorter time intervals in projections of thaw. We developed a simple spatially explicit inventory model to assess the impact of peatland permafrost thaw scenarios on the stocks of C and N as well as GHG fluxes.

The applied permafrost thaw scenarios (SI Appendix, Fig. S2A) assume that, once the temperature threshold for thaw is crossed, the peatlands are affected by active-layer deepening for a period of ∼25 to 75 y (with a mean of 50 y) until the thaw progresses into ice-rich, deeper peat. This time period was calculated based on active-layer deepening of 1 cm per year (estimated from refs. 46, 80) and that the average depth to ice-rich peat from the bottom of the active layer in permafrost peatlands is ∼25 to 75 cm (calculated from data in refs. 72, 81–83). If thaw progresses into the ice-rich core of the permafrost peatland, thermokarst (ground collapse) occurs. Postthaw thermokarst peatlands or lakes were assumed to gradually transition to mature thermokarst systems over 50 to 150 y (a mean of 100 y). This time period of transition into mature thermokarst was based on an average of studies on postthaw chronosequences, which suggested somewhat longer transition times of ∼150 to 200 y (estimated from refs. 13, 18, 54, 61) and remote-sensing studies that showed substantial lake drainage or fen-vegetation infilling in some areas over periods of a few decades (56, 83).

For flux scaling, we separated nonpermafrost and permafrost peatlands from postthaw peatlands. All classes were further separated into minerotrophic and ombrotrophic peatlands, but only if there were statistically significant differences in C accumulation rates or GHG balances. The spatial extent of minerotrophic and ombrotrophic peatlands was scaled from the Canadian Peatland Map (84), as fractions within tundra, boreal, and other biomes (includes temperate, oceanic, mountain, and prairie climate regions; biome distributions from ref. 32; SI Appendix, Table S5).

The C balance of stable peatlands was modeled based on observed long-term apparent C accumulation in the late Holocene (last 2,000 y) from northern (n = 122; ref. 26) and tropical (n = 7; ref. 2) peatlands (Dataset S4). The best model fit was achieved with a logistic model (S-shaped curve) that is able to model growth with saturation at both high and low temperatures (85) (R² = 0.3, Akaike information criterion = 28,822; SI Appendix, Fig. S9B). This model predicted mean potential C accumulation scaled for the extratropical northern hemisphere under present baseline climate is 34 g C⋅m⁻²⋅y⁻¹.

The CO₂-C fluxes during thaw stages were based on a meta-analysis of full year budgets from thawing permafrost in the literature (SI Appendix, section S1.5 and Dataset S1). The net C budget following permafrost thaw was based on chronosequence studies of postthaw permafrost peatlands (13). Old permafrost C is lost following thaw, while increased ecosystem productivity in the young thermokarst (postcollapse) means that the surface peat is gaining C. In the early thaw stages, the loss of old C is much more rapid than the gain of new C. The loss of old permafrost C can be estimated as a function of prethaw C stock (13, 18) (SI Appendix, Table S7). The C loss during the first 100 y after thaw was estimated from the peatland C stock maps using the simplified equation y = 1.1451x^−0.0771, where y is the fraction of prethaw C that is lost in 100 y after thaw and x is the stock of prethaw C in kilograms of C per square meter (R² = 0.93, from 100 y in SI Appendix, Table S12). We scaled the changes in N pools from the C pools based on typical C:N ratios of permafrost peatlands and nonpermafrost peatlands in tundra regions and boreal regions (SI Appendix, Table S6).

All data for estimated CH₄ fluxes were from a recent synthesis of year-round CH₄ fluxes in northern wetlands (86). We used only sites with organic soils and separated nonpermafrost, permafrost, and postthaw sites. We further distinguished the minerotrophic peatlands (swamp, marsh, and fen classes, following the Canadian wetland classification system) from ombrotrophic peatlands (bogs).

For minerotrophic and ombrotrophic permafrost-free peatlands, we used annual N₂O budgets from a synthesis of N₂O fluxes from northern soils (87). Annual/seasonal N₂O data from Arctic peatlands are limited to a single site located in western Russia with discontinuous permafrost. We used published N₂O flux data from this site (41, 88) as N₂O emission estimates for minerotrophic and ombrotrophic (bare and vegetated) permafrost peatlands (Dataset S1). Data on N₂O (46) and CO₂ (47) fluxes from peat mesocosms during simulated permafrost thaw were used to develop a scaling ratio of N₂O release relative to C release (see SI Appendix for more details).

The model of permafrost fraction in peatlands was derived using the method developed in ref. 40, where a relationship between permafrost fractional coverage and MAAT was fitted by minimizing RMSE between MAAT and mapped permafrost fraction in peatlands. The equation used is as follows:

where ERFC is the complementary error function (using the pracma R package).

As in ref. 40, the curve was refitted using “maximum” and “minimum” permafrost fraction to give upper and lower estimates of permafrost fraction, as well as a central estimate. The maximum and minimum extents were derived from the highest and lowest per-pixel estimates of permafrost fraction in the national polygon maps and SoilGrids, respectively. Thus, three different parameter values (for central, upper, and lower curves, respectively) were fitted for μ (1.95, 0.7, and 3.1), σ (7.35, 6.1, and 4.5), and f_max (0.92, 0.96, and 0.86). We assumed that mapped permafrost extent at +0.5 °C global warming (relative to preindustrial levels) was in quasi-equilibrium with the climate of the 1960 to 1990 period, which we also consider to be representative for the permafrost extent in the maps. Global warming stabilization scenarios at 0.5 °C intervals up to a maximum of 6 °C were used for the future projections.

The projected GHG budgets, including CO₂, CH₄, and N₂O fluxes, from the spatial model were used to calculate the future radiative forcing effect. A range of GHG flux scenarios (available in Dataset S3) were exported from the spatial model and used as input in a radiative forcing model (39), with additional parameterization for N₂O and modifications to atmospheric CO₂ lifetimes (89). Separate GHG flux scenarios were calculated for stabilized permafrost conditions at 0.5° increments from 0° to +6 °C global warming stabilization (background concentrations were stable anthropogenic present-day emissions). Separate runs were also done for fluxes resulting from the transient thaw scenarios for each incremental warming, with the added radiative forcing from permafrost thaw calculated from the net changes in GHG fluxes relative to stable baseline peatland GHG balances at present (Fig. 2 and Dataset S1). The net radiative effect of the transient thaw is calculated as the difference between stable scenarios and the transient scenarios. To compare the magnitude of permafrost–peatland thaw emissions to anthropogenic emissions, radiative forcing from anthropogenic emissions together with peatland thaw emissions were compared to anthropogenic emissions alone. The projections for anthropogenic emissions were retrieved from Climate Scoreboard (60) and are computed using the C-ROADS climate policy model (90). For these calculations, we assume that +0.5 °C global warming is consistent with peatland fluxes in 1990 to 2000 (assuming decadal lags in thaw from the 1960 to 1990 climate normal) and that +1 °C warming is consistent with present day.