A Permafrost Implementation in the Simple Carbon-Climate Model Hector

Permafrost, soil that remains below 0 ◦C for two or more years, currently stores more than a fourth of global soil carbon. A warming climate makes this carbon increasingly vulnerable to decomposition and release into the atmosphere in the form of greenhouse gases. The resulting climate feedback can be estimated using Earth system models (ESMs), but the high complexity and computational cost of these models make it challenging to use them for estimating uncertainty, exploring novel scenarios, and coupling with other models. We have added a representation of permafrost to the simple, open-source 5 global carbon-climate model Hector, calibrated to be consistent with both historical data and 21st century ESM projections of permafrost thaw. We include permafrost as a separate land carbon pool that becomes available for decomposition into both CH4 and CO2 once thawed; the thaw rate is controlled by region-specific air temperature increases from a pre-industrial baseline. We found that by 2100 thawed permafrost carbon emissions increased Hector’s atmospheric CO2 concentration by 10-15% and the atmospheric CH4 concentration by 10-20%, depending on the future scenario. This resulted in around 0.5 °C of additional 10 warming over the 21st century. The fraction of thawed permafrost carbon available for decomposition was the most significant parameter controlling the end-of-century temperature change and atmospheric CO2 concentration in the model and became increasingly significant over even longer timescales. The addition of permafrost in Hector provides a basis for the exploration of a suite of science questions, as Hector can be cheaply run over a wide range of parameter values to explore uncertainty and easily coupled with integrated assessment models to explore the economic consequences of warming from this feedback. 15

I enjoyed reading this manuscript and believe it will be of broad interest to the scientific community as well as informative to IAMs. The Discussion section does a great job at documenting the model limitations (as a field ecologist, I appreciate that); I also think the section that compared this study's model outputs with that of other models was beneficial to the reader.
We appreciate the reviewer's interest and overall positive assessment.
I'd be interested to see a few more details pertaining to: (1) How is the "static" (non-labile) C fraction in the permafrost determined? (2) Why is the CH4 emission fraction from thawed permafrost set to 2.3% (please add references or some kind of explanation), and how much effect would a lower or higher fraction have on GHG (provide graphs)? (3) How sensitive is the model to the permafrost pool size (provide graphs of projected GHG under different pool sizes)? 1) We have added the following sentence to better explain where the non-labile fraction comes from: "While in reality turnover times of soil organic carbon fall anywhere along the range from a few days to thousands of years, we group soil decomposition into fast and slow pools, following previous work (Shädel et al. 2014), where the slow fraction decomposes on the order of up to thousands of years and is assumed to be inert for the purpose of this analysis." 2) We revised the relevant sentence to read: "We set the share of CH 4 to be 2.3% of total emissions, a value derived based on expert assessment in Schuur et al. (2014)." 3) We have added the following figure to visualize the sensitivity of the model to the potential ranges of the key parameters in our sensitivity analysis, including the permafrost thaw parameter μ, the non-labile fraction, the methane fraction of heterotrophic respiration, and the initial permafrost carbon fraction. (Note that the exact values on this figure will change in the revised manuscript but we include it here for reference). gmd-2020-377 Minor comments: If possible, add Hugelius et al 2020 (PNAS;https://www.pnas.org/content/117/34/20438) to Table 4, and discuss their findings in light of yours.
As our results do not separate out peatlands, it is not straightforward to directly compare our results to this paper in Table 4, but we have added the following to our discussion addressing the implications of Hugelius et al. 2020 for our results.
"Abrupt thaw is also a key process for permafrost in peatland soils, and a recent analysis estimates an additional 40 Pg of permafrost carbon stored in peat than had been found previously (Hugelius et al. 2020). Based on our sensitivity analysis, increasing the initial permafrost by this amount might translate to a little over an additional 1% increase in overall temperature change by 2100 in Hector." l-39: "models" is there twice gmd-2020-377 Thanks for the catch! Corrected.
l-181: fix the typo in the word "parameter" Corrected.

RC2 comments
This paper discusses the implementation of a simple permafrost carbon module within the Hector simple climate model. The authors talk about uncertainties but the results have very few details of their impact.
We appreciate the reviewer's careful examination of our paper and thoughtful comments here and below.
We have added a figure (new Figure 6 -above) showing the impact of parameter uncertainties on global mean temperature; this output variable is integrative, capturing the effects on other key outcomes as well (permafrost thaw, methane and carbon dioxide emissions). We have also expanded further on uncertainties in the discussion (described further in various responses below).
In particular, the model results have a high permafrost carbon feedback temperature compared to other results in the literature. I think this is fine if the uncertainties are made more high profile throughout the whole document particularly plumes in Figure 3 but also maybe table 3. Q quick thought -do the authors think this high feedback temperature is caused by a relatively large methane contribution?
We agree that the temperature response in this model is high compared to many existing estimates. In our revised manuscript we update to a substantially higher value of the non-labile fraction, 0.72 -based on the more recent and comprehensive Shädel et al. (2014), compared to the previous value of 0.4, which ends up bringing our temperature response down to between 0.1 and 0.2 ℃, depending on the scenario. We also will tune the model parameterization of permafrost thaw versus temperature to CMIP6 data in the revised manuscript, which will provide a more robust foundation for our default configuration of this model. We note that one of Hector's strengths is in its ability to calibrate to different assumptions and explore their consequences, so that our default configuration should only be seen as a starting point for looking at the effects of permafrost in this model.
We additionally add more discussion about the comparison between our temperature estimates and those from previous results throughout the paper. We think this higher temperature impact is driven by both the inclusion of methane, which not all previous studies have done and is responsible for ~25% of our end-of-century warming, as well as our use of a relatively high warming factor parameter. We have also added the following sentence to compare the estimated methane-driven warming to that in previous studies: "The methane contribution to gmd-2020-377 permafrost-driven temperature change estimated by Hector was around 25%, higher than the average of previous estimates (16%) found in Schaefer et al. (2014), but within the range of the studies included in that analysis." I think more details and discussion are required about the physical permafrost change parametrisation. This is a key addition to the model and not explored in much detail. I see it is discussed more in Figure 2 but I think it needs to be compared to something other than Kessler.
This is a good point, and we have expanded Figure 2 (below) and the discussion to include other projections of high latitude temperature change versus permafrost thaw fraction.

Why is a volume fraction not used when considering physical permafrost change? This is much more relate-able to carbon amount.
Because of Hector's structure and global spatial scale, the permafrost component is based only on the movement of carbon by mass, with no explicit representation of depth or area, so this is most accurately thought of as a permafrost mass fraction. We have clarified this in the methods section.

It would be good to have units included for all variables. Also check all acronyms are defined.
In general the paper contains most of the required information but sometimes later than I would like it.
We have carefully gone through the paper and checked signs and added units and acronym definitions in all the places we found they were missing.
Minor comments: Line 3 -permafrost C feedback is hardly ever estimated using ESMs currently.
We have edited this sentence to read: "The resulting climate feedback can be estimated using land surface models, but the high complexity and computational cost of these models make it challenging to use them for estimating uncertainty, exploring novel scenarios, and coupling with other models."

Line 6 -?? ESM permafrost estimate
As we tuned to several different values, we have left the details of these estimates for later in the manuscript in favor of simplicity within the abstract.
Line 10 -0.5 degree feedback temperature is relatively high.
It is, and while diverging from existing estimates is not necessarily cause for concern, some changes we have made to our default parameterization have ended up bringing down this temperature feedback. See response above for details.

Line 30 -JULES -Joint UK Land Environment Simulator
Added.

Line 31 -remove both 'can' s
Thank you for catching the typo! Corrected. gmd-2020-377 Line 35 -these more complex models still have missing/incorrectly parameterized processes and would definitely benefit from uncertainty quantification.
We have added the fact that these models do also benefit from uncertainty quantification to the sentence as follows: "While high complexity models benefit from uncertainty quantification, they require large numbers of inputs and are computationally expensive, making it difficult to do uncertainty analysis directly with these models." Line 39 -remove 'models' Corrected.
Line 48 -The paragraph above suggests that these simple climate models do not have good spatio-temporal resolution but then this line suggests that Hector will be used to evaluate regional impacts.
Unlike an ESM, Hector is not spatially explicit. However, Hector does have the ability to account for some spatial heterogeneity in land surface processes by dividing the global land surface into biomes, each with their own parameters, pools and fluxes. We have adjusted this sentence for clarity. It now reads: "Including a representation of permafrost in this model will allow for the consideration of permafrost in future such analyses with Hector, and, thanks to Hector's ability to represent separate biomes or regions, will be particularly important to evaluating the specific impacts of climate change in high latitudes."

Line 71 -Equation 1 defines a land-to-atmosphere flux but Equation 3 suggests an atmosphere-to-land flux. Line 71 -I think the last sign in Equation 1 should be + and the definition of FL should specifically state that NPP and Respiration act in opposite directions Line 73 -FL is the difference between NPP (uptake) and respiration (loss).
The land flux should be understood as positive into the land, and we have added language to clarify this. We have additionally corrected the language and sign in the F L equation and description to make it clear that NPP and RH act in opposite directions. We have kept the negative sign in Eqn 1 to be consistent with an understanding of this term being defined as positive into the land.

Line 78/79 -what are the 1/4 and 1/50 factors for?
These are the annual fractions of respiration carbon transferred to detritus and soil, respectively. We have replaced these with f rd and f rs and defined these terms in the paragraph following these gmd-2020-377 equations.

Equation 4 has an air temperature change as a power and Equation 5 -has a running mean air temperature as a power? That looks a bit odd?
Equation 4 is merely the canonical Q 10 response of biological processes (in this case, heterotrophic respiration from the land to the atmosphere) to temperature (see for example Davidson and Janssens 2006https://onlinelibrary.wiley.com/doi/abs/10.1111/j. 1365-2486.2005.01065.x). Equation 5 is intended to provide a buffered response of soil temperature (and thus decomposition) to air temperature changes, and we agree that it's a highly empirical and only an approximation. We now note in the discussion that this should be an area of future focus and development for the Hector model.

Line 85 -Where does the 200 years come from? Can this be justified further?
The 200 year value for this smoothing is somewhat arbitrary. We have added this language to the sentence as follows: "Detritus respiration increases with region-specific air temperature change (T[i,t]) while soil respiration increases with the 200-year running mean of air temperature (T 200 [i,t]), a somewhat arbitrary choice of smoothing used in Hector as a proxy for soil temperatures." See also response above; we have added a note about this in the discussion.

Line 99 -Equation 6/7 -DCperm is added to both Cperm and Cthawed? Surely it should be subtracted from one and added to the other?
We have corrected Equation 6 to subtract this term.

Line 99 -What is Fthawed-atm? Is the sign correct?
We have corrected the sign and added the following language to explain this term: "F thawed-atm is the flux of carbon, in Pg C, from the thawed permafrost pool to the atmosphere, including both CO 2 and CH 4 emissions."

Line 103 -shouldn't PHI be a volume fraction?
Phi is a mass fraction of permafrost carbon, as Hector does not directly compute permafrost area or volume. We have added that language for clarity.
Line 109 -Please give more details on Equation 9 and explore its validity.
We have added the following language to the methods section: "The lognormal CDF was chosen for several reasons. Its curvature captures the "activation energy" of permafrost thaw with respect to temperature for low temperature change (left side of the curve), and, more importantly, the "diminishing returns" of permafrost thaw at higher temperatures because the more accessible near-surface permafrost has already thawed by that point. Additionally, its parameters are readily interpretable in terms of the timing of 50% permafrost loss (exp(mu)) and the rate of permafrost loss around the 50% point relative to earlier/later in the process (sigma), which facilitates the use of this framework to emulate global permafrost dynamics in more complex models. Finally, it is naturally bounded between 0 and 1, which is appropriate as a model of the remaining permafrost fraction.
There are a variety of possible choices for this functional form and others can be explored in future model development efforts. Fortunately, the modular design and coding best practices of Hector make it simple to substitute this equation with alternatives." We will also compare this relationship to that derived from CMIP6 data in our revised manuscript.

Line 139 -please state earlier that the 308 Pg C soil carbon is 'non-permafrost' carbon
We have clarified the land pools that are created when permafrost is added, specifying that soil is non-permafrost soil C at the beginning of the permafrost methods section and have specified non-permafrost soil when it is mentioned after this point. We have not checked the model's sensitivity to these parameters as we have focused our sensitivity analysis on parameters we have added to the model, such as the methane fraction and temperature controls on permafrost thaw. In our revised manuscript we will use CMIP6 outputs to derive ranges for these numbers. gmd-2020-377 We have corrected this name to read f CH4 throughout the text.
In our revised manuscript we will update to use CMIP6 data instead of CMIP5. We will update Table 2 to use our own analysis of CMIP6 data in our revised manuscript.

Line 155 -why do we expect mu/sigma to fall within the range of Kessler 2017?
To clarify, we used Kessler 2017 as a range for these two parameters, as the reviewer notes, and estimated the mu and sigma within this range that produced a best-fit result relative to CMIP5 (CMIP6 in the revision). We could of course let the parameters exceed this range, but feel that this nonetheless provided a useful a priori guide for this exploratory work.

Line 159 -this 70 % was not included in the uncertainty range. Any reason?
We have updated our default value of the static fraction to 72% (Schädel et al. 2014) and re-run our analysis. Our range is now defined around this baseline based on the Schädel paper.
Line 161 -uncertainty ranges in Table 1 do not reach the 4.3 % suggested in the text. Why not?
Thank you for catching this. We have updated the table and added the Schädel citation.

Line 178 -check name f_RH_CH4
Corrected to be consistent with previous equations.
Line 179 -how were the parameters sampled from prior distribution? Latin Hypercube?
They were uniformly sampled. We have added language to specify this.

Line 181-184 -please give more details on what these parameters mean and what the approach of LeBauer is.
We have added the following text to this section: "Briefly, the coefficient of variation describes the uncertainty in the parameter (calculated as the parameter variance divided by the mean), the elasticity describes the sensitivity of the model to a relative change in the parameter, and the partial variance synthesizes these two metrics to describe the relative contribution of uncertainty in a parameter to the total predictive uncertainty in the model output (i.e., the parameters that have the highest partial variance are those that are highly uncertain and to which the model is highly sensitive; parameters that are highly uncertain but to which the model is relatively uncertain, and conversely, parameters to which a model is highly sensitive but whose values are known precisely, would both have low partial variance)." And we elaborate on the approach of LeBauer as follows: "We generally follow the approach of Lebauer et al. (2013) which samples from parameter distributions to generate an ensemble of model runs that approximate the posterior distribution of model output that can be used in the sensitivity analysis. Their sensitivity analysis is based on univariate perturbations of each parameter of interest, and the relationship between each parameter and model output is approximated by a natural cubic spline. Their model sensitivity is then based on the derivative of the spline at the parameter median." Line 188 -these 300-400 Pg C are not yet decomposed so comprise the Cthawed pool?I These values include carbon that has been decomposed as well as that which is still in the thawed pool. We have clarified this. The sentence now reads: "In RCP 4.5, 6.0, and 8.5, permafrost losses, including both thawed permafrost and permafrost carbon that has been decomposed and emitted to the atmosphere, reached 300-400 Pg C by 2100 and mostly leveled off after this point." Line 195 -is there any evidence that the ~3 Pg C /year has been found in other models. This is quite high.
We agree this is high compared to, for example, Burke et al. (2017). Since adjusting the default static fraction of thawed permafrost carbon in the model to 0.72, we find a substantially lower flux, closer to 1.5 Pg C/year. We have adjusted the model to use a higher static fraction of thawed permafrost, substantially reducing this temperature effect, and have also added a figure (see response above) to the sensitivity analysis section showing the sensitivity of temperature to the various permafrost controls in the model. We appreciate the tip about standard colors and have adjusted this figure accordingly.
Line 221 -Quite a lot of this carbon remains in the atmosphere (expect ~50 %, ~25% to land/~25% to ocean). Is this a function of the model structure?
We agree that this fraction is high compared to estimates of, e.g., the fraction of anthropogenic emissions that are estimated to remain in the atmosphere (Knorr 2009). Interestingly, however, more recent estimates of the airborne fraction sometimes estimate considerably higher values-see for example Yin et al. (2020). https://www.nature.com/articles/s41467-020-15852-2 Regardless, this potential discrepancy is interesting and will help Hector developers prioritize future efforts and model evaluation exercises. A current focus of these efforts is adding a carbon-tracking feature to the model, which will allow researchers to more clearly diagnose and evaluate the model's carbon flows in the future. We have simplified this figure to only include the coefficient of variation, elasticity, and partial variance in order to more closely follow the revised description of these statistics given in the methods section (see response above) and make it easier to follow. We have also included some of the description of these statistics from the revised methods section into the figure caption as follows: "The coefficient of variation describes the uncertainty in the parameter (parameter variance divided by the mean), the elasticity describes the sensitivity of the model to a relative change in the parameter, and the partial variance synthesizes these two metrics to describe the relative contribution of uncertainty in a parameter to the total predictive uncertainty in the model output." Line 235 -the 30-45% is 'partial variance'? gmd-2020-377 Yes. We have clarified this sentence to reflect that.
Line 275 -this should be encompassed by the parameter uncertainty.
This is an interesting point-but we don't believe it would be included by current parameter uncertainty, because the solution would be to add a completely new parameter controlling the decomposition rate of thawed permafrost (versus non-permafrost soil). In any event, we have removed this paragraph in the revised manuscript.
Line 275 -280 -this relates to the abrupt thaw processes discussed above and the two discussed together.
We have moved this paragraph up in section 4.1 and have connected it better to abrupt thaw.

Line 304 -Chadburn et al assumed the soil and air temperatures were in equilibrium in their analysis.
We appreciate this point of clarification and have added the following language to that sentence in the discussion: "though this analysis used equilibrium temperatures and does not give us information about the potential for these insulation effects to play a role in mitigating transient thaw."

How well does the 200 year temperature term represent the thermal inertia of the permafrost?
It actually is not intended to represent the thermal inertia of permafrost, as this value is only used for estimating respiration from thawed soil. We regret that this was not clear in the paper and have added some clarifying language. See response above for additional comments on this point.

Section 4.1 probably need to mention nutrient limitation.
We appreciate this suggestion and have added a mention of nutrient limitation as follows: "Our results exclude any potential changes in plant productivity as a result of permafrost thaw, including any due to changes in nutrient availability, as the sign of these effects remains highly uncertain." Line 334 -please define/reference GCAM. Added.

CEC1 comments
-The code of the model must be stored in a permanent archive-for example, Zenodo. Github is not a repository acceptable for purposes of long-term storage.
We will add the code to a Zenodo archive to be cited in the revised manuscript.
-You must include the version of the model that you use in the title of the manuscript.
We will revise the title to include the model version number.