The limits to northern peatland carbon stocks

Northern peatlands have been a persistent natural carbon sink since the last glacial maximum. If there were no 10 limits to their growth, carbon accumulation in these ecosystems could offset a large portion of anthropogenic carbon emissions until the end of the present interglacial period. Evaluation of the limits to northern peatland carbon stocks show that northern peatlands will potentially play an important role, second only to the oceans, in reducing the atmospheric carbon dioxide concentration to the level that is typical of interglacial periods if cumulative anthropogenic carbon emissions will be kept below 1000 Pg of carbon. 15

We use the words "peat C addition" to denote the amount of carbon that enter to catotelm.
Peat accumulation is the difference between peat addition and peat decomposition.
The words "peat C addition" are changed to "annual C input to Yes, "potential for growth" and "limits to growth" have the same meaning, and "potential for growth" sound more positive.
However, we would like to keep here a "connotation" to the paper "The limits to peat bog growth" by Clymo (1984)  Here we refer to the article of "According to this notion, northern peatlands were providing a persistent but variable sink for atmospheric carbon (Yu, 2011 (Alexandrov et al., 2016)." -P2,

Introduction 20
The recent compilations of peatland data (Loisel et al., 2014) largely confirm the conventional wisdom of the carbon (C) sink provided by northern peatlands (Loisel et al., 2014;Treat et al., 2019) largely confirm the conventional notion of the carbon (C) sink provided by northern peatlands, namely the peatlands distributed across the northern mid-and high-latitude regions located north of 45°N, since the Last Glacial Maximum (Loisel et al., 2017). According to this notion, northern peatlands were providing a persistent but variable sink for atmospheric carbon (Yu, 2011). The variationsVariations in the 25 sink magnitude are explained by changes in the rate of peatland expansion and in the rate of peat accumulation. In the early Holocene, both the rate of peatland expansion and the rate of carbon accumulation appear to be highest (Yu et al., 2010) as compared to the later Holocene periods. However, during the last 5000 years,Since the area of peatlands remained relatively stable in the late Holocene (Adams and Faure, 1998), and therefore the growth in peat depth could be the major cause of the carbon sink provided by (Adams and Faure, 1998;MacDonald et al., 2006;Yu et al., 2010), the major part of the carbon sink provided by northern peatlands during this period could be attributed to the growth in peat depth, not to the growth of the area occupied by the northern peatlands.
There is, of course, no guarantee that the current interglacial will last for another 20,000 years; howeverThere has been little research, however, on estimating the potential magnitude of the cumulative carbon removal associated with the natural development of peatland ecosystems. Individual peatland development may lead to reduction of the carbon sink potential: 10 the closer the peatland ecosystem is to its steady state, that is, to the equilibrium between organic matter production and decomposition, the lower is the carbon sink magnitude. Therefore, the amount of carbon that northern peatlands could remove from the atmosphere will be less than that estimated above.
The process of reaching equilibrium can be conceptualized as follows, see also (Clymo, 1984;Alexandrov et., 2016). Peat is accumulated due to protection of plant litters in the catotelm, the lower layer of a peat deposit that is permanently saturated 15 with water. The plant litters do not enter the catotelm directly, but instead they first enter the upper layer of the peat deposit, the acrotelm, that is not permanently saturated with water. Despite intense aerobic decomposition of organic matter in the acrotelm, at least a small portion of the organic matter that enters the acrotelm always reaches the catotelm in an accumulating peatland. This is, of course, not true in the case of a degrading peatland, but degrading peatlands do not fall within the scope of this study. 20 Precisely speaking, the organic matter does not reach the catotelm, it is rather "flooded" by elevating groundwater. The rise of groundwater is caused by the rise of the peatland surface that in turn results from accumulation of organic matter. This loop cannot elevate the groundwater infinitely. The maximum height of the water table, and thus the potential peat depth, is determined by the amount of effective rainfall, drainage system density (the length of draining streams per unit area) and the hydraulic conductivity of peat and mineral materials below the peat (Alexandrov et al., 2016). 25 The purpose of our study is to estimate the potential peat depth and carbon stocks over NH area north of 45°C and arrive to conclusion about the cumulative carbon removal associated with the natural development of northern peatlands by the end of the current interglacial. Although it is not completely clear how long the current interglacial will last, the recent attempts to estimate its possible duration lead to conclusion that a glacial inception is unlikely to happen within the next 50,000 years if cumulative carbon emissions exceed 1000 PgC (Berger et al., 2016) . Since the duration of the current interglacial depends on 30 the cumulative carbon emissions, it should also depend on the cumulative carbon removal that may offset the effect of carbon emissions, and therefore our study contributes also to the discussion on whether the Earth System would remain in the present delicately balanced interglacial climate state for an unusually long time.
There has been little research, however, on estimating the potential magnitude of the cumulative carbon removal associated with the natural development of peatland ecosystems. Individual peatland development may lead to reduction of the carbon sink potential: the closer the peatland ecosystem is to its steady state, the lower is the carbon sink magnitude. Therefore, the 5 amount of carbon that northern peatlands could remove from the atmosphere might be less than that estimated above.
The process of reaching equilibrium could be conceptualized as follows. Peat is accumulated due to protection of plant litters in the catotelm, the lower layer of a peat deposit that is permanently saturated with water. The plant litters do not enter the catotelm directly, but instead they first enter to the upper layer of the peat deposit, the so-called acrotelm, that is not permanently saturated with water. Despite intense aerobic decomposition of organic matter in the acrotelm, at least a small 10 portion of the organic matter that enters the acrotelm always reaches the catotelm in an accumulating peatland.Precisely speaking, the organic matter does not reach the catotelm, it is rather "flooded" by elevating groundwater. The elevation of groundwater is caused by the elevation of the peatland surface that in turn results from accumulation of organic matter. This loop cannot elevate the groundwater infinitely. The maximum height of the water table, and thus the potential peat depth, is determined by the amount of effective rainfall, drainage system density and the hydraulic conductivity of peat and mineral 15 materials below the peat.

Model equations
To calculate the potential peat depth, we apply an equation derived (see Methods) derivedSupplement) from the impeded drainage model used in our previous study (Alexandrov et al., 2016) . This equation relates the maximum peat depthheight of 20 the water table above the level of the draining system, hmax, at a given watershed to the fraction of its area covered by peatland, fP,obs, and the average depth to bedrock, and thusg: This allows us (see Supplement) to estimate the potential amount of carbon that could be accumulated in northern peatlands from, based on gridded data of soil properties (Batjes, 2016) and depth to bedrock (Hengl et al., 2014)., the potential average 25 peat depth, pd,max, in a grid cell as The gridded data on soil properties give the fraction of a grid cell covered by peatlands. To estimate the fraction of a watershed covered by peatlands, which is needed for calculating the potential peat depth, one should make an assumption about the peatland distribution within the grid cell. We address the uncertainty associated with making such assumptions by giving three estimates of the potential amount of carbon that could be accumulated in northern peatlands: conservative estimate, non-conservative estimate and less-conservative estimate. The conservative estimate assumes uniform distribution of peatlands over all grid cells, the non-conservative estimate assumes clustered distribution over all grid cells, and the lessconservative estimate is derived using a rule-based algorithm categorizing the grid cells into those where peatland distribution is uniform and those where peatland distribution is clustered (see Methods). 5

Methods
The maximum depth of peat that could be accumulated in a watershed is a function of effective rainfall (the difference between precipitation and evapotranspiration), the density of draining system, and the average height of the watershed above the level of draining system. The particular form of this function, derived by using the impeded drainage model (IDM) approach, implies that the maximum carbon stock in a grid cell, pC,max , can be estimated using the following equations: 10 where g is the average height of the watershed above the level of the draining system, in m; d is the maximum depth of acrotelm, in m; where d is the maximum depth of acrotelm, in m (set at 0.4 m), and then to estimate the maximum carbon stock in the grid 15 cell, pC,max as where c is the bulk carbon density of peat, in gC m -3 ; (set at 58 KgC m -3 ); A is the area of the grid cell in m 2 ; fP,obs is the fraction of the area occupied by peatlands; and hmax is the maximum height of water table above the level of draining system, in m. 20

Input data
The values of g at the cells of 0.1º×0.1º geographic grid (Figure 1) were estimated from the data on depth to bedrock provided by SoilGrids (Hengl et al., 2014) (Hengl et al., 2014) . The use of these data for estimating g in permafrost landscapes is somewhat challenging, because the hydraulic conductivity of permafrost could be as low as that of bedrock under some conditions. Due to this reason, we find it more suitable to use the maximum depth of the active layer for 25 estimating g on these landscapes, for example, by setting g at 2 meters for the regions where mean annual temperature is below -2°C, that is, assuming that the southern boundary of permafrost could be approximated by the -2°C isotherm of mean annual temperature (Riseborough et al., 2008) and that the active layer thickness does not exceed 2 m. The latter is an ad hoc assumption based on the recent discussion of uncertainties in the methods for estimating active layer thickness at regional scale (Mishra et al., 2017).
To determine the present-day peatland extent, we relied on the WISE30sec data set (Batjes, 2016) of soil properties at 30' ' 5 resolution. The data set contains a classification of soil type for each mapping unit, and to diagnose peatland extent we determined the fraction of each 0.1°×0.1° grid cell covered by soils of histosol type (soil code HS in FAO90 classification).
These data allow us to estimate the values that fP,obs may take at the cells of the 0.1°×0.1° geographic grid ( Figure 2) and assume that peatlands occupy athe total area of, 2.86 ×10 6 km 2 , that peatlands occupy in the land north of 45ºN.
This estimate of the peatland area does not go beyond the recent estimates (Yu, 2012) (that fall in the range of 2-4 million 10 km 2 ), but it cannot be easily interpreted as the estimate of the actual peatland area. The estimates of the actual peatland area may vary depending on the criteria that are used to distinguish peatlands from other types of land surface. The minimal depth of the peat layer, which is used to classify a land unit as peatland, is the criterion that affects the estimates of peatland area.
Since the data in soil properties do not allow us to evaluate the actual depth of peat layer, it would be better to interpret them as the area that could be potentially occupied by peatlands under the present climate. (Xu et al., 2018). Since peatland extent 15 is diagnosed by the extent of histosols, 2.86 ×10 6 km 2 should be interpreted as an estimate of the area of peatlands with peat depth exceeding 40 cm (according to FAO definition of histosols).
Besides, the use of regular grid for representing the spatial 2.3 Uncertainty associated with peatlands distribution of peatland area imply large uncertainty in theover a grid cell The gridded data on soil properties (Batjes, 2016) give the fraction of a grid cell covered by peatlands. To estimate of pC,max 20 for a given cell.the fraction of a watershed covered by peatlands, fPW, which is needed for calculating hmax, one should make an assumption about the peatland distribution within the grid cell. This problem couldcan be illustrated with the following example. The fact that 36% of a grid cell is covered by peatlands (fP,obs=0.36) may mean that peatlands cover 36% of each watershed within thisthe grid cell (a conservative interpretationfPW=0.36), or that only 48% of watersheds are occupied by peatlands, (fWP=0.48), and theypeatlands cover 75% (0.48×75=36) of each of these watersheds (a non-conservative 25 interpretationfPW=0.75; fP,obs = fWP ×fPW =0.48×0.75=0.36).
Another illustration of the uncertainty associated with interpretation of fP, obs is provided by Table 1, where the estimates of potential peat carbon density in the central part of peatlands are compared to the values observed at 33 peatland sites (Billings, 1987;Borren et al., 2004;Jones et al., 2009;Robinson, 2006;Turunen et al., 2001;Yu et al., 2009). As it can be seen from Table 1, the conservative estimates of the potential peat carbon density, that is, estimates based on the 30 conservative interpretation of fP,obs, are often lower than the actual peat carbon density at the sites that fall within the cells where fP,obs is less than 20%.
The non-conservative interpretation of fP,obs provides a much higher estimate of the potential carbon stocks in peatlands within latitudinal belt between 45° and 84° N than the conservative estimate: 1258 vs 665 PgC. This large uncertainty cannot be easily reduced by using a 1 km grid, because one cannot expect that each watershed falls within one grid cell. However, moving to finer grid is not the only approach for reducing uncertainty in the spatial distribution of peatlands.We address this uncertainty by giving three estimates of the potential amount of carbon that could be accumulated in northern peatlands: the 5 uniform estimate, the clumped estimate and the conductivity-dependent estimate. The uniform estimate assumes a uniform distribution of peatlands over all grid cells (fPW= fP,obs; fWP=1), the clumped estimate assumes a non-uniform distribution over all grid cells (fPW=0.75;fWP=fP,obs/0.75), and the conductivity-dependent estimate is derived using a rule-based algorithm categorizing the grid cells into those where peatland distribution is uniform and those where peatland distribution is nonuniform. The value of the hydraulic conductivity coefficient, K, calculated from the amount of annual precipitation, potential 10 evapotranspiration, fp.obs, and g (see Supplement) could beis used in this algorithm as an indicator of clusterednon-uniform peatland distribution within a grid cell. If itK is above the typical value, Kc, then one may assumeit can be assumed that peatland occupy fp.obs / fp.est fraction of watersheds and cover fp.est fraction of area of each of these watersheds, where fp.est is set at the value that brings K to Kc. Setting The typical values of hydraulic conductivity vary in a relatively wide range. Due to this reason, we set Kc at 157 m yr -1 ( 15 0.5×10 -5 m s -1 ) the value that leads to the estimate (Figure 3) that could be derived fromof the potential carbon stocks in northern peatlands to that implied by the peat decomposition model thatemployed by Yu (Yu, 2011) employed for estimating actual carbon stocks.. This model suggests that the peat accumulationgrowth of carbon stock in peatlands is limited by the ratio of peatannual C addition rateinput to catotelm to the decay constant. Based on the data from peat cores, the peatannual C addition rateinput to catotelm is estimated at 74.8 TgC yr -1 and decay constant at 0.0000855 yr -1 (Yu, 2011). Thus, the 20 potential carbon stock in northern peatlands could be estimated at 875 PgC (74.8/0.0000855=874,853.8 TgC 875 PgC), and due to uncertainty in the peatannual C addition rateinput to catotelm and decay constant may range from 750 to 1000 PgC (see Supplement). Therefore, we set Kc at the value, namely at 157 m yr -1 ( 0.5×10 -5 m s -1 ), that makes the conductivitydependent estimate of the potential carbon stocks in northern peatlands equal to 875 PgC.
The use of this approach to addressing uncertainty is illustrated by Table 1, where the estimates of potential peat carbon 25 density in the central part of peatlands are compared to the values observed at 33 peatland sites (Billings, 1987;Borren et al., 2004;Jones et al., 2009;Robinson, 2006;Turunen et al., 2001;Yu et al., 2009). As it can be seen from Table 1, the estimates of the potential peat carbon density based on the uniform interpretation of fP,obs (fPW= fP,obs; fWP=1) are often lower than the actual peat carbon density at the sites that fall within the cells where fP,obs is low. For example, the actual peat carbon density at site #30, a raised bog that falls within a cell of which 6% are covered by peatlands, is equal to 214 kgC m -2 , whereas the 30 estimate of the potential peat carbon density based on the uniform interpretation of fP,obs is equal to 65 kgC m -2 . This example shows that in this case assuming a uniform distribution of peatlands could be wrong. The clumped interpretation of fP,obs (fPW=0.75; fWP=0.08) gives much higher value of the potential peat carbon density, 1350 kgC m -2 , that, in its turn, may overestimate the potential peat carbon density at this site if the bog covers less than 75% of the watershed area. The conductivity-dependent interpretation of fP,obs (for Kc =157 m yr -1 ) suggests that the bog covers 53% of the watershed area and its potential peat carbon density is equal to 636 kgC m -2 .

Results
The conductivity-dependent estimates of the potential carbon stocks in the cells of 0.1º×0.1º geographic grid for Kc =157 m 5 yr -1 are shown on Figure 3 (in kilotons of C per square kilometer of the cell area).The sum of the potential carbon stocks for all cells north of 45ºN gives the conductivity-dependent estimate of the potential carbon stock in northern peatlands, which is equal to 875 PgC.
Since northern peatlands have already accumulated 547±74 PgC (Yu, 2011), the conductivity-dependent estimate of their potential carbon stock suggests that the total amount of carbon that they could remove from the atmosphere during the period 10 from present to the end of the current interglacial is limited to 328±74 PgC.
The full range of uncertainty for the estimate of the amount of carbon that northern peatlands may accumulate from the start to the end of the current interglacial could be characterised by the uniform and clumped estimates. The former is equal to 665 PgC, and the latter is equal to 1258 PgC. However, our study shows that neither uniform interpretation nor clumped interpretation of the data on peatland extent is applicable everywhere, and hence the most likely range of uncertainty could 15 be narrower than 665-1258 PgC.

Discussion
The limits to northern peatlands carbon stock were estimated here for the first time in the literature, although the methodology for obtaining such estimate were developed more than 30 years ago by Clymo (1984). We adapted this methodology for use at the Earth system scale based on gridded data (Hengl et al., 2014) representing geomorphological 20 aspects of peat bog growth.
We also characterized the uncertainty in the estimate of the limits to northern peatlands carbon stock induced by sub-grid distribution of peatland. This uncertainty cannot be easily reduced by using a finer grid, because it cannot be expected that each watershed falls within one grid cell. Therefore, we elaborated an approach for reducing uncertainty in the spatial distribution of peatlands that allows us to make a conclusion about the most likely value, 875 PgC, for this estimate. 25 Analyzing the uncertainty in the data on present-day peatland extent goes beyond the scope of this study. Improving the accuracy of these data is a well known task tackled by ISRIC, the International Soil Reference and Information Centre, (Batjes, 2016;Hengl et al., 2014), and by networks of peatland scientists such as C-Peat (Treat et al., 2019) and PeatDataHub (Xu et al., 2018). Hence, it might be more important to update the estimates of potential carbon stocks on a regular basis to keep pace with improvements in the accuracy of the data on present-day peatland extent.
The results of our study suggest that even the conservativeuniform estimate of the potential carbon stocks (665 PgC) is still higher than Gorham's (1991) estimate of 455 PgC in the actual carbon stocks of northern peatlands. Gorham's estimate, based on peat-volume approach (Loisel et al., 2014), is the product of the four numbers: mean depth of peatlands (2.3 m), 5 mean bulk density of peat (112×10 3 g Kg m -3 ), carbon content of its dry mass (0.517), and the area of peatlands (3.42×10 12 m 2 ). Our conservativeuniform estimate of potential carbon stocks implies that the potential mean depth of peat could be as high as 4 m for the same values of mean bulk density of peat and carbon content of its dry mass, and for smaller area of peatlands (2.86×10 12 m 2 ). The conservativeuniform estimate is also higher than the Yu's (2011) estimate of actual carbon stocks, 547±74 PgC, based on the time history approach. However, it is lower than the estimate of the potential carbon stock 10 of 875±125 PgC implied by the model of peat accumulation that Yu employed for estimating actual carbon stocks. This latter estimate of 875±125 PgC could be obtained under the less conservative interpretation of the data on soil properties (see Methods). The map of potential carbon density corresponding to this estimate is shown at Figure 3 (Yu et al., 2010), suggesting that northern peatlands in total would accumulate in the future more carbon than they store now.
The highestclumped estimate of the potential carbon stocks at, 1258 PgC that could be obtained within the range of possible 15 interpretations of the data on soil properties, is beyond the range of uncertainty, 760-1006 PgC, in the estimate of potential carbon stocks that could be derived using the Yu's (2011) model of peat accumulation, which ranges from 760 to 1006 PgC.
(see Supplement). Hence, one could find it is reasonable to agree that the estimate of 875±125 PgC, as obtained from two completely independent methods, is the most expedient estimate of potential carbon stocks in northern peatlands, and that 330±200 PgC is the most expedient estimate of the total amount of carbon that they could remove from the atmosphere 20 during the period from present to the end of the current interglacial.
The estimate of potential carbon stocks, 875±125 PgC, corresponds to the present climate, and therefore assumes that the present climate is typical for the present interglacial period. This assumption, however, might not be relevant to the scenarios of dramatic changes in the Earth system, jeopardizing peatlands development. The recent analysis of mitigation pathways compatible with global warming of 1.5°C above pre-industrial levels (Rogelj et al., 2018) shows that holding the global 25 average temperature increase to well below 2°C is difficult but not impossible. To achieve this goal, cumulative CO2 emissions from the start of 2018 until the time of net zero global emissions must be kept well below 1430 GtCO2, (i.e., 390 PgC), that corresponds to 66 th percentile of transient climate response to cumulative carbon emissions (Rogelj et al., 2018; Table 2.2). Since cumulative CO2 emissions through to year 2017 are estimated at 610 PgC (Le Quéré et al., 2018), 1000 PgC of cumulative carbon emissions, the sum of historical (610 PgC) and the future cumulative emissions compatible with validity of the most expedient estimate of potential carbon stocks in northern peatlands. In brief, if cumulative carbon emissions do not exceed 1000 PgC, the northern peatlands play an important role in global carbon cycle recovery.

Discussion
The limits to northern peatlands carbon stock were estimated here for the first time in the literature, although the methodology for obtaining such estimate were developed more than 30 years ago by Clymo ( 1984). We adapted this 5 methodology for use at the global scale and for taking into account geomorphological aspects of peat bog growth represented by the gridded data on the depth to bedrock (Hengl et al., 2014).
Moreover, this estimate corresponds to the present climate, and therefore assumes that the present climate is somewhat typical for the present interglacial period. This assumption, perhaps, is not relevant to the scenarios of dramatic changes in the Earth system that might take place if cumulative carbon emissions exceed 1000 PgC. But if cumulative carbon emissions 10 would not exceed 1000 PgC, the northern peatlands would play an important role in global carbon cycle recovery.
The ultimate recovery of the global carbon cycle from anthropogenic emissions is a long-term process. (Archer, 2005). The current understanding of this process suggests that oceans absorb the majority of cumulative carbon dioxide emission within several centuries, the minor portion within several thousand years, and the remaining part will be removed through weathering of silicate rocks that may take hundreds of thousands of years (Archer, 2005;Archer and Brovkin, 2008;Brault 15 et al., 2017). In plainother words, the larger the perturbation of the Earth system, the lower the chances that the pre-industrial state will be restored in course of the current interglacial.
Including peatlands in the consideration of global carbon cycle recovery allows us to evaluate the level of the Earth system perturbation that would not last too long to "break" the glacial-interglacial cycle. The results of numerical experiments (see Supplement) performed using an Earth system model of intermediate complexity  imply that keeping 20 cumulative carbon dioxide emissions below 1000 PgC essentially reduces the risk of human intervention of natural glacialinterglacial cycle (Figure 4). The northern peatlands are capable to remove in relevant time frame, that is, over the next [5][6][7][8][9][10][11][12][13][14][15] thousand years, the amount of carbon that ocean won'twill not able to remove, and thus to reduce the atmospheric carbon dioxide concentration to the level that is typical of interglacial periods.

Conclusions 25
Northern peatlands accumulate organic carbon and serve as a slow but persistent land carbon sink since the beginning of the current interglacial. If there were no limits to their growth in the absence of anthropogenic or natural CO2 sources to the atmosphere, they could eventually reduce the atmospheric carbon dioxide concentration to the level at which a next precession-driven decline in the summer insolation in the high northern latitudes would trigger the onset of glaciation.
Our study, however, shows that the cumulative carbon removal associated with the natural development of peatland ecosystems is limited. The most expedient estimate of its potential magnitude, 875±125 PgC, was obtained under the assumption that the present climate is somewhat typical for the current interglacial period. Unless future scenarios of changes in the Earth system would leave no room for northern peatlands, the northern peatlands will play an important role in global carbon cycle recovery from anthropogenic emissions. While studies of this process are now focused on the strength 5 and capacity of the ocean carbon sink, our results open a new perspective for the research on global carbon cycle recovery and on the measures needed to protect the northern peatlands as an important element of the Earth's climate system. Data availability. All data used in this study are available from public databases or literature, cited in the Methods section.
The data produced in course of this work are available from Georgii Alexandrov (g.alexandrov@ifaran.ru) upon request. 10 Author Contributions. All authors contributed to the conception of the work, to data processing and to writing of the paper. G.A.A. drafted the manuscript with inputs from V.A.B., T.K., and Z.Y.
Competing interests. The Authors declare no conflict of interests. 15 Table 1. Potential peat carbon density at the central part of peatland estimated under conservativeuniform (PCD1) and nonconservativeclumped (PCD2) interpretation of fP, obs as compared to the observed peat carbon density (PCD0) at 33 peatland sites (Yu et al., 2009