Diverging responses of high-latitude CO2 and CH4 emissions in idealized climate change scenarios

The present study investigates the response of the high latitude’s carbon cycle to inand decreasing atmospheric greenhouse gas (GHG) concentrations in idealized climate change scenarios. To this end, we use an adapted version of JSBACH – the land-surface component of the Max-Planck-Institute for Meteorology’s Earth system model (MPI-ESM) – that accounts for the organic matter stored in the permafrost affected soils of the high northern latitudes. To force the model, we use different climate scenarios that assume an increase in GHG concentrations, based on the Shared Socioeconomic Pathway 5 and the 5 Representative Concentration Pathway 8.5, until peaks in the years 2025, 2050, 2075 or 2100, respectively. The peaks are followed by a decrease in atmospheric GHGs that returns the concentrations to the levels at the beginning of the 21st century. We show that the soil CO2 emissions exhibit an almost linear dependency on the global mean surface temperatures that are simulated for the different climate scenarios. Here, each degree of warming increases the fluxes by, very roughly, 50% of their initial value, while each degree of cooling decreases them correspondingly. However, the linear dependency does not 10 mean that the processes governing the soil CO2 emissions are fully reversible on short timescales, but rather that two strongly hysteretic factors offset each other – namely the vegetation’s net primary productivity and the availability of formerly frozen soil organic matter. In contrast, the soil methane emissions show almost no increase with rising temperatures and they are consistently lower after than prior to a peak in the GHG concentrations. Here, the fluxes can even become negative and we find that methane emissions will play only a minor role in the northern high latitudes’ contribution to global warming, even 15 when considering the gas’s high global warming potential. Finally, we find that the high-latitude ecosystem acts as a source of atmospheric CO2 rather than a sink, with the net fluxes into the atmosphere increasing substantially with rising atmospheric GHG concentrations. This is very different to scenario simulations with the standard version of the MPI-ESM in which the region continues to take up atmospheric CO2 throughout the entire 21st century, confirming that the omission of permafrostrelated processes and the organic matter stored in the frozen soils leads to a fundamental misrepresentation of the carbon 20 dynamics in the Arctic. 1 https://doi.org/10.5194/tc-2020-164 Preprint. Discussion started: 17 August 2020 c © Author(s) 2020. CC BY 4.0 License.


Introduction
High-latitude terrestrial ecosystems are recognised as an increasingly important factor for the global carbon cycle. On the one hand, global warming is expected to increase the vegetation cover and primary productivity -a trend termed Arctic greening 25 (Keenan and Riley, 2018;Pearson et al., 2013;Zhang et al., 2018), which could significantly increases the terrestrial uptake of atmospheric CO 2 (Qian et al., 2010;McGuire et al., 2018). On the other hand, there are large quantities of effectively inert organic matter stored within the frozen soils of the Northern Hemisphere and a significant fraction of these could become exposed to microbial decomposition in a warmer climate. Areas underlain by permafrost -defined by soil temperatures below the freezing point for at least 2 consecutive years -contain between 1100 -1700 Gt of carbon, the largest fraction of which is 30 stored within the frozen part of the ground (Zimov, 2006a;Tarnocai et al., 2009;Hugelius et al., 2014). With the temperature increase in the high-latitudes being about twice as large as the global average (Stocker et al., 2013), the last decades have already seen substantial changes in the permafrost-affected regions. Regional soil temperatures have increased by up to 2 K and there is a pronounced retreat in the extent of permafrost-affected areas combined with an increase in active layer depth, which leaves large quantities of organic matter vulnerable to decomposition (Biskaborn et al., 2019;Stocker et al., 2013;Et-35 zelmueller et al., 2011;Osterkamp, 2007;Shiklomanov et al., 2010;Frauenfeld, 2004;Wu and Zhang, 2010;Callaghan et al., 2010;Isaksen et al., 2007;Brown and Romanovsky, 2008;. Climate change scenarios project the artic temperatures to increase by between 3 K and 8 K until the end of the 21st century (Stocker et al., 2013). Many modelling studies have investigated the resulting decrease in organic matter stored in the 40 permafrost-affected regions and, for the high emission scenarios -corresponding to a temperature increase of 8 K -, the soils are expected to emit around 120 ± 80 Gt of carbon until the year 2100 Schaefer et al., 2014;McGuire et al., 2018). Increasing temperatures also accelerate the Arctic greening trend and it is highly uncertain at which point the carbon release from thawing soils would surpass the additional carbon uptake by vegetation. However, it is generally assumed that the artic ecosystem will turn from a carbon sink into a carbon source within the 21st century (Schaefer et al., 2011;Schuur 45 et al., 2015). The (net) carbon release will further increase the atmospheric greenhouse gas (GHG) concentrations, leading to a positive feedback. Studies indicate, that this feedback will not only notably accelerate the global warming for high emission scenarios -which result in a near-disappearance of the terrestrial near-surface permafrost -but even for the temperature-target of the Paris Agreement (MacDougall et al., 2012;Burke et al., 2017bBurke et al., , 2018Comyn-Platt et al., 2018). 50 It is exceedingly difficult to estimate the Arctic's contribution to future warming. One issue is the timescale on which the carbon would be released from permafrost-affect soils. While local observations indicate that the change processes, affecting the soil carbon emissions, are locally confined and act on very short timescales, large scale modelling studies suggest that the increase in emissions is likely to occur gradually over a timescale of hundreds of years (Schuur et al., 2015). Another important issue is the fraction of carbon that is released in the form of CH 4 rather than CO 2 . Methane is a much more potent GHG 55 (Stocker et al., 2013), and even a small fraction of formerly frozen carbon that is released as CH 4 would increase the respec-tive global warming potential substantially. Methane is produced during the decomposition under anaerobic conditions and these require soils to be water saturated. Hence, future methane emissions are highly dependant on changes in the sub-surface hydrology in permafrost-affected regions (Olefeldt et al., 2012). It is difficult to represent saturated soils at the typical spatial resolution of present-day Earth system models, making it hard to determine the areas in which the decomposition occurs under 60 anaerobic conditions. Furthermore, the hydrological response to permafrost degradation is very complex and there is some disagreement between land-surface models even as to whether high-latitude soils would in general become drier or wetter in the future (Berg et al., 2017;Andresen et al., 2019). Thus, there are comparatively few studies that use large-scale models to investigate the change in soil methane emissions for future warming scenarios Burke et al., 2012;von Deimling et al., 2012;Koven et al., 2015;Oh et al., 2020). 65 The present study aims to improve our understanding of the arctic ecosystem's importance for the global carbon cycle not only by providing additional estimates of the carbon fluxes under a future warming scenario. More importantly, the study's goal is to provide a better understanding of the processes that govern these fluxes -in particular the soil methane emissionsin permafrost-affected regions. The current anthropogenic GHG emissions make it increasingly likely that temperatures will 70 overshoot any temperature target, before atmospheric GHG concentrations could be stabilized at a desirable level (Geden and Löschel, 2017;Parry et al., 2009;Huntingford and Lowe, 2007;Nusbaumer and Matsumoto, 2008). But while many studies have investigated the response of the arctic ecosystem to increasing GHG concentrations, only a few studies exist that investigate its response to a decrease in concentrations (Boucher et al., 2012;Eliseev et al., 2013) and it is still an open question how the high-latitude carbon cycle responds to overshooting temperatures. Thus, we do not only target the system's response 75 to increasing temperatures, but also during a consecutive temperature decline.
Our investigation is based on simulations with the land-surface component of the MPI-ESM1.2 (Mauritsen et al., 2019), the latest release of the Max-Planck-Institute for Meteorology's Earth system model. However, we could not use the standard JSBACH model, as it includes certain parametrizations that are not well adapted to the specific conditions that are characteristic 80 for high latitudes and neither does it account for freezing and melting of soil water, nor for the methane production in the soil.
In the following we will describe the required modifications to the model, together with a more detailed description of the simulations that were performed in the context of this study (Sec. 2). Section 3 details our findings with respect to the soil CO 2 and CH 4 under in-and decreasing temperatures, while section 4 discusses them in the context of the global carbon cycle.

Model
The changes that were made to JSBACH include the implementation of 3 new modules that represent the formation of inundated areas (sec. 2.1.3) and wetlands (sec. 2.1.3) as well as the soil methane production including the gas-transport in soils (sec. 2.1.4).
Furthermore, we adapted the model's soil physics and carbon cycle to include the processes that are relevant for permafrostaffected regions. In JSBACH, the soil carbon dynamics are simulated by the YASSO model, which calculates the decomposition of organic matter at and below the surface considering five different lability classes -acid-hydrolyzable, water-soluble, ethanol-soluble, non-soluble/non-hydrolyzable and a more recalcitrant humus pool (Liski et al., 2005;Tuomi et al., 2011). The decomposition rates are determined by the standard mass loss parameter, which differs between the lability classes, and two factors that 95 account for the temperature and moisture dependencies of the decomposition process. The standard YASSO model does not consider a vertical distribution of the organic matter within the soil and the decomposition rates depend on the simulated surface temperatures, and precipitation rates. This approach works well in regions in which most of the soil carbon is stored close to the surface, but it is problematic for permafrost-affected regions. The vertical carbon transport in these regions is dominated by very effective processes -cryoturbation (Schuur et al., 2008) -and soils can store organic matter in depth of several meters. Thus, the 100 conditions under which this organic matter decomposes are not well approximated by surface temperatures and precipitation rates. To improve the representation of the carbon cycle in permafrost-affected regions, we implemented a vertical structure of the soil carbon pools and calculate the decomposition rates using depth dependant soil temperature and liquid soil water content. Furthermore, we added a simple parametrization to distribute the carbon inputs according to idealized root-profiles and a scheme to account for the accumulation of organic matter at the top of the soil column and the vertical transport due to 105 bio-and cryoturbation.
Structure of the soil carbon pools JSBACH distinguishes between above and below ground carbon pools, a separation that is -in the standard model -only relevant for the computation of the fuel load required by the model's fire module. However, fresh litter at the surface -such as branches or leaves -has very different properties than organic matter that is encompassed in the soil. To be able to account 110 for these differences, the new structure maintains the separation of above and below ground carbon but introduces a vertical discretization of the below ground carbon pools. As we also maintain the conceptual structure of the lability classes, the new scheme represents soil carbon by 5 lability classes on every model soil layer and 4 above ground carbon pools (note that the humus pool does not exist at the surface, see below).

115
The present model version distinguishes between anoxic and oxic decomposition in the inundated and the non-inundated fractions of the grid box (see below) and the soil carbon pools need to be separated accordingly. Here, we do not simulate the respective pools explicitly. Instead we calculate r cin , the ratio between the carbon concentrations in the inundated and the non-inundated fractions, for each of the soil carbon pools after the decomposition is computed. In the consecutive time step, the soil carbon is distributed between inundated and non-inundated carbon pools according to r cin before the decomposition is 120 calculated. For changes in the inundated area, r cin is modified between two calls of the decomposition routine. This approach 4 https://doi.org/10.5194/tc-2020-164 Preprint. Discussion started: 17 August 2020 c Author(s) 2020. CC BY 4.0 License. allows us to separate oxic and anoxic respiration without having to calculate the entirety of relevant processes -such as land cover changes, disturbances, etc. -for two sets of carbon pools.
For technical reasons we chose to represent the soil carbon pools simultaneously on two different vertical resolutions (soil 125 layers). The coarse layering corresponds to the one used to represent soil temperatures and the hydrological processes while the spacing of the finer layers can be chosen freely. The second structure was implemented because we found JSBACH's standard vertical resolution to be too coarse to properly represent the vertical mixing due to cryoturbation, while it is comparatively expensive to represent all soil processes on a fine grid.

Carbon inputs
In JSBACH, the litter inputs are divided into above and below ground litter fluxes, with 70% of the coarse and 50% of the fine litter entering the above ground pools. We maintain this separation but distribute the below ground litter inputs on the vertical soil layers according to vegetation type specific root profiles. Similarly, the below ground carbon inputs that result from disturbances and land-use change as well as root exudates are distributed according to these profiles. The cumulative root 135 fraction, Y , is described by: with z being the depth below the surface. The parameter β is taken from Jackson et al. (1996) and matched to the plant functional types employed by JSBACH. Furthermore, the cumulative root fraction is scaled by a maximum depth, which is limited by the lower of either the model's prescribed rooting depth or the previous year's maximum thaw depth. The latter is 140 done because JSBACH uses a rooting depth that is fixed in each grid box, but we assume that plants do not extend their roots into the perennially frozen regions of the soil.

Transport
The vertical carbon transport in permafrost-affected regions is dominated by frost heave and freeze-thaw cycles (Schuur et al., 2008). However, cryoturbation involves a variety of complex processes that depend on small scale features of the soil and, even 145 though process models exists (Peterson et al., 2003;Nicolsky et al., 2008), these are not applicable on the scales of land-surface models. Thus, we follow the approach of Koven et al. (2009Koven et al. ( , 2013 and described the vertical mixing of soil organic matter as a diffusive transport: with C being the carbon concentration of the lability class lc, D the diffusion coefficient and z the depth below the surface. Similar to Burke et al. (2017a) we use a constant diffusivity -not varying between grid boxes-to represent bioturbation in regions that are not affected by near-surface permafrost. At the surface we use a diffusivity of 1.5 cm 2 year −1 and for the deeper layers we assume the mixing rates to decline linearly with increasing depth up to a maximum depth of 3 meters or up to the bedrock border. In permafrost regions the mixing rates are much larger and vary based on soil conditions. It is assumed that 155 cryoturbation is more effective in wetter soils and when the freezing during winter and the thawing during spring extends over a long periods -weeks to month -during which the soil repeatedly thaws and refreezes. To account for these effects, we assume a maximum diffusivity of 15 cm 2 year −1 which is scaled by two terms representing the (previous year's mean) saturation of the active layer and the number of days in which temperatures crossed the freezing point. At the surface, diffusivity D [cm 2 year −1 ] is given by: where w atl is the saturation of the active layer, N dc0 the number of days per year in which surface temperatures crossed the freezing point and N dc0,ref a respective reference value which was set to 40 days year −1 . For the depth-dependency of the mixing rates in permafrost-affected regions there are two options included in the scheme. Either a constant diffusivity is 165 assumed throughout the active layer -or until the border with the bedrock -, or the mixing rates are assumed to decline linearly throughout the active layer.
The present model structure separates the organic matter into above and below ground pools and the vertical mixing described above is only applied to the below ground carbon. The organic matter that is deposited above the surface needs to be 170 incorporated into the soil before it can be transported into the deeper layers. The separation between above and below ground litter is a mere conceptual one -used to account for the different properties of the organic matter -and the above ground litter occupies the same physical space as the below ground pools representing the top soil layer. Hence, the transfer of carbon from the above to the below ground pools requires a change in properties rather than in space, and there are two ways by which this can happen. The decomposition at the surface turns a given fraction of the organic matter into humus and with this transfor-175 mation we assume a change in physical properties, that transfers the carbon from the above to the below ground pools (hence there is no above ground humus pool). Furthermore, organic matter builds up at the surface in grid boxes in which the long term carbon input at the surface is larger than the respiration rates. Here, we assume the load of organic matter at the surface to affect its properties, as the latter are largely dependant upon the material's bulk density which is reduced under pressure.
Thus, when the load of organic matter exceeds a given threshold -for the present study we choose ≈ 10 kg m −2 -, the excess 180 material is transferred to the corresponding below ground pools. This corresponds to a surface organic layer with a maximum depth of around 15 cm -averaged over the grid box area-, when assuming a litter density of ≈ 75 kg m −3 , which is well within the range of typical organic layer thickness's (Yi et al., 2009;Lawrence et al., 2008;Johnstone et al., 2010) and very similar to the soil organic layer used in the study of Ekici et al. (2014).
With respect to the decomposition rates, k lc , we follow the same approach as the standard YASSO model in which a labilityclass-specific mass loss parameter, α lc , is multiplied by factors accounting for the temperature and moisture dependencies of decompositiond temp and d mois : For the above ground carbon pools we use the parametrizations of the standard model: ; with β 1 = 0.095 and β 2 = −0.0014, where T surf is the surface temperature [ • C] and P the precipitation rate [m year −1 ]. To account for the different decomposition rates under aerob and anaerob conditions we calculate the moisture dependency in inundated areas as (Kleinen et al., 2019): It should be noted that we assume that only a fraction of the above ground organic matter in inundated areas decomposes under anaerob conditions. As discussed above, the above ground carbon pools in the model occupy the same physical space as the below ground pools representing the top layer of the soil column. In reality, however the litter that falls on top of a fully saturated soil column would still decompose aerobically unless there is standing water on top of the surface. Even then it is 200 highly uncertain how much of the litter decomposes under anaerob or aerob conditions as this depends very much on it's shape and on the depth of the standing water -a twisted branch may be located largely above the water while a straight branch would be fully submerged. In the model we deal with this uncertainty be including the fraction of the above ground organic matter that decomposes anaerobically as an input parameter that can be varied between simulations (see below).

205
For the below ground decomposition rates, we evaluated a variety of functions to represent the moisture and temperature dependencies (Sierra et al., 2015), some of which are included as options in the present version of JSBACH. The goal of this evaluation was to establish a combination of dependencies that changes the carbon dynamics in the non-permafrost-affected regions as little as possible, while preserving the organic matter stored within the perennially frozen ground. For this study, we chose the temperature dependency parametrization of the YASSO model in combination with a simplified version of the 210 moisture limitation function used in the CENTURY ecosystem model (Kelly et al., 2000). The temperature and moisture dependencies, d temp (z), d mois (z) and d mois,inu (z) in depth z (z = s) are given by:

215
where T z is the temperature in depth z and w * z represents the relative saturation of the soil, considering only the liquid water content. Note however, that we do not use the saturation of the soil directly, because certain formulations in the model's soil hydrology module prevent the soil moisture from dropping below a certain threshold or to increase beyond the soil's field capacity. In order to account for this, w * z is not given relative to the soils pore space, but relative to the range between the wilting point and the field capacity. Additionally, we apply a subgrid scale distribution of the soil water in order to determine 220 the inundated grid box fraction (see below). Thus w * z does not correspond to mean saturation of the grid box but to the saturation of the non-inundated fraction. In the inundated fraction soils are fully saturated and d mois,inu (z = s) has a fixed value of 0.32, it is assumed however that decomposition in the inundated areas can only occur when the liquid water content in a soil layer (w liq,z ) is larger than the layer's ice content (w ice,z ), even though it should be noted that in reality microbes do not necessarily require free water in the soil to survive and they can maintain viability for thousands of years within frozen soils (Gilichinsky 225 et al., 2003). However, we assumed the activity under these conditions to be negligible.

Permafrost-physics
The representation of the physical, permafrost-related processes in the soil are largely based on the implementation of Ekici et al. (2014). However, there are certain important differences, which will be described in more detail in the following. Most importantly, we adapted the approach to representing soil organic matter from a pervasive organic-top-soil-layer to explicitly 230 simulating the organic matter at the surface and within each of the vertical soil layers. Furthermore, we adapted the formulations of transpiration and the water limitations of plants to account for perennially frozen soils. It should be noted that the model accounts for the heat generated by decomposition (Khvorostyanov et al., 2008b), even though the effects are negligible in all the simulations.

Soil properties 235
The present model version represents the organic matter at and below the surface explicitly and accounts for respective effects on a given soil property, X soil (z), by aggregating the respective properties of organic, X org (z), and mineral material, X min , according to their volumetric fractions, f org (z) and (1 − f org (z)): (12) where ρ c (z) is the mass concentration of carbon at depth z , r c2b the carbon to biomass ratio and ρ org (z) the dry bulk density of organic matter. The estimates of ρ org (z) vary strongly depending on the quality of organic matter and whether it pertains to litter at the surface or to organic matter that is integrated in the soil (O'Donnell et al., 2009;Ahn et al., 2009;Chojnacky et al., 2009). For the present study, we chose ρ org (s) = 75 kg m −3 for above ground organic matter and ρ org (z = s) = 150 kg m −3 245 for the organic matter below ground. Likewise the properties of the organic matter, X org (z), differ between above and below ground organic matter (Peters-Lidard et al., 1998;Beringer et al., 2001;O'Donnell et al., 2009;Ahn et al., 2009;Chojnacky et al., 2009;Ekici et al., 2014). r c2b was set to 0.5.
This aggregation was applied to all soil properties with the exception of the saturated hydraulic conductivity, for which we 250 follow the approach of the Community Land Model (Oleson et al., 2013). Here, it is assumed that connected flow pathways form, once the fraction of organic matter exceeds a certain threshold. These need to be accounted for in the bulk hydraulic conductivity, k sat (z): where f uncon (z) is the grid box fraction in which no conected pathways exist, k sat,uncon (z) the saturated hydraulic conduc-255 tivity in this fraction and k sat,org (z) the conductivity in the grid box fraction in which pathways form.

Soil and surface hydrology
A given fraction of the water within the soil remains liquid even at sub-zero temperatures. In reality, supercooled water exists in the presence of certain chemicals, such as salts, that lower the freezing temperature, but also because of the absorptive and capillary forces that soil particles exert on the surrounding water. The model does not represent the soil's chemical composition 265 and we only account for the thin film of supercooled water that forms around the soil particles, which can be described by a freezing-point depression (Ekici et al., 2014;Niu and Yang, 2006). However, the liquid water is bound to the soil particles and it is questionable whether it is able to move through the surrounding soil-ice matrix. Thus, in the present model version, we assume the supercooled liquid water in the soil to be immobile. As the vertical movement of water requires flow pathways to be available, percolation of liquid water within the soil is inhibited when more than half of the soil's pore space is occupied by ice.
Additionally, the standard model version assumes lateral drainage from the soil at any given depth, which means that soil ice in deeper layers has very little effect on the saturation of the soil column above. In the present model version, we allow drainage only at the bedrock border, which results in permafrost acting as an effective barrier that strongly impedes drainage.

275
Finally, we changed the conditions controlling infiltration at the surface. In the standard model, infiltration is partly temperature dependant, with no infiltration at the melting point. This condition was removed so that infiltration is controlled purely by the saturation of the near-surface soil and the topography within the grid cell.
In JSBACH, transpiration and the plant's water stress are calculated based on the degree of saturation within the rootzone.

280
However, the respective parametrizations become very problematic in the presence of soil ice because they use a fixed parameter -the maximum rootzone soil moisture -relative to which the degree of saturation is calculated. In reality, the rootzone in permafrost-affected regions is confined to depths above the perennially frozen regions of the soil, while in the standard model, the rootzone can not adapt to the permafrost Similarly, bare soil evaporation is determined by the saturation of the top 6.5 cm of the soil -considering only the liquid water content relative to the entire pore space not to the ice-free pore space. Consequently, evaporation can be reduced substantially when there is ice in the top soil layer, despite enough liquid water being present at the surface. In the present model version we deal with this issue by accounting for the presence of ice and computing the saturation of the rootzone and the top soil layer relative to the ice-free pore space.

Wetlands and inundated areas
In its standard version, JSBACH accounts neither for surface water bodies nor inundated areas and, for the present study, we implemented two schemes that represent different aspects of their formation. Note that in the result section we make no differentiation between wetlands and inundated areas because they have a very similar effect on the carbon cycle, in that they both constitute areas in which soil organic matter decomposes under anaerobic conditions. The first scheme simulates the effect of 295 ponding -the formation of wetlands because water can not infiltrate fast enough and pools at the surface, while the second scheme accounts for inundated areas that form in highly saturated soils, due to low drainage fluxes.
The ARNO model, which is used by JSBACH to determine the infiltration rates, does not account for ponding effects, instead all water arriving on the soil surface is either infiltrated or converted into surface runoff (Dümenil and Todini, 1992;Todini, 300 1996). In the present version of JSBACH, we implemented a WEtland Extent Dynamics (WEED) scheme based on a concept developed for the global hydrology model MPI-HM (Stacke and Hagemann, 2012). WEED adds a water storage to the land surface which intercepts rainfall and snow melt prior to soil infiltration and runoff generation. Based on the storage's surface area fraction f pond and depth h pond , evaporation E pond and outflow R pond are computed as Outflow accounts for topography in form of the outflow lag λ pond computed based on the orographic standard deviation σ oro : resulting in an increased outflow when either the storage contains a large amount of water or the orographic variability in the grid cell is high. Runoff is subdivided into direct infiltration and lateral runoff. The former is diagnosed as the soil moisture 310 saturation deficit of the uppermost soil layer for the wetland covered grid cell fraction and directly added to the soil moisture storage. The latter is further processed into surface runoff and soil infiltration according to the standard soil scheme (Hagemann and Stacke, 2015;Dümenil and Todini, 1992). Runoff is assumed to be zero when temperatures fall below the freezing point.
Considering all these fluxes, the water storage S pond changes according to: Due to the coarse model resolution it is not reasonable to quantify f pond for a given storage state explicitly from highresolved topographical data. Instead, we attribute any change in the wetland's water volume V pond = S pond · f pond · A cell to changes in the wetlands depth and extent using the topographical standard deviation of the grid cell: Thus, any change in surface water is divided equally between water depth and extent if the orographic standard deviation of the grid cell equal a given critical orography standard deviation σ crit . Thus, cells with a high orographic variation exhibit rather deep but small inundated fractions, while flat cells result in very shallow but extensive inundated fractions with a strong seasonality.

325
The WEED scheme is able to represent a realistic wetland distribution with extensive wetlands in the high northern latitudes and tropical rainforest regions. An extensive evaluation of the simulated water bodies is beyond the scope of the current study, but will be presented in an upcoming publication (Stacke et al, in preparation).
To determine the extent of inundation areas dynamically, we use an approach based on the TOPMODEL hydrological 330 framework (Beven and Kirkby, 1979). TOPMODEL employs sub-gridscale topographic information contained in the compound topographic index (CTI) to redistribute the grid-cell mean water table, raising the sub-grid-scale water table in areas of high CTI and lowering it where CTI is low. We employ the CTI index product by Marthews et al. (2015) for the CTI index at a resolution of 15 arcseconds to determine the distribution of CTI values within any particular grid cell and thus determine the fraction of the grid cell where the water tale is at or above the surface. A detailed description of the approach is given by 335 Kleinen et al. (2019).

Gases in the soil
The standard version of JSBACH does not differentiate between aerobic and anaerobic soil respiration and, to be able to determine the methane emissions from saturated soils, we implemented the methane model proposed by Kleinen et al. (2019). These simulations cover the historical period -1850 to 2015 -and a scenario period ranging between the years 2016 and 2100.

360
The present study aims to investigate the high latitude's response to increasing and decreasing atmospheric green-house gas concentrations and, because it is often easier to understand the underlying mechanisms when the effects are large, we investigate a high GHG emission trajectory based on the Shared Socioeconomic Pathway 5 and the Representative Concentration uncertainty and our strategy was not to choose the best estimate for the many of the parameters, but rather vary them within the plausible range to capture the uncertainties that are involved in the respective parametrisations.

380
One key factor, determining the processes in the high latitudes, is the treatment of the soil properties, especially the treatment of the organic matter at and below the surface (Lawrence et al., 2008;Ekici et al., 2015;Jafarov and Schaefer, 2016;Zhu et al., 2019). For the study we chose 5 different configurations, which result in substantially different soil thermal and hydrological properties and, consequently, substantial differences in the simulated sub-surface dynamics (Fig. 2 a-c). In 3 of 385 these configurations we use the entirety of adaptations described above and merely vary the properties assumed for the soil organic matter (Note that in one configuration we also prevented the infiltration at sub-zero temperatures and allowed the movement of supercooled water in the soil to be as close to the standard model as possible). The other two configurations are much more similar to the approach used by Ekici et al. (2014). In the approach, the organic matter only has an impact on the soil properties of the first soil layer, while the lower layers have the properties of mineral soil. The difference between these 2 390 configurations is the properties assumed for the organic matter. Another important factor is the nitrogen limitation in the model.
The changes in the model's hydrology scheme increase the leaching of mineral nitrogen in the high latitudes, which reduces the nitrogen availability substantially. The corresponding nitrogen limitations are much higher than in the standard model which has a drastic impact on the simulated vegetation dynamics and decomposition rates (which are also limited by the nitrogen availability). But rather than re-tune this highly uncertain parametrization we performed an additional set of simulations in 395 which the nitrogen limitations were neglected, capturing the range between highly over-and underestimated nitrogen limitations ( Fig. 2 d-f). Furthermore, a key parameter with respect to the methane emissions is the fraction of above-ground organic matter that decomposes anaerobically in the inundated fraction of the grid box. Here, we performed one set of simulations in which we choose a low value -0.4 -and one set of simulations with a high value -0.8 -for this parameter. Finally, we varied key parameters in in the methane module, most importantly the maximum oxidation velocities, which resulted in two sets of simulations which are largely identical but give very different CO 2 and CH 4 emissions. In total, this gives: n sim = n c,soil · n c,nitro · n c,decomp · n c,CH4 , where n sim is the total number of simulations (40), n c,soil the number of configurations with respect to the soil properties (5), n c,nitro the number of configurations with respect to nitrogen limitations (2), n c,decomp the number of configurations with respect to the fraction of above-ground organic matter that decomposes anaerobically in the inundated fraction (2) and n c,CH4 405 the number of configurations with respect treatment of gases in the soil (2).

Initial carbon pools
Determining the initial soil carbon concentrations is very challenging especially for the northern high latitudes where organic matter was stored in the frozen soils under the cold climate during and since the last glacial period (Zimov et al., 2006;Zimov, 2006b;Schuur et al., 2008). Simulations that target the build up of the soil carbon pools in permafrost-affected regions 410 need to cover the carbon dynamics over a similar period (von Deimling et al., 2018). The respective simulations require many simplified assumptions and, because of the extensive timescale, even small uncertainties propagate into substantial differences between simulated and observed carbon pools. Another strategy is to initialize the simulations with observed soil carbon concentrations . These rely on the spatial extrapolations of thousands of soil profiles (Batjes, 2009(Batjes, , 2016Hugelius et al., 2013) and can be considered much closer to reality than any modelling effort that we are aware of.

415
However, this approach has the disadvantage that the carbon pools are not necessarily consistent with the simulated climate, which can result in unrealistic carbon fluxes at the beginning of a simulation. Furthermore, there is only little information on the quality of the soil organic matter, making it very difficult to separate the carbon into the lability classes used by the model.
Here, we choose a combination of the two approaches to achieve some consistency with both, observed soil carbon pools and the simulated climate.

420
To initialize the soil carbon concentrations, we use the vertically resolved, harmonized soil property values provided by the WISE project (Batjes, 2009(Batjes, , 2016. While the dataset only covers the top 2 meters of the soil column -other datasets provide information up to a depth of 3 meters (Hugelius et al., 2013) -it has the important advantage that it is consistent with the FAO soil units which were used to derive the soil properties for the JSBACH model (Hagemann et al., 2009;Hagemann and Stacke, 425 2015). The dataset provides no information on the quality of the organic matter and we distribute the soil carbon among the lability classes according to the pre-industrial equilibrium distribution that is simulated with the MPI-ESM. However, we do not assume the same distribution in all soil layers and make additional assumptions for different lability classes. The highly labile organic matter has a mass loss parameter that corresponds to a reference decomposition time ranging from a few days to a few years and the respective organic matter decomposes before it can be mixed throughout the soil column. Thus we assume 430 that its vertical profile resembles that of the carbon inputs and distribute the highly labile carbon according to an idealized root profile. In contrast the humus pool has a reference decomposition time of several hundred years, allowing it to be well mixed throughout the soil, and we assume a similar humus concentration in all layers.
435 where C f ast (l) is the concentration of highly labile carbon in layer l and C slow (l) the humus concentration. f f ast,sim and f slow,sim are the respective shares in the total soil carbon as simulated with the MPI-ESM. C obs (l) is the observed carbon concentration in layer l, T C obs the total amount of carbon in a given grid box, y(l) the root-fraction and dz(l) the thickness of layer l.

440
In a first step we calculate C f ast with the above formula but apply the condition that wherever C f ast (l) > C obs (l), the excess in carbon was shifted to the nearest layer in which C f ast (l) < C obs (l). In a second step we calculated C slow , iteratively applying the condition that wherever C slow (l) > (C obs (l) − C f ast (l)) the excess was added to the nearest layer in which C slow (l) < (C obs (l) − C f ast (l)). After C f ast (l) and C slow (l) have been determined, C med (l), the concentration of organic matter with a decomposition timescale of tens of years, is calculated as the difference between C obs (l), and the sum of C f ast (l) 445 and C slow (l): C med (l) = C obs (l) − (C f ast (l) + C slow (l)).
To minimize the inconsistency between observed carbon concentrations and simulated climate conditions, we initialize the experiments with the observation-based, present-day carbon pools but start the simulations in the year 1850 at the end of the pre-industrial period. In the high northern latitudes this allows the carbon concentrations within the (simulated) active layer to 450 adapt to the simulated climate conditions during the historical period, while the perennially frozen regions of the soil conserve the observed carbon concentrations. As the soil thermal and hydrological dynamics vary depending the treatment of the soil properties, this initialisation approach results in substantially different soil carbon pools at the end of the spin-up period.

Results
At the beginning of the 21st century, regions that are affected by near-surface permafrost, here defined as featuring perennially 455 frozen soils within the top 3 meters of the soil column (Andresen et al., 2019), contain roughly between 800 and 1400 Gt of carbon (Fig. 2 a, Tab. 1). 375 to 620 Gt(C) of these are located within the active layer, where they are exposed to microbial decomposition, and the resulting soil CO 2 emissions range between roughly 3.0 and 4.5 Gt(C) year −1 . In most simulations, these emissions are not fully balanced by the soil's carbon uptake, resulting in net fluxes of between -0.1 and 0.6 Gt(C) year −1 meaning that, most likely, the high-latitude soils no longer act as a carbon sink at the beginning of the 21st century. This is 460 also the case for the terrestrial ecosystem as a whole, that is when also accounting for changes in vegetation biomass (Fig.   3 a). Almost all simulations estimate the ecosystem's carbon flux into the atmosphere to be positive at the beginning of the 21st century -on average around 0.1 Gt(C) -and none of the simulations predicts the transition from carbon sink to source to occur after the year 2023. The net carbon emissions increase substantially with rising GHG concentrations, and even for the temperature target of the Paris Agreement -global mean surface temperatures about 1.5 K above pre-industrial levels -the 465 simulated net fluxes increase by a factor of 5 -6, while for a temperature rise of 2.5 K the net emissions would be over 10 times larger than at the beginning of the century. Here, the ecosystem emissions exhibit a non-linear and strongly hysteretic dependency on the simulated surface temperatures and the fluxes into the atmosphere are substantially lower after than before the GHG peaks in 2050, 2075 and 2100. The ecosystem's net carbon flux is largely determined by CO 2 exchange between the land -soil and vegetation -and the atmosphere, while methane emissions contribute very little to the overall carbon flux (Fig. 3 b,c).

Soil CO 2 flux and carbon uptake
The CO 2 emissions from permafrost-affected soils are very sensitive to changes in the atmospheric GHG concentration making a substantial increase in soil CO 2 fluxes very likely, should 21st-century GHG emissions follow the SSP5-RCP8.5 scenario.
The fluxes depend on a number of factors that are affected by the atmospheric GHG concentrations, most importantly the 475 changing climate conditions. As a result, the soil emissions exhibit a non-linear dependency on the atmospheric CO 2 concentrations (not shown) but an almost linear dependency on the simulated surface temperatures (Fig. 3 b) -where, very roughly, each degree of warming increases the soil CO 2 fluxes by 50%, relative to the emissions at the beginning of the century. For the temperature target of the Paris Agreement, the soil emissions increase by about 25% to 40% and if GHG concentrations follow SSP5-RCP8.5 until the year 2100, the soil CO 2 fluxes potentially increase by over 150%, resulting in fluxes of roughly 6.5 to 480 13.5 Gt(C) year −1 . The (almost) linear temperature dependency of the soil CO 2 fluxes is also valid for decreasing temperatures and, when the GHG forcing is reversed, the soil CO 2 fluxes decrease on a trajectory very similar to the increase prior to the GHG peak. However, this does not necessarily mean that the main processes governing the changes in soil CO 2 emissions are fully reversible on a decadal to centennial timescale (see below).

485
One reason for the strong increase in soil CO 2 fluxes with rising temperatures is the degradation of near-surface permafrost and the corresponding increase in active layer depth. For a GHG peak in 2100 over 80% of the near-surface permafrost disappear (Fig. 4 a), exposing large amounts of organic matter to conditions which allow for microbial decomposition. The rise in soil respiration, resulting from the increased exposure of formerly frozen carbon, reduces the soil organic matter substantially and for a GHG peak in 2100 the total amount of carbon stored in permafrost-affected soils could be reduced by close to 15% 490 (Fig. 4 b). This corresponds to a loss of roughly 150 ± 50 Gt(C), which is at the higher end of previous estimates Schaefer et al., 2014). Because the soils remain net carbon emitters even after the GHG forcing is reversed, the soil carbon pools continue to decline after the GHG peak, resulting in a carbon loss of up to 20%. When temperatures have returned to the level of the beginning of the century the overall loss in soil carbon amounts to 120 -240 Gt(C), a large fraction of which is irreversible -at least on short timescales -because it pertains to carbon stored within the frozen soils.

d).
Here, the amount of labile active layer carbon starts to decrease even before the GHG peak in 2100 is reached, while the soil CO 2 fluxes continue to increase. Furthermore, the labile carbon in the active layer shows a strongly hysteretic behaviour and after the GHG peak there is substantially less labile organic matter in the active layer than prior to the GHG peak, while there is slightly more stable organic matter. This indicates that, in addition to the overall carbon loss, the permafrost affect soils undergo major compositional changes which should lead to a reduction of the soil CO 2 fluxes.

505
Another important driver of the soil CO 2 fluxes is the carbon input into the soil, which is largely dependant upon the vegetation's net primary productivity (NPP). The NPP, in turn, depends on the atmospheric GHG concentrations directly, as this determines the CO 2 uptake by leaves, and indirectly, through the resulting climate conditions -namely surface temperatures and water availability -and the vegetation distribution -characterised by the type of vegetation and the vegetation cover.

510
Overall, the changes in climate conditions make the high latitudes much more habitable for plants. The surface temperatures in the Arctic increase about twice as fast as global mean temperatures (Fig. 5 a). This does not only extend the growing season but also increases the water availability for plants, because higher soil temperatures cause the soil ice to melt earlier and refreeze later during the year. Together with the increase in precipitation (Fig. 1 d), this raises the plant-available water by as much as 55% (Fig. 5 b). The changes in climate conditions also increase the vegetation cover and facilitate a (relative) shift from grasses 515 and shrubs towards more productive trees (Fig. 5 c). In combination with the direct effect of CO 2 fertilization, the changes in climate and vegetation increase the NPP in permafrost-affected regions substantially, roughly doubling the productivity for the GHG peak in 2100 (Fig. 5 d).
This increase in NPP corresponds to an increase in carbon input into the soil of up to 4 Gt(C) year −1 and, while GHG 520 concentrations increase, about half of the increase in soil CO 2 emissions is balanced by the increase in primary productivity.
After the GHG peak, the NPP is consistently larger than before the peak, mostly because the tree cover remains very high, resulting in substantially larger carbon inputs than prior to the peak. This balances the reduced amount of labile carbon in the active layer, explaining why soil CO 2 emissions can even be larger when temperatures have returned to the levels of the beginning of the century despite there being up to 30% less labile carbon and about 20% less highly labile carbon in the 525 active layer (Tab. 1). Thus, the predominant absence of hysteresis in the simulated soil CO 2 emissions does not mean that the governing processes are fully reversible on short timescales, but that it is the result of two strongly hysteretic factors offsetting each other -before the GHG peak the large CO 2 fluxes are supported by the deepening of the active layer, while it is largely the increase in NPP that drives the post-peak CO 2 fluxes. The larger carbon uptake by plants following the GHG peak also explains the hysteresis of the ecosystem's net carbon flux -with similar soil CO 2 fluxes and a substantially larger NPP the flux 530 into the atmosphere is consistently smaller after than prior to the GHG peak.

CH 4
The methane emissions from permafrost-affected soils behave very differently than the soil CO 2 fluxes. Most importantly, the soils' CH 4 emissions are 3 orders of magnitude smaller than the respective CO 2 fluxes, indicating that methane will play only a minor role in the northern high latitudes' contribution to global warming, even when considering the respective difference 535 in global warming potential. At the beginning of the 21st century the simulated net CH 4 emissions from high latitude soils amount to roughly 9 Mt(C) year −1 -or about 9 Tg(CH 4 ) year −1 . With a global warming potential of 28 times that of CO 2 (Stocker et al., 2013), this corresponds to a CO 2 flux of 0.25 Gt(C) year −1 . Here it should be noted that the spread in the simulated methane fluxes is substantial, but even the largest net CH 4 emission of any of the simulations is around 34 Mt(C) year −1 , which has the warming potential of a CO 2 flux of 0.9 Gt(C) year −1 .

540
One reason for the low soil CH 4 fluxes is the temperature dependency that determines the ratio of CH 4 and CO 2 which are being produced during the anaerobic decomposition of organic matter. For the low temperatures that are characteristic for the high northern latitudes, only a small fraction -on average around 10% -of the organic matter is converted into methane, even though this rate can increase to 50% for more moderate temperatures. Furthermore, the area in which organic matter 545 decomposes anaerobically is comparatively small. The vast majority of all inundated areas are only seasonally flooded so that anaerobic conditions in the soil exist only temporarily, predominantly during late spring and early summer (Fig. 2 g-i). Thus, while there are regions in western Siberia where up to 40% of the surface are inundated during the snow melt season, the average inundated fraction in permafrost-affected regions ranges roughly between 4% and 6%, which increases to about 10% to 12% when only considering the period of April -June. On one hand this means that the amount of organic matter that is 550 decomposed under anaerobic conditions is roughly an order of magnitude smaller than the amount decomposed under aerobic conditions -around 0.2 Gt(C) year −1 compared to 4 Gt(C) year −1 . On the other hand, it means that the largest fraction of the high-latitude soils produces no methane but actually takes up atmospheric CH 4 , oxidising it to CO 2 .
The way the CH 4 emissions react to the changes in atmospheric GHG concentrations also differs substantially from the 555 CO 2 fluxes (Fig. 3 c). There is almost no increase with rising GHG levels and when global mean surface temperatures rise beyond 2.5 K above pre-industrial levels, the emissions even start to decrease. This is very different to the results of previous modelling studies, who found a strong, positive connection between the 21st century temperature rise and methane emissions (Khvorostyanov et al., 2008a;Burke et al., 2012;von Deimling et al., 2012;Lawrence et al., 2015). In large parts, the behaviour of the simulated methane emissions is a result of the dependency of the net CH 4 fluxes on the at-560 mospheric methane concentrations, as the former are determined by the CH 4 gradient between the soil and the atmosphere.
The SSP5-RCP8.5 scenario predicts the atmospheric methane concentrations to double by the end of the 21st century (Fig. 1   b). Consequently, the methane concentrations in the soil have to increase similarly merely to maintain constant CH 4 fluxes.
The same is true for the CO 2 fluxes, however there is an important difference. When the soil-atmosphere-CO 2 gradient decreases due to increasing atmospheric GHG concentrations, CO 2 rapidly accumulates in the soil until an equilibrium with the atmospheric concentrations is reached that allows to respire the CO 2 generated by decomposition. In contrast, it is much more difficult for the CH 4 concentrations to build up within the soils because a large fraction of the methane is constantly converted into CO 2 in the oxygen-rich layers near the surface. Additionally, larger atmospheric methane concentrations increase the CH 4 flux into non-saturated soils that do not produce methane. For atmospheric methane concentrations that correspond to the temperature target of the Paris Agreement, this could potentially lead to the average CH 4 fluxes from permafrost-affected 570 soils becoming negative. For larger atmospheric GHG concentrations, the high northern latitudes could continue to act as a net methane sink even while rising temperatures increase the CH 4 production within the soil.
However, the atmospheric CH 4 concentration can not explain why the methane fluxes are substantially smaller following the GHG peak and why the high latitude soils may remain a methane sink when the forcing is fully reversed and atmospheric is decomposed under anaerobic conditions. On one hand rising temperatures increase evapotranspiration in summer which decreases the inundated areas after the spring snow melt. On the other hand, ice within the deeper layers of the soil acts as a barrier and the soils drain much more readily when this barrier is melted. Rising GHG concentrations also increase precipitation rates by up to 25% ( Fig. 1 d), partly balancing the negative effect that higher temperatures have on the extent of inundated areas. The combined effect of increasing temperatures and precipitation rates is a slight drying of the soils, which has also 590 been found by other models (Andresen et al., 2019). However, the drying of the soil has very little impact on the overall extent of inundated areas (Fig. 6 a). Rather than shrink, the spatial distribution of inundated areas adapts to the changes in climate conditions -with their extent decreasing in lower and increasing in higher latitudes (Fig. 7). Furthermore, the liquid soil water content in regions that feature large wetland areas increases substantially (Fig. 6 b), while the overall water content (including soil ice) declines only slightly (not shown).

595
With a similar extent in the inundated area and increasing temperatures, the soil methane production mainly benefits from the changes in climate. However, the oxic soil respiration in the adjacent non-inundated areas increases even more because, in addition to the effects of rising temperatures, the increase in liquid water content reduces the moisture limitations on the oxic decomposition rates. Furthermore, the vertical distribution of organic matter in the soil changes in a way that also increases the CO 2 production relative to that of CH 4 . The largest fraction of the carbon inputs occurs above or close to the surface. Thus the increase in NPP, resulting from rising GHG concentrations, primarily increases the carbon concentrations at the top of the soil column ( Fig. 8 a, b). At the surface the oxic decomposition rates are much larger than those under anoxic conditions, while this difference is less pronounced deeper within the soil. Consequently, the (relative) increase in organic matter at the surface further increases the difference between the oxic and anoxic respiration rates. When the forcing is reversed, the factors that 605 determine the difference between the decomposition rates do not return to their state prior to the GHG peak -the distribution of inundated areas remains very different (Fig. 7), the liquid water content in the soil remains much higher (Fig. 6 b) and the soil carbon profile still exhibits higher concentrations of organic matter closer to the surface (Fig. 8 c). Consequently, the difference between oxic and anoxic decomposition rates are also larger after than prior to the GHG peak.

610
The difference in the decomposition rates is highly relevant for the soil methane production. The largest fraction of the anaerobic decomposition takes place in seasonally saturated soils, which means that, for a given period, the organic matter in the respective areas decomposes under aerobic conditions. Consequently, the soil methane production depends on both, the anoxic and the oxic decomposition rates, as the latter determine how much organic matter is decomposed during the drier month, hence how much organic matter is available at the onset of inundation.

615
The inundated area is largest during spring and early summer -April to July -with a peak in May and June followed by a sharp decline due to a strong increase in evapotranspiration (Fig. 9 a). The productivity, on the other hand, peaks in July and the decomposition rates and the litter flux even later. Thus, a large fraction of the carbon input into the soils occurs when conditions favour oxic over anoxic decomposition. Increasing temperatures strengthen this effect because the extent of inundated soils in-620 creases during spring -due to larger snow melt fluxes -while it decreases during summer, owing to higher evapotranspiration and drainage rates ( Fig. 9 b). Most importantly, the increase in oxic decomposition rates is far larger than the increase in anoxic decomposition rates. Thus when GHG concentrations increase, the largest fraction of the additional organic matter -resulting from the increased productivity and the deepening of the active layer -is being respired under oxic conditions in summer and fall when the extent of inundated areas is relatively small. Additionally, higher soil temperatures lead to a larger fraction of 625 the fresh litter being decomposed aerobically during winter, resulting in a lower soil carbon availability during and after the following snow melt season. As a result the relative increase in oxic soil respiration in regions that feature large wetland areas is roughly twice as large as the relative increase in methane production (Fig. 6 c,d).
When the GHG forcing is reversed, there is less organic matter within the active layer, especially during summer and fall, but 630 the decomposition rates and the litter flux remain higher than prior to the peak (Fig. 9 c). In case of the aerobic soil respiration, the increase in decomposition rates is so large that the soil respiration during winter and spring is actually larger than prior to the GHG peak, despite the reduced availability of organic matter. This partly balances the reduction in soil respiration during summer and fall and the annually averaged aerobic soil respiration following the GHG peak is only slightly smaller than prior to the peak. In contrast, the anoxic decomposition rates increase just enough to balance the reduced soil carbon availability dur-ing the winter and spring months, while there is substantially less anoxic soil respiration during summer and fall. In case of the anoxic respiration, the effect of the lower soil carbon availability is even more severe because more organic matter -including the litter input -has been respired aerobically during the previous dry period. Thus, because the anaerobic decomposition rates are much lower than the aerobic rates, while both anoxic and oxic respiration (partly) depend on the same carbon pools, the relative changes in oxic and anoxic respiration are very different. In regions that feature large extents of inundated areas, the 640 oxic respiration is only a few percent lower after than before the GHG peak, while the soil methane production is reduced by up to 30% (Fig. 6 c,d).
Still, even the reduced methane production can not fully explain the difference in the net CH 4 emissions before and after the GHG peak. Another important factor is the methane-transport in the soil and how this is affected by the changes in the 645 vertical soil carbon profile. There are two major pathways by which methane is transported towards the surface, one being the diffusive transport through the soil, the other being plant-mediated transport (Note that the model also simulates ebulition, but the respective fluxes can be neglected). Even in saturated soils, the layers near the surface have a high oxygen content and the better part of the methane that diffuses upwards is oxidised within these layers. Consequently, the largest methane fluxes at the surface do not result from vertical diffusion, but from the methane release by (vascular) plants, whose roots absorb methane 650 within the soil. Here, the changes in the vertical soil carbon profile alter the fraction of methane that is transported towards the surface by a given transport mechanism. When the share of organic matter increases close to the surface and decreases in the deeper layers, a smaller fraction of the methane produced in the soil can be absorbed by roots, substantially reducing the respective emissions at the surface. In a addition there is a 5% decrease in cover fraction of grasses, which are the most effective gas transporters (not shown). As a result, the methane emissions by plants decrease by up to 60% when the GHG 655 forcing is fully reversed, which is a reduction roughly twice as large as the relative decrease in the soil methane production ( Fig. 10 a). In contrast, the oxidation rates differ comparatively little before and after the GHG peak ( Fig. 10 b), because the lower methane concentrations in the soil lead to a larger uptake of atmospheric CH 4 . Thus, with a reduced methane production and a larger CH 4 uptake the soil's net emissions decrease substantially, potentially turning the permafrost-affected regions into a net CH 4 sink.

Conclusion
One of the most important factors in the permafrost-carbon-climate feedback is the fraction of soil carbon that is released in the form of CH 4 , which is tightly connected to the question of whether the high latitudes will become wetter or drier in the future (Schuur et al., 2015). Here, many land surface models indicate a drying of the high latitudes (Anderson, 2019) which will most likely constrain the decomposition under anaerobic conditions (Oberbauer et al., 2007;Olefeldt et al., 2012;Elberling et al., 665 2013; Schaedel et al., 2016;Lawrence et al., 2015). In parts, this is confirmed by results of our study and while we did not find a decrease in the extent of areas with saturated soils, we also did not find a significant expansion of the inundated areas in high latitudes despite the pronounced increase in precipitation resulting from the SSP5-RCP8.5 scenario. Furthermore, we could show that the high latitude methane fluxes are strongly limited by the increase in oxic decomposition because most of the soils are only seasonally saturated and the availability of organic matter depends on the respiration during the drier month of the 670 year. But most importantly, we found the methane oxidation in the soil to be the dominant constraint on the soil CH 4 emissions.
Even for present conditions, only a quarter of the methane produced in permafrost-affected soils was actually emitted at the surface. With the atmospheric CH 4 concentration increasing and the vertical methane transport by vascular plants decreasing, this fraction could be reduced to less than 10 % in the future. Because of these limitations, there was not a single yearfrom 20000 years of simulation (500 years for 40 ensemble members) -in which the methane emissions from permafrost-675 affected soils exceeded 55 Mt(C) year −1 . Considering the global warming potential of methane, this corresponds to a CO 2 flux of about 1.5 Gt(C) year −1 , which is an order of magnitude smaller than the largest simulated CO 2 flux -about 14 Gt(C) year −1 .
Thus, our results indicate that the soil methane fluxes in permafrost-affected regions do not constitute an important contributor to the climate-carbon feedback. In contrast, the soil CO 2 emissions are so large that the (terrestrial) arctic ecosystem acts as 680 a source for atmospheric CO 2 , with net emissions of up to 0.4 Gt(C) year −1 at the beginning and up to 3.9 Gt(C) year −1 at the end of the 21st century. This is very different to scenario simulations with the standard version of the MPI-ESM1.2 in which the region continues to take up atmospheric CO 2 throughout the entire 21st century, with the uptake increasing from around 0.005 to around 0.015 Gt(C) year −1 (not shown). These differences confirm that the non-consideration of permafrost-related processes and the organic matter stored in the frozen soils leads to a fundamental misrepresentation of the carbon dynamics in 685 the Arctic.
Here, the differences in simulated net emissions in the high northern latitudes -between the standard and the modified JS-BACH version -are so large that they become highly relevant for the global carbon budget. For the present-day (2007 -2016) the terrestrial carbon sink has been estimated to about 3 Gt(C) year −1 (Quéré et al., 2018) and net emissions from permafrost-690 affected soils of 0.4 Gt(C) year −1 would reduce the land's capacity to take up carbon by more than 10%. It should be noted that the net emissions of 0.4 Gt(C) year −1 may even be a conservative estimate as our study considered only regions that still featured near-surface permafrost in the year 2000, while the area in which the high latitude soils act as a source of atmospheric CO 2 was substantially larger. With the advancing degradation of near-surface permafrost, the respective CO 2 fluxes will become even more important for future carbon budgets and already by 2050 the net emissions from permafrost-affected soils could be on par with (present-day) global land use change emissions.
Despite their importance, the processes governing the carbon dynamics in permafrost-affected regions are not fully taken into account in the present generation of Earth-system models. Substantial advances have been made within the last decade -many land surface models now include some physical and biogeochemical permafrost processes (McGuire et al., 2016;Chadburn 700 et al., 2017) -, but not one of the models that participated in CMIP6 included an adequate representation of the soil physics in high latitudes, while simulating (interactive) vegetation dynamics as well as the carbon and nitrogen cycle. Consequently, model based studies, at present, can merely provide qualitative answers to the question how permafrost-thaw may contribute to global warming (Schuur et al., 2015). Even when the most relevant processes are included, there are large uncertainties regarding the respective parametrizations as well as the initial-and boundary conditions used by the model -many of which 705 are only poorly constrained by observations. In the present study we used an ensemble of simulations in which we varied key parameters within the uncertainty range and while the simulations agreed on the system's relative response to increasing and decreasing GHG concentrations, the spread in the absolute values was substantial between the ensemble members, e.g. the simulated carbon pools at the beginning of the 21st century ranged between 800 and 1400 Gt while the soil methane emissions ranged between 0.5 and 34 Mt(C) year −1 . To be able quantify the impact of permafrost degradation on the climate system, 710 these uncertainties need to be reduced substantially.
Code and data availability. The primary data is available via the German Climate Computing Center's long-term archive for documentation data (https://cera-www.dkrz.de/........, to be specified before publication). The model, scripts used in the analysis and other supplementary information that may be useful in reproducing the authors' work are archived by the Max Planck Institute for Meteorology and can be 715 obtained by contacting publications@mpimet.mpg.de.   concentrations. b) Atmospheric CH4 concentrations. c) Global mean surface temperature relative to the pre-industrial temperature. Note that the model is not forced by surface temperatures directly, but by atmospheric temperatures at a height of roughly 30 m and the surface incoming long-and short-wave radiative fluxes. d) Precipitation rates, averaged over the latitudinal band between 60 • and 90 • North. Grey lines show the forcing according to the SSP5-RCP8.5 scenario. The coloured lines show the forcing-pathways that are used to reverse the forcing to the state at the beginning of the 21st century -after an assumed peak in the year 2025 (blue), 2050 (green), 2075 (yellow) and 2100 (red). In case of temperature (and the surface radiative fluxes) and precipitation rates, the forcing was derived from CMIP6 scenario simulations with the fully coupled MPI-ESM. Shown is the 20-year moving average.