Estimated effect of the permafrost carbon feedback on the zero emissions commitment to climate change

Zero Emissions Commitment (ZEC), the expected change in global temperature following the cessation of CO2 emissions has recently been assessed by the Zero Emissions Commitment Model Intercomparison Project (ZECMIP). ZECMIP concluded that the component of ZEC from CO2 emissions will likely be close to zero in the decades following the cessation of emissions. However, of the 18 Earth system models that participated in ZECMIP only two included a representation of the permafrost carbon feedback to climate change. To better assess the potential impact of permafrost carbon decay on ZEC a 5 series of perturbed parameter experiments are here conducted with an Earth system model of intermediate complexity. The experiment suggest that the permafrost carbon cycle feedback will directly add 0.06 [0.02 to 0.14]◦C to the benchmark ZEC value assesses 50 years after 1000 PgC of CO2 has been emitted to the atmosphere. An additional 0.04 [0 to 0.06]◦C is likely to been added relative to the benchmark ZEC value from the thaw-lag effect unaccounted for in the ZECMIP experiment design. Overall we assess that the permafrost carbon feedback is unlikely to change the assessment that ZEC is close to zero on decadal 10 timescales, however the feedback is expected to become more important over the coming centuries.

The soils of the northern hemisphere permafrost region are estimated to contain between 1100 and 1500 PgC of organic matter (Hugelius et al., 2014), about half of which is held in the perennially frozen zone of these soils (Hugelius et al., 2014). As climate warms and permafrost thaws organic matter in permafrost affected soils is exposed to increased periods of time where temperature is above freezing, and hence to enhanced rate of decay releasing CH 4 and CO 2 to the atmosphere (e.g. Schuur et al., 2015). A recent informal (non-CMIP) model intercomparison exercise quantifying the permafrost carbon 30 feedback estimated a release of carbon of between 74 to 652 PgC by year 2300 under the high end Representative Concentration Pathway 8.5 scenario, with substantially lower release or even gain of soil carbon under mitigation scenarios (McGuire et al., 2018). Thus the permafrost carbon feedback to climate change has the potential the affect the value of ZEC in a fashion that was poorly quantified by ZECMIP.
Uncertainty in projections from Earth system models can be classified into three components: 1) structural uncertainty, 2) 35 parameter uncertainty, and 3) scenario uncertainty (e.g. MacDougall and Knutti, 2016). Structural uncertainty is created from the discrepancy between the system the model is intended to represent and the system the model actually describes (Smith, 2007;Eyring et al., 2016). Examining different models of the same system with standardized forcings through model intercomparison projects such as ZECMIP are the principle way of quantifying structural uncertainty in climate sciences Eyring et al. (2016). Parameter uncertainty is uncertainty about the value model parameter should take on (Smith, 2007). As model 40 parameters are sometimes quantities measurable in the natural world parameter uncertainty is some cases equivalent to measurement uncertainty. In other cases parameters represent an amalgam of natural processes, in such cases defining parameter uncertainty becomes more ambiguous (Smith, 2007). Parameter uncertainty can be quantified with perturbed parameter experiments, wherein ensembles of model variants with parameter values selected from defined probability distribution functions are run under the same experiment conditions (e.g. Forest et al., 2002). Several such experiments have been conducted to assess 45 uncertainty in the permafrost carbon feedback (Schneider von Deimling et al., 2015;MacDougall and Knutti, 2016;Gasser et al., 2018). Scenario uncertainty is created by uncertainty about what humans will do in the future, and is well explored by the coordinated scenario framework of CMIP (Eyring et al., 2016;O'Neill et al., 2017).
For the permafrost carbon feedback uncertainty ranges derived from structural uncertainty and parameter uncertainty assess- In addition to 16 out of the 18 ZECMIP models having no permafrost carbon module, the experimental design of ZECMIP is ill designed to quantify the permafrost carbon feedback. The top-tier idealized ZECMIP experiments branches from the 55 idealized 1pctCO2 experiment where atmospheric CO 2 concentration rises at 1% a year compounded leading to a quadrupling of CO 2 concentration in 140 years (Jones et al., 2019;Eyring et al., 2016). Following this protocol CO 2 and hence global temperatures rise much faster than in the historical trajectory, and since it takes time for permafrost soil to thaw and organic mater within them to decay, the experimental protocol will tend to underestimate release of carbon from permafrost soils (MacDougall, 2019). Here we will use a perturbed parameter ensemble approach to estimate the contribution to CO 2 -only ZEC from the release of carbon from permafrost soils, following the ZECMIP protocol. We will also conduct an experiment following a more realistic CO 2 emission trajectory in order to quantity the thaw-lag effect from the high emission rates of the ZECMIP protocol.  (Weaver et al., 2001). The version of the model used here (version 2.9pf) has full representation of the oceanic and terrestrial carbon cycles (Schmittner et al., 2008;Meissner et al., 2003). The oceanic carbon cycle has representations of ocean carbonate chemistry (Weaver et al., 2001), phytoplankton-zooplankton-detritus ocean biology scheme (Schmittner et al., 2008), and 70 interaction between ocean sediments and alkalinity (Archer, 1996). The terrestrial component is composed of the Top-down Representation of Interactive Foliage and Flora Including Dynamics (Triffid) dynamic vegetation model (Cox et al., 2001;Meissner et al., 2003), a multi-layer representation soil respiration (MacDougall et al., 2012), and a permafrost carbon module (MacDougall and Knutti, 2016).
The terrestrial subsurface of the model is composed of 14 layers, reaching a total depth of 250 m (Avis et al., 2011). The top 8 75 layers (10 m) are active and the hydraulic cycle and deeper layers are impermeable bedrock (Avis et al., 2011). The freeze-thaw physics of the soil accounts for the effect of soil valence forces on freezing point and the frozen and unfrozen fraction of the soil water is calculated using equations the minimize Gibbs free energy (Avis, 2012). The top 6 layers of the model (3.35 m) are active in the carbon cycle. Carbon is assigned to soil layers from Triffid based on the root density in each soil layer, with remaining dead plant matter added to the top soil layer (MacDougall et al., 2012). Root density varies by plant function type 80 and the temperature of the soil layer (roots do not grow in frozen soil) (MacDougall et al., 2012). In model grid cells where permafrost exists (where soil layers have been below 0 • C for two or more consecutive years) a diffusion based cryoturbation scheme is used to redistribute soil carbon in the soil column. The scheme was originally developed by Koven et al. (2009)  In the version of the UVic ESCM used here (MacDougall and Knutti, 2016) permafrost carbon is a separate carbon pool.
Permafrost carbon is created when carbon is advected across the permafrost table by the cryoturbation scheme and can only be destroyed by being respired into CO 2 . The pool is characterized by a decay rate constant (κ p ), a fraction of the pool that 90 is available for decay (available fraction, A f ), and a passive pool transformation rate (κ tf ), which is the rate at which the passive permafrost carbon becomes part of the available fraction. The available fraction is essentially the combined size of the fast and slow carbon pools as conceptualized in incubation experiments (Schädel et al., 2014;MacDougall and Knutti, 2016). This scheme accounts for the large fraction of permafrost carbon that is very resistant to decay (Schädel et al., 2014), while still allowing the pool to decay over millennial time periods (MacDougall and Knutti, 2016). A fourth parameter, the saturation 95 factor (S) from the cryoturbation scheme allows the size of the permafrost carbon pool to be tuned (MacDougall and Knutti, 2016).
The model experiments here use the permafrost carbon variant of the UVic ESCM 2.9 detailed in MacDougall and Knutti (2016). A newer version of the UVic ESCM (version 2.10) is now available (Mengis et al., 2020). We use the older version of the model to allow for the use of legacy code and legacy model spin-ups from MacDougall and Knutti (2016). Note that the 100 terrestrial component of UVic ESCM 2.10 was taken from the version developed for MacDougall and Knutti (2016), and thus the terrestrial components of the model versions are virtually identical (Mengis et al., 2020).
By changing the flow of outgoing longwave radiation to space as a function of global surface temperature anomaly the climate sensitivity of the UVic ESCM can be altered (Zickfeld et al., 2009). Similarly by changing the meridional diffusivity of the atmosphere within the model, arctic amplification can also be altered (Fyke et al., 2014).

Perturbed Parameter Experiments
To assess the uncertainty in the strength of the permafrost carbon cycle feedback to climate change MacDougall and Knutti (2016) generated 250 variants of the UVic ESCM by perturbing six model parameters. Four of these parameters control the size and susceptibility to decay of the permafrost carbon pool and two (climate sensitivity and arctic amplification) are physical climate parameters. The four permafrost carbon parameters are: 1) the permafrost carbon decay constant; 2) the available 110 fraction; 3) the passive pool transformation rate; and 4) the permafrost carbon saturation factor -which controls the size of the permafrost carbon pool. The permafrost carbon decay constant controls how fast available permafrost carbon can decay given the temperature and moisture of the soil. The available fraction is the fraction of permafrost carbon that is allowed to decay, effectively the fraction of permafrost carbon that is unprotected or weakly protected from decay. The passive pool transformation rate is the rate at which highly protected permafrost carbon becomes weakly protected. The Probability Density 115 Functions (PDFs) for the permafrost carbon decay constant and the available fraction were taken from the meta-analysis of permafrost carbon incubation experiments conducted by Schädel et al. (2014). The passive pool transformation rate is constrained primarily by the non-existence of a remnant mid-latitude permafrost carbon pool from the last glacial maximum, yielding an estimated value of 0.25×10 −10 to 4×10 −10 s −1 , with a best guess of 1×10 −10 s −1 (Trumbore, 2000;MacDougall and Knutti, 2016). Fortunately, the analysis of MacDougall and Knutti (2016) showed that the passive pool transformation rate 120 has only a weak effect of the permafrost carbon feedback on decadal and centennial time scales. The estimated uncertainty in the size of the permafrost carbon pool was taken from Hugelius et al. (2014). For each model variant 5000 year model spin-ups were conducted with year 1850 CE radiative forcing to bring the permafrost carbon into equilibrium, representing an investment of 1.25 million model-years of simulation time. We have re-used these model spin-ups for the present study.
MacDougall and Knutti (2016) also perturbed two physical model parameters, climate sensitivity and arctic amplification.

125
These parameter do not affect the model spin-up as both effect deviations from the pre-industrial climate, thus they can be changed for the present study. To my knowledge there has been no major update in the uncertainty range of arctic amplification PDFs as previous papers (e.g. Olson et al., 2012;MacDougall and Knutti, 2016), a product of two Normal Inverse Gaussian functions. To get new parameter values for the climate sensitivity PDF a Monte-Carlo method was used to fit the function to the distribution outlined by Sherwood et al. (2020): a 5th to 95th percentile range of 2.3 to 4.7 • C a 66% range of 2.6 to 3.9 • C and a median value of 3.0 • C. The new parameter values for the PDF are given in Appendix A.

Model Experiments
To quantify the effect of the permafrost carbon feedback on ZEC we have three key questions: 1) How much warming will the permafrost carbon feedback add to ZEC? 2) What is the magnitude of the thaw-lag permafrost effect from using the 1pctCO2 experiment to quantify ZEC? And 3) How sensitive is the permafrost carbon feedback contribution to ZEC to total CO 2 emitted 140 before cessation of emissions?
To answer the first question we would ideally compare simulations with and without permafrost carbon that are otherwise identical. In the UVic ESCM framework we can create a version without permafrost carbon by setting the cryoturbation diffusion parameter to zero during model spin-up. Without cryoturbation there will be no permafrost carbon pool and active layer carbon pool will also be reduced in size. However, the presence of carbon in soils subtly changes soil thermal and   SSP4-6.0 is the lowest emission SSP that reaches 1000 PgC whilst maintaining near-peak CO 2 emissions. Thus using SSP4-6.0 maximizes the time need to reach 1000 PgC and therefore is optimal for assessing the thaw-lag effect. Under SSP4-6.0 CO 2 emissions 1000 PgC is reached in year 2067, allowing permafrost the time to thaw and the organic matter within it to decay.      c. d.  change the expected value of ZEC of decal time-scales. Thus, the overall conclusion that ZEC will be close to zero in the decades following cessation of emissions remains unchanged.

220
Here we found that the permafrost carbon feedback contribution to ZEC was insensitive to cumulative CO 2 emissions, despite a larger release of CO 2 to the atmosphere from permafrost soils in the A3 (2000 PgC) experiment. The linear relationship between cumulative emissions of CO 2 and global temperature change is generated by an atmosphere-ocean phenomena (Mac-Dougall, 2017), thus the atmosphere-ocean system does not distinguish between emissions from the terrestrial biosphere and fossil fuel emissions (Simmons and Matthews, 2016). Therefore this insensitivity of temperature to carbon released from per-225 mafrost soils appears anomalous. However, for intermediate complexity models like the UVic ESCM the cumulative emissions of CO 2 versus global temperature change curve is only approximate linear and the change in temperature with a unit of CO 2 emitted declines at high cumulative emission totals (MacDougall and Friedlingstein, 2015). Thus the emissions of CO 2 from permafrost soil is less effective at warming after 2000 PgC of CO 2 has been emitted than when 1000 PgC has been emitted in the UVic ESCM. However, it has been shown the full Earth system models do not have non-linear cumulative emissions 230 of CO 2 versus global temperature change curves (Tokarska et al., 2016). Therefore the insensitivity of the permafrost carbon feedback contribution to ZEC to cumulative CO 2 emissions found here should be treated with caution.
To date no Earth system model (McGuire et al., 2018) accounts for abrupt thaw processes in permafrost systems. These processes including thermokarst production, active hill slope erosion, and coastal erosion, could accelerate thaw processes by 40% over the next few centuries (Turetsky et al., 2020). Additionally the 2.9pf version of the UVic ESCM (used here) and  methane production scheme has recently been added to a newly developed thread of the model which preliminarily suggests a warming effect from CH 4 production from permafrost soils of 0 to 0.24 • C, depending on parameter values and scenario followed (Nzotungicimpaye, 2021). Such values are consistent with expert assessment of (Schuur et al., 2015). Accounting for these processes will likely increase the estimated effect of the permafrost carbon cycle feedback on ZEC and therefore the 240 effect of the permafrost carbon feedback on ZEC should be reassessed when these processes are better accounted for in Earth system models.
UVic ESCM 2.10 was one of the two models that participated in ZECMIP that included a permafrost carbon scheme (the other was CESM). The ZEC 50 years after CO 2 emissions cease for the A1 experiment (1000 PgC (Mengis et al., 2020), with the newer version of the model having substantially improved ocean dynamics and a state-of-the-art representation of ocean biogeochemistry (Mengis et al., 2020). Ocean heat and carbon 250 uptake are two of the processes the determine the value of ZEC (MacDougall et al., 2020), therefore it is not unexpected that differences in the representation of the ocean would change the ZEC value of a model.

Conclusions
Here though the effect of abrupt thaw remains unaccounted for and the feedback is of greater concern over longer timeframes. The functional form of the climate sensitivity PDF was taken to be the product of two Normal Inverse Gaussian functions 265 following Olson et al. (2012). The Normal Inverse Gaussian function is: where µ is location, α is tail heaviness, β is an asymmetry parameter, δ is a scale parameter and K 1 is a modified Bessel function of the third kind (Olson et al., 2012). The Python scipy.stats software package was used to compute the PDF. shown in Table A1 below.