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

Significance Over many millennia, northern peatlands have accumulated large amounts of carbon and nitrogen, thus cooling the global climate. Over shorter timescales, peatland disturbances can trigger losses of peat and release of greenhouses gases. Despite their importance to the global climate, peatlands remain poorly mapped, and the vulnerability of permafrost peatlands to warming is uncertain. This study compiles over 7,000 field observations to present a data-driven map of northern peatlands and their carbon and nitrogen stocks. We use these maps to model the impact of permafrost thaw on peatlands and find that warming will likely shift the greenhouse gas balance of northern peatlands. At present, peatlands cool the climate, but anthropogenic warming can shift them into a net source of warming.

used for this purpose. The peatland extent (organic soil extent) in WISE30sec dataset is from the HWSD and can be traced back to the original FAO soil map of the World (17). The regional soil maps were combined and harmonized for this study. Such maps are available for the geographic regions Canada, Europe (except Italy), Mongolia, USA and Russia. The soil maps used were the Northern and Mid-latitude soils maps (14) and the NCSCDv2 (10,18). We refer to these references for details about how the maps were made. These two databases were projected to equal area projections and merged, preserving NCSCDv2 where they overlapped. The dataset was checked for consistency and adjusted so that the various mapped coverages added up to 100% (case by case solution applied, whenever the source of error was not evident water cover was adjusted, but preserving the relative coverage between other classes). Sliver polygons where removed by identifying those polygons that had an area <5000 m 2 and a circularity index value of <0.01. The extent of mapped peatlands and permafrost in peatlands in the maps was evaluated against individual peat cores and local-scale high resolution maps (see supplement).

S1.4.1 Spatial upscaling of peat depth
Upscaling of point observations to areal coverage can be made either based on thematic mean upscaling (e.g. mean values for different regions, soil units or landscape units) or some form of digital mapping method. The peat depth dataset is highly spatially clustered and exhibits characteristics that makes it unsuitable for analyses with predictive interpolation methods that rely on semivariogram modelling, such as kriging or co-kriging. This conclusion was based on extensive examination of the properties of the data, including semivariograms at multiple distances, which consistently showed a lack of local-scale autocorrelation. While Geographically Weighted Regression (GWR) can work well on spatially clustered data it can be sensitive to limited local-scale spatial autocorrelation in the data. Machine learning methods have been shown to produce relatively robust digital maps with a wide range of input data, even with spatially clustered or noisy data (19). Random Forest Machine Learning (RFML) has been identified as a suitable machine learning method for the spatial interpolation of environmental variables and several studies have shown its successful and often superior application in digital soil mapping (20)(21)(22), including peatland depth and SOC stocks (23,24) . RFML is a tree based machine learning method that can use an ensemble of bootstrapped trees for regression estimates (25). In RFML each regression tree depends on a random selection of environmental predictors at nodes. The ensemble of all trees are then averaged for prediction. We selected 1000 trees (ntree). A great advantage of tree based machine learning methods is their ability to cope with non-linear relationships between the predicted variable and covariates representing the environmental space (26). This makes them particularly suitable for regional to global mapping approaches. Modelling was performed using R statistical software (2), following standard procedures as outlined in ref (27). The RFML was applied to the 12 environmental predictor variables (Table S2) and a subset of the peat depth cores for which it was possible to extract corresponding values of the environmental variables. The model was then trained using 10-fold cross-validation with five repetitions to develop a stable model using the caret package (28). The best performing model was automatically selected based on the lowest RMSE (142.2 cm) corresponding to an R 2 0.304 at a mtry value of 3. Predictor variables reflecting soil parent material texture, mean summer temperature and pixel peatland cover have the highest variable importance in predicting peat depth at circumpolar scale (Fig. S4). To reduce an observed regression to mean effect in the original RFML model, we applied bias correction using residual rotation on the model, by adapting an existing approach to raster maps (29). The effect on observed against predicted peat depth in values of the raster map of peat depth is presented in figure S5. The residual rotation decrease the MSE from 9984 to 9252 and improved Lin's concordance correlation coefficient from 0.767 to 0.868, which is a measure of accuracy and precision of agreement to a 45• line (30). Because of the very high variability in peat depth, and problems with some imprecise geolocations of point data, the spatial resolution of the analyses was reduced to 10 km pixels (from an original 5 km resolution). Data users should note that the maps are generalized and meant to show mean peat depths regionally. As with all generalized and broads-scale maps, it is not appropriate, or good scientific practice, to relate information from a single site or peat core to the generalized 10*10 km pixel values. There is sometimes large differences between individual peatland cores and the mean mapped depths, which is evident from a map of residuals between observed and predicted peat depths (Fig. S6).

S1.4.2 Calculating peat depth, volume and stocks of C and N
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 modelled peat depths were used to estimate peat organic C and N stocks (kg C or N per m 2 ) using linear relationships formulated based on the peat core data (Fig. S3). The linear functions to estimate peat C and N stocks from peat depth were done separately for permafrost-free and permafrost peatlands, respectively. The linear models were fit on log-transformed data using Major Axis (MA) linear models (31) (considered robust with limited data) with intercepts set to zero (Fig. S3, see supplement section S2 for more details). The estimated C and N stocks were then used to calculate total C and N mass per pixel.

S1.4.3 Uncertainties in estimates of peat depth, volume and stocks of C and N
In addition to the 10-fold cross validation of the RFML map we also performed a consistent uncertainty estimate of peat depth, volume and stocks of C and N against the full point datasets for depth, C stock and N stock. For each variable, we calculated the root mean square error (RMSE) between observed and modelled values for all observations. We also calculated trimmed RMSE (based on 5th/95th percentiles of residuals), R 2 between modelled and observed values and deviation of the mean in modelled and observed. Given the very high variability in the point data (also over short distances), the trimmed RMSE was considered as the most robust estimate of error and this is reported in the main text and used to calculate uncertainty ranges for upscaled peat volume, C mass and N mass in Pg.

S1.5 Scaling C and N balances and projecting permafrost thaw
The baseline C and N balances, including GHGs (CO2, CH4 and N2O) of peatlands 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 to estimate constrain long-term net C budgets. Paleo-observations have the advantage of integrating long time periods and, thus, can be used to constrain long-term C budgets in a way that flux observation timeseries cannot. Another advantage of the paleo-observations is that they account for the total C and N losses including both atmospheric and aquatic losses, while most of the flux studies report just one of the two components. Syntheses of flux measurement data were used to estimate GHG balances and lateral losses (32) for shorter time intervals (such as post-thaw thermokarst) and to help attribute longterm changes in bulk C stocks to gaseous or lateral (hydrological) fluxes. We used a spatial model to assess the impact of peatland permafrost thaw scenarios on the stocks of C and N as well as GHG fluxes. The values, approaches, and data sources are summarized in Figure 2 and Dataset S1 of the supplement, but also described below. Permafrost peatland hydrology is locally complex and variable, with thaw causing both drying and wetting depending on local peatland morphology (33,34). We note that such complexities cannot yet be mechanistically accounted for; it would require couple process models with access to extremely high-resolution maps of peatland, soil and ground-ice properties. Here, we use a simplified conceptual model to represent the main stages of thaw and represent spatial and temporal uncertainty by statistical uncertainty range assessment and probability distributions. The applied permafrost thaw scenarios (Fig 2a) assume that once the temperature threshold for thaw is crossed, the peatlands are affected by active layer deepening for a period of ca. 25-75 years (mean 50 years) 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. (35,36) and that the average depth to ice rich peat from the bottom of the active layer in permafrost peatlands is ca. 25 to 75 cm (calculated from refs. (9,37,38). If thaw progresses into the ice-rich core of the permafrost peatland, thermokarst (ground collapse) occurs. Post-thaw thermokarst peatlands or lakes were assumed to gradually transition to mature thermokarst systems over 50-150 years (mean 100 years). This time period of transition into mature thermokarst was based on an average of studies on post-thaw chronosequences which suggested somewhat longer transition times of ca 150-200 yrs (estimated from refs. (39)(40)(41)(42)) and remote-sensing studies that showed relatively rapid lake drainage or fen-vegetation infilling in some areas ((43, 44). Due to the large imprecision of these time-transition estimates, we generated wide probability distributions for these time phases which span a wide range of possible occurrences, this means that the active layer deepening phase could, in principle, last anywhere from 1 to 150 years and the young thermokarst phase from could occur already at year 1 after thaw and proceed up to 250 years ( Fig  S11). The uncertainty from the time ranges were propagated into the estimate of upscaled flux uncertainties. Also, during these stages there is variability in soil drainage conditions (wetness) which, in turn affects the GHG-balance (45). We cannot explicitly account for this fine-resolution variability, but we propagate the added uncertainty from this into the estimate of upscaled flux uncertainties. This is based on ref. (46) who in their thaw experiment included active layer deepening under both dry and wet conditions, from their data, we find that the range between wet and dry condition adds +24% to the flux uncertainty range. The last phase, mature thermokarst, transitioned over several centuries to stable peatlands (assuming with an equal coverage of minerotrophic and ombrotrophic peatlands in mature peatlands). The time constraint on this transitional period is very poorly known and in this study it was assumed to exceed 300 years. For the calculations of transitional flux rates it was set to 750 years (500-1000 year range following ref. (41)), but this assumption has a limited effect the results in the main text or figures since the projections are not reported beyond 300 years. In the scaling of fluxes, we separated non-permafrost and permafrost peatlands from post-thaw peatlands. All classes were further separated into minerotrophic and ombrotrophic peatlands, but only if there were statistically significant differences in the respective C accumulation rates or GHG balances. The mapping of minerotrophic and ombrotrophic peatlands is challenging as there are no global or hemispheric maps of these properties. But there is a regional-scale map over Canada (47) and we assumed that the general patterns of peatland morphology across biomes in Canada were broadly representative for the whole study area. The spatial extent of minerotrophic and ombrotrophic peatlands was scaled from the Canadian Peatland Map (47), as fractions within Tundra, Boreal and other biomes (includes Temperate, Oceanic, Mountain and Prairie regions; biome distributions from ref. (48). This data is summarized in table S5.

S1.5.1 Calculations of C and N balances
The CO2 fluxes and lateral fluxes of C and N in peatlands were assessed using a combination of approaches. In stable peatlands (with or without permafrost), the baseline apparent long-term C accumulation was modelled based on MAAT, while lateral C losses as dissolved organic carbon (DOC) in aquatic systems were estimated based on literature (Dataset S1). The CO2-C flux was calculated as the residual between long-term C accumulation, lateral C flux and CH4-C flux. In the permafrost thaw scenarios, CO2-C fluxes are based on a literature meta-analyses of measured annual fluxes from thawing permafrost (see table S4) while the net loss of old peatland C was modelled from observations from paleo-reconstructions in post-thaw chronosequences (see section 1.5.1.2).

S1.5.1.1 Estimating long-term apparent C accumulation rates
The C balance of stable peatlands was modelled based on observed long-term apparent C accumulation in the late Holocene (last 2000 years) from northern (n=122; ref. (5)) and tropical (n=7; ref. (49)) peatlands (data summarized in table S10). We used the late Holocene rates because longer time-records are characterized by variability caused by non-climatic factors and we were mainly interested in establishing relationships between C accumulation and present climate. There was no significant difference in the mean (or median) C accumulation of minerotrophic and ombrotrophic peatlands (ANOVA, Mann-Whitney test; p>0.05); nor were there a difference in the slope of the relationships between MAAT and C accumulation of minerotrophic and ombrotrophic peatlands (ANCOVA; p>0.05). There was insufficient data to estimate whether there was a difference in C accumulation between non-permafrost and permafrost peatlands. Thus, for the subsequent C accumulation analyses, all peatland types were grouped together. There was a significant linear relationship between C accumulation and MAAT (R2=0.27; p<0.01; AIC=29240; Fig. S9a). However, such a linear model unrealistically predicted negative C accumulation at MAAT below -18°C, where peatlands are known to exist. Based on our data points, the observed climate envelope of peatland sites has MAATs from -20° to +27° C. The best non-linear model fit is achieved with a logistic model (S-shaped curve) that is able to model growth with saturation at both high and low temperatures(50) (AIC=28822; Fig. S9b). This model predicted a C accumulation of ca. 5 g C m -2 yr -1 at MAAT -18°C, where the linear model becomes negative. Modelled C accumulation approaches zero at -40° C, but never becomes negative. At very high temperatures, it saturates below 64 g C m -2 yr -1 . The predicted long-term mean C accumulation scaled for the extratropical northern hemisphere under present baseline climate is 34 g C m -2 yr -1 . This model represents the best fit to the available data, but it is important to point out that the relationship is tenuous, as each site has been affected by multiple factors over the past 2000 years. To assess the validity of our approach for estimating CO2-C as a residual from long-term paleo-data, we compare it to data from the literature. For non-permafrost peatlands we estimate residual CO2-C fluxes of 36 and 61 g C m -2 yr -1 , respectively for ombrotrophic and minerotrophic peatlands while the numbers for permafrost peatlands are 27 and 34 g C m -2 yr -1 (based on MAAT of +6° C and -2° C, respectively). These numbers are similar to typical values reported from the literature (summarized in Dataset S1). We note that there is a potential bias in our estimates since the long-term C accumulation rate accounts for how natural fire dynamics have influenced peatland C, while this is not included to short-term GHG flux measurements.

S1.5.1.2 Net loss of C and N after permafrost thaw from chronosequences
The net C budget following permafrost thaw was based on chronosequence studies of post-thaw permafrost peatlands (39,41). Old permafrost C is lost following thaw, while increased ecosystem productivity in the young thermokarst (post-collapse) means that the surface peat is gaining C. In the early thaw stages the loss of old C is more rapid than the gain of new C, but there are differences between peatland types and sites. A compilation of five permafrost peatland thaw chronosequences (three epigenetic, two syngenetic, compiled in figure S7 in Turetsky et al. (2020) (41), show that epigenetic permafrost peatlands loose less deep C than syngeneic peatlands. There is an average net C loss of 19% of the total initial C storage after 100 yrs across all five chronosequences. Jones et al.
(2017) (39) show that the loss of old permafrost C can be estimated as a function of pre-thaw C storage, but these results are based on syngenetic peatlands only. Assuming that the difference between epigenetic and syngenetic sites scale equally across sites with different initial peat C storage, we adapted the curves describing loss as a function of pre-thaw peat C storage from Jones et al to the mean of the five chronosequences. Table S7 of the supplement shows this adapted relationship between pre-thaw C and fraction of reminaing post-thaw C as a function of time extracted from figure 6 in Jones et al (39). The C loss during the first 100 years 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 pre-thaw C that is lost in 100 years after thaw and x is the storage of pre-thaw C in kg C m -2 (R² = 0.93, from 100 years in Tab. S12). We scaled the changes in N pools from the C pools based on typical C:N ratios of permafrost peatlands and non-permafrost peatlands in tundra regions and boreal regions (table S6). We assumed these net peat C and N losses occur during the active layer deepening and young thermokarst stages and attributed the residual between our GHGs budgets (Tab. S4) and the predicted thaw-loss from chronosequences to lateral losses as dissolved or particulate organic matter in fluvial/aquatic systems. We considered a fluvial pathway more likely than GHG losses. To reconcile the C balance, annual C losses of several kg C yr -1 are needed (see ref. (45)) and our meta-analyses of GHG fluxes reveal no records of post-thaw C fluxes of that magnitude. However, we note that the observational network for monitoring thawing permafrost peatlands is very sparse, and rapid thawpulses of C could occur unobserved (41). Our approach here differed from that of ref. (41), who assumed that all the post-thaw losses were gaseous (CO2 or CH4). In our study, CO2-C flux for the short-term post thaw stages were based on literature syntheses, but the mature collapse scar peatland and stable post-thaw peatlands were based on C accumulation modelled from projected future stabilized MAAT.

S1.5.1.3 Estimates of CH4 fluxes
All data for estimated CH4 fluxes were from a recent synthesis of year-round CH4 fluxes in northern wetlands (51). We used only sites with organic soils (equals peatlands). We separated nonpermafrost, permafrost and post-thaw sites. We further distinguished the minerotrophic peatlands (Swamp, Marsh and fen classes, following the Canadian wetland classification system) from ombrotrophic peatlands (bogs). We tested if minerotrophic and ombrotrophic sites had significantly different median, or mean, CH4 (Mann-Whitney test, ANOVA; p<0.05). For both non-permafrost and permafrost peatlands, there were significantly higher median/mean fluxes from minerotrophic peatlands. The Post-thaw landforms did not show significant differences between ombrotrophic and minerotrophic peatlands and no further distinction was made in the modelling. The annual CH4 flux was not correlated to MAAT or MAP in non-permafrost, permafrost or post-thaw sites (Pearson linear correlation, R 2 <0.1; p>0.05) so no attempt was made to model CH4 flux as a function of climate.

S1.5.1.4 Estimates of N2O fluxes
For minerotrophic and ombrotrophic permafrost-free peatlands, we used annual N2O budgets from a synthesis of N2O fluxes from northern soils (52). Annual/seasonal N2O data from Arctic peatlands are sparse: to our knowledge, they are limited to a single site located in Western Russia with discontinuous permafrost ("Seida", 67°03'N, 62°57'E, 100m a.s.l.). We used published N2O flux data from this site (53,54) as N2O emission estimates for minerotrophic and ombrotrophic (bare and vegetated) permafrost peatlands (Dataset S1). During the initial stage of gradual deepening of the active layer we modelled an increase in post-thaw N2O emissions caused by N2O production from N forms released with peat decomposition, and to a smaller extent, release of trapped N2O (36). Data on N2O (36) and CO2 (46) fluxes from peat mesocosms during simulated permafrost thaw were used to develop a scaling ratio of N versus C release. First, we determined the theoretical amount of mobilized N, based on the C:N ratio with C lost as CO2. From this theoretical amount of N (including all N forms such as N stored in the microbial biomass and N lost as leaching or N2) we determined the percentage of N2O -N lost based on measured N2O fluxes. We used the obtained N2O -N loss values of 7.38% (bare peat) and 0.25% (vegetated peat) to calculate post-thaw N2O fluxes relative to the estimated CO2 fluxes. Bare peat is conservatively assumed to cover 3% of the peatland surface (55).

S1.5.2 Model of permafrost fraction in peatlands
The model of permafrost fraction in peatlands was derived using the method developed in ref. (56), where a relationship between permafrost fractional coverage and MAAT was fitted, by minimizing root mean squared error. For this study, the relationship was re-fitted using the fraction of permafrost in peatlands rather than the overall landscape fraction. Because the upper layers of peat often insulate the ground in summer and keep the soil cool the permafrost fraction in peatlands is usually higher than the overall landscape permafrost fraction because. The relationship between MAAT and peatland permafrost extent was represented by an error function with two parameters (μ and σ), based on ref. (57). For this study, the peatland maps did not reach 100% permafrost coverage even at very low MAAT. To model this, we introduced an additional overall scaling, fmax, representing the maximum permafrost coverage. This is consistent with small fractions of non-permafrost peatlands occurring also in very cold climates, usually caused by proximity to streams or lakes (37) (58). This gives the equation used as: where ERFC is the complementary error function (using the pracma R package). As in ref (56), the curve was re-fitted 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 (with no MAAT correction applied to SoilGrids which allow larger spread of potential ecosystem-protected permafrost). 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 fmax (0.92, 0.96 and 0.86). The μ values were smaller than in ref. (56) indicating that, indeed, peatlands allow permafrost to persist in warmer climates than is typical for landscapes dominated by mineral soils. The air temperature map for present day and future scenarios was as in ref. (56). We assumed that mapped permafrost extent at +0.5ºC global warming (relative to preindustrial levels) was in quasiequilibrium with the climate of the 1960 to 1990 period, which we also consider to be representative for the mapped permafrost extent in the maps of peatland permafrost fraction in ca. year 1990-2000 (approximate age of the map input data, which allows a lag time in peatland permafrost extent response to climate of several decades, consistent with ref. (38)). Global warming stabilization scenarios at 0.5ºC intervals up to a maximum of 6ºC (representing a very high-end RCP 8.5 scenario) were used for the future projections.

S1.5.3 Modelling radiative forcing
The projected GHG budgets, including CO2, CH4, and N2O fluxes, from the spatial model were used to calculate the future radiative forcing effect. A range of GHG flux scenarios (available in appendix Dataset S3) were exported from the spatial model and used as input in a radiative forcing model (59) with additional parameterization for N2O and modifications to atmospheric CO2 lifetimes (60). Separate GHG flux scenarios were calculated for stabilized permafrost conditions at 0.5 degree 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 (Fig. 2, Dataset S1). The net radiative effect of the transient thaw is calculated as the different 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. These emissions scenarios were retrieved from Climate Scoreboard (61) and are computed using the C-ROADS climate policy model (62). For these calculations we assume that +0.5°C global warming is consistent with peatland fluxes in 1990-2000 (assuming decadal lags in thaw response relative to the 1960-1990 climate normal) and that +1°C warming is consistent with present day.

S2: Maps of peatland and permafrost extent and assessments against local scale data
For the northern peatland region there are two, fully independent, sources of soil maps with substantially higher spatial resolution than the global maps in WISE30sec: (i) a collection of harmonized national soil survey maps created by soil scientists from field data, topographic maps aerial photographs and satellite remote sensing data (14) (further adapted for this study, see methods) and (ii) the SoilGrids250 m dataset (13), a digital soil map created from a combination of environmental variables and soil profile data using machine learning algorithms. These maps are based on entirely different input data, where made using different methods and at different times by different soil science experts. Because of this, we consider them to be independent and partially complementary, also acknowledging that combining the two maps may yield more robust scaling of peatland extent and permafrost extent in peatlands. There is not a complete overlap in these maps (Fig. S1). In SoilGrids, 95% of the northern extratropical peatland extent occurs within the region that is covered by the national soil maps There is a remarkable consistency in total mapped peatland extent between SoilGrids and the regional polygon maps (table S1). However, there are differences in the spatial distribution of peatlands and very substantial difference in the areas of mapped permafrost extent (Fig. S1). Because we are discussion soil maps we will from now on use USDA Soil taxonomy terminology (63). Here, nonpermafrost organic soils are called Histosol and permafrost affected organic soils are called Histels. In the national soil maps, the combined cover of Histels and Histosols in the northern permafrost region (as defined by ref (64) (64) and 97% within the permafrost extent of another map based on global temperatures and topography (57). Permafrost in peatlands is known to be less sensitive to thaw than mineral soil permafrost (due to the insulating properties of peat (66)). The mapped extent in SoilGrids is thus not entirely unrealistic, but certainly warrants further investigation and ground-truthing against local scale datasets. Ref. (67) reviewed the extent of peatlands within the Siberian Yedoma region and found limited evidence of peat formation there. On the other hand, ref. (68) performed a meta-analyses of local vs regional soil maps and found evidence of underestimation of permafrost peatlands in the regional soil maps in e.g. Arctic Canada.
To assess which maps are more accurate, we evaluated them against local scale observations of peatland and permafrost extent, including the >7000 peat cores (point-to-pixel comparisons) in our database as well as high-resolution and extensively ground-truthed local scale maps which span over several pixels or polygons in the maps. In general, SoilGrida maps a more dispersed peatland distribution, while the national soil maps map denser peat coverage in core regions known for extensive peat formation (e.g. the West Siberian lowlands and the Hudson Bay lowlands) (Fig. S1).
Overlaying the point observations of peat cores (n >7000) on the maps suggest that SoilGrids better captures sites that are located in areas of low peat density while the regional polygon soil maps seem to better describe regions of very high peatland density (Fig. S2). More than 2000 of the sampled peat cores are in locations that show zero peatland cover in the national soil maps (or WISE30sec). Such distribution patterns can partly be explained by the differences in scale between the datasets. Minimum mapping units in national soil surveys will have prevented mapping of small isolated peat complexes which can be mapped by SoilGrids250m. At the same time, SoilGrids is based on machine learning with no described validation and expert input of peatland distributions. This leads us to reason that the maps are partly complementary and that using a per-pixel mean of the products may be the most robust option. For ground-truthing of permafrost extent (fraction of mapped organic soils that are Histels), only the Cstock peat cores (n=ca 700) have data on peatland type. We compared the permafrost extent in the maps against applying MAAT thresholds but find that both maps outperform any MAAT threshold applied. This analyses, and the scientific literature, shows that Histels may occur at MAAT >+1°C, but it is rare and usually occurs in areas where the permafrost persist as relict from previous colder climates (66,69). Because we are interested in broad-scale patterns and less in isolated spots or relict permafrost, we apply a threshold where no Histels are mapped at MAAT >+1°C. This also addresses unrealistic distributions of Histels in some regions from SoilGrid250m. This threshold eliminates some problems with SoilGrids. For instance, there is a definite overestimation of permafrost extent in some regions within the original SoilGrids datasets. This is exemplified by extensive Histel occurrences in central Scandinavia which are known to be permafrost free.
Comparisons, like the ones described above, that uses individual point observations to evaluate coverage in larger aggregated polygons or pixels will always suffer from a scale miss-match (70). We are also able to evaluate the maps against local/regional maps of high accuracy which span over several pixels or polygons (Table S3). The results from the map comparisons are not fully conclusive, but they do show that a mean between the two maps performs better than either of the individual maps in most cases. Taken together these ground-truthing analyses suggest that the maps are to some degree complementary and that a mean of the two maps is the most robust scaling available.

S3: Results on peat depth, stocks of OC and TN
The peat depth, and stocks of peat C and N are characterized by very high spatial variability (Fig. 1c,  d). In all regions where there is a high data-density, there is also very substantial variability in both peat depth and peat C storage. The mean depth in all cores is 234±175 cm (range 40-2,000 cm; n=7,111), the mean peat C storage is 106±66 kg C m -2 (range 0.4-593 kg C m -2 ; n= 782) and the mean peat N storage is 3.9±2.4 kg N m -2 (ranges from 0.3-11.5 kg N m -2 ; n=105). Sites with data on OC and TN stocks are somewhat shallower. This is because the spatial distribution of peat cores is uneven, particularly for sites with OC and N data which are clustered in northern regions (Fig. 1c) with shallower peat. We used linear models to estimate peat C and N storage from peat depth for permafrost-free and permafrost peatlands, respectively. We developed a separate function to estimate peat C storage from peat depth for permafrost-free and permafrost peatlands, respectively (Fig. S3ab). The slope of the linear relationship between peat depth and peat C storage is significantly flatter for permafrost sites than for other peatland types (ANCOVA, p<0.05), which is the reason two separate models were used instead of joining the data. This log-transformed data using Major Axis (MA) linear models (31); intercepts set to zero; n= 334 and 446; slopes 0.85011 and 0.88263; R 2 0.77 and 0.57; both p<0.0001). The low R 2 for scaling of C storage from depth in permafrost peatland is likely caused by the natural variability in ground-ice affecting the C density of frozen peat.
Although the input data is limited and variable, we also developed separate function to estimate peat N storage from peat depth for permafrost-free and permafrost peatlands, respectively (Fig S3c,d; logtransformed data using Major Axis linear models; intercepts set to zero; n= 16 and 85; slopes 0.18162 and 0.31196; R 2 0.35 and 0.42; p<0.05 and p<0.0001). The large spread in this data is reflected in the broad uncertainty ranges of scaled peatland N stocks reported in the main text.

S4: Analyses of observed lateral losses from thawing peatlands in relation to projected chronosequence peat losses
In this sections we outline a brief analyses of how our projections of potential lateral peat C fluxes inferred from chronosequnces of permafrost peatland thaw compare to observations of lateral flux from the field. Datasets that allow quantification of DOC/POC losses from areas with a strong signal of thawing peat (but limited interference of other permafrost thaw signals) are rare. We have found two examples, one from Scotty Creek in Western Canada and one from the West Siberian Lowlands, and they show very different results. The projected average net lateral losses for a +1° C (consistent with the warming signal we could detect with contemporary data) warming scenario from our spatial model is 11 kg C m -2 for the thawed peatland, spread over 100 years. We do not define how this C loss is distributed over the 100 year time-period, but we assume that losses are higher in the earlier period and should be in the range of 100-200 g C m -2 yr -1 . The Scotty Creek catchment in the discontinuous permafrost zone of western Boreal Canada is a well studied site where data on peatland thaw over time, GHG flux data, chronosequence studies and lateral flux data is available. The data from this site does not support our projections as there is little sign of lateral peat C losses despite ongoing and past thaw. Since 1970, the areal coverage of permafrost peat plateaus has declined from 35 km2 to <20 km (total catchment area is 139 km2), a net areal loss of 15% (71). If all these peatlands lost C at the rates projected, it would result in a DOC export from the thawing peatlands alone at the catchment outlet of 15-30 g C m -2 yr -1 . But the DOC export from the catchment is observed to be only ~1.5 g C m -2 yr -1 (72). There is very limited POC export and also radiocarbon data indicating that most of the DOC is modern (72). We thus conclude that at Scotty Creek, there is presently no large lateral flux component that can explain the post-thaw losses. But there are also no large gaseous losses at thawing sites (47), while the chronosequence at this site shows substantial long-term losses which we are unable to constrain with the available flux data.
There is data to support large lateral losses from the West Siberian Lowlands (57). In a transect from no permafrost to continuous permafrost across small peatland dominated catchments. We do not have data to constrain the rate of peatland permafrost thaw in these catchments over recent decade. Even though this region is dominated by peatlands, it is unlikely that the average catchments in that area have lost as much as 15% of the full catchment area as peatlands alone; we assume that they were lower or comparable to those at Scotty Creek along the southern fringes of permafrost. Serikova et al. (57) show a distinct peak in DOC fluxes in the sporadic permafrost zone, 15-20 g C m -2 yr -1 , declining to 6-8 g g C m -2 yr -1 in the continuous permafrost zone. We thus conclude, that in these streams the observational data supports lateral fluxes of the magnitude projected from the limited chronosequence data. There is also strong indirect support for substantial lateral looses of permafrost DOC from peatlands in the permafrost thaw-signal detected in the large Siberian Arctic rivers. Wild et al. (73) make use of time series data, including C isotopic data, to show that the DOC export attributable to permafrost thaw is by far the strongest in the Ob river (with limited permafrost cover but much thawing peatland area), but rather small in the Kolyma river (with much permafrost cover but very little peatland cover). Figure S1. Maps of fractional peatland cover, subdivided into Histosols and histels, in in the SoilGrids dataset (panels a and b) and in the regional polygon soil maps (panels c and d). White land areas in c and d are not covered by the regional polygon soil maps. In the regional soil maps, the combined cover of Histels and Histotols in the northern permafrost region (as defined by Brown et al., 2002) is 2.1 M km2, similar to the extent of Histels in SoilGrids. In SoilGrids, 95% of the northern extratropical peatland extent occurs within the region that is covered by the regional soil maps. Figure S2. Percentiles of the peat core sites plotted against fraction of peatland cover in the different map products.     Figure S10. Curves showing the predicted permafrost fraction as a function of mean annual air temperature (MAAT). The original curve describing mineral soil permafrost (in pink) has been adapted to show predicted permafrost distribution in peatlands. The dotted lines of the green curves represent the spread between the lowest and highest mapped permafrost fractions in the national soil maps and SoilGrids products. Note that the peatland curve predicts ecosystem-protected permafrost occurring at MAAT >0°. The green curve includes a multiplicative factor which allows the curve to asymptote below 100% permafrost at cold air temperatures. This maximizes fit to the map data and implies that even in very cold places, some permafrost-free peatlands exist (due to e.g. local thermokarst or hydrological controls).  Table S1. Table summarizing estimated peatland extent for various geographic regions based on regional polygon soil maps, ref (14) modified in this study), SoilGrids (13), WISE30sec (12) and published peat inventories based on national or regional statistics from tables (74)(75)(76) Table S3. Summary of estimated coverage of peatlands in local study areas compared to SoilGrids (SG) (13), national/regional harmonized soil survey maps (NM) (14) and the mean between the two map products (which is what was used in this paper).  Table S4. Summed of annual fluxes of greenhouse gases. Includes annual stable baseline fluxes from all peatlands at present day after temperature stabilization at +1.5 C above pre-industrial temperatures as well as transient annual fluxes from post-thaw peatlands, including active layer deepening and collapse scars. The standard error (SE) is calculated based on the coefficient of variation from the flux data used for scaling (table S4), for N20 it is the propagated error from N20 and C fluxes since N2O loss is scaled based on C loss.   Table S7. The loss of old permafrost C can be estimated as a function of pre-thaw C stock. Table  shows the relationship between pre-thaw C (initial permafrost peatland C stock) and fraction of prethaw C as a function of time. The fraction of C remaining at 100 yrs is based on a set of five chronosequences (2 syngenetic and 3 epigenetic peatland complexes) from ref (41) (Figure S7), the time dynamics and shifts with C stocks are based on figure 6 from ref. (39) , assuming that the dynamics scale linearly after including additional chronosequences sites.