Quantifying the Terrestrial Carbon Feedback to Anthropogenic Carbon Emission

The surface warming response to carbon emission is dependent on feedbacks operating in both the physical climate and carbon cycle systems, with physical climate feedbacks quantified via linearly combinable climate feedback terms, λclimate in watt per square meter per kelvin. However, land carbon feedbacks are often quantified using a two‐parameter description, with separate cumulative carbon uptake responses to surface warming, γL in petagram of carbon per kelvin, and rising atmospheric CO2 concentration, βL in petagram of carbon per parts per million. Converting the γL and βL responses to an overall terrestrial carbon feedback parameter, λcarbon in watt per square meter per kelvin, has remained problematic, with λcarbon affected by significant nonlinear interactions between carbon‐climate and carbon‐concentration responses and a nonlinear relation between atmospheric CO2 and subsequent radiative forcing. This study presents new relationships quantifying how the overall steady state terrestrial carbon feedback to anthropogenic emission, λcarbon, is dependent on the terrestrial carbon responses to rising CO2 and temperature, βL, and γL, and the physical climate feedback, λclimate. Nonlinear interactions between βL and γL responses to carbon emission are quantified via a three‐parameter description of the land carbon sensitivities to rising CO2 and temperature. Numerical vegetation model output supports the new relationships, revealing an emerging sensitivity of land carbon feedback to climate feedback of ∂λcarbon/∂λclimate ~ 0.3. The results highlight that terrestrial carbon feedback and physical climate feedback cannot be considered in isolation: Additional surface warming from stronger climate feedback is automatically compounded by reduced cooling from terrestrial carbon feedback, meanwhile around half the uncertainty in terrestrial carbon feedback originates from uncertainty in the physical climate feedback.


Introduction
The surface warming response to anthropogenic carbon emission is dependent on feedbacks operating both in the physical climate system (e.g., Knutti et al., 2017) and in the biogeochemical cycling of carbon (e.g., Friedlingstein et al., 2006). Feedbacks operating in the physical climate system include the Planck, water vapor-lapse rate, cloud, and snow and sea ice albedo feedbacks (e.g., IPCC, 2013). These individual feedbacks are quantified in terms of their climate feedback responses, λ in watt per square meter per kelvin, which are added together to find the total physical climate feedback, λ climate in watt per square meter per kelvin (IPCC, 2013;Knutti et al., 2017) The land carbon system responds to rising global temperatures and CO 2 levels via a number of feedback mechanisms (e.g., Friedlingstein et al., 2006;IPCC, 2013). The Net Primary Productivity (NPP) of land ecosystems, removing CO 2 from the atmosphere into the land system, is thought to increase with rising atmospheric CO 2 levels through CO 2 fertilization (e.g., Alexandrov et al., 2003). NPP is also sensitive to global mean temperatures, including via NPP sensitivity to other factors that are themselves linked to changes in global temperature such as the hydrological cycle. The rate of microbial respiration of soil carbon, returning land carbon to the atmosphere, is thought to increase with global mean temperature due to increased metabolic rate (e.g., Friedlingstein et al., 2006). The strengths of all these sensitivities are highly uncertain globally (e.g., Arora et al., 2013;Friedlingstein et al., 2006;Gregory et al., 2009;IPCC, 2013), in part due to uncertainty in how other factors affect the land carbon system such as nutrient availability. Other land carbon feedback processes occur over long timescales, such as how permafrost thawing releases locked carbon to the atmosphere (e.g., MacDougall & Knutti, 2016;Schuur et al., 2015).
Carbon cycle feedbacks, excluding long timescale responses such as permafrost, are often quantified in terms of two sensitivity terms representing the land and ocean carbon cycle responses to rising atmospheric CO 2 and temperature (e.g., Arora et al., 2013;Friedlingstein et al., 2006;Gregory et al., 2009). For the land carbon cycle, the carbon-climate feedback expresses the sensitivity of cumulative land carbon uptake to rising global mean surface temperature, γ L in petagram of carbon per kelvin, while the carbon-concentration feedback expresses the sensitivity of cumulative land carbon uptake to rising atmospheric CO 2 , β L in petagram of carbon per parts per million (e.g., Arora et al., 2013;Friedlingstein et al., 2006;Gregory et al., 2009).
There are difficulties in expressing the overall terrestrial carbon feedback to rising CO 2 and temperature as a λ carbon term in watt per square meter per kelvin, such that terrestrial carbon feedbacks can be easily compared to and combined with physical climate feedbacks (e.g., Arora et al., 2013 ;Gregory et al., 2009). First, the carbon-climate and carbon-concentration feedbacks are interdependent, such that nonlinear interactions altering the effective values of γ L and β L significantly affect the terrestrial carbon response to a scenario with both rising CO 2 and temperature (Arora et al., 2013;Gregory et al., 2009). Second, γ L and β L are known to be time evolving and path dependent, such that their values at any given time depend on history of the temperature and CO 2 (Arora et al., 2013). Third, the land carbon-concentration and carbonclimate feedback terms, γ L and β L , calculate the cumulative carbon uptake by the land system in petagram of carbon, and not the radiative forcing from the change in atmospheric CO 2 due to land carbon uptake in watt per square metre. Gregory et al. (2009) convert γ L and β L into a land carbon feedback in watt per square metre per kelvin by 1. Assuming that the carbon uptake by the land system causes an equal and opposite carbon loss by the atmosphere, and 2. Assuming this carbon loss by the atmosphere (changing atmospheric CO 2 in parts per million) relates linearly to the carbon feedback impact on radiative forcing in watt per square meter.
However, neither of these assumptions holds when isolating the land carbon feedback. First, although the carbon taken up by the land system does initially originate from the atmosphere, as atmospheric CO 2 is being removed by the land system this inevitably leads to outgassing of CO 2 from the ocean to the atmosphere across the air-sea interface through chemical exchange (e.g., Zeebe & Wolf-Gladrow, 2001). Therefore, the net reduction in atmospheric carbon due to land carbon uptake is less than the amount of carbon taken up by the land because, due to the induced air-sea exchange, the carbon removed into the land originates from both the atmosphere and ocean (e.g., Goodwin et al., 2008). Note that this effect does not impact the total land and ocean carbon feedback as analyzed by Gregory et al. (2009), but only the extraction of the land component. Second, radiative forcing in watt per square meter is related approximately logarithmically to the change in atmospheric CO 2 (Myhre et al., 2013), such that the same reduction in atmospheric carbon (Pg C or ppm) has approximately half the radiative forcing impact (W m −2 ) if background atmospheric CO 2 is doubled. Goodwin et al. (2019) recently identified how the magnitudes of terrestrial carbon uptake and surface warming since the preindustrial can be used to calculate the overall land carbon feedback, λ carbon in watt per square meter per kelvin, finding λ carbon = 0.31 ± 0.09 W m −2 K −1 for the present day based on observational reconstructions. However, it is not known how this single land carbon feedback term, λ carbon , relates to the more established land carbon-climate γ L , and land carbon-concentration β L , parameters that are typically evaluated in coupled model simulations (e.g., Arora et al., 2013;Friedlingstein et al., 2006;Gregory et al., 2009).
This study identifies how the steady state terrestrial carbon feedback (λ carbon in W m −2 K −1 ) following anthropogenic carbon emission is related to the land carbon-climate, land carbon-concentration, and physical climate feedbacks (γ L , β L , and λ climate respectively). The relationships account for both the subsequent air-sea exchange of CO 2 due to land carbon uptake and the logarithmic relationship between atmospheric CO 2 and radiative forcing. The link between λ carbon and λ climate is first explored for small perturbations, and then the impact of nonlinear interactions between carbon-climate and carbon-concentration responses are quantified, identifying a relationship for λ carbon for large carbon emission perturbations. Using these new relationships for λ carbon , this study then shows how the amplification of anthropogenic warming due to terrestrial carbon feedback is dependent on both the equilibrium climate sensitivity (ECS, in K) and carbon emission size, with more likelihood of warming amplification by the terrestrial carbon system when ECS and emission size are large. Note that this steady state analysis does not consider slow land carbon responses with centennial and millennial timescales, such as the permafrost carbon feedback (e.g., MacDougall & Knutti, 2016;Schuur et al., 2015).

Warming Response in the Absence of Terrestrial Carbon Feedback
A pulse of CO 2 initially emitted into the atmosphere will eventually partition between the atmosphere, ocean, and land systems. The total carbon emitted, δI em (t) in petagram of carbon, will at any time t be equal to the sum of the increases in carbon inventories within the atmosphere, δI atm (t), ocean, δI ocean (t), and land, δI land (t), carbon systems, The question is, how will the partitioning of the carbon emission between the atmosphere, ocean, and land systems evolve over time? First, consider an atmosphere-ocean only system at an initial steady state, with no carbon exchanges allowed with the land system such that δI land = 0 in (1) over all time t. The atmosphereocean system is then perturbed by an instantaneous pulse of carbon emission at time t 0 , δI em . At the initial moment of the emission pulse all of the emitted carbon enters the atmosphere, and the increase in atmospheric carbon is therefore equal to the emission pulse size, δI em (t 0 ) = δI atm (t 0 ). Subsequently, this rise in atmospheric CO 2 inevitably leads to a flux of carbon into the ocean due to chemical exchange across the air-sea interface and the emitted carbon is now partitioned between the atmosphere and ocean, I em (t) = δI atm (t)+δI ocean (t).
Over many centuries, the system reaches a new steady state with the carbon emission partitioned between the atmosphere and ocean (Goodwin et al., 2007). Once CO 2 crosses the air-sea interface, it combines with water and dissociates forming three chemical species (e.g., Zeebe & Wolf-Gladrow, 2001) comprising dissolved inorganic carbon (DIC): an uncharged form consisting of aqueous CO 2 and carbonic acid (CO 2 * ), a single-charged bicarbonate ion form (HCO 3 − ), and a double-charged carbonate ion form (CO 3 2− ). Air-sea exchange of CO 2 is determined by the CO 2 * component of DIC, and at the preindustrial chemical state, the approximate ratios of CO 2 * :HCO 3 − :CO 3 2− are around 1:100:10. However, as more CO 2 dissolves in the ocean seawater becomes more acidic and the relative fraction of DIC composed of CO 2 * increases, while the relative fraction composed of CO 3 2− decreases.
Due to this nonlinear response ocean carbonate chemistry, the fraction of the emitted carbon that remains in the atmosphere depends on the emission size: as more carbon is emitted the ocean becomes more acidic and less soluble to further CO 2 . With no land carbon response, the change in atmospheric CO 2 over multicentury timescales, t = t cent , once the emitted carbon becomes chemically partitioned between the atmosphere and ocean is related to the cumulative carbon added to the air-sea system through carbon emission, δI em in petagram of carbon, via (Goodwin et al., 2007(Goodwin et al., , 2008(Goodwin et al., , 2009, using the notation δln x = ln(x+δx) − lnx, and where I B is the preindustrial atmosphere-ocean buffered

10.1029/2019EF001258
Earth's Future carbon inventory of around 3,451 ± 96 Pg C in the Climate Model Intercomparison Project phase 5 (CMIP5) models evaluated in Williams et al. (2017). The buffered carbon inventory I B represents the amount of CO 2 and DIC that is capable of redistributing between the atmosphere and ocean in the atmospheric CO 2 , ocean CO 2 * , and ocean CO 3 2− pools and excludes the ocean DIC stored in the HCO 3 − pool: , where V is the volume of the ocean (Goodwin et al., 2009). Equation (2) holds for carbon emissions up to δI em~5 ,000 Pg C, because the value of I B can be assumed constant since the increases in I atmos and V [CO 2 * ] as more carbon is emitted into the system are opposed by a decrease in V[CO 3 2− ] (Goodwin et al., 2007(Goodwin et al., , 2009. The impacts of the CaCO 3 system on atmospheric CO 2 acting over multimillennial timescales (e.g., Archer, 2005;Goodwin & Ridgwell, 2010) are ignored in this study, which focusses on a century timescale response.
The logarithmic term in equation (2) expresses the impact of nonlinear ocean carbonate chemistry on the air-sea partitioning of carbon emitted into the air-sea system: The fraction of emitted carbon remaining in the atmosphere increases with the cumulative carbon emission size, while the fraction of emitted carbon taken up by the ocean decreases with emission size, due to the decreasing solubility of CO 2 in seawater as the ocean becomes more acidic (see also, e.g., Zeebe & Wolf-Gladrow, 2001). Thus, the nature of ocean carbonate chemistry implies that a constant sensitivity of ocean carbon uptake, δI ocean , to atmospheric CO 2 , via either δCO 2 or δI atm , cannot be defined: the sensitivity is itself dependent on the atmospheric CO 2 level. Therefore, this study refrains from attempting to define an ocean carbon-concentration feedback strength (e.g., see Friedlingstein et al., 2006), β ocean , in units of Pg C ocean uptake per unit ppm increase in atmospheric CO 2 . Instead, ocean carbonate chemistry is utilized to express the sensitivity of ocean carbon uptake to atmospheric CO 2 via (2). Note that the nonlinear ocean carbonate chemistry does not affect the ability to define a land carbon-concentration feedback, β L .
The rise in CO 2 from carbon emission in (2) induces a radiative forcing (Myhre et al., 2013), which in turn induces a surface warming (Williams et al., 2012). The radiative forcing from the increase in CO 2 due to carbon emission, δR em CO2 , is related to the increase in the log of atmospheric CO 2 , δR em CO2 ¼ aδln CO 2 , where a=5.35 ± 0.27 W m −2 is the CO 2 radiative forcing coefficient (Myhre et al., 2013), making δR em CO2 linearly related to carbon emission (2) (Goodwin et al., 2009). This radiative forcing from carbon emission induces surface warming until the radiative forcing is balanced by an increase in outgoing radiation from elevated surface temperatures, λ climate δT 0 in watt per square meter, via, where δT 0 is the steady state temperature change from carbon emission in the absence of terrestrial carbon response in K and λ climate is the climate feedback in watt per square meter per kelvin. The climate feedback is formally defined as the sensitivity of Earth's net radiation balance to changes in surface temperature, λ climate = − δR feedback /δT, where δR feedback is the change in net downward energy flux at the top of the atmosphere due to the change in global mean surface temperature, δT. Note that the sign convention adopted here is such that λ climate is positive, because a surface warming (δT > 0) causes a net upward radiation flux (R feedback < 0). The next section considers how this relationship between warming and emissions (3) is altered by the presence of a terrestrial carbon system for small perturbations.

Impact of Terrestrial Carbon Feedback for Small Perturbation
Section 3.1 finds a relationship to calculate steady state λ carbon following a small carbon emission, in terms of λ climate , β L , and γ L . Section 3.2 tests this relationship using numerical model simulations.

Theory
Now consider an atmosphere-ocean-land system at an initial steady state then perturbed by a CO 2 emission, with no perturbations to other sources of radiative forcing. Once the system reaches a new steady state, the carbon emission will be partitioned between the atmosphere, ocean, and land systems. The component of emitted carbon that remains in the atmosphere will increase atmospheric CO 2 and induce a radiative forcing that causes a rise in surface temperatures.

Earth's Future
Terrestrial carbon storage, I ter in petagram of carbon, is sensitive to changes in both atmospheric CO 2 levels and climate, with global mean surface temperature commonly used to represent the level of climate change (e.g., Friedlingstein et al., 2006). At steady state, a small perturbation in terrestrial carbon storage, δI ter in petagram of carbon, is related to small perturbations in atmospheric CO 2 , δCO 2 in parts per million, and global mean surface temperature change, δT in K, via, The empirically determined carbon-concentration feedback, β L ¼ ∂Iter ∂CO2 T in petagram of carbon per parts per million and carbon-climate feedback, γ L ¼ ∂Iter ∂T CO2 in petagram of carbon per kelvin, represent the cumulative terrestrial carbon uptake sensitivities to atmospheric CO 2 (at constant preindustrial temperature) and global mean surface warming (at constant preindustrial CO 2 ) respectively, following the framework set out in Friedlingstein et al. (2006). Note that β L may change with background temperature and γ L may change with background CO 2 , leading to significant nonlinearities between the carbon-climate and carbonconcentration feedbacks (Arora et al., 2013;Gregory et al., 2009). Therefore, (4) is only strictly applicable to small perturbations.
Relative to the case in the absence of terrestrial carbon feedback (equations (2) and (3)), this change in terrestrial carbon storage (4) will alter the steady state rise in atmospheric CO 2 (1) and so also alter the radiative forcing from atmospheric CO 2 and the global mean surface warming (3). For a hypothetical atmosphereland only system, with no coupled ocean, an increase in land carbon storage would cause and equal and opposite decrease in atmospheric carbon storage. However, for a coupled atmosphere-ocean-land carbon system, an increase in land carbon storage leads to, and is balanced by the sum of, decreases in both the atmosphere and ocean carbon storage. Initially, the additional carbon stored in the land system comes from the atmosphere, but over time this decrease in atmospheric CO 2 then causes an inevitable ocean outgassing due to air-sea chemical exchange. By similarity to equation (2), we find that the change in atmospheric CO 2 over multicentury timescales due to an initial uptake of carbon by the terrestrial system, after accounting for subsequent air-sea gas exchange, is given by (Goodwin et al., 2008), where the carbon added to the air-sea system in (2) due to emission, δI em , is replaced here by the carbon added to the air-sea system by terrestrial carbon uptake, −δI ter , noting the minus sign arises because an increase in terrestrial carbon storage removes carbon from the air-sea system.
When terrestrial carbon uptake is considered in the context of a coupled atmosphere-ocean-land carbon system peterbed by anthropogenic carbon emission, the atmosphere-ocean relationship for the total log CO 2 change at steady state, equation (2), is modified to contain an additional term representing how the total carbon added to the air-sea system now has contributions from both carbon emission, δI em , and terrestrial carbon uptake, −δI ter , (Goodwin et al., 2007(Goodwin et al., , 2008(Goodwin et al., , 2009(Goodwin et al., , 2015, giving This relationship calculates the long-term atmospheric CO 2 change in response to carbon emission and terrestrial carbon uptake, accounting for the inevitable air-sea gas exchange over many centuries through the I B terms.

Earth's Future
The rise in surface temperature from carbon emission in the presence of terrestrial carbon uptake, δT in K, is then given by this total radiative forcing accounting for both the emissions and terrestrial carbon response, noting the identity δlnx = δx/x, The climate feedback may be considered in terms of the change in Earth's radiation balance from physical climate system induced changes per unit increase in surface temperature: λ climate = − δR feedback /δT. By similarity, we may define the terrestrial carbon feedback in terms of the change in Earth's energy balance from terrestrial carbon system induced changes in atmospheric CO 2 per unit surface warming (Goodwin et al., 2019): λ carbon ¼ −δR ter CO2 =δT in watt per square meter per kelvin. Substituting δCO 2 = (λ climate CO 2 /a)δT from (8) into (4) reveals δI ter = β L (λ climate CO 2 /a)δT+γ L δT and then dividing both sides by δT gives δI ter /δT = β L (λ climate CO 2 /a)+γ L . Finally, multiplying both sides by a/I B , to express δI ter in terms of δR ter CO2 using (7), δR ter CO2 ¼ − a=I B ð ÞδI ter , reveals how steady state terrestrial carbon feedback, λ carbon , is related to λ climate , β L , and γ L , This relationship (9) solves for the terrestrial carbon feedback following carbon emission once atmosphereocean-land carbon partitioning has reached a steady state, and temperatures have stabilized with respect to the elevated atmospheric CO 2 . Equation (9) predicts that for given land carbon-concentration and carbonclimate responses to an infinitesimal carbon emission, the carbon feedback, λ Carbon ¼ −δR ter CO2 =δT, is linearly related to climate feedback, λ climate . Using the mean and standard deviation values of β L (0.92 ± 0.44 Pg C ppm −1 ) and γ L (−58.4 ± 28.5 Pg C K −1 ), from the CMIP5 models evaluated by Arora et al. (2013) following a 4 × CO 2 experiment, equation (9) predicts: λ carbon = (0.30 ± 0.14)λ climate +(−0.09 ± 0.04), assuming normal error propagation and adopting I B = 3451 ± 96 Pg C, a = 5.35 ± 0.27 W m −2 , and CO 2 = 1120 ppm.
Diagnosing λ carbon from land carbon uptake, δI ter , and surface warming, δT, using equation (9) rests on two assumptions: the use of the buffered carbon inventory I B to calculate δlnCO 2 and the use of the radiative forcing coefficient, a, to calculate the radiative forcing from δlnCO 2 . The discrepancy in δln CO 2 as predicted using I B via equations (2) or (6) remains under 3% for carbon perturbation up to the approximate magnitude of the entire land carbon reservoir, δI ter~2 ,000 Pg C, when compared to multicentury numerical simulations with explicit representations of ocean carbonate chemistry (Goodwin et al., 2007). Once partitioned between the atmosphere and ocean, a δI ter~2 ,000 Pg C magnitude perturbation would change atmospheric CO 2 by around δlnCO 2~2 ,000/I B~0 .6. The discrepancy in δR CO2 when using δR CO2 = aδlnCO 2 , with a = 5.35 W m −2 , is around 5% when compared to calculations containing second-order terms (δR CO2 = 5.32δlnCO 2 +0.26[δlnCO 2 ] 2 : Byrne & Goldblatt, 2014). Thus, the two assumptions in equation (9) are valid for plausible magnitude land carbon perturbations.
Utilizing (9) and (7) in (8) then relates steady state surface warming to cumulative carbon emission in the presence of terrestrial carbon responses to rising CO 2 and temperature, Inspecting equation (10) shows that λ carbon is directly comparable to, and linearly combinable with, physical climate feedbacks evaluated in watt per square meter per kelvin such as the water vapor-lapse rate and cloud feedbacks that make up λ climate (IPCC, 2013;Knutti et al., 2017).

Comparison of Theory to Numerical Model Output
This section tests the prediction from (9) that λ carbon is linearly related to λ climate under a fixed CO 2 perturbation using numerical Dynamic Global Vegetation Model (DGVM) output from Pugh et al. (2018) and output from an efficient Earth system model (Goodwin, 2016(Goodwin, , 2018 (Figure 1a, black dots), where λ climate is diagnosed here from the model temperature response to CO 2 using a = 5.35 W m −2 : λ climate = aΔlnCO 2 /ΔT.
Separately, an ensemble of 6,270 observation-constrained simulations of the efficient Warming Acidification and Sea level Projector (WASP) Earth system model (Appendix A: Goodwin, 2016Goodwin, , 2018Goodwin et al., 2019) is integrated to steady state following a 1-year increase in CO 2 from 280 to 850 ppm. The ensemble of 6,270 simulations are generated from the Monte Carlo combined with history matching approach set out in Goodwin et al. (2018), using the WASP model configuration of Goodwin (2018). In this configuration the value of λ climate is allowed to vary over multiple response timescales linked to the different timescales of climate feedback processes (Goodwin, 2018). For example, there is an instantaneous contribution to λ climate from the Planck feedback, while contributions from the fast cloud response and water vapor-lapse rate response occur over order 10 days linked to the residence time of water vapor in the atmosphere, and contributions to λ climate from the way that changes in sea surface warming pattern alter the cloud response occur over decades.
The Monte Carlo combined with history matching method generates 6,270 observation-consistent simulations in the following way. First, an initial ensemble of 10 million simulations is generated with varying model parameter values, using the parameter input distributions of Goodwin (2018). The model parameters for climate feedback from different processes are varied to span the ranges evaluated in CMIP5 models (Goodwin, 2018). Also, the initial ensemble sensitivities of terrestrial NPP and soil carbon residence time to global temperature and CO 2 are varied to span the range of sensitivities seen the in the C4MIP models analyzed by Freidlingstein et al. (2006-see Figure 3 therein). These 10 million initial simulations are integrated from preindustrial to present day and evaluated for observational consistency against observational reconstructions of surface warming (Hansen et al., 2012;IPCC, 2013;Morice et al., 2012;Smith et al., 2008;Vose et al., 2012), ocean heat content (Balmaseda et al., 2013;Cheng et al., 2017;Giese & Ray, 2011;Good et al., 2013;Levitus et al., 2012;Smith et al., 2015), and ocean and terrestrial carbon uptake (IPCC, 2013;le Quéré et al., 2018) after Goodwin et al. (2018). The observation consistency test of Goodwin (2018 -see Table 2  A total of 6,273 simulations pass the updated observation-consistency test. Three of these simulations are excluded as nonphysical, since their values of λ climate become negative on long timescales, leaving a final ensemble of 6,270 observation-consistent simulations (Goodwin et al., 2019). This final ensemble therefore contains ranges of terrestrial carbon sensitivities to temperature and CO 2 that agree with both the analyzed sensitivities of the C4MIP models (Friedlingstein et al., 2006) and observational reconstructions of cumulative carbon uptake (le Quéré et al., 2018).
Each of the 6,270 observation-consistent ensemble members are reinitialized to preindustrial conditions, and forced with a 1-year step function increase in CO 2 from 280 to 850 ppm. Each ensemble member is then integrated for 500 years to reach a new steady state, without any imposed noise in the surface temperature. The values of ΔI ter and λ climate in the efficient model simulations are diagnosed at the end of the 500-year simulations to represent the new steady state reached. The observation-consistent ensemble of efficient model simulations shows a similar increasing trend in steady state ΔI ter at high λ climate (Figure 1a, blue

10.1029/2019EF001258
Earth's Future transparent dots) to the DGVM simulations (Figure 1a, black dots), but with greater variation reflecting the greater extent of parameter space explored.

Results From Model Output
Next, for the 22 DGVM and 6,270 efficient model simulations, λ carbon is calculated from ΔI ter and ΔT using equation (9): . The emergent linear link between steady state carbon feedback and climate feedback predicted from theory (equation (9)) is identified on both the DGVM ensemble (Figure 1b, black) and an efficient model ensemble (Figure 1, blue). In the DGVM simulations, with identical carbon cycle configurations, over 90% of the variance in λ carbon is explained by the variation in λ climate (Figure 2b, black: R 2 = 0.96). The efficient model ensemble contain variation in the carbon cycle model parameter values (Goodwin, 2018) and so will have variation in the effective values of β L and γ L between ensemble members. Despite this variation, around half of the variance in λ carbon is explained by the variation in λ climate (Figure 2b, blue: R 2 = 0.49). This demonstrates the robustness of the emergent link identified between terrestrial carbon and physical climate feedback, λ carbon and λ climate (equation (9)). The sensitivity of λ carbon to λ climate of~0.3 in the DGVM and efficient model ensembles (Figure 1b) is consistent with the sensitivity predicted using equation (9) for the CMIP5 model values of β L and γ L . Note that the analysis here is for steady state λ carbon , but that for transient cases λ carbon will vary over time as δI ter and δT vary (Goodwin et al., 2019).

Impact of Terrestrial Carbon Feedback for Large Perturbation
Nonlinear terms will affect the terrestrial carbon uptake response for large emission sizes, leading to errors when applying β L and γ L using equation (4) (e.g., Arora et al., 2013;Gregory et al., 2009). The question is, how will carbon feedback, λ carbon , alter for large emission perturbations due to these nonlinear terms compared with the expected value for small perturbations, equation (9)? Here instead of representing the sensitivity of the terrestrial carbon cycle to rising CO 2 and temperature via the carbon climate and carbon CO 2 feedback parameters, β L ¼ ∂Iter ∂CO2 T and γ L ¼ ∂Iter ∂T CO2 , respectively, the terrestrial carbon system is characterized in terms of empirical feedback parameters for aspects of the carbon system that allow nonlinear interactions between carbon cycle responses to temperature and CO 2 to be considered.

Earth's Future
First, consider a simple two box representation of the terrestrial carbon cycle coupled to an atmosphere (Figure 2), where the total terrestrial carbon storage, I ter , is the sum of the soil carbon reservoir, I soil in petagram of carbon, and the vegetation carbon reservoir, I veg in petagram of carbon: I ter = I veg +I soil . The vegetation carbon pool has an incoming carbon flux from the atmosphere due to NPP, F NPP in petagram of carbon per year. There is then a flux from the vegetation carbon pool into the soil carbon pool due to leaf litter, F leaflitter in petagram of carbon per year, and a flux from the soil carbon pool into the atmosphere due to soil carbon respiration, F respiration in petagram of carbon per year (Figure 2). At steady state the leaf litter and soil carbon respiration carbon fluxes must equal NPP, F NPP = F Leaflitter = F respiration , and note that a subscript 0 is used to denote the value of a quantity at the initial steady state (Figure 2).
Next, consider the carbon fluxes in the terrestrial carbon system to be sensitive to atmospheric CO 2 and temperature via the following parameters ( Figure 2): 1. A dimensionless CO 2 fertilization coefficient (Alexandrov et al., 2003) representing the fractional change in NPP flux per unit log change in CO 2 , β CO2 , such that at constant temperature F NPP = F NPP,0 (1 +β CO2 δlnCO 2 ); 2. A coefficient representing the fractional change in NPP per unit change in global mean surface temperature, c NPP in K −1 , such that at constant CO 2 F NPP = F NPP,0 (1+c NPP δT); and 3. A coefficient representing the fractional change in soil carbon residence time, τ soil in year, per unit change in global mean surface temperature, c soil in K −1 , such that τ soil = τ soil,0 (1+c soil δT) where F respiration = I soil /τ soil .
Adopting this representation for the CO 2 and T dependences of carbon fluxes within the land carbon system (Figure 2) allows the steady state terrestrial carbon storage be expressed in terms of the log CO 2 and warming perturbations, δlnCO 2 and δT, the initial NPP, F NPP,0 , and the initial residence timescales of carbon in the vegetation and soil carbon pools, τ veg,0 and τ soil,0 respectively, via (Appendix B) Note that this equation (11) solves for the steady state terrestrial carbon storage, I ter , for defined values of CO 2 and T, and so time dependencies are not shown. However, if the terrestrial carbon reservoir responds more quickly than slowly evolving changes in temperature or CO 2 , then (11) still applies and the terms in I ter , δlnCO 2 , and δT can be considered time dependent.
Substituting steady state relationships for δlnCO 2 = (λ climate CO 2 /a)δT and λ carbon ¼ −δR ter CO2 =δT ¼ aδI ter ð Þ= I B δT ð Þ, from (8) and (9), respectively, into equation (11) results in a second-order polynomial equation for λ carbon in δT (Appendix B), This relationship for steady state terrestrial carbon feedback, λ carbon , preserves nonlinear interactions between carbon-concentration and carbon-temperature responses, equation (12). It is noted that additional nonlinearities affecting λ carbon may exist that are not captured in (12), if the values of the coefficients β CO2 , c NPP , and c soil change with perturbation size.
By inspecting the leading order terms in (12), and comparing to (9), we can express the carbon-concentration and carbon-climate feedbacks in terms of the alternative carbon-system parameters, β CO2 , c NPP , and c soil : β L ≈ Iter;0 CO2 β CO2 and γ L ≈ I ter,0 c NPP +I soil,0 c soil .

Earth's Future
In equation (12) the term in δT is much larger than the term in δT 2 for reasonable temperature changes and parameter values. Therefore, the sensitivity of λ carbon to surface warming, at constant λ climate , from the nonlinear interactions between terrestrial carbon-climate and carbon-concentration responses simplifies to Noting that c soil and c NPP are likely negative (see Friedlingstein et al., 2006- Figure 3 therein), equation (13) therefore predicts that there will be a near-linear decrease in λ carbon as the perturbation in δT is increased for a given value of λ climate . Example carbon sensitivity values suggests a linear reduction in λ carbon with increasing temperature anomaly of order ∂λ carbon /∂T~− 0.015 W m −2 K −2 : using β CO2~0 .45, c NPP~− 0.04 K −1 , and c soil~− 0.02 K −1 (each within the range of the Earth system models analyzed in Friedlingstein et al., 2006-Figure 3 therein), along with λ climate = 1.2 W m −2 K −1 , a = 5.35 W m −2 (Myhre et al., 2013), I B = 3,451 Pg C (Williams et al., 2017), I ter,0 = 2,000 Pg C, and I soil,0 = 1,500 Pg C. Note that a positive λ carbon implies that terrestrial carbon feedback is negative, reducing surface warming, and so from (13) terrestrial carbon feedback is expected to become a less negative feedback (or even a positive feedback) with increasing temperature anomaly.

Comparison of Theory to Numerical Model Output
This section tests the prediction from (13), which λ carbon linearly reduces with perturbation size δT for a given value of λ climate , using published numerical DGVM output (Pugh et al., 2018) and output from the observation-consistent ensemble of 6,270 efficient Earth system model (Goodwin, 2016(Goodwin, , 2018 simulations (Figure 3). Pugh et al. (2018) integrate multiple DVGMs to steady state under CO 2 forced climate scenarios with warming of ΔT = 1, 2, 3, 4, and 5 K (Pugh et al., 2018-"DGVM ensemble" therein, using HYLAND, SDGVM, ORCHIDEE, TRIFFID, and LPJ models) and also integrate a DGVM within a coupled Earth system model at multiple CO 2 forcing scenarios achieving different levels of warming (the HadCM3LC simulations in Pugh et al., 2018). In these DVGM simulations (Pugh et al., 2018), the magnitude of steady state cumulative land carbon uptake, ΔI ter , initially increases with ΔT for CO 2 only forcing (Figure 3a, dots, diamonds, and dashed lines). However, the rate of increase in ΔI ter per unit additional surface warming reduces for all DVGMs, with some models showing a rate of decrease per unit additional warming as total warming nears 5 K (Figure 3a).

Descriptions of Model Output
The 6,270 efficient WASP model simulations are reinitialized to a preindustrial steady state and perturbed this time with carbon emissions scenarios that interactively restores atmospheric CO 2 to produce a range of specified surface warming targets, of ΔT = 1, 2, 3, 4, and 5 K (for description of the restoring method used in the WASP model see Nichols et al., 2018). Again, the simulations are integrated without imposed noise in the surface temperature (Goodwin, 2018). For each warming target, all 6,270 simulations are integrated for 500 years until a new steady state is reached. Global mean surface temperature anomaly is stabilized to within ±0.02 K of the desired target in at least 99% of the 6,270 simulations. This efficient model ensemble shows a similar pattern of change in steady state ΔI ter with increasing steady state ΔT to the DVGMs (Figure 3a, compare blue solid line and shading to black and color dots, grey diamonds, and associated dashed lines):

Results From Model Output
Terrestrial carbon feedback, λ carbon , is diagnosed from the model output ( Figure 3a The reduction in λ carbon for larger perturbations in δT in the DGVM and efficient model simulations (Figure 3) shows that the nonlinear interactions between carbon-climate and carbon-concentration feedbacks (equation (9)) are significant for large carbon emissions, in agreement with previous studies (e.g., Arora et al., 2013;Gregory et al., 2009), implying that the standard β L and γ L representation of the terrestrial carbon feedback may lead to significant error.
The agreement between the predicted near-linear decrease in λ carbon with increasing δT from equations (12) and (13) and the behavior of the ensemble of DGVM simulations (Figure 3) indicates that nonlinearities between terrestrial carbon-climate and carbon-concentration feedbacks may be captured by

10.1029/2019EF001258
Earth's Future considering a relatively simple representation of the terrestrial carbon cycle (Figure 2 and Appendix B). This simple representation includes three empirically determined sensitivities: a dimensionless CO 2 fertilization sensitivity of NPP, β CO2 ; a temperature sensitivity of NPP, c NPP in K −1 ; and a temperature sensitivity of the soil carbon respiration timescale, c soil in K −1 (equations (B4) and (B5)). Note that additional nonlinearities may become significant for the real terrestrial carbon cycle that are not captured in the DGVM simulations.

Implications for Gain in Surface Warming From Terrestrial Carbon Feedback
A gain factor for land carbon feedback on surface warming, G L , may be defined as the ratio of warming in presence of terrestrial carbon feedback divided by warming in the absence of terrestrial carbon feedback (Gregory et al., 2009). G L is expressed in terms of either the emission size and change in terrestrial carbon reservoir or the carbon and climate feedbacks, as where eqns. (9) and (12) show how λ carbon relates to λ climate for infinitesimal and finite carbon emission perturbations, respectively. For infinitesimal emission perturbations, where nonlinear interactions between CO 2 and T responses of terrestrial carbon cycle can be ignored, the gain G L becomes β L is likely positive, since NPP increases with rising CO 2 due to CO 2 fertilization. However, γ L is negative in many models (Arora et al., 2013), since soil carbon storage decreases with rising T as soil carbon residence Figure 4. Steady state gain in surface warming due to land carbon feedback varies with equilibrium climate sensitivity for a given CO 2 stabilization. The gain G L as a function of equilibrium climate sensitivity is calculated from equation (14) for different relationships between λ climate and λ carbon and for different β L and γ L values. For the efficient WASP model ensemble (blue solid line) and the DGVM TRIFFID model ensemble of Pugh et al. (2018) (black solid line), G L is calculated as a function of ECS using the relationships between λ climate and λ carbon identified in Figure 1b (equation (14)). G L is calculated as a function of ECS using values of β L and γ L identified by Arora et al. (2013) for the CMIP5 ensemble mean (orange solid line) and individual CMIP5 models (color dashed lines). Equation (14) is then applied assuming climate is stabilized with approximately present-day CO 2 levels of 410 ppm in (equation (14)).

10.1029/2019EF001258
Earth's Future time reduces (Friedlingstein et al., 2006). This means that the gain factor for terrestrial carbon feedback, G L , will increase at higher ECS or lower λ climate (14) (Figure 1b).
By inspecting how λ carbon changes for different values of λ climate for the DGVM and efficient model ensemble output in Figure 1b and converting λ climate into ECS, we can see from equation (14a) that the gain G L <1 for ECS ≤ 5.8 K (Figure 4, solid black and blue lines), such that terrestrial carbon feedback reduces surface warming from anthropogenic carbon emissions. However, the gain switches to G L >1 (equation (14a)) for ECS ≥ 5.9 K, such that the terrestrial carbon feedback increases surface warming from carbon emissions (Figure 4). For larger perturbation sizes (Figure 3 and equation (13)), the nonlinear interactions cause the carbon feedback to become less negative (or more positive), and so the ECS value above which terrestrial feedbacks switch from damping to amplifying anthropogenic warming would decrease. Other model ensembles may yield different results. The values of β L and γ L from the CMIP5 models analyzed by Arora et al. (2013), when applied to equation (14b) and considering perturbations that stabilize CO 2 at approximately present-day levels (CO 2 = 410 ppm), suggest considerable uncertainty in the value of ECS above which gain transitions from damping gain, G L < 1, to amplifying gain, G L > 1 (Figure 4, solid orange line and color dashed lines). The CMIP5 multimodel mean values of β L and γ L (Arora et al., 2013) suggest land carbon feedbacks amplify steady state warming for ECS above 4.5 K (G L >1) and dampen steady state warming below ECS of 4.5 K (G L <1), for CO 2 stabilization at 410 ppm ( Figure 4, orange solid line). However, the ECS values above which land carbon feedback amplifies steady state warming range from a low as 2.4 to above 10 K ( Figure 4, dashed color lines). Note that the values of β L and γ L for the CMIP5 models analysed by Arora et al. (2013) may be scenario or time dependent. Therefore, the variation in G L with ECS calculated in Figure 4 should be considered illustrative, for the CMIP5 values of β L and γ L , and not precise predictions of what would occur in the terrestrial components of the CMIP5 models if run to steady state with CO 2 levels of 410 ppm.

Discussion
Two of the most significant sources of uncertainty in the sensitivity of warming to anthropogenic carbon emission arise from uncertainties in the strength of feedbacks operating in the physical climate system (e.g., IPCC, 2013; Knutti et al., 2017) and the land carbon system (e.g., Arora et al., 2013;Friedlingstein et al., 2006). This study shows how, when the land carbon system reaches a new steady state following carbon emission, the strength of these physical climate and terrestrial carbon feedbacks is linked via theoretical relationships (equations (9) and (12)) and in numerical model simulations (Figures 1 and 3).
This identified link between terrestrial carbon and physical climate feedbacks implies that the impact on surface warming of the two systems should not be considered in isolation. First, when calculating surface warming from carbon emissions for a different climate feedback, one must also consider the impact on terrestrial carbon feedback. For example, an increase in the expected surface warming due to stronger than expected cloud feedback would be compounded by the subsequent reduction in the damping of surface warming from terrestrial carbon feedback (equations (9), (10), (12), and (14) and Figure 1). Second, a significant component of the uncertainty terrestrial carbon feedback arises from the uncertainty in physical climate feedback. In a large observation-constrained ensemble of many thousands of simulations containing significant variation in the carbon cycle responses to rising CO 2 and temperature (Goodwin, 2018;Goodwin et al., 2018), around half the uncertainty in steady state terrestrial carbon feedback (W m −2 K −1 ) originates from uncertainty in physical climate feedback (Figure 1b, blue: R 2 = 0.49).
In the present transient state, the terrestrial carbon feedback appears to be robustly negative with the terrestrial carbon cycle absorbing anthropogenic CO 2 from the atmosphere (le Quéré et al., 2018) and thus reducing anthropogenic warming from carbon emissions. However, the present transient state is also characterized by a lag between rising atmospheric CO 2 and rising surface temperatures, because the transient climate response is lower than the equilibrium climate sensitivity (IPCC, 2013;Knutti et al., 2017). As the terrestrial carbon cycle likely has opposing sensitivities to rising CO 2 and rising temperature (Friedlingstein et al., 2006;Gregory et al., 2009), the future response of the terrestrial carbon uptake depends crucially on the climate sensitivity determining the relative increases in CO 2 and temperature following carbon emission (equations (9) and (12)).

Earth's Future
The analysis presented here implies that we may not simply assume that the terrestrial carbon feedback will remain robustly negative at steady state, in line with previous studies finding that the current land carbon sink may become a carbon source over the 21st century or beyond (e.g., Cox et al., 2000;Friedlingstein et al., 2006). There is an increased likelihood that terrestrial carbon feedback will turn positive, enhancing future anthropogenic warming (equation (12)), either at high climate sensitivity (Figure 1 and equations (9) and (14)) or for large warming perturbations caused by increased anthropogenic emissions (Figure 3 and equations (12) and (13)).

Appendix A: Model Simulations and Analysis of Model Output
This appendix provides details of how the efficient model simulations are performed and how the both DVGM and efficient model output is analyzed.
To generate the WASP simulations, the observational constraints of Goodwin (2018) are applied, updated here with a consistency test for land carbon uptake based on the Global Carbon Budget (le Quéré et al., 2018). All consistency tests and ranges remain as in Goodwin (2018), and based on observational constraints on surface warming, ocean heat content, and carbon fluxes (see Goodwin, 2018- The cumulative residual land carbon uptake (excluding carbon emitted from land use change) from preindustrial to 2017 is considered observation consistent if the simulation lies between 96 and 331 Pg C, replacing the equivalent terrestrial carbon uptake range in Goodwin (2018- Table 2 therein). All other ranges for observational constraints are as in Goodwin (2018).
When converting simulated terrestrial carbon uptake, ΔI ter , into terrestrial carbon feedback, λ carbon , using equation (9) the values of a and I B must be known (Figures 1 and 3). To analyze the DGVM simulations of Pugh et al. (2018), values of a = 5.35 W m −2 and I B = 3,451 Pg C are used and λ climate is diagnosed using λ climate = aΔlnCO 2 /ΔT. For the efficient WASP model simulations the values of a and I B are considered individually for each ensemble member. The values of λ climate evolve over time in the efficient model (Goodwin, 2018), so the values at the end of the 500 year simulations are used (Figures 1b and 3).

Appendix B: Steady State Carbon Uptake in the Idealized Terrestrial Carbon System B.
This appendix provides the derivation of cumulative terrestrial carbon uptake, ΔI ter , following carbon emission including nonlinear interaction between the land carbon responses to rising CO 2 and surface warming, ΔT.
Consider an idealized system consisting of an atmosphere containing CO 2 coupled to a vegetation carbon reservoir and a soil carbon reservoir (Figure 2). The flux of carbon from the atmosphere to the vegetation pool is due to NP, F NPP in petagram of carbon per year. The flux of carbon from the vegetation to the soil carbon pool due to leaf litter, F LeafLitter in petagram of carbon per year, is equal to the vegetation carbon inventory, I veg in petagram of carbon, divided by the vegetation carbon residence timescale, τ veg in years, The flux of carbon from the soil to the atmosphere due to respiration, F Respiration , in petagram of carbon per year, is similarly equal to the soil carbon inventory, I soil in petagram of carbon, divided by the soil carbon residence timescale, τ soil in year 10.1029/2019EF001258 Earth's Future These three fluxes must be equal both at the initial steady state, F NPP,0 = F Leaflitter,0 = F respiration,0 (where subscript 0 is used to denote the initial conditions) and once the system reaches a new steady state after perturbation.
At steady state the terrestrial carbon storage, I ter in petagram of carbon, is given by the sum of the vegetation and soil carbon reservoirs, also written in terms of the residence timescales using (B1) and (B2) such that the initial steady state is written I ter,0 = I veg,0 +I soil,0 = F NPP,0 (τ veg,0 +τ soil,0 ).
Next, we assume that F NPP varies from the initial steady state with both the log change in atmospheric CO 2 , due to CO 2 fertilization (e.g., Alexandrov et al., 2003), and the global mean surface temperature (e.g., Friedlingstein et al., 2006) via where β CO2 is the empirically determined dimensionless CO 2 fertilization coefficient relating the sensitivity of NPP to the log change in CO 2 (Alexandrov et al., 2003) and c NPP is the empirically determined fractional sensitivity of NPP to global mean surface temperature in per kelvin.
Soil carbon residence time is also known to be sensitive to global mean surface temperature due to temperature effects on microbial respiration (e.g., Friedlingstein et al., 2006). Here we assume an idealized relationship where c soil is the empirically determined fractional sensitivity of soil carbon residence time to global mean surface temperature in per kelvin.