the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Timescale dependence of airborne fraction and underlying climate–carboncycle feedbacks for weak perturbations in CMIP5 models
Guilherme L. Torres Mendonça
Julia Pongratz
Christian H. Reick
The response of the global climate–carboncycle system to anthropogenic perturbations happens differently at different timescales. The unravelling of the memory structure underlying this timescale dependence is a major challenge in climate research. Recently the widely applied α–β–γ framework proposed by Friedlingstein et al. (2003) to quantify climate–carboncycle feedbacks has been generalized to account also for such internal memory. By means of this generalized framework, we investigate the timescale dependence of the airborne fraction for a set of Earth system models that participated in CMIP5 (Coupled Model Intercomparison Project Phase 5). The analysis is based on published simulation data from C^{4}MIPtype (Coupled Climate–Carbon Cycle Model Intercomparison) experiments with these models. Independently of the considered scenario, the proposed generalization describes at global scale the reaction of the climate–carbon system to sufficiently weak perturbations. One prediction from this theory is how the timescaleresolved airborne fraction depends on the underlying feedbacks between climate and the carbon cycle. These feedbacks are expressed as timescaleresolved functions depending solely on analogues of the α, β, and γ sensitivities, introduced in the generalized framework as linear response functions. In this way a feedbackdependent quantity (airborne fraction) is predicted from feedbackindependent quantities (the sensitivities). This is the key relation underlying our study. As a preparatory step, we demonstrate the predictive power of the generalized framework exemplarily for simulations with the Max Planck Institute (MPI) Earth System Model. The whole approach turns out to be valid for perturbations of up to an about 100 ppm CO_{2} rise above the preindustrial level; beyond this value the response becomes nonlinear. By means of the generalized framework we then derive the timescale dependence of the airborne fraction from the underlying climate–carboncycle feedbacks for an ensemble of CMIP5 models. Our analysis reveals that for all studied CMIP5 models (1) the total climate–carboncycle feedback is negative at all investigated timescales, (2) the airborne fraction generally decreases for increasing timescales, and (3) the land biogeochemical feedback dominates the model spread in the airborne fraction at all these timescales. Qualitatively similar results were previously found by employing the original α–β–γ framework to particular perturbation scenarios, but our study demonstrates that, although obtained from particular scenario simulations, they are characteristics of the coupled climate–carboncycle system as such, valid at all considered timescales. These more general conclusions are obtained by accounting for the internal memory of the system as encoded in the generalized sensitivities, which in contrast to the original α, β, and γ are scenarioindependent.
 Article
(7250 KB)  Fulltext XML
 BibTeX
 EndNote
The global carbon cycle plays a key role in determining the sensitivity of Earth's climate to anthropogenic emissions from fossilfuel burning, cement production, and landuse change. The increase in atmospheric CO_{2} concentrations driven by those emissions is considered the main radiative forcing driving climate change (see e.g. Gulev et al., 2021, Fig. 2.10). But to infer the resulting changes in atmospheric CO_{2} it is not sufficient to consider anthropogenic emissions alone; also the response of the global carbon cycle to these emissions must be taken into account. In reaction to rising emissions, land and ocean store increasing amounts of carbon, on average 56 % of emissions, a number that stayed surprisingly constant over the last 6 decades (Canadell et al., 2021, Fig. 5.7). This storage fraction strongly depends on environmental conditions: in years with a positive phase of the El Niño–Southern Oscillation (ENSO) it can be as low as 20 %, while in negative phases it may rise up to 75 % (Canadell et al., 2021, Fig. 5.7). One must thus suspect that with rising atmospheric CO_{2} the resulting climate change will affect this fraction of emissions stored away with potentially large consequences for the amount of anthropogenic CO_{2} remaining in the atmosphere. In this way, the carbon cycle may either slow or accelerate climatic changes resulting from anthropogenic emissions. Understanding the dynamics of this cycle, and in particular how it responds to perturbations from emissions and interferes with climate, thus constitutes an essential step in advancing climate research (Marotzke et al., 2017).
To gain insight into this combined dynamics of carbon cycle and climate, one must in particular study climate–carboncycle feedbacks. Such feedbacks arise from the already mentioned reaction of the global carbon cycle to a change in atmospheric CO_{2} that may amplify or counteract the initial change. For example, a rise in CO_{2} concentrations caused by fossilfuel emissions drives CO_{2} into the ocean by increasing the difference between the CO_{2} partial pressure in the atmosphere and that in surface waters (Takahashi et al., 2009). This flux of CO_{2} into the ocean reduces the initial increase in atmospheric CO_{2}. On the other hand, increasing atmospheric CO_{2} leads by the greenhouse effect to a rise in air temperatures. The warmer climate, in turn, leads to an acceleration of soil microbial activity and thereby to an increase in soil respiration rates (Raich and Potter, 1995). The resulting enhanced land CO_{2} flux into the atmosphere leads to even more CO_{2} in the atmosphere so that the initial increase in atmospheric CO_{2} is amplified. In principle, the global response of the whole carbon cycle to emissions may be quantified by summing the contributions from all such feedback mechanisms.
Depending on the speed at which the various feedbacks unfold, climate change may develop differently. Generally, the dynamics of the coupled climate–carboncycle system arising in response to perturbations depends on the spectrum of internal timescales of the various processes involved in the response. For instance, the speed at which global climate is warming in reaction to anthropogenic emissions depends largely on the rate at which the oceans can take up heat, and this rate – actually an inverse timescale – is determined by the temporal characteristics of the internal ocean dynamics, like the rate of mixing between upper and deep ocean and the speed at which heat is redistributed by ocean currents. Similar remarks apply to the uptake and redistribution of CO_{2} by the oceans. An indication of the involved timescales is obtained by noting that most of the carbon taken up by the ocean since preindustrial times still resides in its upper layers (74 % in the top 700 m; Frölicher et al., 2015). For land, the temporal characteristics of the response of the carbon cycle to rising CO_{2} is determined by the turnover times of the various biogeochemical processes by which CO_{2} is sequestered and once more decomposed in the various ecosystems. It is well known that in particular our incomplete knowledge of the internal timescales of the land carbon cycle is severely limiting our ability to predict the development of the land carbon sink (ToddBrown et al., 2013; Friend et al., 2014; Exbrayat et al., 2014; Koven et al., 2015; He et al., 2016; Yan et al., 2017; Zhou et al., 2018). To improve the understanding of the dynamics of coupled climate–carboncycle system one thus needs to investigate together with the feedbacks also the issue of internal timescales.
Concerning the analysis of feedbacks, a large step forward was the seminal work by Friedlingstein et al. (2003), who proposed a mathematical framework to quantify their contributions to the response. The basic idea underlying their α–β–γ framework is that at global scale one can identify two types of climate–carboncycle feedback: a “biogeochemical feedback”, which arises through the direct effect of changes in atmospheric CO_{2} concentrations on global carbon stocks, and a “radiative feedback”, which affects the carbon cycle indirectly from the change in climate that arises from the radiative forcing of CO_{2} perturbations. Key elements of this framework to quantify the two types of feedback – also known as carbonconcentration and carbon–climate feedback (Arora et al., 2013) – are the β and γ sensitivities that, respectively, characterize the response of stored land and/or ocean carbon to changes in CO_{2} and in climate. As in studies of the physical system by means of “pattern scaling” (e.g. Mitchell, 2003), climate is in this framework collectively represented by the single quantity temperature – whose sensitivity to CO_{2} perturbations is quantified by α. The framework of Friedlingstein et al. (2003) led to important insights into the role of climate–carboncycle feedbacks for climate change and stimulated a vast amount of research in the field (Friedlingstein et al., 2006; Gregory et al., 2009; Arneth et al., 2010; Zickfeld et al., 2011; Boer and Arora, 2013; Arora et al., 2013; Schwinger et al., 2014; Friedlingstein et al., 2014; Friedlingstein, 2015; Adloff et al., 2018; Williams et al., 2019; Goodwin et al., 2019; Jones and Friedlingstein, 2020).
Quantitative results on the size of global climate–carboncycle feedbacks were particularly obtained as part of the Coupled Climate–Carbon Cycle Model Intercomparison (C^{4}MIP) project (Friedlingstein et al., 2013–2016; Arora et al., 2020) from Earth system simulations tailored for the application of the α–β–γ framework. In terms of radiative forcing, when raising atmospheric CO_{2} at a fixed rate of 1 % per year from its preindustrial value to 4 times this value, the negative biogeochemical feedback is more than 4 times stronger than the positive radiative feedback (Canadell et al., 2021). This results in a net land and ocean carbon sequestration over the whole simulation period of between 33 % and 50 % across the participating models (Arora et al., 2020). Most of the spread in these numbers arises from differences in the representation of land carbon cycling in the various models, in particular from the internal timescales assumed for the turnover of vegetation and soils (Arora et al., 2020).
The size of these feedbacks depends on the considered timescale (e.g. Enting, 2022). In the original α–β–γ framework, this timescale dependence shows up only implicitly (see discussion in Sect. 2), so results from this framework are specific to the considered perturbation scenario (Gregory et al., 2009; Torres Mendonça et al., 2021b). This limitation is overcome by the recently proposed generalization of the α–β–γ framework (Heimann, 2014; Rubino et al., 2016; Enting and Clisby, 2019; Enting, 2022) that instead explicitly quantifies the timescale dependence of the climate–carboncycle feedbacks independently of the scenario. Here generalized sensitivities α, β, and γ are introduced as timedependent linear response functions, where the term “linear” indicates that they specify the response only to linear order in a Volterra expansion of the response into the perturbation (see Torres Mendonça et al., 2021b); practically this means that this approach applies only as long as the perturbations are sufficiently weak. As will be discussed below, in principle this generalization gives – at a globally aggregated level – a complete description of the linear dynamics of the coupled climate and carbon cycle in terms of the responses and feedbacks at different timescales.
In the present study we employ this generalized framework to study the role of feedbacks and their timescale dependence for airborne fraction. Airborne fraction is classically defined as the fraction of emitted CO_{2} that stays in the atmosphere after accounting for the induced land and ocean uptake. Accordingly, airborne fraction quantifies the amount of emissions that effectively contributes to climatic change and is therefore a key quantity of climate research (Oeschger and Heimann, 1983). It is closely related to the climate–carboncycle feedbacks because, as discussed above, the reduction in atmospheric CO_{2} caused by ocean and land uptake depends itself on the changes in atmospheric CO_{2} induced by the emissions. Because of its importance, much effort has been put into quantifying and investigating the airborne fraction (e.g. Canadell et al., 2007; Raupach et al., 2008; Archer et al., 2009; Gregory et al., 2009; Gloor et al., 2010; Jones et al., 2013; Le Quéré et al., 2009; Raupach et al., 2014; Bennedsen et al., 2019). The remarkable constancy of the airborne fraction over the last decades – a consequence of the abovementioned constancy of the land and ocean carbon storage fraction – indicates that land and ocean uptake has not saturated yet, so their response is still linear. But as with the original α, β, and γ sensitivities, this classical definition of airborne fraction also neglects that its value must depend on the internal timescales at which the land and ocean carbon cycles react to emissions (see e.g. Enting, 2007). Addressing this timescale dependence, Raupach (2013) argues that the observed constancy of the airborne fraction is actually only a reflection of the linearity of the response of land and ocean carbon reservoirs together with the exponential nature of the historical rise in anthropogenic emissions. Thus, as anthropogenic emissions cease to behave exponentially – for instance as a consequence of a cut in emissions – the airborne fraction is expected to deviate from its historical value. Hence, as in the case of the feedbacks quantified by the α–β–γ framework, the airborne fraction in its standard definition cannot be seen as an invariant property of the carbon cycle alone but only as a metric that depends on the considered perturbation scenario. This is different when tackling airborne fraction by means of the generalized α–β–γ framework. As demonstrated by Rubino et al. (2016) and Enting and Clisby (2019) when studying the prehistorical and historical development of airborne fraction, by such an approach not only α, β, and γ but also airborne fraction can be generalized to a dynamic quantity that describes the response of the coupled climate–carbon system to emissions at its various internal timescales for any sufficiently weak perturbation scenario and is thus a true property of the coupled climate–carbon system itself.
To gain further insight into the timescale dependence of the airborne fraction and in particular how this timescale dependence emerges from the underlying feedbacks, in the present study we analyse by means of the generalized α–β–γ framework published simulations of different Earth system models. More precisely, we analyse C^{4}MIPtype simulation results from models participating in the fifth phase of the Coupled Model Intercomparison Project (CMIP5; see Taylor et al., 2012). These C^{4}MIPtype simulations were designed to separately quantify the biogeochemical and radiative feedbacks from the original α, β, and γ values. We use these simulations to obtain for the respective models the linear response functions of the generalized framework necessary to study the timescale dependence of airborne fraction. For their robust recovery, we employ the response function identification method that we developed for this purpose in Torres Mendonça et al. (2021a).
Most of our present study relies on a single theoretical result of the generalized α–β–γ framework, namely the relation expressing the timescaleresolved airborne fraction exclusively by the underlying feedbacks (see Eq. 15 below). In this formula the feedbacks are represented by timescaledependent functions that in turn are fully determined by the generalized α, β, and γ sensitivities. By this relation the generalized framework makes a strong claim, namely that airborne fraction (a quantity fully determined by feedbacks) can be predicted from the generalized sensitivities (quantities that are independent of feedbacks). The validity of this relation fully relies on the assumption that at climate timescales the overall feedback is dominated by the pair of biogeochemical and radiative feedbacks, which with respect to the latter includes the assumption that nearsurface temperature is a good measure to represent also all other climate aspects, in particular those related to the eminently important hydrological cycle. As a consequence, in order to employ this relation in our study to derive the timescale dependence of airborne fraction via the underlying feedbacks from the generalized sensitivities, we first demonstrate the predictive power of this relation when applied to Earth system simulations. We will do this exemplarily for the Max Planck Institute Earth System Model (MPIESM) by performing additional simulations not available for the other CMIP5 models. This demonstration is also interesting on its own because so far theoretical inferences from the generalized α–β–γ framework have not been tested, although this is a prerequisite for its faithful application.
Overall, our analysis of the simulation results from the considered set of CMIP5 models will show that one can understand the dynamics of the airborne fraction from the behaviour of the climate–carboncycle feedbacks and that it is possible to pinpoint the particular feedback that dominates the observed model spread in the airborne fraction at different timescales. Moreover, it will become clear that by applying the generalized α–β–γ such results are scenarioindependent, although they are obtained from simulations performed for particular scenarios.
The outline of the paper is as follows. In the next section we introduce the generalized α–β–γ framework that underlies our whole analysis. Then we demonstrate in Sect. 3 its predictive power exemplarily for MPIESM. This part of the study serves also a second methodological purpose: here we develop and test through the example of MPIESM all the necessary numerical details to derive from simulation data via the generalized sensitivities α, β, and γ the timescaleresolved airborne fraction also for other models in the subsequent Sect. 4. This section then contains the main part of our study where we investigate the timescale dependence of airborne fraction and the underlying feedbacks for an ensemble of CMIP5 models. In the final Sects. 5, 6, and 7 we summarize and discuss our results. To keep the paper more readable, the extensive technical parts of our investigations are presented in the appendices.
To prepare for our investigation of the timescale dependence of the airborne fraction and the underlying feedbacks, we introduce here the generalized α–β–γ framework (Heimann, 2014; Rubino et al., 2016; Enting and Clisby, 2019; Enting, 2022).
To introduce this framework, we start from the carbon balance equation that couples the different subsystems of the global carbon cycle:
where ΔC_{A}, ΔC_{L}, and ΔC_{O} are the differences in global carbon content in the atmosphere, land, and ocean with respect to preindustrial times (formally denoted here and below as t=0 in all equations), and E(t) is the timedependent flux of (anthropogenic) carbon emissions. Following the framework of Friedlingstein et al. (2003), the land and ocean carbon changes in Eq. (1) are assumed to depend only on the atmospheric CO_{2} concentration (driving the biogeochemical response) and on temperature (driving the radiative response). Considering these changes as separate responses to CO_{2} and temperature perturbations, they can for the weak perturbations assumed here be approximated as the linear term of a Volterra expansion around the preindustrial state (see Torres Mendonça et al., 2021b) and thus must also combine linearly for this order of approximation to give
Here Δc and ΔT are the changes in CO_{2} concentration and in land (L) and ocean (O) global mean nearsurface air temperature with respect to the preindustrial state, while ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{L}\right)}$, ${\mathit{\chi}}_{\mathit{\gamma}}^{\left(\text{L}\right)}$, ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{O}\right)}$, and ${\mathit{\chi}}_{\mathit{\gamma}}^{\left(\text{O}\right)}$ are the linear response functions that generalize the land and ocean sensitivities β^{(L)}, γ^{(L)}, β^{(O)}, and γ^{(O)} of the original α–β–γ framework. Additionally, the temperature can in the same way be understood as a response to weak CO_{2} perturbations:
where ${\mathit{\chi}}_{\mathit{\alpha}}^{\left(\text{L}\right)}$ and ${\mathit{\chi}}_{\mathit{\alpha}}^{\left(\text{O}\right)}$ are the linear response functions that generalize the land and ocean sensitivities α^{(L)} and α^{(O)}. Following Torres Mendonça et al. (2021b), the response functions that generalize the α, β, and γ sensitivities are referred to as generalized sensitivities.
As discussed above, atmospheric CO_{2} depends not only on the emissions perturbing the system but also on the responses of land and ocean CO_{2} exchanges induced by them. How land and ocean carbon and thus also the exchange fluxes react to such perturbations is characterized by the generalized sensitivities just introduced. With their help, a compact equation relating the change in atmospheric carbon content to the perturbing emissions including all feedbacks is obtained by employing the last two equations to eliminate temperature in Eqs. (2) and (3), inserting these into the carbon balance Eq. (1) and using the known relation C_{A}(t)=kc(t), k=2.12 Pg C ppm^{−1} (CO_{2}) (Ciais et al., 2013, p. 417), between atmospheric CO_{2} concentration and carbon content. After applying a Laplace transform to the resulting integrodifferential equation one obtains
where the tilde denotes Laplacetransformed quantities. Applying the Laplace transform has the technical advantage that linear integral equations turn into linear algebraic equations for the transformed quantities that are much easier to handle – e.g. one can directly solve the last equation for $\mathrm{\Delta}{\stackrel{\mathrm{\u0303}}{C}}_{\text{A}}\left(p\right)$ (see Eq. 7 below). The other, somewhat challenging consequence is that the interpretation of Laplacetransformed quantities is less intuitive because they are not functions of time but of the rates p, whose inverses $\mathrm{1}/p$ should be understood as timescales. Similar Laplacetransformed equations were derived in Enting (2007), Rubino et al. (2016), and Enting and Clisby (2019).
In the absence of feedbacks the atmospheric change in carbon content $\mathrm{\Delta}{\stackrel{\mathrm{\u0303}}{C}}_{\text{A}}\left(p\right)$ would be fully determined by the cumulated emissions $\stackrel{\mathrm{\u0303}}{E}\left(p\right)/p$ (this is the Laplace transform of ∫E(s)ds, $E(t\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}=\phantom{\rule{0.125em}{0ex}}\phantom{\rule{0.125em}{0ex}}\mathrm{1850})\phantom{\rule{0.125em}{0ex}}=\phantom{\rule{0.125em}{0ex}}\mathrm{0}$) in Eq. (6). Thus the feedbacks enter by the second righthandside term. This additional contribution to $\mathrm{\Delta}{\stackrel{\mathrm{\u0303}}{C}}_{\text{A}}\left(p\right)$ that depends on the response itself characterizes the total climate–carboncycle feedback (Roe, 2009). This gets even clearer by rewriting Eq. (6) analogously to the feedback equations of the original α–β–γ framework (Roe, 2009; Gregory et al., 2009; Adloff et al., 2018): setting
with
one obtains from Eq. (6) by termwise identification
In this way the full information on how atmospheric carbon $\mathrm{\Delta}{\stackrel{\mathrm{\u0303}}{C}}_{\text{A}}\left(p\right)$ (and thus atmospheric CO_{2}) responds to a trajectory of emissions $\stackrel{\mathrm{\u0303}}{E}\left(p\right)$ is contained in the function $\stackrel{\mathrm{\u0303}}{f}\left(p\right)$, which we call, following the terminology of Roe (2009) and Adloff et al. (2018), the total feedback function. As also explained there, from the point of view of feedback analysis $\mathrm{1}/(\mathrm{1}\stackrel{\mathrm{\u0303}}{f}(p\left)\right)$ is the gain of the system: depending on whether $\mathrm{1}/(\mathrm{1}\stackrel{\mathrm{\u0303}}{f}(p\left)\right)$ is larger or smaller than 1, the inclusion of the feedbacks amplifies or dampens the response of atmospheric CO_{2} in comparison to a reference system that would simply store the cumulated emissions without reaction in land and ocean fluxes^{1}. In other words, depending on whether $\mathrm{1}/(\mathrm{1}\stackrel{\mathrm{\u0303}}{f}(p\left)\right)$ is larger or smaller than 1, the inclusion of the feedbacks results in CO_{2} fluxes into or out of the atmosphere in addition to the emissions.
The total feedback function $\stackrel{\mathrm{\u0303}}{f}\left(p\right)$ is defined in Eq. (8) as the sum of the feedback functions ${\stackrel{\mathrm{\u0303}}{f}}_{\mathit{\beta}}^{\left(\text{L}\right)}\left(p\right)$, ${\stackrel{\mathrm{\u0303}}{f}}_{\mathit{\gamma}\mathit{\alpha}}^{\left(\text{L}\right)}\left(p\right)$, ${\stackrel{\mathrm{\u0303}}{f}}_{\mathit{\beta}}^{\left(\text{O}\right)}\left(p\right)$, and ${\stackrel{\mathrm{\u0303}}{f}}_{\mathit{\gamma}\mathit{\alpha}}^{\left(\text{O}\right)}\left(p\right)$. These functions quantify for land and ocean the biogeochemical $\left({\stackrel{\mathrm{\u0303}}{f}}_{\mathit{\beta}}^{\left(\text{L}\right)}\text{and}{\stackrel{\mathrm{\u0303}}{f}}_{\mathit{\beta}}^{\left(\text{O}\right)}\right)$ and the radiative $\left({\stackrel{\mathrm{\u0303}}{f}}_{\mathit{\gamma}\mathit{\alpha}}^{\left(\text{L}\right)}\text{and}{\stackrel{\mathrm{\u0303}}{f}}_{\mathit{\gamma}\mathit{\alpha}}^{\left(\text{O}\right)}\right)$ feedback, also known as carbonconcentration and carbon–climate feedback (Arora et al., 2013). The feedback functions, in turn, are determined from the generalized sensitivities via Eqs. (9) and (10).
Concerning the timescale dependence it is important to note that in Eq. (6) all terms depend on the same rate p, which means that at each timescale $\mathrm{1}/p$ the response in atmospheric carbon to the forcing by emissions is fully determined by the properties of the system at that very timescale alone. This independent behaviour at different timescales is a consequence of the assumption that the forcing is sufficiently weak so that the system behaviour is already well approximated by the linear term in the Volterra expansions of the response in the perturbations (Eqs. 2 to 5) when Laplace transformed; taking the Volterra expansion to higher order would introduce terms involving mixed timescales (see e.g. Schetzen, 2010, Eqs. 2.1, 2.2).
Such independent behaviour is also the reason for the identical structure of the Laplacetransformed formulas of the generalized α–β–γ framework and those of the original framework in the time domain (Heimann, 2014), which gets obvious by comparing the relation between atmospheric carbon and emissions (Eqs. 7–10) with its analogue from the original framework (Gregory et al., 2009; Adloff et al., 2018; Jones and Friedlingstein, 2020):
where the t argument emphasizes the time dependence of α, β, and γ (Adloff et al., 2018). But despite this striking similarity, these are fundamentally different formulations: while Eq. (11) is a diagnostic way of writing the response of atmospheric carbon to emissions by means of sensitivities that generally differ for different scenarios, Eqs. (7)–(10) predict this response for any (weak) emissions scenario by means of unique system properties – the generalized sensitivities, which are completely independent of the scenario. It is this predictive power that we test in the next section (see more below).
Note also that the timescale dependence of feedbacks cannot be obtained from the original α–β–γ framework, even if one computes the α, β, and γ sensitivities underlying its feedback quantification as timedependent quantities α(t), β(t), and γ(t) as above. To understand this one must realize that, in contrast to the original α–β–γ framework where feedbacks are quantified for a particular scenario, in the generalized framework feedback strengths are internal properties of the climate–carbon system, consistent with the viewpoint that these strengths depend only on internal system characteristics (such as the sensitivity of soil microbial activity to changes in temperature or the sensitivity of plant photosynthesis to changes in CO_{2} concentrations). If one understands feedbacks in this way, it gets clear that by calculating the timedependent α(t), β(t), and γ(t) sensitivities and combining them to quantify feedbacks one is obtaining only implicit information on the timescaledependent feedback strengths because the combined values of these sensitivities reflect not internal system feedbacks alone but also the external forcing scenario (e.g. Gregory et al., 2009; Boer and Arora, 2013; Arora et al., 2013; Torres Mendonça et al., 2021b). Accordingly, from the time dependence of those sensitivities one cannot infer the timescale dependence of the feedback strengths. In the generalized framework contributions from forcing and feedback are disentangled so that the timescale dependence of the climate–carboncycle feedbacks is instead explicitly quantified. This more general quantification of feedback strengths, which arises when considering the internal memory of the climate–carbon system, may be even more clearly understood by noting that the α, β, and γ sensitivities can be predicted by their generalized counterparts for any weak perturbation scenario (see Sect. 4.1).
Our main topic in this study is the timescale dependence of the airborne fraction. As explained in the following, by the generalized framework this timescale dependence can be fully traced back to that of the feedback functions. In its standard definition (e.g. Oeschger and Heimann, 1983; Raupach, 2013), the airborne fraction AF(t) is specified by
where the lefthand side is the rate at which emitted carbon accumulates in the atmosphere. Airborne fraction obtained its name because this accumulation rate can also be viewed as the fraction of the emitted carbon flux that remains airborne. As already discussed in the Introduction, despite its constancy over the last decades, AF(t) cannot be seen as a property of the climate–carbon system itself but only as a metric dependent on the emissions scenario E(t). But following an analogous strategy as for the α, β, and γ sensitivities, a scenarioindependent generalized airborne fraction – denoted by A(t) – can be obtained by expanding $\mathrm{d}{C}_{\text{A}}/\mathrm{d}t$ in a Volterra series up to linear order in the emissions E(t) around the preindustrial state ($\mathrm{d}{C}_{\text{A}}/\mathrm{d}t=\mathrm{0}$). Following Enting (2007), the standard definition (12) of airborne fraction thereby generalizes to
Compared to Eq. (12), this generalized response formula accounts not only for the effect of emissions at time instant t but also for their effect during their whole past history. Accordingly, Eq. (13) accounts for the memory of the carbon cycle, and having introduced the generalized airborne fraction A(t) as the kernel of the linear term of a Volterra expansion about the preindustrial state, A(t) is for weak perturbations a property of the system itself independent of the emissions scenario E(t). This generalized airborne fraction is a generalization not only of the standard airborne fraction AF(t) but also of the cumulative airborne fraction^{2} CAF(t).
To relate the generalized airborne fraction to the feedbacks, Eq. (13) is first Laplace transformed. Using then $\mathrm{d}{C}_{\text{A}}/\mathrm{d}t=\mathrm{d}\mathrm{\Delta}{C}_{\text{A}}/\mathrm{d}t$ and noting that a standard property of the Laplace transform ℒ of the time derivative of ΔC_{A} is that $\mathcal{L}\mathit{\{}\mathrm{d}\mathrm{\Delta}{C}_{\text{A}}/\mathrm{d}t\mathit{\}}=p\mathcal{L}\mathit{\left\{}\mathrm{\Delta}{C}_{\text{A}}\right(t\left)\mathit{\right\}}$ for ${lim}_{t\to {\mathrm{0}}^{+}}\mathrm{\Delta}{C}_{\text{A}}\left(t\right)=\mathrm{0}$ (where “+” indicates the onesided limit of t approaching zero from positive t values), one obtains
By comparing this with Eq. (7), one finds that the generalized airborne fraction is identical to what was above called gain:
This is the key relation underlying our study. It demonstrates that at each timescale the generalized airborne fraction is fully determined by the values of the total feedback function $\stackrel{\mathrm{\u0303}}{f}\left(p\right)$ at that very timescale $\mathrm{1}/p$. An analogous relation was obtained by Gregory et al. (2009) and Jones and Friedlingstein (2020) employing the original α–β–γ framework and by Rubino et al. (2016) and Enting and Clisby (2019) in this generalized form.
To follow our subsequent investigation of airborne fraction it is important to note that the generalized α–β–γ framework is more than only a way to describe the coupled climate–carbon system at global scale: actually it is a theory about this system with predictive power. Basic to this whole framework is the generalization of the original α, β, and γ sensitivities to response functions. Already by this first step some predictive power is gained because once these generalized sensitivities are known, the response to any sufficiently weak CO_{2} perturbation scenario can be predicted (Torres Mendonça et al., 2021b). But more important for the present study is another type of predictive power of the generalized framework that arises by describing the climate–carbon system in terms of assumed key feedbacks: by Eq. (15) the airborne fraction is via Eqs. (8) to (10) fully determined by the generalized sensitivities ${\stackrel{\mathrm{\u0303}}{\mathit{\chi}}}_{\mathit{\alpha}}$, ${\stackrel{\mathrm{\u0303}}{\mathit{\chi}}}_{\mathit{\beta}}$, and ${\stackrel{\mathrm{\u0303}}{\mathit{\chi}}}_{\mathit{\gamma}}$. These characterize the responses of subsystems to specific perturbations and have at first sight nothing to do with feedbacks. The feedbacks come about only by their combined action as described by the generalized framework. This is particularly obvious when considering the airborne fraction: as explained in the Introduction, this quantity embodies by its very nature the effect of all the ruling feedbacks. Hence, our key Eq. (15) predicts a quantity shaped by feedbacks – the airborne fraction – from quantities that are independent of feedbacks – the generalized sensitivities. Thereby naturally the question arises of how good such a prediction will be. This question will be answered in the next section exemplarily for simulations performed with the MPIESM. Finding a good predictability will justify deriving the airborne fraction by Eq. (15) merely from the knowledge of the generalized sensitivities also for other CMIP5 Earth system models in the main part of this study.
The present section prepares for the main investigation of our study (next section). This involves two issues. The first was already shortly addressed at the end of the previous section, namely that we have to demonstrate the predictive power of Eq. (15) before we can reliably use it to calculate the generalized airborne fraction $\stackrel{\mathrm{\u0303}}{A}\left(p\right)$. We demonstrate this by first determining $\stackrel{\mathrm{\u0303}}{A}\left(p\right)$ directly by its definition (13) from simulated atmospheric CO_{2} and then comparing it with its prediction obtained from Eq. (15) via the generalized sensitivities ${\stackrel{\mathrm{\u0303}}{\mathit{\chi}}}_{\mathit{\alpha}}$, ${\stackrel{\mathrm{\u0303}}{\mathit{\chi}}}_{\mathit{\beta}}$, and ${\stackrel{\mathrm{\u0303}}{\mathit{\chi}}}_{\mathit{\gamma}}$. While the sensitivities necessary for the prediction can in principle be calculated for all considered CMIP5 models from published simulation results (see technical issues below), to obtain $\stackrel{\mathrm{\u0303}}{A}\left(p\right)$ directly from simulated atmospheric CO_{2} one needs additional simulations. Therefore we restrict our demonstration to MPIESM for which we perform these additional simulations. It should be noted that this demonstration is also interesting on its own because the validity of inferences from the generalized framework has so far never been demonstrated. For conciseness, in the following we call the airborne fraction computed by application of its definition (13) from simulated atmospheric CO_{2} the true airborne fraction, while the airborne fraction calculated via Eq. (15) of the generalized framework from simulated generalized sensitivities the predicted airborne fraction.
The second preparatory issue tackled in this section concerns the technical aspects of the calculation of the generalized sensitivities from simulation data. As explained in Torres Mendonça et al. (2021a), this calculation is generally not trivial. There are two main complications.
 i.
Noise. The problem of deriving response functions such as the generalized sensitivities from perturbation experiments data is mathematically “illposed”. In practice this means that by trying to solve it by classical numerical methods one obtains sensitivities corrupted by the noise in the data.
 ii.
Nonlinearities. The generalized framework is a linear approach. Therefore, when recovering the generalized sensitivities, one has to make sure that the response data are not contaminated by strong nonlinear contributions that would otherwise hinder the recovery.
To derive the generalized sensitivities we employ our recently developed response function identification (RFI) method (Torres Mendonça et al., 2021a, b). This method recovers the generalized sensitivity from a single realization of an arbitrary perturbation experiment. In particular the RFI method addresses complication (i): it filters out the noise in the recovered generalized sensitivity by means of Tikhonov–Phillips regularization (Phillips, 1962; Tikhonov, 1963), with the regularization parameter determined via the discrepancy principle (Morozov, 1966) from an estimate of the noise level in the data (obtained from the associated control simulation).
To address complication (ii), we employ following Torres Mendonça et al. (2021b) two additional procedures: first, we pretransform the data by different techniques to try to account for known nonlinearities in the response for which we want to derive the generalized sensitivity (e.g. the response of land carbon to changes in CO_{2} concentrations characterized by ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{L}\right)}\left(t\right)$ – first term in Eq. 2); second, by means of additional perturbation experiments we estimate the maximum perturbation strength limiting the extent of the linear regime of that response. By accounting for known nonlinearities in the response, the first procedure allows one to take data at higher perturbation strengths and thus higher signaltonoise ratio, which leads to a better recovery of the generalized sensitivity. Therefore another technical aspect of the present section is to determine for each generalized sensitivity the preprocessing technique that gives the best recovery. The second procedure assures that the taken data contain no strong nonlinearities that could hinder the recovery. This second procedure serves also a different purpose: to estimate the range of perturbation strengths for which the generalized framework as a whole is applicable. Since this range is generally different for the different involved responses (i.e. for each term on the righthand side of Eqs. 2–5), the linear regime of the generalized framework as a whole is determined by the smallest of the maximum perturbation strengths limiting the linear regime of those responses separately. Since for employing these procedures additional experiments are needed, we restrict our analysis to the MPIESM – for which we perform those experiments – and assume in the next section that the preprocessing techniques and ranges of linearity identified for MPIESM apply also to other CMIP5 models.
To demonstrate the predictive power of the generalized framework, all these technical issues must be tackled before we can invoke Eq. (15) to reliably compute the predicted airborne fraction. We tackle them by following the procedures given in Torres Mendonça et al. (2021b). Since these technical parts of our study reveal no further scientific insight, we have relegated their rather lengthy description to Appendix A. The obtained results concerning the size of the linear regime and the best preprocessing technique are summarized in Table A2.
3.1 Determining the true airborne fraction from simulated atmospheric CO_{2}
As explained above, to demonstrate that indeed the timescaleresolved airborne fraction $\stackrel{\mathrm{\u0303}}{A}\left(p\right)$ is reliably predicted by Eq. (15) of the generalized framework, we compare it with the true airborne fraction calculated by application of its defining Eq. (13). This section explains how to obtain this true airborne fraction from a simulation with prognostic atmospheric CO_{2}, i.e. from an emissiondriven simulation.
From given time series for atmospheric carbon content C_{A}(t) and emissions E(t) one could in principle obtain $\stackrel{\mathrm{\u0303}}{A}\left(p\right)$ by solving the defining Eq. (13) by means of our RFI method for A(t) followed by a Laplace transform. But to proceed in this way one had first to calculate $\mathrm{d}{C}_{\text{A}}/\mathrm{d}t$ from C_{A}(t) (compare Eq. 13), which introduces numerical noise that deteriorates the quality of recovery of A(t) (Torres Mendonça et al., 2021a). Therefore we proceed differently. To linear order a Volterra expansion of C_{A} into the perturbing emissions gives
which defines another response function^{3} χ_{ζ}(t). A Laplace transform then gives (Enting, 1990)
Note that, in contrast to Eq. (14), no p factor shows up in Eq. (17). This is essentially because Eq. (13), from which Eq. (14) is obtained, has on the lefthand side the time derivative of the lefthand side of Eq. (16) (from which Eq. 17 is obtained), and the Laplace transform of this time derivative is $p\mathrm{\Delta}{\stackrel{\mathrm{\u0303}}{C}}_{\text{A}}\left(p\right)$ (see text introducing Eq. 14).
By comparing now Eq. (17) with the Laplacetransformed definition of the generalized airborne fraction (14), one thus obtains (as also noted by Enting and Clisby, 2019)
By these considerations, the true airborne fraction $\stackrel{\mathrm{\u0303}}{A}\left(p\right)$ can be determined from emissiondriven simulations in three steps: first, solve Eq. (16) by our RFI method for χ_{ζ}(t) using simulation data for ΔC_{A}(t) and E(t) (no numerical derivative needed). Second, Laplacetransform the recovered χ_{ζ}(t) to obtain ${\stackrel{\mathrm{\u0303}}{\mathit{\chi}}}_{\mathit{\zeta}}\left(p\right)$. Finally, apply Eq. (18) to determine $\stackrel{\mathrm{\u0303}}{A}\left(p\right)$ from ${\stackrel{\mathrm{\u0303}}{\mathit{\chi}}}_{\mathit{\zeta}}\left(p\right)$.
For our demonstration of predictive power we performed impulsetype emissiondriven experiments with MPIESM and obtained $\stackrel{\mathrm{\u0303}}{A}\left(p\right)$ from the resulting simulation data following these three steps (see Appendices B and C for details). The resulting true $\stackrel{\mathrm{\u0303}}{A}\left(p\right)$ is plotted in Fig. 1.
3.2 Determining the predicted airborne fraction from generalized α–β–γ sensitivities
The second step towards the demonstration of the predictive power of the generalized framework is to calculate the predicted airborne fraction by application of Eq. (15) from the generalized sensitivities. For a proper comparison with the true airborne fraction from above, this predicted airborne fraction is obtained from simulations with the same model (MPIESM) as the true airborne fraction although from different simulation experiments.
To predict airborne fraction via the total feedback function $\stackrel{\mathrm{\u0303}}{f}\left(p\right)$ from Eq. (15) one needs to know the generalized sensitivities (see Eqs. 8 to 10). These we obtain from two standard C^{4}MIPtype experiments performed with MPIESM that are published in the international CMIP5 repository (see Appendix A for details). These C^{4}MIPtype simulations were tailored for separate determination of the α, β, and γ sensitivities of the original framework (Taylor et al., 2012; Arora et al., 2013) but are similarly suited for separate determination of their generalized counterparts (Torres Mendonça et al., 2021b). In both of the experiments atmospheric CO_{2} concentration is prescribed to rise from its preindustrial level by 1 % per year, but for separate identification of the different sensitivities this rising CO_{2} is made to act differently in the two simulations: in the radiatively coupled (“rad”) simulation the CO_{2} rise affects only the atmospheric radiation code, while in the biogeochemically coupled (“bgc”) simulation the rise in CO_{2} affects only biogeochemical processes (ocean pCO_{2}, leaf CO_{2}); in both simulations for the respective other aspect CO_{2} stays at its preindustrial level.
To determine the generalized sensitivities from the simulation data we once more invoke our RFI method (Torres Mendonça et al., 2021a, b). In the bgc simulation, climate change is largely suppressed so that the changes in ocean and land carbon are to first order determined only by the rising CO_{2}. The rather small change in land temperature in this simulation arises by various indirect effects, among them by a reduction in transpiration cooling because of the closing of plant stomata under higher CO_{2} (Arora et al., 2013). Ignoring this comparably small temperature rise, one can assume ΔT=0 in Eqs. (2) and (3) so that only the integrals over the rising CO_{2} remain in these equations. These are the equations that we solve by means of the RFI method for the ocean and land χ_{β}(t) sensitivities using the data for the rising CO_{2} and the stored land and ocean carbon from the bgc simulation. The generalized α and γ sensitivities are obtained from the rad simulation. In this simulation the effect of rising CO_{2} on the carbon chemistry is missing; i.e. stored land carbon and ocean carbon change only because of changing climate, collectively represented by temperature in the α–β–γ framework. Hence in Eqs. (2) and (3) one can now drop the integrals over rising CO_{2} so that only the integrals over rising temperature remain; these reduced equations are then solved by the RFI method for the ocean and land χ_{γ}(t) under the integral. Finally, because the α sensitivities measure the direct response of rising CO_{2} via its greenhouse effect on land and ocean temperatures, these sensitivities are determined from the rad simulation as well. The respective equations to be solved for the χ_{α}(t) sensitivities are (4) and (5). Note that the sensitivities need not be derived from bgc and rad simulations: also either of these two together with a “fully coupled” experiment – where both the radiation and the biogeochemical code are affected by CO_{2} changes – suffices (see e.g. Arora et al., 2020).
Actually, as already explained in the introduction to this section, to obtain linear response functions reliably by the RFI method, additional preparatory effort is needed concerning selection of a preprocessing technique and checks assuring that the underlying linearity assumption is valid for the simulation data used. For those purposes we performed additional rad and bgc simulations with MPIESM for a variety of different CO_{2} forcing scenarios (see Table A1 in Appendix A). Based on this preparatory analysis we then obtain the generalized sensitivities of MPIESM in the time domain (see Appendix A).
The final steps to obtain the generalized airborne fraction as predicted by Eq. (15) from the generalized framework are to Laplacetransform the obtained generalized sensitivities (done analytically; see Torres Mendonça et al., 2021a), calculate from Eqs. (8)–(10) the total feedback function f(p), and then finally obtain through our key Eq. (15) the predicted timescaleresolved airborne fraction. The result is seen in Fig. 1.
Please note that the way of deriving here the predicted generalized airborne fraction for MPIESM is exactly how we derive it also for the other CMIP5 models in the next section, except that the additional preparatory analysis and checks cannot be performed because of the lack of the necessary additional simulations. Accordingly, we will assume that the size of the linear regime obtained for MPIESM applies also to these other models and will preprocess also their data by the technique identified to be best for MPIESM (see summary of linear regime and best preprocessing technique for each sensitivity in Table A2).
3.3 Demonstration of predictive power by comparing predicted with true airborne fraction
So far in this section, the generalized airborne fraction $\stackrel{\mathrm{\u0303}}{A}\left(p\right)$ has been derived for MPIESM in two ways: first directly from simulated atmospheric CO_{2} (“true” airborne fraction) and then by employing Eq. (15) of the generalized framework (“predicted” airborne fraction). The results are plotted in Fig. 1. The two curves differ by up to around 5 % for timescales below 5 years but match well for the longer timescales up to 100 years^{4} (mean difference less than 1 %). To judge from the closeness of the two curves how well airborne fraction is predicted one should note that their coincidence for $\mathrm{1}/p\to \mathrm{0}$ is not a hint for good predictive power but merely a hint of the reliability of the numerics by which the curves were obtained: from carbon conservation it follows that χ_{ζ}(0)=1 (see Appendix B) so that
where the second equality follows from the initial value theorem of Laplace transforms (e.g. Beerends et al., 2003, p. 292). Hence the two curves match at short timescales for theoretical reasons and not because of the quality of the prediction. More insight into the behaviour of $\stackrel{\mathrm{\u0303}}{A}\left(p\right)$ will be given in Sect. 4 when discussing its dependence on the climate–carboncycle feedbacks calculated for the CMIP5 models.
The discrepancy between the two estimates of the airborne fraction observed at timescales shorter than 5 years is expected from two types of error that might have affected the results. The first type affects the predicted airborne fraction and arises from the illposedness of the deconvolution problem that must be solved to derive the generalized sensitivities employed in Eq. (15). This illposedness obscures information at short timescales and therefore deteriorates the recovery of the sensitivities at those scales (see Torres Mendonça et al., 2021a). The second type of error affects the true airborne fraction and arises from the fact that χ_{ζ}(t), from which then $\stackrel{\mathrm{\u0303}}{A}\left(p\right)$ was derived via Eq. (18), was not obtained from a perfect impulse experiment (see details in Appendix B). Although we derived χ_{ζ}(t) in Appendix B enforcing its known value ${\mathit{\chi}}_{\mathit{\zeta}}(t=\mathrm{0})=\mathrm{1}$ as a numerical constraint to partially account for this error, the recovery of χ_{ζ}(t) at short timescales might still not be fully correct. Despite these discrepancies, as timescales lower towards $\mathrm{1}/p=\mathrm{0.01}$ years the agreement improves once more, in line with the theoretical expectation discussed above.
Overall, the close agreement between the two estimates of the airborne fraction demonstrates the predictive power of the generalized α–β–γ framework: since $\stackrel{\mathrm{\u0303}}{A}\left(p\right)$ encodes all information needed to predict atmospheric carbon response to any (weak) emission scenario – accounting for all climate–carboncycle feedbacks – this close agreement shows that, at least for MPIESM, the generalized framework correctly predicts at global scale the linear dynamics of the coupled climate–carbon system. In addition, because the two curves were obtained from very different simulations, their agreement adds confidence that the numerical methods employed can be trusted. These two results suggest that the generalized framework and our numerical methods are appropriate to confidently predict the airborne fraction and the underlying climate–carboncycle feedbacks by Eq. (15) of the generalized framework from the concentrationdriven CMIP5 experiments, as we do in the next section.
In the present section we extend our analysis of the timescale dependence of airborne fraction to the set of CMIP5 models listed in Table 1. In particular we study the importance of the different climate–carbon feedbacks for this timescale dependence.
Wu and Xin (2015b)Wu and Xin (2015c)Wu and Xin (2015d)Wu and Xin (2015a)Lindsay (2013b)Lindsay (2013c)Lindsay (2013d)Lindsay (2013a)Dunne et al. (2014a)Dunne et al. (2014b)Dunne et al. (2014c)Dunne et al. (2014d)Liddicoat et al. (2014a)Liddicoat et al. (2014b)Jones et al. (2014)Webb et al. (2014)IPSL (2011)IPSL (2011)IPSL (2011)Caubel et al. (2016)JAMSTEC et al. (2015b)JAMSTEC et al. (2015a)Torres Mendonca et al. (2023)Torres Mendonca et al. (2023)Torres Mendonca et al. (2023)Giorgetta et al. (2012)Tjiputra et al. (2012a)Tjiputra et al. (2012b)Bentsen et al. (2011)Tjiputra et al. (2012c)The whole investigation is based on the calculation of the generalized airborne fraction by means of Eq. (15), whose power to predict airborne fraction from the generalized sensitivities has been demonstrated exemplarily for MPIESM in the previous section. As part of this demonstration, methods to derive the necessary generalized sensitivities from MPIESM standard C^{4}MIP simulations had to be developed (see Appendix A). They are applied here to calculate also for those other CMIP5 models the generalized sensitivities from published 1 % bgc and 1 % rad simulation data.
4.1 Generalized sensitivities of CMIP5 models
In this subsection we present our results for the generalized sensitivities of the considered CMIP5 models. The robustness of the recovered sensitivities depends on the quality of the data (Torres Mendonça et al., 2021a, b) and on how appropriate it is to apply the numerical techniques selected in Appendix A for MPIESM to the other CMIP5 models as well. In principle this robustness should be examined for these other CMIP5 models by means of additional simulations, as was done for the MPIESM in Torres Mendonça et al. (2021b) and in Appendix A. But such simulations are not available. Nevertheless to get an idea of the quality of the recovered generalized sensitivities we use them to predict the time dependence of the standard α, β, and γ sensitivities for published 1 % simulations and compare them with their values obtained directly from the simulation data.
We start by discussing the identified generalized sensitivities. In Fig. 2 we show the ocean sensitivities in the first row and the land sensitivities in the second row. Figure 2a and d show the ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{O}\right)}\left(t\right)$ and ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{L}\right)}\left(t\right)$ sensitivities. The plotted vertical lines indicate the end of that part of the time series that is used to derive the sensitivities according to the techniques summarized in Table A2. As seen, these two sensitivities are for almost all models positive at all times analysed. This is because an increase in atmospheric CO_{2} concentrations results in an increase in land and ocean carbon stocks: for land this positive response is a consequence of the CO_{2} fertilization effect, which increases plant productivity and vegetation growth, while in the oceans it is a consequence of the increase in the difference between atmospheric and oceanic CO_{2} partial pressure, leading to a positive input flux of CO_{2} into the ocean (Arora et al., 2013; Friedlingstein et al., 2006). The surprising negative values of ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{L}\right)}\left(t\right)$ for the HadGEM2ES after 70 years are most likely a consequence of nonlinearities in the response of this model: by the very nature of the biogeochemical response it can be shown that ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{L}\right)}\left(t\right)$ should decrease monotonically to zero (see Torres Mendonça et al., 2021b, Appendix C), so the negative values must be an artefact of the numerical recovery caused either by the nonlinearity of the response being stronger than expected from MPIESM or by deterioration from noise (Torres Mendonça et al., 2021a). Noting that the order of magnitude of the estimated signaltonoise ratio in the HadGEM2ES response is equal to or larger than that in the response of the other models (not shown), the negative values are most likely not caused by noise but are related to nonlinearities. This result suggests that to reliably derive ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{L}\right)}\left(t\right)$ for this model and to minimize the influence of nonlinearities, one should take data until a perturbation strength smaller than that assumed following our investigation with MPIESM (Appendix A).
There is a close agreement between the obtained generalized ocean sensitivities ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{O}\right)}\left(t\right)$ in Fig. 2a: in most models ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{O}\right)}\left(t\right)$ decays rapidly at a similar pace in the first years. Such overall agreement is in contrast to the results for the land sensitivities ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{L}\right)}\left(t\right)$ in Fig. 2d, which spread largely for the different models. In particular the ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{L}\right)}\left(t\right)$ sensitivities behave differently for the NorESM1ME and CESM1BGC models, which have very small values at all times. This behaviour may be explained by noting that these models account for the coupling between the nitrogen and carbon cycles, which reduces the strength of CO_{2} fertilization because of nitrogen limitation (Zaehle et al., 2010; Arora et al., 2013). Excluding NorESM1ME and CESM1BGC, ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{L}\right)}\left(t\right)$ for the other models agrees better at shorter than at longer times.
Figure 2b and e present the results for the ${\mathit{\chi}}_{\mathit{\gamma}}^{\left(\text{O}\right)}\left(t\right)$ and ${\mathit{\chi}}_{\mathit{\gamma}}^{\left(\text{L}\right)}\left(t\right)$ sensitivities. Both generalized sensitivities are negative for all times. This is because globally the land and ocean lose carbon to the atmosphere when only the climatic effect of CO_{2} is taken into account. As clarified by Arora et al. (2013), this loss is explained by noting that rising temperatures over land result in an increase in heterotrophic soil respiration and almost everywhere to a decrease in net primary production (NPP), while over the oceans rising temperatures result primarily in a decrease in CO_{2} solubility and thus degassing.
In contrast to ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{O}\right)}\left(t\right)$, for ${\mathit{\chi}}_{\mathit{\gamma}}^{\left(\text{O}\right)}\left(t\right)$ the model spread is large over the whole time range (compare Fig. 2a and b). Particularly different is the behaviour of the sensitivities for MIROCESM and GFDLESM2M: while for all other models the magnitude of ${\mathit{\chi}}_{\mathit{\gamma}}^{\left(\text{O}\right)}\left(t\right)$ decreases rapidly in the first years, for these two models only a slow decrease resulting from strong contributions of long timescales in the generalized sensitivities (not shown) is observed. In fact the decrease is so slow that it looks as if they had constant values throughout the whole time range, but this is a misimpression induced by the scale of the plot.
The magnitude of ${\mathit{\chi}}_{\mathit{\gamma}}^{\left(\text{L}\right)}\left(t\right)$ is – at least at short times – much larger than that of ${\mathit{\chi}}_{\mathit{\gamma}}^{\left(\text{O}\right)}\left(t\right)$ for almost all models except for NorESM1ME and CESM1BGC. As explained by Arora et al. (2013), in these two models the coupling between the nitrogen and carbon cycles weakens not only the biogeochemical but also the radiative response because temperaturedriven nitrogen remineralization enhances plant productivity, which counteracts the parallel carbon loss from the enhanced soil respiration in the warmer climate (see also Melillo et al., 2002; Thornton et al., 2009). Analogously to ${\mathit{\chi}}_{\mathit{\gamma}}^{\left(\text{O}\right)}\left(t\right)$ in MIROCESM and GFDLESM2M, because of strong contributions from long timescales, the generalized sensitivity ${\mathit{\chi}}_{\mathit{\gamma}}^{\left(\text{L}\right)}\left(t\right)$ decays in the two other models NorESM1ME and CESM1BGC so slowly (see Fig. 2e) that it seems as if it had a constant value throughout the time series, which is also only a misimpression due to the scale. Overall there is a better model agreement for ${\mathit{\chi}}_{\mathit{\gamma}}^{\left(\text{L}\right)}\left(t\right)$ at long rather than at short times.
Finally, Fig. 2c and f present the results for the ${\mathit{\chi}}_{\mathit{\alpha}}^{\left(\text{O}\right)}\left(t\right)$ and ${\mathit{\chi}}_{\mathit{\alpha}}^{\left(\text{L}\right)}\left(t\right)$ sensitivities that characterize the response of ocean and land temperature to atmospheric CO_{2} perturbations. For both ${\mathit{\chi}}_{\mathit{\alpha}}^{\left(\text{O}\right)}\left(t\right)$ and ${\mathit{\chi}}_{\mathit{\alpha}}^{\left(\text{L}\right)}\left(t\right)$ there is a relatively good agreement among models. A larger spread is nevertheless found for values at short times, for which the recovery is less robust due to illposedness of the deconvolution problem that must be solved to recover the generalized sensitivities (Torres Mendonça et al., 2021a). We note also that the values of both sensitivities are closely related: this is expected from the wellknown fact that by various mechanisms the ratio $\mathrm{\Delta}{T}_{\text{L}}/\mathrm{\Delta}{T}_{\text{O}}=:a$ of land to ocean temperature is around 1.4–1.7 (Lee et al., 2021, Sect. 4.5.1.1.1; Eyring et al., 2021, Fig. 3.2b). By their definition in the Laplace domain ${\stackrel{\mathrm{\u0303}}{\mathit{\chi}}}_{\mathit{\alpha}}^{\left(\text{O}\right)}\left(p\right)=\mathrm{\Delta}{\stackrel{\mathrm{\u0303}}{T}}_{\text{O}}\left(p\right)/\mathrm{\Delta}\stackrel{\mathrm{\u0303}}{c}\left(p\right)$ and ${\stackrel{\mathrm{\u0303}}{\mathit{\chi}}}_{\mathit{\alpha}}^{\left(\text{L}\right)}\left(p\right)=\mathrm{\Delta}{\stackrel{\mathrm{\u0303}}{T}}_{\text{L}}\left(p\right)/\mathrm{\Delta}\stackrel{\mathrm{\u0303}}{c}\left(p\right)$, it follows that ${\stackrel{\mathrm{\u0303}}{\mathit{\chi}}}_{\mathit{\alpha}}^{\left(\text{L}\right)}\left(p\right)=a\mathrm{\Delta}{\stackrel{\mathrm{\u0303}}{T}}_{\text{O}}\left(p\right)/\mathrm{\Delta}\stackrel{\mathrm{\u0303}}{c}\left(p\right)=a{\stackrel{\mathrm{\u0303}}{\mathit{\chi}}}_{\mathit{\alpha}}^{\left(\text{O}\right)}\left(p\right)$ so that these two generalized sensitivities should differ by about that factor. Interestingly, this is indeed seen for most models (see inset in Fig. 2f) but only for the first decades – over the last years a large spread arises.
We now turn to the analysis of the plausibility of the recovered generalized sensitivities by means of the prediction of the original α, β, and γ sensitivities for standard C^{4}MIP 1 % simulations. For this analysis we first employ the recovered generalized sensitivities to predict the original α, β, and γ sensitivities and then compare the results to the actual α, β, and γ sensitivities obtained directly from data. That the original α, β, and γ sensitivities can in principle be predicted from the generalized sensitivities may be seen by noting that
where X stands for land (L) or ocean (O), the superscripts “bgc” and “rad” indicate data taken from bgc and rad simulations (see Table A1), and quantities prefixed by Δ stand for a difference with respect to preindustrial times. The first equalities are the definitions of the standard α, β, and γ sensitivities (Friedlingstein et al., 2003). The second equalities are obtained by inserting Eqs. (2)–(5) while accounting for the specific setup of the bgc and rad simulations: to a good approximation (Friedlingstein et al., 2003) temperature remains unchanged in the bgc simulations, while the CO_{2} rise acts only via climate on land and ocean carbon storage in the rad simulations. In the following analysis we compare the prediction from the generalized sensitivities (second equalities) to the actual α, β, and γ values computed by their definition directly from the simulation data (first equalities). The whole comparison is performed for the 1 % rad and the 1 % bgc simulations (Table A1).
The results of these calculations are shown in Fig. 3. We plot the α, β, and γ sensitivities as a function of time for the first 30 years of the 1 % simulations so that CO_{2} forcing strengths are within the estimated linear regime of the generalized framework (around 94 ppm; see Appendix A). Because of natural climate variability there is some uncertainty in the choice of the values of preindustrial temperature and carbon stocks needed to calculate the differences ΔT_{X}(t) and ΔC_{X}(t) in the definitions of the sensitivities (see Eqs. 20–22). In Fig. 3 we used the variability from the associated control simulations to estimate the resulting uncertainty range for the dataderived sensitivities (shaded area in the figures) – details are found in Appendix E. As seen in the figure, for γ(t) the uncertainty range sometimes gets extremely large; this happens when the ΔT_{X}(t) found in the denominator comes close to zero (see Eq. E5). In principle, such an uncertainty from initial values enters also the denominator of the predicted γ sensitivities, but we do not display this to keep the figures simple, and the interpretation of this comparison would not change.
The results for MPIESMLR can be considered a reference for the achievable agreement between dataderived and predicted sensitivities because for MPIESMLR the generalized sensitivities were obtained in a qualitycontrolled way by means of additional simulations (see Appendix A) not available for the other models. In view of this achievable agreement, the predicted α and γ sensitivities excellently match for all models the respective dataderived sensitivities, judged by noting that the predicted sensitivities stay within the amplitude range of interannual variability that is by principle not predictable by the linear response methods employed here because they predict the ensemble mean instead of the individual system development (see the discussion in Torres Mendonça et al., 2021a). In contrast, for all models the predicted β^{(O)}(t) is – for times longer than 15 years – systematically too high. Additionally, the predicted β^{(L)}(t) has a slope and/or offset that systematically differs from those of the respective dataderived β^{(L)}(t).
Such systematic deviations in β^{(O)}(t) and β^{(L)}(t) are largely a result of the effect of nonlinearities in the respective carbon responses $\mathrm{\Delta}{C}_{X}^{\text{bgc}}\left(t\right)$: in fact, if the predicted sensitivities are calculated by using not the untransformed forcing Δc(t) under the integral in Eq. (20) but rather the transformed forcings ${c}_{\text{PI}}\mathrm{ln}\left(c\right(t)/{c}_{\text{PI}})$ and ΔNPP(t) – these are used in the derivation of the generalized sensitivities ${\mathit{\chi}}_{\mathit{\beta}}^{\left(X\right)}$ to account for nonlinearities (see Appendices A2 and A5) – the quality of agreement for β^{(O)}(t) and β^{(L)}(t) considerably improves (not shown). We nevertheless chose to perform the predictions with the untransformed forcing Δc(t) to conform with the standard formulation of the generalized framework (compare Eqs. 2–5). Despite the encountered deviations in β^{(X)}(t), in all cases the magnitude and tendency of the predicted sensitivities match those of the dataderived sensitivities.
Overall, we consider the results of this comparison as sufficiently convincing to add confidence to the validity of the recovered generalized sensitivities (Fig. 2) that underlie the predicted α, β, and γ sensitivities in Fig. 3. We note also that the predicted sensitivities, in contrast to the often noisy dataderived sensitivities, are typically well defined because, as mentioned above, the generalized sensitivities predict the response not in noisy individual realizations but in a smooth ensemble mean. Nevertheless, it should be kept in mind that this comparison is a mere plausibility check because essentially the same data used to predict the α, β, and γ sensitivities were also used to derive the generalized sensitivities.
4.2 Additivity of responses
Before the main question of this study on the role of feedbacks for airborne fraction can finally be addressed in the next section, another preparatory step is necessary. Key to investigate this question will be Eq. (15) from the generalized framework to predict the generalized airborne fraction via the feedback functions from the generalized sensitivities as already explained in Sect. 3.2. In applying this relation to the data from the different CMIP5 models one must realize that the accuracy of such predictions depends on two aspects: (1) the quality of the numerical recovery of the generalized sensitivities (see previous subsection) and (2) the validity of the assumption underlying the generalized framework that for weak perturbations the carbon response to CO_{2} is determined by the sum of the biogeochemical and radiative responses; this assumption of additivity is implicit to Eqs. (2)–(3), where the first term represents the biogeochemical response and the second (after insertion of Eqs. 4–5) the radiative response. Ideally, for each model one should fully check these two aspects with the aid of additional simulations, as we did for the MPIESM in Torres Mendonça et al. (2021b) and in Appendix A. Unfortunately, such additional simulations are at present not available for other models, so a full check is not possible. Nevertheless, since all CMIP5 models provide in addition to the 1 % rad and bgc simulations also a 1 % fully coupled simulation (a 1 % simulation where both the biogeochemical and the radiative effects of CO_{2} are active), one can at least check whether the biogeochemical and radiative responses are indeed additive for a certain range of perturbation strengths. The rationale underlying this check is that the 1 % rad and bgc simulations separately give the values for the two righthandside terms in Eqs. (2)–(3), while the 1 % fully coupled simulation gives the lefthand sides.
To check additivity, we plot in Fig. 4 for each of these models the response in carbon storage for land, ocean, and global (land plus ocean) carbon from the 1 % fully coupled experiment along with the sum of the responses from the 1 % bgc and 1 % rad experiments. If additivity holds for a certain range of perturbation strengths, then within this range these two curves (1 % fully coupled and the sum of 1 % bgc and 1 % rad) must agree. As seen, for all models there is indeed agreement at least within the estimated range of linearity (94 ppm CO_{2} rise; see Appendix A) for land, ocean, and global carbon stock changes, with larger discrepancies for land and global carbon in MIROCESM and NorESM1ME. For those two models, the generalized framework may not fully describe the linear dynamics of the system in the fully coupled setup. But overall, within the linear regime one may say that the biogeochemical and radiative responses are approximately additive in the CMIP5 models. This result, together with the evidence for the overall plausibility of the generalized sensitivities obtained in the last subsection, gives some confidence that our numerical methods and the framework as a whole are describing in reasonable approximation the linear dynamics of the carbon cycle in these models.
4.3 Climate–carboncycle feedbacks and airborne fraction
In this section we tackle the main question of our study, namely how the climate–carboncycle feedbacks shape the timescale dependence of the generalized airborne fraction. From here on we take for granted that by the methods presented in the previous sections Eq. (15) indeed reliably predicts the generalized airborne fraction for the considered CMIP5 models. We thus proceed to estimate the feedbacks and the airborne fraction for each model via Eqs. (9), (10), and (15). The results are shown in Fig. 5.
In Fig. 5a, one sees that for almost all CMIP5 models the timescaledependent airborne fraction decreases as the timescale $\mathrm{1}/p$ increases, all starting at 1 for short timescales and spreading from 0.56 to 0.75 at a timescale of 10 years and from 0.26 to 0.5 at a timescale of 100 years. The only exception is HadGEM2ES, whose timescaledependent airborne fraction once more increases at long timescales, which, as will be seen below, is related to a reduction in the magnitude of its land biogeochemical feedback. That the airborne fraction has a value of 1 at short timescales is not a property of the models but follows from its definition (compare Eq. 19). But this is not so for the decrease in airborne fraction at long timescales: as will become clear below this behaviour is a consequence of the biogeochemical feedback being stronger than the radiative feedback in the considered CMIP5 models, so this decrease is presumably a genuine property of the Earth system.
How the airborne fraction changes in the timescale domain is determined by the climate–carboncycle feedbacks. As seen in Fig. 5b, for $\mathrm{1}/p\to \mathrm{0}$ all feedback functions approach zero, which in view of Eq. (15) is consistent with $\stackrel{\mathrm{\u0303}}{A}\left(p\right)\to \mathrm{1}$ for the airborne fraction. Besides the mathematical reasons explained in connection with Eq. (19), this behaviour can intuitively be understood by noting that at short timescales the ocean and land carbon cycles have not yet reacted to emissions, so at these scales $\stackrel{\mathrm{\u0303}}{A}\left(p\right)\to \mathrm{1}$ and thus atmospheric CO_{2} simply follows emissions (Eq. 14). To understand how the airborne fraction behaves as the timescale increases, one must look at the behaviour of the separate feedback functions. For longer timescales, the feedback functions ${\stackrel{\mathrm{\u0303}}{f}}_{\mathit{\gamma}\mathit{\alpha}}^{\left(\text{L}\right)}$ and ${\stackrel{\mathrm{\u0303}}{f}}_{\mathit{\gamma}\mathit{\alpha}}^{\left(\text{O}\right)}$ that quantify the radiative feedbacks become increasingly positive, while the feedback functions ${\stackrel{\mathrm{\u0303}}{f}}_{\mathit{\beta}}^{\left(\text{L}\right)}$ and ${\stackrel{\mathrm{\u0303}}{f}}_{\mathit{\beta}}^{\left(\text{O}\right)}$ that quantify the biogeochemical feedbacks become in general increasingly negative. The sign of these feedbacks is in agreement with current process understanding (Friedlingstein et al., 2006; Gregory et al., 2009; Arora et al., 2013, 2020). But that these feedbacks are either positive or negative for all timescales is a nontrivial result that could not be obtained by the original α–β–γ framework because there the timescaledependent feedback strengths show up only combined with the external forcing that enters the quantification by the underlying α, β, and γ sensitivities (see discussion in Sect. 2). The observed uniformity of the sign of the feedbacks is explained by the fact that almost all generalized sensitivities in Fig. 2 are either always negative or always positive: since the feedback functions are proportional either to the Laplace transform of χ_{β} or to the product of the Laplace transforms of χ_{γ} and χ_{α} (Eqs. 9 and 10), one sees following the lines of the argument for $\stackrel{\mathrm{\u0303}}{A}\left(p\right)\ge \mathrm{0}$ in Appendix D that by the positivity or negativity of these response functions in time also their Laplace transforms and the associated feedback functions must have a unique sign at all timescales. Interestingly, for increasing timescales in almost all models – except for the HadGEM2ES at long timescales – the sum (land plus ocean) of the biogeochemical feedbacks becomes increasingly larger than the sum of the radiative feedbacks. As a result, in all these models the total feedback function $\stackrel{\mathrm{\u0303}}{f}\left(p\right)$ gets increasingly negative (not shown) and by Eq. (15) the airborne fraction always decreases. In the HadGEM2ES, at long timescales the magnitude of the land biogeochemical feedback starts to decrease, thereby reducing the magnitude of the (negative) total feedback function and as a consequence increasing the airborne fraction.
In the mean over all models, the land biogeochemical feedback is at all timescales longer than the ocean biogeochemical feedback: at a timescale of 10 years, it is 1.4 times larger, and at a timescale of 100 years, it is 1.8 times larger. The picture is qualitatively similar for the radiative feedback: at a timescale of 10 years, the land feedback is, despite its small value of 0.03, orders of magnitude larger than the almost absent ocean feedback, and at a timescale of 100 years, the land feedback is 7.4 times larger than its ocean counterpart. Aggregating land and ocean, the mean of the biogeochemical feedback is 22 times larger than the radiative feedback at a timescale of 10 years and 5.6 times larger at a timescale of 100 years. These results are in particular at short timescales in contrast to previous estimates (Gregory et al., 2009; Arora et al., 2013) using Friedlingstein's framework, which suggested that the biogeochemical feedback is about 4 times larger than the radiative feedback. One must note though that our and previous estimates are not entirely comparable: while previous estimates were made for a particular scenario, our estimate is valid for any scenario. In addition, here only the linear regime is considered, so the saturation of the land and ocean carbon sinks (which reduces the values of β^{(L)} and β^{(O)} when higher perturbation strengths are considered) is not taken into account.
By Fig. 5b the model spread is for the land feedbacks much larger than that for the ocean feedbacks, as expected from previous studies (Gregory et al., 2009; Arora et al., 2013; Friedlingstein et al., 2014). Because of the nonlinear relationship between $\stackrel{\mathrm{\u0303}}{A}\left(p\right)$ and the feedback functions (see Eq. 15), it is not immediately clear how the model spread in the feedbacks propagates to the airborne fraction. Assuming small, independent spreads, this propagation may be computed as
where σ^{2} denotes the spread (variance) for each quantity. This follows by expanding $\stackrel{\mathrm{\u0303}}{A}\left(p\right)$ around the mean into the feedback functions, assuming small spreads in the functions so that only linear terms are kept, and then calculating the variance of the result assuming independent spreads so that covariance terms are ignored (Roe, 2009; Barlow, 1989, p. 55). Figure 6a shows the terms on the righthand side of Eq. (23), the resulting approximation of ${\mathit{\sigma}}_{\stackrel{\mathrm{\u0303}}{A}}^{\mathrm{2}}\left(p\right)$ (sum of those terms), and the true spread in the airborne fraction. As seen, the true variance in the airborne fraction follows closely the component arising from the land biogeochemical feedback, with a slightly larger discrepancy at timescales above 30 years, indicating that on longer timescales other feedbacks' spread becomes relatively more important. This indicates that most of the model spread in the airborne fraction arises from the spread in the land biogeochemical feedback. This result agrees with that obtained in a recent study (Jones and Friedlingstein, 2020) that also evaluated how the model spread in the airborne fraction is affected by the spread in the climate–carboncycle feedbacks but employing Friedlingstein's original framework for the analysis.
An even clearer view about the impact of the different feedbacks on the airborne fraction may be gained by artificially changing the values of these feedbacks to study hypothetical situations and then evaluating the resulting change in the airborne fraction. For instance, one can illustrate how strongly the model spread in the airborne fraction depends on the spread in the land biogeochemical feedback by recalculating the statistics of the airborne fraction taking ${\stackrel{\mathrm{\u0303}}{f}}_{\mathit{\beta}}^{\left(\text{L}\right)}$ for all models equal to its model mean. As shown in Fig. 6b, it turns out that if the exact values of only this feedback function were known and equal to the model mean (with all other feedback spreads remaining the same), the spread in the airborne fraction would reduce by about 82 % at a timescale of 10 years and 61 % at a timescale of 100 years. This once more makes obvious that the main reason for the model spread in the airborne fraction is the large spread in the land biogeochemical feedback.
The dynamics of the global carbon cycle can be understood in terms of feedbacks arising via the land and ocean carbon cycle when atmospheric CO_{2} is perturbed. To disentangle and separately quantify those feedbacks, Friedlingstein et al. (2003) developed the α–β–γ framework. Although this framework gives insight into the main effects of atmospheric CO_{2} perturbations in the global carbon cycle, by not accounting for the internal timescales of the system, the resulting quantification of feedbacks is valid only for a particular perturbation scenario and time period. Such limitations are overcome by employing the recently proposed generalization of this framework (Heimann, 2014; Rubino et al., 2016; Enting and Clisby, 2019; Enting, 2022). By assuming weak perturbations and accounting for the memory of the carbon cycle, the generalized α–β–γ framework quantifies feedbacks independently of the perturbation scenario at different timescales. As a result of the generalization, this framework gives in principle a complete description of the linear dynamics of the global carbon cycle in terms of the underlying feedbacks. But so far its applicability to study the dynamics of the climate–carbon system had not been systematically investigated.
Here, we employed this generalized framework to study the timescale dependence of the climate–carboncycle feedbacks and the associated airborne fraction for an ensemble of CMIP5 models. In Sect. 3, we have shown for the example of MPIESM that the generalized sensitivities identified from concentrationdriven simulations correctly predict via Eqs. (7)–(8) and (15) the generalized airborne fraction derived from emissiondriven simulations. This demonstrates that the generalized α–β–γ framework has indeed predictive power.
Based on experience with MPIESM, we quantified in Sect. 4 the timescaledependent airborne fraction and feedbacks for various other CMIP5 models. As can be seen from Eq. (14), the timescaledependent airborne fraction quantifies the fraction of emissions staying in the atmosphere at a particular timescale. At short timescales this fraction is known to be 1; i.e. atmospheric carbon simply follows emissions because feedbacks that could change it need sufficient time to react. We found that for almost all models, the airborne fraction strictly decreases towards long timescales. This decrease is slow: even at a timescale of 100 years the airborne fraction has dropped down only to values ranging from 0.26 to 0.5, meaning that even a century after CO_{2} emissions happened a considerable amount is still airborne, which is consistent with results from impulse experiments (Archer et al., 2009).
Considering global carbon, in the model mean the biogeochemical feedback was found to be 22 times larger than the radiative feedback at a 10year timescale and 5.6 times larger at a 100year timescale. This result suggests that at least over shorter timescales the difference between these feedbacks may be even greater than previously thought (Gregory et al., 2009; Arora et al., 2013). Nevertheless, one must note that here only the linear perturbation regime is considered, so a possible saturation of the land and ocean carbon sinks at high CO_{2} (that would reduce the values of β^{(L)} and β^{(O)} compared to the case where no saturation is present) is not reached.
The influence of the model spread of the different feedbacks on the airborne fraction was also investigated. It was found that the spread in the airborne fraction arises mostly from the spread in the land biogeochemical feedback, especially for timescales below 30 years. By considering the hypothetical case where this particular feedback would be equal to the model mean, we found that the spread in the airborne fraction would decrease by 82 % at a 10year timescale and by 61 % at a 100year timescale, which demonstrates even more clearly that it is indeed the land biogeochemical feedback that causes the spread in airborne fraction between the different models.
While the generalized framework was shown here to reasonably describe the linear dynamics of the global carbon cycle in MPIESM, the results obtained for the other CMIP5 models depend on two basic assumptions. The first is that the generalized sensitivities in the CMIP5 models are recoverable with sufficient quality by the same numerical approaches that were appropriate to recover the sensitivities in MPIESM. This involves the assumption that for all other considered CMIP5 models the linear perturbation regime is of similar extent to that found for MPIESM for the different response variables invoked to recover the sensitivities. This might not be the case – and is probably not for ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{L}\right)}$ in the HadGEM2ES (see discussion of Fig. 2d). The second basic assumption is that the generalized framework itself correctly describes the dynamics of the global carbon cycle in those models. This involves the assumption that the biogeochemical and radiative responses are additive – which was confirmed to a good degree of approximation within the linear regime (see Fig. 4) – but also that the whole carbon cycle dynamics can be described in terms of the responses to atmospheric CO_{2} and temperature alone (see Eqs. 2–5). Ideally, each of these assumptions should be carefully investigated for each model separately by aid of additional simulations, similarly to what was done here for the case of MPIESM.
With these cautionary remarks in mind, our conclusion that the spread in the airborne fraction arises mostly from the spread in the land biogeochemical feedback corroborates the recent finding by Jones and Friedlingstein (2020), who performed a similar analysis employing Friedlingstein's original α–β–γ framework. This agreement adds confidence to the results obtained here and suggests that research on climate–carboncycle feedbacks should focus especially on understanding the land biogeochemical feedback, since it is the largest source of model uncertainty in the airborne fraction in our and in their analysis. But additionally, the conclusion drawn here is valid not only for the particular perturbation scenario underlying the CMIP5 protocol but for all scenarios within the linear regime. To this extent, this study at least partially answers the question raised by Jones and Friedlingstein (2020) about the behaviour of climate–carboncycle feedbacks in scenarios with different characteristics: as long as we stay within the linear regime, there is no need for calculating separate feedback metrics for the different scenarios; the feedback functions of the generalized framework describe this behaviour for all scenarios at once.
Estimates of the timescaleresolved airborne fraction, by means of Eq. (18), and of timescaledependent generalized sensitivities had already been attempted (Enting, 2007, 2022; Rubino et al., 2016; Enting and Clisby, 2019). These attempts, focused on the observed carbon cycle, were based on ad hoc assumptions on the internal memory in terms of the analytical structure of the underlying response functions and values of internal timescales. Such assumptions could be circumvented in the present study by employing our RFI method (Torres Mendonça et al., 2021a, b) that instead derives the internal timescale spectra of the system by fully accounting for the illposedness of the underlying inverse problem. But it should be noted that our approach is tailored to simulation data, and whether it may be applicable to observation data needs to be seen. One of the involved challenges for such an application is that our setting here is limited to an idealized case where CO_{2} is the only forcing, while in observations one would have to account for further perturbations such as nonCO_{2} greenhouse gases, land use change, and aerosols.
As explained in Sect. 2, the timescaleresolved airborne fraction can be understood as a generalization of the airborne fraction in its standard definition (defined as a ratio of atmospheric CO_{2} fluxes to emission fluxes; see Eq. 12). It is well known that the nearconstancy of the value of the standard airborne fraction is a result of the linearity of response in combination with the exponential character of the increase in historic emissions (Raupach, 2013). Accordingly, once emissions get nonexponential – as it must be if future emissions are significantly reduced – the standard airborne fraction must deviate from its constant value (Raupach, 2013). Thus, standard airborne fraction cannot be seen as an invariant property of the climate–carbon system. In contrast, the generalized airborne fraction investigated here describes the response of atmospheric CO_{2} to any emissions scenario and is therefore indeed an invariant property of the system. As such, the generalized airborne fraction may be employed to predict the standard airborne fraction by Eq. (13) for any emissions scenario as long as emissions are sufficiently weak (see Appendix F for an exemplary demonstration). This is the case, for example, for emission scenarios with low or even negative future emissions, as investigated in Jones et al. (2016). For such scenarios the generalized airborne fraction – once it has been determined for one or more models – could be used as an emulator to see without expensive Earth system simulations how atmospheric CO_{2} develops; and a similar approach could be applied to changes in the land and ocean carbon reservoirs by determining the appropriate response functions.
Besides investigating the timescale dependence of airborne fraction, our study also demonstrated for MPIESM the predictive power of the generalized framework (see Sect. 3). This demonstration indicates that this framework may be invoked also to tackle other problems involving timedependent components of the climate–carbon system. Directly related to the airborne fraction is the transient climate response to cumulative emissions (TCRE), which quantifies the change in global mean temperature in response to cumulative emissions (Ciais et al., 2013). As recognized by Gregory et al. (2009) and Jones and Friedlingstein (2020), the TCRE can be investigated using the standard climate–carbon feedback framework by means of airborne fraction and the global temperature sensitivity to CO_{2}. We believe this extends to the generalized framework investigated here – note that although we take separate temperature sensitivities for land and ocean, sensitivities defined for a single global temperature can be easily obtained from them (see Appendix G). Jones and Friedlingstein (2020) in particular quantified the contribution of the spread in each sensitivity in Friedlingstein's framework to the spread in TCRE, finding in their multimodel ensemble that the spread in the β^{(L)} sensitivity has the second largest contribution, smaller only than that from the spread in the α sensitivity (see their Fig. 4). It would be interesting to check this finding in the generalized framework employed here, where the additional complications arising from the scenario dependence of Friedlingstein's framework are absent. Further, the generalized framework could be used to study how the TCRE is determined by contributions from the different climate–carboncycle feedbacks at different timescales. Such analysis could lead to a better understanding of the dynamics behind this metric.
Furthermore, the generalized framework may be invoked to investigate the contribution of the different feedbacks to committed climate change, where one is interested in understanding the behaviour of the system once atmospheric CO_{2} stabilizes or emissions cease (Wetherald et al., 2001; Meehl et al., 2005; Wigley, 2005; Plattner et al., 2008; Mauritsen and Pincus, 2017; MacDougall et al., 2020). Also, since the generalized framework gives a consistent formalism for quantifying climate–carboncycle feedbacks, it may be possible to apply it to study climate–carboncycle feedbacks and physical climate feedbacks in a consistent unified framework (Gregory et al., 2009; Williams et al., 2019; Goodwin et al., 2019).
One aspect emphasized throughout this study is that the generalized framework is valid only for weak perturbations. In fact, we have found in application to the MPIESM that the linear regime extends only up to about 100 ppm atmospheric CO_{2} increase beyond the preindustrial level, i.e. up to a value of about 380 ppm. This perturbation level was reached already around 2005 (Dlugokencky and Tans, 2023), so this linear framework cannot be employed to study future climate change. But gaining understanding of the largescale dynamics and the underlying memory structure of the coupled climate–carboncycle system should be easier in the linear regime where complications from nonlinearities that are expected to get increasingly important during future climate change are absent. Given the promising results obtained here and the potential applications outlined above, we believe that the generalized framework is thus the right tool for such investigations. But in principle the Volterra expansion underlying the generalized framework can be extended beyond the linear term (Ruelle, 1998; Lucarini, 2009), so one could also think of a nonlinear generalized feedback formalism applicable to nearfuture climate change (Roe, 2009). For such a research program it would be useful to have simulations with a better signaltonoise ratio to recover the response functions. In the present study we used published 1 % simulations from C^{4}MIP, but as we showed in Torres Mendonça et al. (2021a, b) simulations forced by a stepfunction CO_{2} rise would be more suitable for their recovery. Therefore it could be an idea for a future C^{4}MIP protocol to switch to such simulations.
Finally it may be noted that our study is an example for the application of linear response theory – known from statistical mechanics (Kubo et al., 2012) – to a dissipative system (Lucarini et al., 2014), namely to the global carbon cycle. Recently it has been noted that the RFI method underlying our recovery of the generalized sensitivities is limited in scope because it assumes the absence of an imaginary part of the eigenvalues of the evolution operator (Santos Gutiérrez and Lucarini, 2022). In the case that some eigenvalues have a nonzero imaginary part, the system contains internal oscillatory modes. But this is hardly believable to be true for the carbon cycle at the timescales from 1 year up to a century considered in the present study. While the intraannual oscillations seen, for example, in atmospheric CO_{2} are caused externally by the effect of the seasonal changes in insolation on climate and photosynthesis, ecological communities in the ocean and on land may in principle be capable of oscillatory behaviour in their population dynamics at multiannual scales (see e.g. Murray, 1993). But there is no evidence that such processes may be relevant at global scale. In contrast to the carbon cycle, the climate system is known to have internal oscillatory modes (e.g. ENSO and North Atlantic Oscillation). Currently it is unclear how one could account for nonzero imaginary parts of the eigenvalues in the recovery of linear response functions from data; the RFI method gains part of its simplicity from the particular assumption of noncomplex eigenvalues. To this extent, the carbon cycle seems to be that part of the Earth system to which linear response theory may be applied most easily.
This appendix complements the results from Torres Mendonça et al. (2021b) to derive for the MPIESM all generalized sensitivities. The results from this appendix are needed for (i) testing the predictive power of the generalized α–β–γ framework in Sect. 3, (ii) identifying the best data preprocessing techniques to optimally recover the generalized sensitivities for the investigation of the feedbacks and of airborne fraction in CMIP5 models in Sect. 4, and (iii) obtaining an estimate of the linear regime for which the generalized α–β–γ framework is valid.
Since the land carbon sensitivities ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{L}\right)}$ and ${\mathit{\chi}}_{\mathit{\gamma}}^{\left(\text{L}\right)}$ were already derived in Torres Mendonça et al. (2021b), here we derive the remaining generalized sensitivities for ocean carbon and land and ocean temperature ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{O}\right)}$, ${\mathit{\chi}}_{\mathit{\gamma}}^{\left(\text{O}\right)}$, ${\mathit{\chi}}_{\mathit{\alpha}}^{\left(\text{L}\right)}$, and ${\mathit{\chi}}_{\mathit{\alpha}}^{\left(\text{O}\right)}$ (see Sect. 2). We recover these sensitivities by the RFI method (Torres Mendonça et al., 2021a) using data taken from standard C^{4}MIP 1 % experiments (see Table A1). To obtain the generalized sensitivities with the best possible quality and also estimate the linear regime for which the generalized α–β–γ framework is valid, following Torres Mendonça et al. (2021b) we proceed in three steps:

Select a technique to pretransform the data to account for known nonlinearities in the response. Accounting for these nonlinearities allows for recovering the generalized sensitivity from experiments with higher perturbation strengths and thus higher signaltonoise ratio, which improves the quality of the results.

Determine the maximum forcing strength for which no strong nonlinearities are present in the response. This gives the best tradeoff between signaltonoise ratio and nonlinearity for a particular pretransformed response data, thereby further improving the quality of the recovery.

Calculate the linear regime of the response, i.e. the range of forcing strengths for which the generalized sensitivity can be used to predict the response of the system. By analysing this linear regime for all generalized sensitivities we will be able to determine the overall linear regime for which the generalized α–β–γ framework as a whole can predict the dynamics of the coupled global carbon cycle accounting for all climate–carboncycle feedbacks.
The final result of this analysis is summarized in Table A2 in Sect. A6. The results for the best pretransformation technique and maximum forcing strength for the recovery of the MPIESM sensitivities are used to derive the generalized sensitivities for all CMIP5 models in Sect. 4.
A1 Procedure to analyse the recovery of the generalized sensitivities
To perform the analysis described in the three steps above we employ a simple procedure introduced in Torres Mendonça et al. (2021b). The idea behind the procedure can be understood as follows. First, using the RFI method (Torres Mendonça et al., 2021a), we recover the generalized sensitivity taking pretransformed data from increasingly longer time periods of the 1 % experiment. For each time period, we employ the recovered generalized sensitivity to predict the response of additional 0.5 % and 0.75 % experiments covering that same time period (for a description of the experiments see Table A1). Then, we measure the quality of the recovery of the generalized sensitivity by the quality of the prediction of the responses for these additional simulations. The quality of the prediction is quantified by the (dimensionless) prediction error:
where χ⋆Δf^{k} gives the predicted values, ⋆ stands in short for the convolution operation, ΔY^{k} and Δf^{k} are the response and the perturbation in experiment “k”, and χ is the response function recovered from the 1 % experiment. Because in the 1 % experiment the perturbation strength increases with time, also the signaltonoise ratio of the data increases. On the other hand, higher perturbation strengths increase nonlinearities. This results in the following tradeoff: while a higher signaltonoise ratio results in a recovery with better quality, larger nonlinearities deteriorate the quality of the recovery. From this tradeoff, by analysing the prediction error (A1) for different perturbation strengths of the 1 % experiment and also for different data pretransformation techniques, one can (1) select the best pretransformation technique for the recovery of the response function, (2) determine the maximum forcing strength for which the response function can be optimally recovered, and (3) estimate the linear regime of the response. For a more detailed description of this procedure please refer to Torres Mendonça et al. (2021b).
A2 Generalized sensitivity ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{O}\right)}$
Pretransformation techniques to recover ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{O}\right)}$
Similarly to Torres Mendonça et al. (2021b), we recover ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{O}\right)}\left(t\right)$ employing the RFI method in combination with two techniques. The first technique consists of simply taking Δc(t) as perturbation and deriving ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{O}\right)}\left(t\right)$ from
where $\mathrm{\Delta}{C}_{\mathrm{O}}^{\mathrm{bgc}}\left(t\right)$ is the ocean carbon response in the bgc setup (see Table A1). Because we directly take Δc(t) for the recovery, we call this the notransform technique.
Since Eq. (A2) is the equation employed in the generalized framework to describe the ocean biogeochemical response (compare Eq. 3), from this notransform technique we will also obtain an estimate of the range of CO_{2} perturbation strengths for which this response can be considered linear. This estimate is needed to address the question of what the linear regime is for which the generalized α–β–γ framework is valid in the MPIESM.
In the second technique, we consider the logarithm of c as perturbation and derive ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{O}\right)}\left(t\right)$ from
where c_{PI} is the preindustrial value for atmospheric CO_{2}. Because we take instead of c its logarithm, we call this the logtransform technique.
This logtransform technique is inspired by the fact that nonlinearities in the ocean carbon uptake come mainly from the dissolution of CO_{2} in ocean surface water (see e.g. Joos and Bruno, 1996). Because with accumulation of CO_{2} in upper layers of the ocean the ability for further uptake of CO_{2} decreases (Hooß et al., 2001), we use a logarithmic representation for the perturbation to try to explicitly describe the nonlinearity between CO_{2} concentration and the carbon flux into the ocean. As explained in Torres Mendonça et al. (2021b, Eqs. (16), (18), and (19) in Sect. 4.1), employing Eq. (A3) has the advantage that ${\mathit{\chi}}_{\mathit{\beta},\mathrm{ln}}^{\left(\text{O}\right)}\left(t\right)={\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{O}\right)}\left(t\right)$; i.e. by deriving ${\mathit{\chi}}_{\mathit{\beta},\mathrm{ln}}^{\left(\text{O}\right)}\left(t\right)$ from Eq. (A3) one obtains also the desired ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{O}\right)}\left(t\right)$.
For both techniques ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{O}\right)}\left(t\right)$ is derived enforcing monotonicity by the RFI algorithm (see Fig. 1 in Torres Mendonça et al., 2021a).
Recovery of ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{O}\right)}$ and linear regime of the biogeochemical response of ocean carbon
Using the pretransformation techniques described above, in the following we recover the generalized sensitivity ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{O}\right)}$ and evaluate the quality of the results.
We start by deriving ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{O}\right)}\left(t\right)$ from Eq. (A2) (notransform technique) taking data from the 1 % bgc experiment. Following the procedure described in Sect. A1, we employ Eq. (A2) to predict the response from the 0.5 % and 0.75 % bgc experiments (see Table A1). Figure A1a shows the resulting prediction error (see Eq. A1) for the increasing time period of the data used to obtain ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{O}\right)}\left(t\right)$ – note that on the x axis the CO_{2} forcing strength in the 1 % bgc experiment at the end of the period is used. Similarly to the biogeochemical response of land carbon in Torres Mendonça et al. (2021b), minima are found for final forcing strengths below 100 ppm (about 94 ppm for the 0.75 % and 58 ppm for the 0.5 % bgc experiments), indicating the presence of strong nonlinearities for forcing strengths above this value. Because both minima are “flat”, i.e. the error does not change substantially for values around the minima, we take as an estimate of the linear regime forcing the strengths below the highest minimum – about 94 ppm in the 0.75 % curve.
To try to improve the quality of the recovery, ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{O}\right)}\left(t\right)$ was derived as well using the logtransform technique (Eq. A3). To see whether the quality of the recovery indeed improves, one must check if Eq. (A3) indeed accounts for some of the nonlinearities in the response. If this is the case, then Eq. (A3) should predict the 0.5 % and 0.75 % bgc responses better than Eq. (A2). To check this, we recovered ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{O}\right)}\left(t\right)$ by this logtransform technique and then employed the recovered ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{O}\right)}\left(t\right)$ in Eq. (A3) to predict these responses. The resulting prediction error is shown in Fig. A1b. Overall the error is substantially smaller than in Fig. A1a. Minima are still found but now at values between 100 and 200 ppm and with smaller optimal errors. In contrast to Fig. A1a, after the minima the error increases only slightly for increasing final forcing strength. Therefore, Eq. (A3) seems to indeed account for some of the nonlinearities in the response. Hence, with this approach one can derive ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{O}\right)}\left(t\right)$ taking data from the 1 % bgc experiment until larger perturbation strengths and therefore higher signaltonoise ratios than with the notransform technique, making it in principle possible to recover ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{O}\right)}\left(t\right)$ with better quality.
But despite the overall reduction in the prediction error, Fig. A1b still shows a slightly increasing trend in the prediction error after the minima, indicating the presence of possibly nonnegligible nonlinearities in the response. Therefore we check whether the logtransform technique indeed gives a better recovery for ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{O}\right)}\left(t\right)$ by comparing the recovery from this technique with that from the notransform technique. Figure A1c shows the response function recovered with each technique taking data from the 1 % bgc experiment until the optimal final forcing strength (we take 94 ppm or 30 years of the 1 % bgc experiment for the notransform technique and 138 ppm or 50 years of the 1 % bgc experiment for the logtransform technique; the values correspond to the minima found in the 0.75 % curve, respectively, in Figs. A1a and b) and the response function recovered with the logtransform technique but taking data until the maximal final forcing strength in the 1 % bgc experiment, i.e. the whole time series (since for those forcing strengths the error in Fig. A1b is still small). As seen, all recovered response functions are very similar. This similarity is most likely due to the combination of an overall high signaltonoise ratio of the ocean biogeochemical response (as indicated by the small errors even at small perturbation strengths in panels a and b) and only small contributions from nonlinearities when employing the logtransform technique. This suggests that for the MPIESM ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{O}\right)}\left(t\right)$ can be equally well recovered by either of the two techniques. Despite this result, we selected the logtransform technique for deriving ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{O}\right)}\left(t\right)$ in Sects. 3 and 4 because this technique accounts for nonlinearities that may cause a more significant error in the recovery for other models. But because panel (b) indicates that nonnegligible nonlinearities may still remain in the response even when employing this logtransform (because of the slightly increasing trend), we choose conservatively to take data only until the optimal forcing strength (138 ppm or 50 years of the 1 % bgc experiment) and not the full time series. As suggested by the small differences between the recoveries for the optimal and maximal forcing strengths in panel (c), taking data only until this optimal forcing strength does not significantly hinder the recovery of timescales on the order of the length of the full time series.
To obtain evidence that the recovered ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{O}\right)}\left(t\right)$ indeed characterizes the biogeochemical response of the ocean carbon of the MPIESM to any temporal development of sufficiently weak atmospheric CO_{2} perturbations, we show how well it predicts the model response for additional CO_{2} perturbation experiments (see Table A1). The prediction is performed by convoluting the recovery of ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{O}\right)}\left(t\right)$ selected above with the different CO_{2} perturbation scenarios first via Eq. (A2), which is the equation employed in the generalized framework (compare with Eq. 3 in Sect. 2), and then via Eq. (A3) to see if by this equation the prediction is indeed further improved for all experiments. The result for the first case is shown in Fig. A1d. To better visualize the linear regime region where the approximation works, we show only the first 30 years. The responses for the 1.5 % and 2 % bgc experiments are marked with dots where the forcing strength exceeds 94 ppm, the maximal forcing strength corresponding to the linear regime estimated from Fig. A1a. As seen, the recovered ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{O}\right)}\left(t\right)$ predicts almost perfectly the response from experiments whose forcing strengths are within the linear regime estimate, i.e. the 0.5 %, 0.75 %, and 1.1× CO_{2} bgc experiments. It also predicts with reasonable quality of agreement the response from the 1.5 % and 2 % bgc experiments for forcing strengths within the estimated linear regime. As forcing strengths exceed the linear regime, the prediction starts to deviate from the model response. Likewise, the prediction fails for the whole 2× CO_{2} bgc response because its forcing strength is larger than the linear regime estimate throughout the whole experiment. These results therefore suggest that ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{O}\right)}\left(t\right)$ indeed characterizes the biogeochemical response of ocean carbon for any temporal development of atmospheric CO_{2}, as long as its strength is within the linear regime.
Figure A1e shows the prediction employing Eq. (A3). As seen, now the predictive power of ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{O}\right)}\left(t\right)$ improves considerably, extending to the whole time series for all experiments. Some discrepancies are nevertheless encountered, especially for the last years of the 2× CO_{2} bgc response. Such discrepancies are probably related to the relatively limited amount of information that data from 1 % experiments provide to recover response functions, as observed for the land carbon in Torres Mendonça et al. (2021b).
Despite these discrepancies, overall Fig. A1d and e show that the recovered ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{O}\right)}\left(t\right)$ predicts the response from different experiments with reasonable quality up to certain perturbation strengths. Equation (A3) – employed for the prediction in Fig. A1e – demonstrates successful prediction of the response of the model up to higher perturbation strengths than Eq. (A2). Therefore, Eq. (A3) is probably more appropriate when the aim is to predict the separate biogeochemical response of the model to prescribed atmospheric CO_{2} perturbations. On the other hand, Eq. (A2) is the equation employed in the generalized α–β–γ framework (first term on the righthand side of Eq. 3). Therefore, when estimating the linear regime for the application of the generalized α–β–γ framework one has to consider for the ocean biogeochemical response the estimate obtained from Eq. (A2), i.e. the linear regime estimate indicated in panel (a).
The conclusions from this subsection therefore suggest that the best pretransformation technique to derive ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{O}\right)}\left(t\right)$ is the logtransform technique (Eq. A3), taking data until 138 ppm or 50 years of the 1 % bgc experiment, which corresponds to the first minimum for the prediction error in Fig. A1b. In addition, the linear regime for the ocean biogeochemical response in the generalized α–β–γ framework (first term on the righthand side of Eq. 3) is about 94 ppm, as estimated from Fig. A1a.
A3 Generalized sensitivity ${\mathit{\chi}}_{\mathit{\gamma}}^{\left(\text{O}\right)}$
Technique to recover ${\mathit{\chi}}_{\mathit{\gamma}}^{\left(\text{O}\right)}$
In this subsection we recover ${\mathit{\chi}}_{\mathit{\gamma}}^{\left(\text{O}\right)}\left(t\right)$. Since, as will be seen, no strong nonlinearities are present in the response, ${\mathit{\chi}}_{\mathit{\gamma}}^{\left(\text{O}\right)}\left(t\right)$ is well recovered without any pretransformations from
where $\mathrm{\Delta}{C}_{\text{O}}^{\text{rad}}\left(t\right)$ is the ocean carbon response in the rad setup. Because the quality of the recovery depends on the quality of the estimate of the noise in the data, and this estimate may be improved by enforcing monotonicity (see Torres Mendonça et al., 2021a), we assume that ${\mathit{\chi}}_{\mathit{\gamma}}^{\left(\text{O}\right)}$ is monotonic and recover it enforcing monotonicity by the RFI method. An indication of the quality of the recovery will be obtained below by checking how well the recovered ${\mathit{\chi}}_{\mathit{\gamma}}^{\left(\text{O}\right)}$ predicts the model response in different perturbation experiments.
Recovery of ${\mathit{\chi}}_{\mathit{\gamma}}^{\left(\text{O}\right)}$ and linear regime of the radiative response of ocean carbon
We start by recovering ${\mathit{\chi}}_{\mathit{\gamma}}^{\left(\text{O}\right)}\left(t\right)$ from the 1 % rad experiment. To evaluate the quality of the recovery, following the procedure from Sect. A1 we employed the recovered ${\mathit{\chi}}_{\mathit{\gamma}}^{\left(\text{O}\right)}\left(t\right)$ in Eq. (A4) to predict the 0.5 % and 0.75 % rad responses. The prediction error is plotted in Fig. A2a as a function of the final forcing strength in the 1 % rad experiment. As seen, the error decreases with increasing final forcing strength, indicating that no strong nonlinearities are present in the response. As a result, the response can be considered linear for the whole range of forcing strengths in the 1 % rad experiment (temperatures up to around 4 K). Hence we choose to recover ${\mathit{\chi}}_{\mathit{\gamma}}^{\left(\text{O}\right)}\left(t\right)$ for the investigation in the main text by taking data from the full 1 % rad experiment.
The resulting ${\mathit{\chi}}_{\mathit{\gamma}}^{\left(\text{O}\right)}\left(t\right)$ recovered in this way is shown in Fig. A2b. The negativity of ${\mathit{\chi}}_{\mathit{\gamma}}^{\left(\text{O}\right)}\left(t\right)$ reflects the fact that as temperatures rise, globally the ocean loses carbon to the atmosphere. This is consistent with the results shown in Fig. A2c, which shows the model responses from the different rad experiments and the prediction from the recovered ${\mathit{\chi}}_{\mathit{\gamma}}^{\left(\text{O}\right)}\left(t\right)$. Because in these experiments only the radiative effect of CO_{2} is active (while the biogeochemical effect is switched off), temperatures rise, leading to the mentioned global loss of ocean carbon to the atmosphere reflected by the negative responses. The quality of agreement between the model responses and the predictions from ${\mathit{\chi}}_{\mathit{\gamma}}^{\left(\text{O}\right)}\left(t\right)$ suggests that indeed the recovered ${\mathit{\chi}}_{\mathit{\gamma}}^{\left(\text{O}\right)}\left(t\right)$ characterizes the radiative response of ocean carbon not only for the few perturbation experiments considered here but to any temporal development of weak CO_{2} perturbations and is therefore a true property of the MPIESM itself.
In summary, since the response can be considered linear over the whole 1 % rad experiment, we chose to derive ${\mathit{\chi}}_{\mathit{\gamma}}^{\left(\text{O}\right)}\left(t\right)$ in Sects. 3 and 4 from Eq. (A4) taking data from the full experiment.
A4 Generalized sensitivities ${\mathit{\chi}}_{\mathit{\alpha}}^{\left(\text{L}\right)}$ and ${\mathit{\chi}}_{\mathit{\alpha}}^{\left(\text{O}\right)}$
Pretransformation techniques to recover ${\mathit{\chi}}_{\mathit{\alpha}}^{\left(\text{L}\right)}$ and ${\mathit{\chi}}_{\mathit{\alpha}}^{\left(\text{O}\right)}$
Following Sect. A2 we recover ${\mathit{\chi}}_{\mathit{\alpha}}^{\left(\text{L}\right)}\left(t\right)$ and ${\mathit{\chi}}_{\mathit{\alpha}}^{\left(\text{O}\right)}\left(t\right)$ by employing two different techniques. In the first technique we recover the sensitivities from Eqs. (4) and (5), i.e. without any pretransformation. From this first technique we will also obtain an estimate of the linear regime for the temperature responses as described in the generalized α–β–γ framework.
Because CO_{2} radiative forcing is known to increase logarithmically with atmospheric CO_{2} concentration (Myhre et al., 1998), we also derive ${\mathit{\chi}}_{\mathit{\alpha}}^{\left(\text{L}\right)}\left(t\right)$ and ${\mathit{\chi}}_{\mathit{\alpha}}^{\left(\text{O}\right)}\left(t\right)$ using a logarithmic pretransformation:
As explained in Torres Mendonça et al. (2021b, Sect. 4), the resulting ${\mathit{\chi}}_{\mathit{\alpha},\mathrm{ln}}^{\left(\text{L}\right)}\left(t\right)={\mathit{\chi}}_{\mathit{\alpha}}^{\left(\text{L}\right)}\left(t\right)$ and ${\mathit{\chi}}_{\mathit{\alpha},\mathrm{ln}}^{\left(\text{O}\right)}\left(t\right)={\mathit{\chi}}_{\mathit{\alpha}}^{\left(\text{O}\right)}\left(t\right)$.
Recovery of ${\mathit{\chi}}_{\mathit{\alpha}}^{\left(\text{L}\right)}$, ${\mathit{\chi}}_{\mathit{\alpha}}^{\left(\text{O}\right)}$, and linear regime of temperature responses
We start by recovering ${\mathit{\chi}}_{\mathit{\alpha}}^{\left(\text{L}\right)}\left(t\right)$ and ${\mathit{\chi}}_{\mathit{\alpha}}^{\left(\text{O}\right)}\left(t\right)$ without any pretransformation (Eqs. 4 and 5) taking data from the 1 % rad experiment. As in the previous subsections, to evaluate the quality of the recoveries we employed them once more in Eqs. (4) and (5) to predict the responses from the 0.5 % and 0.75 % rad experiments. The prediction error is plotted in Figs. A3a and A4a. As seen, in both figures the error decreases until about 279 ppm, presenting after that in Fig. A3a a clear increase and in Fig. A4a a slight decrease followed by an increase for the 0.5 % rad response and an approximately constant behaviour for the 0.75 % rad response. Therefore we conservatively estimate the linear regime as perturbation strengths below 279 ppm for both land and ocean temperature responses.
To assess whether the logtransform technique (Eqs. A5 and A6) indeed improves the recovery of ${\mathit{\chi}}_{\mathit{\alpha}}^{\left(\text{L}\right)}\left(t\right)$ and ${\mathit{\chi}}_{\mathit{\alpha}}^{\left(\text{O}\right)}\left(t\right)$, analogously to Sect. A2 we check whether Eqs. (A5) and (A6) indeed account for nonlinearities in the response. If this is the case, then Eqs. (A5) and (A6) should predict the 0.5 % and 0.75 % rad responses better than Eqs. (4) and (5). To evaluate this, we recovered ${\mathit{\chi}}_{\mathit{\alpha}}^{\left(\text{L}\right)}\left(t\right)$ and ${\mathit{\chi}}_{\mathit{\alpha}}^{\left(\text{O}\right)}\left(t\right)$ from Eqs. (A5) and (A6) taking data once more from the 1 % rad experiment and then employed the recovered response functions in Eqs. (A5) and (A6) to predict the 0.5 % and 0.75 % rad responses. The resulting prediction error is shown in Figs. A3b and A4b. As seen in both figures, now the prediction error clearly decreases for increasing final forcing strength, indicating that no strong nonlinearities are present in the response. This clear decreasing trend in the prediction error indicates that the logtransform technique (Eqs. A5 and A6) is indeed more appropriate to recover ${\mathit{\chi}}_{\mathit{\alpha}}^{\left(\text{L}\right)}\left(t\right)$ and ${\mathit{\chi}}_{\mathit{\alpha}}^{\left(\text{O}\right)}\left(t\right)$.
Therefore in Figs. A3c and A4c we show the response functions recovered with the logtransform technique. Plotted are recoveries by the RFI method both enforcing and not enforcing monotonicity. As seen, the recoveries obtained without enforcing monotonicity are approximately monotonic, so differences between the two types of recovery are small. Hence in the next sections we choose to derive ${\mathit{\chi}}_{\mathit{\alpha}}^{\left(\text{L}\right)}\left(t\right)$ and ${\mathit{\chi}}_{\mathit{\alpha}}^{\left(\text{O}\right)}\left(t\right)$ employing the logtransform and enforcing monotonicity by the RFI method. As already mentioned above (see also Torres Mendonça et al., 2021a), the advantage of enforcing monotonicity is that one may improve the noise level estimate used in the RFI method and thereby improve the quality of the recovery.
To obtain evidence that the recovered ${\mathit{\chi}}_{\mathit{\alpha}}^{\left(\text{L}\right)}\left(t\right)$ and ${\mathit{\chi}}_{\mathit{\alpha}}^{\left(\text{O}\right)}\left(t\right)$ indeed characterize the land and ocean temperature responses to any temporal development of weak CO_{2} perturbations, we show how well they predict additional experiments. Following Sect. A2 we first show the predictive power of ${\mathit{\chi}}_{\mathit{\alpha}}^{\left(\text{L}\right)}\left(t\right)$ and ${\mathit{\chi}}_{\mathit{\alpha}}^{\left(\text{O}\right)}\left(t\right)$ when employed in Eqs. (4) and (5), i.e. without pretransformations. In Figs. A3d and A4d we plot the predictions by Eqs. (4) and (5) taking the recoveries selected above (logtransform technique enforcing monotonicity). The responses for the 0.75 %, 1.5 %, and 2 % rad experiments are marked with dots where the forcing strength exceeds 279 ppm – the maximal forcing strength corresponding to the linear regime estimated from Figs. A3a and A4a. As seen, the recovered response functions predict the model responses with some overestimation but still with reasonable quality of agreement^{5} for forcing strengths within the estimated linear regime, i.e. for the entire 1.1× CO_{2} and 0.5 % rad experiments and for the 0.75 %, 1.5 %, and 2 % rad experiments for values preceding the dots. For these three latter experiments predictions start to strongly deviate from the responses as forcing strengths exceed the estimated linear regime. Likewise, the predictions fail basically for the whole 2× CO_{2} rad response because its forcing strength is larger than the linear regime estimate throughout the whole experiment.
As expected from the known logarithmic relationship between radiative forcing and CO_{2} concentration, the predictive power of ${\mathit{\chi}}_{\mathit{\alpha}}^{\left(\text{L}\right)}\left(t\right)$ and ${\mathit{\chi}}_{\mathit{\alpha}}^{\left(\text{O}\right)}\left(t\right)$ improves when employing Eqs. (A5) and (A6) instead of Eqs. (4) and (5) (see Figs. A3e and A4e). Almost all responses are well predicted for the whole time series. The exceptions are the last years of the 1.5 % and 2 % rad responses, which indicates a deviation from linearity for those levels of perturbation strength.
Torres Mendonça et al. (2021b)Torres Mendonça et al. (2021b)As in Sect. A2, overall the prediction results suggest that indeed the recovered response functions characterize the land and ocean temperature responses to any temporal development of atmospheric CO_{2} up to certain perturbation strengths. While Eqs. (A5) and (A6) are more appropriate when the aim is simply to predict model responses because of their extended predictive power, Eqs. (4) and (5) are the equations actually employed in the generalized α–β–γ framework, so for the application of the framework one must consider the linear regime of the temperature responses estimated from Figs. A3a and A4a.
In summary, the conclusions from this subsection suggest that the best approach to derive ${\mathit{\chi}}_{\mathit{\alpha}}^{\left(\text{L}\right)}\left(t\right)$ and ${\mathit{\chi}}_{\mathit{\alpha}}^{\left(\text{O}\right)}\left(t\right)$ is the logtransform technique (Eqs. A5 and A6), taking data from the full 1 % rad experiment. The linear regime for the land and ocean temperature responses in the generalized α–β–γ framework (Eqs. 5 and 4) is about 279 ppm, as estimated from Figs. A3a and A4a.
A5 Selection of techniques to recover ${\mathit{\chi}}_{\mathit{\gamma}}^{\left(\text{L}\right)}\left(t\right)$ and ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{L}\right)}\left(t\right)$
To recover ${\mathit{\chi}}_{\mathit{\gamma}}^{\left(\text{L}\right)}\left(t\right)$ and ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{L}\right)}\left(t\right)$ we use the conclusions from Torres Mendonça et al. (2021b). Because the radiative response of land carbon can be considered linear for the whole range of perturbation strengths in the 1 % rad experiment (Torres Mendonça et al., 2021b, Sect. 3), ${\mathit{\chi}}_{\mathit{\gamma}}^{\left(\text{L}\right)}\left(t\right)$ is derived taking data without any pretransformation from the full 1 % rad experiment by employing
where $\mathrm{\Delta}{C}_{\text{L}}^{\text{rad}}\left(t\right)$ is the land carbon response in the rad setup. Because the obtained ${\mathit{\chi}}_{\mathit{\gamma}}^{\left(\text{L}\right)}\left(t\right)$ is monotonic (without enforcing monotonicity), we assume monotonicity when deriving ${\mathit{\chi}}_{\mathit{\gamma}}^{\left(\text{L}\right)}\left(t\right)$ for all other CMIP5 models so that all generalized sensitivities are derived enforcing monotonicity by the RFI method (as mentioned above and explained in Sect. 3.5 of Torres Mendonça et al., 2021a, this assumption may further improve the quality of the recovery).
To derive ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{L}\right)}\left(t\right)$, conclusions from Sect. 4 of Torres Mendonça et al. (2021b) suggest that the best approach is the NPPtransform technique analysed there: first derive the response function χ_{NPP}(t) from
and then transform the obtained χ_{NPP}(t) into the desired ${\mathit{\chi}}_{\mathit{\beta}}^{\left(\text{L}\right)}\left(t\right)$ via
where c_{PI} is the preindustrial value for atmospheric CO_{2}. Because the prediction error plot for χ_{NPP}(t) presents minima (see Fig. 6c in Torres Mendonça et al., 2021b), following the reasoning from Sect. A2 we choose conservatively to take data only until the first minimum, which corresponds to 279 ppm or 70 years of the 1 % bgc experiment.
A6 Summary of best techniques for recovery of the generalized sensitivities
With the results presented in the preceding subsections we can now summarize the best techniques identified to recover the generalized sensitivities for our study. A general summary of the identified techniques is given in Table A2. The table indicates which is (1) the best pretransformation technique to derive each generalized sensitivity, (2) the respective equation used for the derivation, (3) the optimal time series length taken for recovery using the best pretransformation technique, and (4) the linear regime of the respective response in the generalized α–β–γ framework (i.e. the linear regime for each term on the righthand side of Eqs. 2–5).
Because the linear regime for the biogeochemical response of land and ocean carbon is restricted to forcing strengths even smaller than that for temperature responses (4) and (5) – and obviously also smaller than that for the radiative responses of land and ocean carbon, which are linear for the whole 1 % rad experiment, as seen in Torres Mendonça et al. (2021b) and Sect. A3 – the applicability of the generalized α–β–γ framework as a whole is limited by the linear regime of the biogeochemical responses, which is estimated as forcing strengths up to about 94 ppm.
With this subsection we complete the recovery of all generalized sensitivities for the MPIESM. The approaches selected here are employed to recover the generalized sensitivities for all CMIP5 models in the main text.
In this appendix we explain in detail how we derived the airborne fraction $\stackrel{\mathrm{\u0303}}{A}\left(p\right)$ from emissiondriven simulations for the demonstration of Sect. 3. For the derivation we followed the three steps described in Sect. 3.1.
In the first step, we recovered χ_{ζ}(t) taking data from an impulseemission experiment. In this experiment, starting from a standard preindustrial control run (see “esmControl” in Taylor et al., 2012), we perturbed the system by a small impulse in emissions of 100 Pg C yr^{−1} (as in Joos et al., 2013) during the first year. Therefore we name this experiment Impulse^{100}. The advantage of using data from this type of experiment is that the impulse response is from a practical point of view directly the desired linear response function χ_{ζ}(t), so errors from the illposedness of Eq. (16) (Torres Mendonça et al., 2021a) can be avoided. Nevertheless, the disadvantage is that the small impulse perturbation strength leads to a response with poor signaltonoise ratio, which deteriorates the recovery of χ_{ζ}(t). Therefore, instead of taking data from only one realization of the experiment, we performed an ensemble of five realizations with initial conditions taken 100 years apart from one another (as in Lembo et al., 2020) and then took the data from the ensemble average of the response. This procedure is in agreement with linear response theory because strictly the linear response function characterizes only the ensemble average of the response (e.g. Torres Mendonça et al., 2021a). Taking the data from the ensemble average we employed the RFI method to recover χ_{ζ}(t). Although in principle a special method is not needed to recover χ_{ζ}(t) because the impulse response is directly the linear response function, by employing the RFI method we obtain not only χ_{ζ}(t) but also the spectrum of internal timescales of the response (Torres Mendonça et al., 2021a), from which the Laplace transform ${\stackrel{\mathrm{\u0303}}{\mathit{\chi}}}_{\mathit{\zeta}}\left(p\right)$ in Eq. (18) can be analytically calculated.
The resulting impulse response after taking the ensemble average and the fit by the recovered χ_{ζ}(t) are shown in Fig. B1a. One sees that the experiment is not a perfect impulse experiment because the impulse extends even beyond the first year. This may be related to internal interpolations in the model when computing the emissions within the time interval of 1 year. Such a problem leads to an error in the estimation of χ_{ζ}(0), which must be 1 by conservation of mass: for an impulse in emissions E(t)=aδ(t), it must be that at t=0 the impulse response $\mathrm{\Delta}{c}_{\mathit{\delta}}\left(\mathrm{0}\right)=a{\mathit{\chi}}_{\mathit{\zeta}}\left(\mathrm{0}\right)=a$ because land carbon and ocean carbon have no time to react; hence χ_{ζ}(0)=1. To avoid this error, we recalculated χ_{ζ}(t) using the same regularization parameter obtained from the RFI method but employing a Lagrange multiplier to account for this constraint (see Appendix C). The result for both recoveries is shown in Fig. B1b. The response functions are almost identical except for χ_{ζ}(0), which is corrected by enforcing the constraint χ_{ζ}(0)=1.
To make sure that the impulse response is within the linear regime and therefore that the recovery of χ_{ζ}(t) is not spoiled by nonlinearities, in line with the procedure in Sect. A, we check the quality of the recovered χ_{ζ}(t) by employing it to predict additional emissiondriven experiments via Eq. (16). For the additional experiments we chose stepemission experiments where starting from the control run the system is perturbed by abrupt, constant emissions of 1.43, 2.86, and 5.71 Pg C yr^{−1}, which by the end of the simulation result in cumulative emissions of 200, 400, and 800 Pg C. For this reason we name these experiments Step^{200}, Step^{400}, and Step^{800}, respectively. The quality of the prediction is shown in Fig. B1c. As seen, the obtained χ_{ζ}(t) can predict almost all responses with reasonable accuracy for the whole time series. The only exception is over the last years of the Step^{800} experiment. The discrepancy encountered there is nevertheless in agreement with results from Torres Mendonça et al. (2021b) and from Appendix A that show that strong nonlinear contributions to the biogeochemical response of land and ocean carbon start to arise at about 94 ppm. As a consequence, the identified discrepancy is not a result of nonlinearities spoiling the recovery of χ_{ζ}(t) but rather a result of nonlinearities that arise in the response of the Step^{800} experiment, which can therefore not be completely predicted by χ_{ζ}(t).
Hence, overall these results suggest that the recovery of χ_{ζ}(t) is not spoiled by nonlinearities and is thus a good candidate to be used to compute the airborne fraction. Therefore χ_{ζ}(t) was Laplacetransformed and the airborne fraction was finally derived by applying Eq. (18).
Further evidence of the reliability of our numerics is obtained by examining the agreement of the resulting airborne fraction (Fig. 1) with theoretical expectations: $\stackrel{\mathrm{\u0303}}{A}\left(p\right)$ indeed converges to 1 for $\mathrm{1}/p\to \mathrm{0}$, as expected by Eq. (19), and $\mathrm{0}\le \stackrel{\mathrm{\u0303}}{A}\left(p\right)\le \mathrm{1}$ for all timescales $\mathrm{1}/p$, as expected by the considerations of Appendix D.
This appendix complements Appendix B to show how χ_{ζ}(t) is calculated by the RFI method enforcing in addition the theoretically exact relation χ_{ζ}(0)=1. The RFI method recovers a response function χ(t) employing the Tikhonov–Phillips regularization, with the regularization parameter λ determined from an estimate of the noise in the data (Torres Mendonça et al., 2021a). But the RFI method does not account for the desired constraint χ_{ζ}(0)=1. Therefore to derive χ_{ζ}(t) enforcing this constraint we proceed in two steps. First, we derive χ_{ζ}(t) by the RFI method without the constraint, obtaining thereby the value of λ. Then, we use the obtained λ to derive χ_{ζ}(t) employing the same regularization procedure from the RFI method but accounting in addition for the constraint χ_{ζ}(0)=1 by the method of Lagrange multipliers (Sundaram et al., 1996; Selesnick, 2013).
After obtaining λ by the RFI method in the first step, the derivation of χ_{ζ}(t) in the second step proceeds as follows. As in the RFI method, we assume χ_{ζ}(t) to be given by
where q(τ) is the spectrum of internal timescales. By plugging Eq. (C1) into Eq. (16) and prescribing a distribution of timescales τ, the problem of finding χ_{ζ}(t) boils down to finding the spectrum q(τ). Once q(τ) is derived, χ_{ζ}(t) follows from Eq. (C1). For more details please refer to Torres Mendonça et al. (2021a).
To understand how the constraint can be enforced, one has to consider
Equation (C2) can be discretized using the prescribed distribution of timescales equally spaced at a logarithmic scale between maximum and minimum values τ_{max} and τ_{min} (Torres Mendonça et al., 2021a). This discretization gives
where M is the number of q_{j} terms and $\mathrm{\Delta}{\mathrm{log}}_{\mathrm{10}}\mathit{\tau}:=({\mathrm{log}}_{\mathrm{10}}{\mathit{\tau}}_{\mathrm{max}}{\mathrm{log}}_{\mathrm{10}}{\mathit{\tau}}_{\mathrm{min}})/M$ is the spacing between the timescales. Following Torres Mendonça et al. (2021a, Sect. 4.2), we take M=30, τ_{min}=0.1, and τ_{max}=10^{5}. Using Eq. (C3) one can write the desired constraint χ_{ζ}(0)=1 in a discrete formulation as
where C is the row matrix $\mathbf{C}:=[\mathrm{1},\mathrm{1},\mathrm{\dots},\mathrm{1}]\mathrm{\Delta}{\mathrm{log}}_{\mathrm{10}}\mathit{\tau}$.
Knowing how to discretely account for the desired constraint the spectrum q can now be derived. The procedure consists of minimizing the standard cost function employed in the Tikhonov–Phillips regularization (see e.g. Hansen, 2010, p. 60), which is also done in the RFI method but now subject to the constraint (C4); i.e.
where $\left\right\cdot \left\right$ denotes the Euclidean norm, q_{λ} is the spectrum vector recovered by regularization, ΔY is the response data vector, A is the matrix obtained after discretization (Torres Mendonça et al., 2021a, Sect. 3.2), and λ is the regularization parameter. The solution to Eq. (C5) can be obtained by the method of Lagrange multipliers (Selesnick, 2013). We first define the Lagrangian:
where μ is the Lagrange multiplier. Taking $\partial \mathcal{L}/\partial {\mathit{q}}_{\mathit{\lambda}}=\mathbf{0}$ and $\partial \mathcal{L}/\partial \mathit{\mu}=\mathrm{0}$ and solving the resulting system for q_{λ} leads to
where I is the identity matrix. Since λ was already obtained in the first step, Eq. (C7) gives a discrete approximation to the spectrum q(τ). Plugging this spectrum in a discrete form of Eq. (C1) (for discretization details see Torres Mendonça et al., 2021a, Appendix B), we obtain χ_{ζ}(t).
The response function χ_{ζ}(t) defined by Eq. (16) is in terms of Laplace transforms closely related to the generalized airborne fraction $\stackrel{\mathrm{\u0303}}{A}\left(p\right)$ defined by Eq. (18). In this appendix it is shown that if χ_{ζ}(t) is nonnegative and monotonically decreasing for all times t (as suggested by the results in Appendix B), then $\mathrm{0}\le \stackrel{\mathrm{\u0303}}{A}\left(p\right)\le \mathrm{1}$ for all timescales $\mathrm{1}/p$, as claimed in Appendix B. This follows from the separate proofs of the two inequalities involved in this claim:

Show ${\mathit{\chi}}_{\mathit{\zeta}}\left(t\right)\ge \mathrm{0}\Rightarrow \stackrel{\mathrm{\u0303}}{A}\left(p\right)\ge \mathrm{0}$: by Eq. (18) one immediately finds
$$\begin{array}{}\text{(D1)}& \stackrel{\mathrm{\u0303}}{A}\left(p\right)=p\underset{\mathrm{0}}{\overset{\mathrm{\infty}}{\int}}{\mathit{\chi}}_{\mathit{\zeta}}\left(t\right){e}^{pt}\mathrm{d}t\ge \mathrm{0},\end{array}$$where for the inequality χ_{ζ}(t)≥0 was used.

Show $\mathrm{d}{\mathit{\chi}}_{\mathit{\zeta}}\left(t\right)/\mathrm{d}t\le \mathrm{0}\Rightarrow \stackrel{\mathrm{\u0303}}{A}\left(p\right)\le \mathrm{1}$: as a prerequisite for this proof one needs to know that
$$\begin{array}{}\text{(D2)}& \stackrel{\mathrm{\u0303}}{A}\left(p\right)=\stackrel{\mathrm{\u0303}}{{\mathit{\chi}}_{\mathit{\zeta}}^{\prime}}\left(p\right)+\mathrm{1},\end{array}$$where the apostrophe notation is used for time derivatives. This relation follows by (i) noting that by the general rules of Laplace transforms $\stackrel{\mathrm{\u0303}}{{\mathit{\chi}}_{\mathit{\zeta}}^{\prime}}\left(p\right)=p\stackrel{\mathrm{\u0303}}{{\mathit{\chi}}_{\mathit{\zeta}}}\left(p\right){\mathit{\chi}}_{\mathit{\zeta}}\left(\mathrm{0}\right)$, (ii) using χ_{ζ}(0)=1 (see Appendix B), and (iii) inserting $\stackrel{\mathrm{\u0303}}{A}\left(p\right)=p\stackrel{\mathrm{\u0303}}{{\mathit{\chi}}_{\mathit{\zeta}}}\left(p\right)$ (see Eq. 18). The rest of the proof is obtained from noting that
$$\begin{array}{}\text{(D3)}& \stackrel{\mathrm{\u0303}}{{\mathit{\chi}}_{\mathit{\zeta}}^{\prime}}\left(p\right)=\underset{\mathrm{0}}{\overset{\mathrm{\infty}}{\int}}{\mathit{\chi}}_{\mathit{\zeta}}^{\prime}\left(t\right){e}^{pt}\mathrm{d}t\le \mathrm{0}\end{array}$$because ${\mathit{\chi}}_{\mathit{\zeta}}^{\prime}\left(t\right)\le \mathrm{0}$ is assumed. Using this in Eq. (D2) completes the proof.
When calculating the α–β–γ sensitivities in Fig. 3, we take into account the uncertainty in the choice of the initial values in Eqs. (20)–(22):
where the initial values ${C}_{X}^{\mathrm{0}}$, ${T}_{X}^{\mathrm{0}}$, and c^{0} are taken as the mean from the control simulation with uncertainty (standard deviation) ${\mathit{\sigma}}_{{C}_{X}^{\mathrm{0}}}$, ${\mathit{\sigma}}_{{T}_{X}^{\mathrm{0}}}$, and ${\mathit{\sigma}}_{{c}^{\mathrm{0}}}$ = 0 (because in the considered simulations atmospheric CO_{2} is prescribed). Assuming small, independent uncertainties, they propagate to the sensitivities, following Barlow (1989):
In the present appendix we demonstrate exemplarily for MPIESM that from the generalized airborne fraction A(t) one can indeed predict, for sufficiently weak emissions, the airborne fraction in its standard definition AF(t), as claimed in Sect. 6. For this demonstration we take data from two perturbation experiments performed with the MPIESM: the Step^{400} experiment – an emissiondriven experiment where starting from the control run a constant amount of 2.86 Pg C yr^{−1} is emitted every year until cumulative emissions reach 400 Pg C (see Appendix B) – and the 1 % fully coupled experiment – a concentrationdriven experiment where starting from the control run the atmospheric CO_{2} concentration is increased by 1 % every year (see Table A1). The demonstration is carried out by first computing the true standard airborne fraction AF(t) via its definition (12) and then comparing it with its prediction obtained from the generalized airborne fraction A(t) via Eqs. (13) and (12).
To compute the true standard airborne fraction AF(t) we directly evaluate the defining Eq. (12) by using $\mathrm{d}{C}_{\text{A}}/\mathrm{d}t$ and E(t) as obtained from the two experiments. The accumulation rate $\mathrm{d}{C}_{\text{A}}/\mathrm{d}t$ is for both experiments numerically calculated from the atmospheric carbon content C_{A}(t). Concerning the emissions E(t), since in the Step^{400} experiment they are prescribed, for this experiment we take them as constant and equal to 2.86 Pg C yr^{−1}. For the 1 % fully coupled experiment the situation is different: here atmospheric CO_{2} concentrations – not emissions – are prescribed; therefore for this experiment we infer the emissions that would be needed to produce the respective changes in atmospheric CO_{2} from the evolution of carbon content in the atmosphere, land, and ocean via the carbon balance Eq. (1). Results for the true AF(t) are shown in Fig. F1.
To predict now the standard airborne fraction AF(t) from the generalized airborne fraction A(t) we proceed as follows. First we compute A(t) by Laplace inverting relation (D2) so that
where χ_{ζ}(t) was already obtained in Appendix B and we take for the numerics the Kronecker delta δ_{i} instead of δ(t). We then plug the resulting A(t) and the emissions E(t) – obtained as described above – into Eq. (13) to calculate $\mathrm{d}{C}_{\text{A}}/\mathrm{d}t$. The standard airborne fraction AF(t) is then finally predicted by plugging the resulting $\mathrm{d}{C}_{\text{A}}/\mathrm{d}t$, once more together with the emissions E(t), now into Eq. (12).
The results are compared in Fig. F1. As seen in Fig. F1a, because in the Step^{400} experiment atmospheric CO_{2} changes are below our estimated linear regime of 94 ppm (see Table A2), there the predicted standard airborne fraction fits the true standard airborne fraction over the whole simulation period. The large variability in the true airborne fraction arises from the illposedness (e.g. Torres Mendonça et al., 2021a) of the numerical differentiation needed to compute $\mathrm{d}{C}_{\text{A}}/\mathrm{d}t$, which substantially amplifies the noise from the already noisy response of C_{A}(t) that results from the small emissions forcing strength of this experiment.
Such large variability in AF(t) is not present in the 1 % fully coupled experiment, as shown by Fig. F1b. Here, the larger signaltonoise ratio of $\mathrm{d}{C}_{\text{A}}/\mathrm{d}t$ allows for a better quality of fit of the prediction from the generalized airborne fraction. But since in this experiment changes in atmospheric CO_{2} get larger than our linear regime estimate, the prediction works only for that first part of the time series where atmospheric CO_{2} concentrations are sufficiently small.
In this appendix we show how from α and γ sensitivities defined using separate land and ocean temperatures one can compute their analogues defined by means of a single global temperature, as claimed in the Outlook (Sect. 7). Global temperature in a model is obtained by
where A_{i} and ΔT_{i} are the area and temperature of grid box i, i∈L and i∈O indicate sum over grid boxes on land and ocean, ΔT_{L} and ΔT_{O} are land and ocean temperatures, and F_{L} and F_{O} are the fractions of global area occupied by land and ocean.
Using a single global temperature, α is defined by
Using separate land and ocean temperatures one defines
Plugging Eqs. (G1), (G3), and (G4) into (G2) gives
Taking a single global temperature, γ is defined by
where X denotes the carbon response over land (L) or ocean (O). Using separate land and ocean temperatures, γ can be defined by
Inserting Eqs. (G3), (G4), (G6), and (G1) into (G7) gives
Since the Laplacetransformed formulation of the generalized framework is completely analogous to that of the original α–β–γ framework, Eqs. (G5) and (G8) extend straightforwardly to the respective generalized sensitivities:
The scripts and data used to produce the results in this paper can be found at https://hdl.handle.net/21.11116/0000000CF6A27 (Torres Mendonca et al., 2023).
GLTM led the study, performed the analysis, and wrote the first draft. All authors contributed to the conception and scientific content of the study. The final manuscript was jointly prepared by GLTM and CHR.
The contact author has declared that none of the authors has any competing interests.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
We would like to thank Matteo Puglini, Peter Landschützer, Ian Enting, and Vivek Arora for their helpful comments on the manuscript and Ludwig Lierhammer for preparing the CMIP5 data used in the study. Guilherme L. Torres Mendonça was supported by the Max Planck Institute for Meteorology and by research grant no. 2023/045795 from the São Paulo Research Foundation (FAPESP).
Guilherme L. Torres Mendonça was supported by the Max Planck Institute for Meteorology and by research grant no. 2023/045795 from the São Paulo Research Foundation (FAPESP).
The article processing charges for this openaccess publication were covered by the Max Planck Society.
This paper was edited by Anja Rammig and reviewed by Ian Enting and Vivek Arora.
Adloff, M., Reick, C. H., and Claussen, M.: Earth system model simulations show different feedback strengths of the terrestrial carbon cycle under glacial and interglacial conditions, Earth Syst. Dynam., 9, 413–425, https://doi.org/10.5194/esd94132018, 2018. a, b, c, d, e
Archer, D., Eby, M, Brovkin, V., Ridgwell, A., Cao, L., Mikolajewicz, U., Caldeira, K., Matsumoto, K., Munhoven, G., Montenegro, A., and Tokos, K.: Atmospheric lifetime of fossil fuel carbon dioxide, Annual Rev. Earth Planet. Sci., 37, 117–134, 2009. a, b
Arneth, A., Harrison, S., Zaehle, S., Tsigaridis, K., Menon, S., Bartlein, P., Feichter, J., Korhola, A., Kulmala, M., O'Donnell, D., Schurgers, G., Sorvari, S., and Vesala, T.: Terrestrial biogeochemical feedbacks in the climate system, Nat. Geosci., 3, 525–532, 2010. a
Arora, V., Boer, G., Friedlingstein, P., Eby, M., Jones, C, Christian, J., Bonan, G., Bopp, L., Brovkin, V., Cadule, P., Hajima, T., Ilyina, T., Lindsay, K., Tjiputra, J., and Wu, T.: Carbon–concentration and carbon–climate feedbacks in CMIP5 Earth system models, J. Climate, 26, 5289–5314, 2013. a, b, c, d, e, f, g, h, i, j, k, l, m, n
Arora, V. K., Katavouta, A., Williams, R. G., Jones, C. D., Brovkin, V., Friedlingstein, P., Schwinger, J., Bopp, L., Boucher, O., Cadule, P., Chamberlain, M. A., Christian, J. R., Delire, C., Fisher, R. A., Hajima, T., Ilyina, T., Joetzjer, E., Kawamiya, M., Koven, C. D., Krasting, J. P., Law, R. M., Lawrence, D. M., Lenton, A., Lindsay, K., Pongratz, J., Raddatz, T., Séférian, R., Tachiiri, K., Tjiputra, J. F., Wiltshire, A., Wu, T., and Ziehn, T.: Carbon–concentration and carbon–climate feedbacks in CMIP6 models and their comparison to CMIP5 models, Biogeosciences, 17, 4173–4222, https://doi.org/10.5194/bg1741732020, 2020. a, b, c, d, e
Barlow, R. J.: Statistics: A Guide to the Use of Statistical Methods in the Physical Sciences, The Manchester Physics Series, John Wiley & Sons, ISBN 0471922943, 1989. a, b
Beerends, R. J., ter Morsche, H. G., Van den Berg, J., and Van de Vrie, E.: Fourier and Laplace transforms, Cambridge University Press, 1st edn., ISBN 9780511673122, 2003. a
Bennedsen, M., Hillebrand, E., and Koopman, S. J.: Trend analysis of the airborne fraction and sink rate of anthropogenically released CO_{2}, Biogeosciences, 16, 3651–3663, https://doi.org/10.5194/bg1636512019, 2019. a
Bentsen, M., Bethke, I., Debernard, J., Drange, H., Heinze, C., Iversen, T., Kirkevåg, A., Seland, y., and Tjiputra, J.: cmip5 output1 NCC NorESM1M piControl, served by ESGF, World Data Center for Climate (WDCC) at DKRZ [data set], https://doi.org/10.1594/WDCC/CMIP5.NCCNMpc, 2011. a
Boer, G. and Arora, V.: Feedbacks in emissiondriven and concentrationdriven global carbon budgets, J. climate, 26, 3326–3341, 2013. a, b
Canadell, J., Monteiro, P., Costa, M., da Cunha, L. C., Cox, P., Eliseev, A., Henson, S., Ishii, M., Jaccard, S., Koven, C., Lohila, A., Patra, P., Piao, S., Rogelj, J., Syampungani, S., Zaehle, S., and Zickfeld, K.: Global Carbon and other Biogeochemical Cycles and Feedbacks, in: Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: MassonDelmotte, V., Zhai, P., Pirani, A., Connors, S., Péan, C., Berger, S., Caud, N., Chen, Y., Goldfarb, L., Gomis, M., Huang, M., Leitzell, K., Lonnoy, E., Matthews, J., Maycock, T., Waterfield, T., Yelekçi, O., Yu, R., and Zhou, B., 673–816, Cambridge University Press, https://doi.org/10.1017/9781009157896.007, 2021. a, b, c
Canadell, J. G., Le Quéré, C., Raupach, M. R., Field, C. B., Buitenhuis, E. T., Ciais, P., Conway, T. J., Gillett, N. P., Houghton, R., and Marland, G.: Contributions to accelerating atmospheric CO_{2} growth from economic activity, carbon intensity, and efficiency of natural sinks, P. Natl. Acad. Sci. USA, 104, 18866–18870, 2007. a
Caubel, A., Denvil, S., Foujols, M. A., Marti, O., Dufresne, J.L., Bopp, L., Cadule, P., Ethé, C., Idelkadi, A., Mancip, M., Masson, S., Mignot, J., Ionela, M., Balkanski, Y., Bekki, S., Bony, S., Braconnot, P., Brockman, P., Codron, F., Cozic, A., Cugnet, D., Fairhead, L., Fichefet, T., Flavoni, S., Guez, L., Guilyardi, E., Hourdin, F., Ghattas, J., Kageyama, M., Khodri, M., Labetoulle, S., Lefebvre, M.P., Levy, C., Li, L., Lott, F., Madec, G., Marchand, M., Meurdesoif, Y., Rio, C., Schulz, M., Swingedouw, D., Szopa, S., Viovy, N., and Vuichard, N.: IPSLCM5ALR model output prepared for CMIP5 1pctCO2 experiment, served by ESGF, World Data Center for Climate (WDCC) at DKRZ [data set], https://doi.org/10.1594/WDCC/CMIP5.IPILc1, 2016. a
Ciais, P., Sabine, C., Bala, G., Bopp, L., Brovkin, V., Canadell, J., Chhabra, A., DeFries, R., Galloway, J., Heimann, M., Jones, C., Le Quéré, C., Myneni, R., Piao, S., and Thornton, P.: Carbon and Other Biogeochemical Cycles, in: Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: Stocker, T., Qin, D., Plattner, G.K., Tignor, M., Allen, S., Boschung, J., Nauels, A., Xia, Y., Bex, V., and Midgley, P., 465–570, Cambridge University Press, ISBN 9781107661820, 2013. a, b
Dlugokencky, E. and Tans, P.: Trends in atmospheric carbon dioxide, National Oceanic and Atmospheric Administration, Global Monitoring Laboratory (NOAA/GML), http://www.gml.noaa.gov/gmd/ccgg/trends/global.html (last access: 26 January 2023), 2023. a
Dunne, J., John, J., Shevliakova, E., Stouffer, R., Griffies, S., Malyshev, S., Milly, P., Sentman, L., Adcroft, A., Cooke, W., Dunne, K., Hallberg, R., Harrison, M., Krasting, J., Levy, H., Phillips, P., Samuels, B., Spelman, M., Winton, M., Wittenberg, A., and Zadeh, N.: NOAA GFDL GFDLESM2M, esmFdbk1 experiment output for CMIP5 AR5, served by ESGF, World Data Center for Climate (WDCC) at DKRZ [data set], https://doi.org/10.1594/WDCC/CMIP5.NFCBe1, 2014a. a
Dunne, J., John, J., Shevliakova, E., Stouffer, R., Griffies, S., Malyshev, S., Milly, P., Sentman, L., Adcroft, A., Cooke, W., Dunne, K., Hallberg, R., Harrison, M., Krasting, J., Levy, H., Phillips, P., Samuels, B., Spelman, M., Winton, M., Wittenberg, A., and Zadeh, N.: NOAA GFDL GFDLESM2M, esmFixClim1 experiment output for CMIP5 AR5, served by ESGF, World Data Center for Climate (WDCC) at DKRZ [data set], https://doi.org/10.1594/WDCC/CMIP5.NGEMx1, 2014b. a
Dunne, J., John, J., Shevliakova, E., Stouffer, R., Griffies, S., Malyshev, S., Milly, P., Sentman, L., Adcroft, A., Cooke, W., Dunne, K., Hallberg, R., Harrison, M., Krasting, J., Levy, H., Phillips, P., Samuels, B., Spelman, M., Winton, M., Wittenberg, A., and Zadeh, N.: NOAA GFDL GFDLESM2M, piControl experiment output for CMIP5 AR5, served by ESGF, World Data Center for Climate (WDCC) at DKRZ [data set], https://doi.org/10.1594/WDCC/CMIP5.NGEMpc, 2014c. a
Dunne, J., John, J., Shevliakova, E., Stouffer, R., Griffies, S., Malyshev, S., Milly, P., Sentman, L., Adcroft, A., Cooke, W., Dunne, K., Hallberg, R., Harrison, M., Krasting, J., Levy, H., Phillips, P., Samuels, B., Spelman, M., Winton, M., Wittenberg, A., and Zadeh, N.: NOAA GFDL GFDLESM2M, 1pctCO2 experiment output for CMIP5 AR5, served by ESGF, World Data Center for Climate (WDCC) at DKRZ [data set], https://doi.org/10.1594/WDCC/CMIP5.NGEMc1, 2014d. a
Enting, I.: Ambiguities in the calibration of carbon cycle models, Inverse Problems, 6, 5, https://doi.org/10.1088/02665611/6/5/001, 1990. a
Enting, I.: Laplace transform analysis of the carbon cycle, Environ. Model. Softw., 22, 1488–1497, 2007. a, b, c, d
Enting, I. and Clisby, N.: Estimates of climatic influence on the carbon cycle, Earth Syst. Dynam. Discuss. [preprint], https://doi.org/10.5194/esd201941, 2019. a, b, c, d, e, f, g, h
Enting, I. G.: Response function analysis of carbon dioxide and climate using the PadéLaplace technique, AIMS Geosciences, 8, 346–365, 2022. a, b, c, d, e
Exbrayat, J.F., Pitman, A. J., and Abramowitz, G.: Response of microbial decomposition to spinup explains CMIP5 soil carbon range until 2100, Geosci. Model Dev., 7, 2683–2692, https://doi.org/10.5194/gmd726832014, 2014. a
Eyring, V., Gillett, N., Achuta Rao, K., Barimalala, R., Barreiro Parrillo, M., Bellouin, N., Cassou, C., Durack, P., Kosaka, Y., McGregor, S., Min, S., Morgenstern, O., and Sun, Y.: Human Influence on the Climate System, Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, 423–552, https://doi.org/10.1017/9781009157896.005, 2021. a
Friedlingstein, P.: Carbon cycle feedbacks and future climate change, Philos. T. Roy. Soc. A, 373, 20140421, https://doi.org/10.1098/rsta.2014.0421, 2015. a
Friedlingstein, P., Dufresne, J.L., Cox, P., and Rayner, P.: How positive is the feedback between climate change and the carbon cycle?, Tellus B, 55, 692–700, 2003. a, b, c, d, e, f, g, h
Friedlingstein, P., Cox, P., Betts, R., Bopp, L. von Bloh, W., Brovkin, V., Cadule, P., Doney, S., Eby, M., Fung, I., Bala, G., John, J., Jones, C., Joos, F., Kato, T., Kawamiya, M., Knorr, W., Lindsay, K., Matthews, H., Raddatz, T., Rayner, P., Reick, C., Roeckner, E., Schnitzler, K.G., Schnur, R., Strassmann, K., Weaver, A. J., Yoshikawa C., and Zeng, N.: Climate–carbon cycle feedback analysis: results from the C4MIP model intercomparison, J. Climate, 19, 3337–3353, 2006. a, b, c
Friedlingstein, P., Jones, C., and Arora, V.: Special collection “ClimateCarbon Interactions in the CMIP5 Earth System Models (C4MIP)” J. Climate, https://journals.ametsoc.org/collection/C4MIP (last access: 10 April 2024), 2013–2016. a
Friedlingstein, P., Meinshausen, M., Arora, V. K., Jones, C. D., Anav, A., Liddicoat, S. K., and Knutti, R.: Uncertainties in CMIP5 climate projections due to carbon cycle feedbacks, J. Climate, 27, 511–526, 2014. a, b
Friend, A. D., Lucht, W., Rademacher, T. T., Keribin, R., Betts, R., Cadule, P., Ciais, P., Clark, D. B., Dankers, R., Falloon, P. D., Ito, A., Kahana, R., Kleidon, A., Lomas, M. R., Nishina, K., Ostberg, S., Pavlick, R., Peylin, P., Schaphoff, S., Vuichard, N., Warszawski, L., Wiltshire, A., and Woodward, F. I.: Carbon residence time dominates uncertainty in terrestrial vegetation responses to future climate and atmospheric CO_{2}, P. Natl. Acad. Sci. USA, 111, 3280–3285, https://doi.org/10.1073/pnas.1222477110, 2014. a
Frölicher, T. L., Sarmiento, J. L., Paynter, D. J., Dunne, J. P., Krasting, J. P., and Winton, M.: Dominance of the Southern Ocean in anthropogenic carbon and heat uptake in CMIP5 models, J. Climate, 28, 862–886, 2015. a
Giorgetta, M., Jungclaus, J., Reick, C., Legutke, S., Brovkin, V., Crueger, T., Esch, M., Fieg, K., Glushak, K., Gayler, V., Haak, H., Hollweg, H.D., Kinne, S., Kornblueh, L., Matei, D., Mauritsen, T., Mikolajewicz, U., Müller, W., Notz, D., Raddatz, T., Rast, S., Roeckner, E., Salzmann, M., Schmidt, H., Schnur, R., Segschneider, J., Six, K., Stockhause, M., Wegner, J., Widmann, H., Wieners, K.H., Claussen, M., Marotzke, J., and Stevens, B.: CMIP5 simulations of the Max Planck Institute for Meteorology (MPIM) based on the MPIESMLR model: The 1pctCO2 experiment, served by ESGF, World Data Center for Climate (WDCC) at DKRZ [data set], https://doi.org/10.1594/WDCC/CMIP5.MXELc1, 2012. a
Gloor, M., Sarmiento, J. L., and Gruber, N.: What can be learned about carbon cycle climate feedbacks from the CO_{2} airborne fraction?, Atmos. Chem. Phys., 10, 7739–7751, https://doi.org/10.5194/acp1077392010, 2010. a
Goodwin, P., Williams, R. G., Roussenov, V. M., and Katavouta, A.: Climate sensitivity from both physical and carbon cycle feedbacks, Geophys. Res. Lett., 46, 7554–7564, 2019. a, b
Gregory, J. M., Jones, C., Cadule, P., and Friedlingstein, P.: Quantifying carbon cycle feedbacks, J. Climate, 22, 5232–5250, 2009. a, b, c, d, e, f, g, h, i, j, k, l, m, n
Gulev, S., Thorne, S. K., Ahn, P. W., Dentener, J., Domingues, F. J., Gerland, C. M., Gong, D., Kaufman, D. S., Nnamchi, H. C., Quaas, J., Rivera, J. A., Sathyendranath, S., Smith, S. L., Trewin, B., von Schuckmann, K., and Vose, R. S.: Changing state of the climate system, in: Climate change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: MassonDelmotte, V., Zhai, P., Pirani, A., Connors, S., Péan, C., Berger, S., Caud, N., Chen, Y., Goldfarb, L., Gomis, M., Huang, M., Leitzell, K., Lonnoy, E., Matthews, J., Maycock, T., Waterfield, T., Yelekçi, O., Yu, R., and Zhou, B., Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, 287–422, https://doi.org/10.1017/9781009157896.004, 2021. a
Hansen, P. C.: Discrete inverse problems: insight and algorithms, SIAM, 7, 213 pp., ISBN 9780898716962, 2010. a
He, Y., Trumbore, S. E., Torn, M. S., Harden, J. W., Vaughn, L. J. S., Allison, S. D., and Randerson, J. T.: Radiocarbon constraints imply reduced carbon uptake by soils during the 21st century, Science, 353, 1419–1424, https://doi.org/10.1126/science.aad4273, 2016. a
Heimann, M.: Carbon Cycle – Climate Models – Theory, Talk given at the Earth System Research Partnership (ESRP) workshop of the Max Planck Society, Weimar, Germany, 21–23 May, 2014, 2014. a, b, c, d
Hooß, G., Voss, R., Hasselmann, K., MaierReimer, E., and Joos, F.: A nonlinear impulse response model of the coupled carbon cycleclimate system (NICCS), Clim. Dynam., 18, 189–202, 2001. a
IPSL: IPSLCM5ALR model output prepared for CMIP5, served by ESGF, World Data Center for Climate (WDCC) at DKRZ [data set], http://cerawww.dkrz.de/WDCC/CMIP5/Compact.jsp?acronym=IPIL (last access: 10 April 2024), 2011. a, b, c
JAMSTEC, AORI, and NIES: MIROCESM model output prepared for CMIP5 1pctCO2, served by ESGF, World Data Center for Climate (WDCC) at DKRZ [data set], https://doi.org/10.1594/WDCC/CMIP5.MIMEc1, 2015a. a
JAMSTEC, AORI, and NIES: MIROCESM model output prepared for CMIP5 piControl, served by ESGF, World Data Center for Climate (WDCC) at DKRZ [data set], https://doi.org/10.1594/WDCC/CMIP5.MIMEpc, 2015b. a
Jones, C., Robertson, E., Arora, V., Friedlingstein, P., Shevliakova, E., Bopp, L., Brovkin, V., Hajima, T., Kato, E., Kawamiya, M., Liddicoat, S., Lindsay, K., Reick, C., Roelandt, C., Segschneider, J., and Tjiputra, J.: Twentyfirstcentury compatible CO_{2} emissions and airborne fraction simulated by CMIP5 earth system models under four representative concentration pathways, J. Climate, 26, 4398–4413, 2013. a
Jones, C., Hughes, J., Jones, G., Christidis, N., Lott, F., Sellar, A., Webb, M., BodasSalcedo, A., Tsushima, Y., and Martin, G.: HadGEM2ES model output prepared for CMIP5 piControl, served by ESGF, World Data Center for Climate (WDCC) at DKRZ [data set], https://doi.org/10.1594/WDCC/CMIP5.MOGEpc, 2014. a
Jones, C. D. and Friedlingstein, P.: Quantifying processlevel uncertainty contributions to TCRE and Carbon Budgets for meeting Paris Agreement climate targets, Environ. Res. Lett., 15, 7, https://doi.org/10.1088/17489326/ab858a, 2020. a, b, c, d, e, f, g, h
Jones, C., Ciais, P., Davis, S, Friedlingstein, P., Gasser, T, Peters, G., Rogelj, J., van Vuuren, D., Canadell, J., Cowie, A., Jackson, R., Jonas, M., Kriegler, E., Littleton, E., Lowe, J., Milne, J., Shrestha, G., Smith, P., Torvanger, A., and Wiltshire, A.: Simulating the Earth system response to negative emissions, Environ. Res. Lett., 11, 095012, https://doi.org/10.1088/17489326/11/9/095012, 2016. a
Joos, F. and Bruno, M.: Pulse response functions are costefficient tools to model the link between carbon emissions, atmospheric CO_{2} and global warming, Phys. Chem. Earth, 21, 471–476, 1996. a
Joos, F., Roth, R., Fuglestvedt, J. S., Peters, G. P., Enting, I. G., von Bloh, W., Brovkin, V., Burke, E. J., Eby, M., Edwards, N. R., Friedrich, T., Frölicher, T. L., Halloran, P. R., Holden, P. B., Jones, C., Kleinen, T., Mackenzie, F. T., Matsumoto, K., Meinshausen, M., Plattner, G.K., Reisinger, A., Segschneider, J., Shaffer, G., Steinacher, M., Strassmann, K., Tanaka, K., Timmermann, A., and Weaver, A. J.: Carbon dioxide and climate impulse response functions for the computation of greenhouse gas metrics: a multimodel analysis, Atmos. Chem. Phys., 13, 2793–2825, https://doi.org/10.5194/acp1327932013, 2013. a, b
Koven, C. D., Chambers, J. Q., Georgiou, K., Knox, R., NegronJuarez, R., Riley, W. J., Arora, V. K., Brovkin, V., Friedlingstein, P., and Jones, C. D.: Controls on terrestrial carbon feedbacks by productivity versus turnover in the CMIP5 Earth System Models, Biogeosciences, 12, 5211–5228, https://doi.org/10.5194/bg1252112015, 2015. a
Kubo, R., Toda, M., and Hashitsume, N.: Statistical physics II: nonequilibrium statistical mechanics, vol. 31, Springer Science & Business Media, https://doi.org/10.1007/9783642967016, 2012. a
Le Quéré, C., Raupach, M., Canadell, J., Marland, G., Bopp, L., Ciais, P., Conway, T., Doney, S., Feely, R., Foster, P., Friedlingstein, P., Gurney, K., Houghton, R., House, J., Huntingford, C., Levy, P., Lomas, M., Majkut, J., Metzl, N., Ometto, J., Peters, G., Prentice, I., Randerson, J., Running, S., Sarmiento, J., Schuster, U., Sitch, S., Takahashi, T., Viovy, N., van der Werf, G., and Woodward, F.: Trends in the sources and sinks of carbon dioxide, Nat. Geosci., 2, 831–836, 2009. a
Lee, J.Y., Marotzke, J., Bala, G., Cao, L., Corti, S., Dunne, J., Engelbrecht, F., Fischer, E., Fyfe, J., Jones, C., Maycock, A., Mutemi, J., Ndiaye, O., Panickal, S., and Zhou, T.: Future Global Climate: ScenarioBased Projections and NearTerm Information, Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, 553–672, https://doi.org/10.1017/9781009157896.006, 2021. a, b
Lembo, V., Lucarini, V., and Ragone, F.: Beyond forcing scenarios: predicting climate change through response operators in a coupled general circulation model, Sci. Rep., 10, 8668, https://doi.org/10.1038/s41598020652972, 2020. a
Liddicoat, S., Jones, C., and Hughes, J.: HadGEM2ES model output prepared for CMIP5 esmFdbk1, served by ESGF, World Data Center for Climate (WDCC) at DKRZ [data set], https://doi.org/10.1594/WDCC/CMIP5.MOGEe1, 2014a. a
Liddicoat, S., Jones, C., and Hughes, J.: HadGEM2ES model output prepared for CMIP5 esmFixClim1, served by ESGF, World Data Center for Climate (WDCC) at DKRZ [data set], https://doi.org/10.1594/WDCC/CMIP5.MOGEx1, 2014b. a
Lindsay, K.: CESM1BGC model output prepared for CMIP5 1 percent per year increasing CO2, served by ESGF, World Data Center for Climate (WDCC) at DKRZ [data set], https://doi.org/10.1594/WDCC/CMIP5.NFCBc1, 2013a. a
Lindsay, K.: CESM1BGC model output prepared for CMIP5 ESM feedback 1, served by ESGF, World Data Center for Climate (WDCC) at DKRZ [data set], https://doi.org/10.1594/WDCC/CMIP5.NFCBe1, 2013b. a
Lindsay, K.: CESM1BGC model output prepared for CMIP5 ESM fixed climate 1, served by ESGF, World Data Center for Climate (WDCC) at DKRZ [data set], https://doi.org/10.1594/WDCC/CMIP5.NFCBx1, 2013c. a
Lindsay, K.: CESM1BGC model output prepared for CMIP5 preindustrial control, served by ESGF, World Data Center for Climate (WDCC) at DKRZ [data set], https://doi.org/10.1594/WDCC/CMIP5.NFCBpc, 2013d. a
Lucarini, V.: Evidence of dispersion relations for the nonlinear response of the Lorenz 63 system, J. Stat. Phys., 134, 381–400, 2009. a
Lucarini, V., Blender, R., Herbert, C., Ragone, F., Pascale, S., and Wouters, J.: Mathematical and physical ideas for climate science, Rev. Geophys., 52, 809–859, https://doi.org/10.1002/2013RG000446, 2014. a
MacDougall, A. H., Frölicher, T. L., Jones, C. D., Rogelj, J., Matthews, H. D., Zickfeld, K., Arora, V. K., Barrett, N. J., Brovkin, V., Burger, F. A., Eby, M., Eliseev, A. V., Hajima, T., Holden, P. B., JeltschThömmes, A., Koven, C., Mengis, N., Menviel, L., Michou, M., Mokhov, I. I., Oka, A., Schwinger, J., Séférian, R., Shaffer, G., Sokolov, A., Tachiiri, K., Tjiputra , J., Wiltshire, A., and Ziehn, T.: Is there warming in the pipeline? A multimodel analysis of the Zero Emissions Commitment from CO_{2}, Biogeosciences, 17, 2987–3016, https://doi.org/10.5194/bg1729872020, 2020. a
Marotzke, J., Jakob, C., Bony, S. Dirmeyer, P. O'Gorman, P., Hawkins, E., PerkinsKirkpatrick, S., Le Quéré, C., Nowicki, S., Paulavets, K.,Seneviratne, S., Stevens, B., and Tuma, M.: Climate research must sharpen its view, Nat. Clim. Change, 7, 89–91, 2017. a
Mauritsen, T. and Pincus, R.: Committed warming inferred from observations, Nat. Clim. Change, 7, 652–655, https://doi.org/10.1038/nclimate3357, 2017. a
Meehl, G. A., Washington, W. M., Collins, W. D., Arblaster, J. M., Hu, A., Buja, L. E., Strand, W. G., and Teng, H.: How much more global warming and sea level rise?, Science, 307, 1769–1772, 2005. a
Melillo, J., Steudler, P., Aber, J., Newkirk, K., Lux, H., Bowles, F., Catricala, C., Magill, A., Ahrens, T., and Morrisseau, S.: Soil warming and carboncycle feedbacks to the climate system, Science, 298, 2173–2176, 2002. a
Mitchell, T. D.: Pattern scaling: an examination of the accuracy of the technique for describing future climates, Clim. Change, 60, 217–242, 2003. a
Morozov, V. A.: On the solution of functional equations by the method of regularization, in: Doklady Akademii Nauk, Russian Academy of Sciences, 167, 510–512, 1966. a
Murray, J. D.: Mathematical Biology. I. An introduction, Springer, Berlin, Heidelberg, Berlin Heidelberg, 3rd Edn., https://doi.org/10.1007/9783662085424, 1993. a
Myhre, G., Highwood, E. J., Shine, K. P., and Stordal, F.: New estimates of radiative forcing due to well mixed greenhouse gases, Geophys. Res. Lett., 25, 2715–2718, 1998. a
Oeschger, H. and Heimann, M.: Uncertainties of predictions of future atmospheric CO_{2} concentrations, J. Geophys. Res.Oceans, 88, 1258–1262, https://doi.org/10.1029/JC088iC02p01258, 1983. a, b
Phillips, D. L.: A technique for the numerical solution of certain integral equations of the first kind, J. ACM, 9, 84–97, 1962. a
Plattner, G.K. Knutti,R., Joos,F. , Stocker,T. F., von Bloh, W., Brovkin, V. , Cameron, D. , Driesschaert, E. Dutkiewicz, S., Eby, M., Edwards, N., Fichefet, T., Hargreaves, J., Jones, C., Loutre, M., Matthews, H., Mouchet, A. , Müller, S., Nawrath, S., Price, A., Sokolov, A., Strassmann, K., and Weaver, A.: Longterm climate commitments projected with climate–carbon cycle models, J. Climate, 21, 2721–2751, 2008. a
Raich, J. W. and Potter, C. S.: Global patterns of carbon dioxide emissions from soils, Global Biogeochem. Cy., 9, 23–36, 1995. a
Raupach, M. R.: The exponential eigenmodes of the carbonclimate system, and their implications for ratios of responses to forcings, Earth Syst. Dynam., 4, 31–49, https://doi.org/10.5194/esd4312013, 2013. a, b, c, d, e
Raupach, M. R., Canadell, J. G., and Le Quéré, C.: Anthropogenic and biophysical contributions to increasing atmospheric CO_{2} growth rate and airborne fraction, Biogeosciences, 5, 1601–1613, https://doi.org/10.5194/bg516012008, 2008. a
Raupach, M. R., Gloor, M., Sarmiento, J. L., Canadell, J. G., Frölicher, T. L., Gasser, T., Houghton, R. A., Le Quéré, C., and Trudinger, C. M.: The declining uptake rate of atmospheric CO_{2} by land and ocean sinks, Biogeosciences, 11, 3453–3475, https://doi.org/10.5194/bg1134532014, 2014. a
Roe, G.: Feedbacks, timescales, and seeing red, Annu. Rev. Earth Planet. Sci., 37, 93–115, 2009. a, b, c, d, e
Rubino, M., Etheridge, D. M., Trudinger, C. M., Allison, C. E., Rayner, P. J., Enting, I., Mulvaney, R., Steele, L. P., Langenfelds, R. L., Sturges, W. T., Curran, M. A. J., and Smith, A. M.: Low atmospheric CO_{2} levels during the Little Ice Age due to coolinginduced terrestrial uptake, Nat. Geosci., 9, 691–694, 2016. a, b, c, d, e, f, g
Ruelle, D.: Nonequilibrium statistical mechanics near equilibrium: computing higherorder terms, Nonlinearity, 11, 5, https://doi.org/10.1088/09517715/11/1/002, 1998. a
Santos Gutiérrez, M. and Lucarini, V.: On some aspects of the response to stochastic and deterministic forcings, J. Phys. AMath. Gen., 55, 425002, https://doi.org/10.1088/17518121/ac90fd, 2022. a
Schetzen, M.: Nonlinear system modelling and analysis from the Volterra and Wiener perspective, in: Blockoriented Nonlinear System Identification, 13–24, Springer, https://doi.org/10.1007/9781849965132_2, 2010. a
Schwinger, J., Tjiputra, J., Heinze, C., Bopp, L., Christian, J., Gehlen, M., Ilyina, T., Jones, C., SalasMélia, D., Segschneider, J., Séférian, R., and Totterdell, I.: Nonlinearity of ocean carbon cycle feedbacks in CMIP5 earth system models, J. Climate, 27, 3869–3888, 2014. a
Selesnick, I.: Least squares with examples in signal processing, Connexions, 4, 1–25, 2013. a, b
Sundaram, R. K.: A first course in optimization theory, Cambridge university press, https://doi.org/10.1017/CBO9780511804526, 1996. a
Takahashi, T., Sutherland, S., Wanninkhof, R., Sweeney, C., Feely, R. Chipman, D., Hales, B., Friederich, G., Chavez, F., Sabine, C., Watson, A., Bakker, D., Schuster, U., Metzl, N., YoshikawaInoue, H., Ishii, M., Midorikawa, T., Nojiri, Y., Körtzinger, A., Steinhoff, T., Hoppema, M., Olafsson, J., Arnarson, T., Tilbrook, B., Johannessen, T., Olsen, A., Bellerby, R., Wong, C., Delille, B., Bates, N., and de Baar, H.: Climatological mean and decadal change in surface ocean pCO_{2}, and net sea–air CO_{2} flux over the global oceans, DeepSea Res. Pt II, 56, 554–577, 2009. a
Taylor, K. E., Stouffer, R. J., and Meehl, G. A.: An overview of CMIP5 and the experiment design, B. Am. Meteor. Soc., 93, 485–498, 2012. a, b, c
Thornton, P. E., Doney, S. C., Lindsay, K., Moore, J. K., Mahowald, N., Randerson, J. T., Fung, I., Lamarque, J.F., Feddema, J. J., and Lee, Y.H.: Carbonnitrogen interactions regulate climatecarbon cycle feedbacks: results from an atmosphereocean general circulation model, Biogeosciences, 6, 2099–2120, https://doi.org/10.5194/bg620992009, 2009. a
Tikhonov, A. N.: Solution of incorrectly formulated problems and the regularization method, in: Dokl. Akad. Nauk., 151, 1035–1038, 1963. a
Tjiputra, J., Roelandt, C., Bentsen, M., Lawrence, D., Lorentzen, T., Schwinger, J., Seland, Ø., and Heinze, C.: cmip5 output1 NCC NorESM1ME esmFdbk1, served by ESGF, World Data Center for Climate (WDCC) at DKRZ [data set], https://doi.org/10.1594/WDCC/CMIP5.NCCNEe1, 2012a. a
Tjiputra, J., Roelandt, C., Bentsen, M., Lawrence, D., Lorentzen, T., Schwinger, J., Seland, Ø., and Heinze, C.: cmip5 output1 NCC NorESM1ME esmFixClim1, served by ESGF, World Data Center for Climate (WDCC) at DKRZ [data set], https://doi.org/10.1594/WDCC/CMIP5.NCCNEx1, 2012b. a
Tjiputra, J., Roelandt, C., Bentsen, M., Lawrence, D., Schwinger, J., Seland, Ø., and Heinze, C.: CMIP5 output1 NCC NorESM1ME 1pctCO2, served by ESGF, World Data Center for Climate (WDCC) at DKRZ [data set], https://doi.org/10.1594/WDCC/CMIP5.NCCNEc1, 2012c. a
ToddBrown, K. E. O., Randerson, J. T., Post, W. M., Hoffman, F. M., Tarnocai, C., Schuur, E. A. G., and Allison, S. D.: Causes of variation in soil carbon simulations from CMIP5 Earth system models and comparison with observations, Biogeosciences, 10, 1717–1736, https://doi.org/10.5194/bg1017172013, 2013. a
Torres Mendonca, G. L., Reick, C. H., and Pongratz, J.: Supplementary material for “Timescale dependence of airborne fraction and underlying climatecarbon feedbacks for weak perturbations in CMIP5 models”, MPG Publication Repository – MPG. PuRe [data set], https://hdl.handle.net/21.11116/0000000CF6A27 (last access: 11 April 2024), 2023. a, b, c, d
Torres Mendonça, G. L., Pongratz, J., and Reick, C. H.: Identification of linear response functions from arbitrary perturbation experiments in the presence of noise – Part 1: Method development and toy model demonstration, Nonlin. Processes Geophys., 28, 501–532, https://doi.org/10.5194/npg285012021, 2021a. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z, aa, ab, ac, ad, ae, af, ag, ah
Torres Mendonça, G. L., Pongratz, J., and Reick, C. H.: Identification of linear response functions from arbitrary perturbation experiments in the presence of noise – Part 2: Application to the land carbon cycle in the MPI Earth System Model, Nonlin. Processes Geophys., 28, 533–564, https://doi.org/10.5194/npg285332021, 2021b. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z, aa, ab, ac, ad, ae, af, ag, ah, ai, aj, ak, al
Webb, M., Williams, J., Andrews, T., BodasSalcedo, A., Tsushima, Y., Hughes, J., and Jones, C.: HadGEM2ES model output prepared for CMIP5 1pctCO2, served by ESGF, World Data Center for Climate (WDCC) at DKRZ [data set], https://doi.org/10.1594/WDCC/CMIP5.MOGEc1, 2014. a
Wetherald, R. T., Stouffer, R. J., and Dixon, K. W.: Committed warming and its implications for climate change, Geophys. Res. Lett., 28, 1535–1538, 2001. a
Wigley, T. M.: The climate change commitment, Science, 307, 1766–1769, 2005. a
Williams, R. G., Katavouta, A., and Goodwin, P.: Carboncycle feedbacks operating in the climate system, Current Climate Change Reports, 5, 282–295, 2019. a, b
Wu, T. and Xin, X.: bcccsm11 model output prepared for CMIP5 1pctCO2 experiment, served by ESGF, World Data Center for Climate (WDCC) at DKRZ [data set], https://doi.org/10.1594/WDCC/CMIP5.BCB1c1, 2015a. a
Wu, T. and Xin, X.: bcccsm11 model output prepared for CMIP5 esmFdbk1 experiment, served by ESGF, World Data Center for Climate (WDCC) at DKRZ [data set], https://doi.org/10.1594/WDCC/CMIP5.BCB1e1, 2015b. a
Wu, T. and Xin, X.: bcccsm11 model output prepared for CMIP5 esmFixClim1 experiment, served by ESGF, World Data Center for Climate (WDCC) at DKRZ [data set], https://doi.org/10.1594/WDCC/CMIP5.BCB1x1, 2015c. a
Wu, T. and Xin, X.: bcccsm11 model output prepared for CMIP5 piControl experiment, served by ESGF, World Data Center for Climate (WDCC) at DKRZ [data set], https://doi.org/10.1594/WDCC/CMIP5.BCB1pc, 2015d. a
Yan, Y., Zhou, X., Jiang, L., and Luo, Y.: Effects of carbon turnover time on terrestrial ecosystem carbon storage, Biogeosciences, 14, 5441–5454, https://doi.org/10.5194/bg1454412017, 2017. a
Zaehle, S., Friend, A., Friedlingstein, P., Dentener, F., Peylin, P., and Schulz, M.: Carbon and nitrogen cycle dynamics in the OCN land surface model: 2. Role of the nitrogen cycle in the historical terrestrial carbon balance, Global Biogeochem. Cy., 24, GB1006, https://doi.org/10.1029/2009GB003522, 2010. a
Zhou, S., Liang, J., Lu, X., Li, Q., Jiang, L., Zhang, Y., Schwalm, C. R., Fisher, J. B., Tjiputra, J., Sitch, S., Ahlström, A., Huntzinger, D. N., Huang, Y., Wang, G., and Luo, Y.: Sources of Uncertainty in Modeled Land Carbon Storage within and across Three MIPs: Diagnosis with Three New Techniques, J. Climate, 31, 2833–2851, https://doi.org/10.1175/JCLID170357.1, 2018. a
Zickfeld, K., Eby, M., Matthews, H. D., Schmittner, A., and Weaver, A. J.: Nonlinearity of carbon cycle feedbacks, J. Climate, 24, 4255–4275, 2011. a
Note that our reference system is different from that used in Friedlingstein et al. (2003) to quantify the feedbacks. While for our reference system atmospheric CO_{2} is fully determined by CO_{2} emissions without any response in the land and ocean carbon fluxes, their reference system includes in addition the biogeochemical response so that they quantify the radiative feedback alone. Therefore our Eq. (8) differs from their Eq. (9) not only because we investigate a generalization of their framework but also because of the different reference system. This difference is also discussed in Gregory et al. (2009).
This follows by deriving Eq. (14) as described in the main text and then inverting it to obtain $\mathrm{\Delta}{C}_{\text{A}}\left(t\right)={\int}_{\mathrm{0}}^{t}A(ts){\int}_{\mathrm{0}}^{s}E\left({s}^{\prime}\right)\mathrm{d}{s}^{\prime}\mathrm{d}s$. This can be viewed as a generalization of the equation defining the cumulative airborne fraction $\mathrm{\Delta}{C}_{\text{A}}\left(t\right)=:\mathrm{CAF}\left(t\right){\int}_{\mathrm{0}}^{t}E\left(s\right)\mathrm{d}s$ (e.g. Raupach, 2013) by accounting for the effect of the whole past history of cumulative emissions.
We thank the reviewer Ian Enting for making us aware that χ_{ζ}(t) has been widely studied in connection with the socalled “global warming potential” (see e.g. Joos et al., 2013).
Note that for our analysis we calculated $\stackrel{\mathrm{\u0303}}{A}\left(p\right)$ only up to a timescale of 100 years because of the restricted size of the linear regime: $\stackrel{\mathrm{\u0303}}{A}\left(p\right)$ can be predicted only until those timescales where the recovery of the generalized sensitivities is reliable; a reliable recovery, in turn, can only be obtained when the underlying data are unaffected by strong nonlinearities that appear at larger perturbation strengths (Torres Mendonça et al., 2021a). Most restrictive in this respect is χ_{β}(t), whose recovery had for this reason to be based on only the first few decades of available simulation data (70 years for land and 50 years for ocean; see Table A2). Nevertheless, experience with MPIESM in Torres Mendonça et al. (2021b, Fig. 8a and b) and in Appendix A2 (Fig. A1e) suggests that the recovery of generalized sensitivities can be relied upon even for timescales longer than those involved in their recovery – even for the whole 140 years of available data. But since we have not tested this for models other than MPIESM, we decided to consider our recovery of $\stackrel{\mathrm{\u0303}}{A}\left(p\right)$ reliable only up to a 100year timescale.
Note that linear response functions characterize only the ensemble average of the response (see Torres Mendonça et al., 2021a, Sect. 2), so they cannot predict highfrequency temperature variations arising from the internal variability in the system.
 Abstract
 Introduction
 Generalized α–β–γ framework
 Predictive power of the generalized α–β–γ framework
 Timescale dependence of climate–carboncycle feedbacks and airborne fraction for weak perturbations in CMIP5 models
 Conclusions
 Discussion
 Outlook
 Appendix A: Calculation of generalized sensitivities for the MPIESM
 Appendix B: Calculation of the timescaledependent airborne fraction from emissiondriven simulations
 Appendix C: Calculation of χ_{ζ}(t) enforcing the constraint χ_{ζ}(0)=1
 Appendix D: Nonnegative and monotonically decreasing χ_{ζ}(t) implies $\mathrm{0}\le \stackrel{\mathrm{\u0303}}{A}\left(p\right)\le \mathrm{1}$
 Appendix E: Calculation of uncertainty range in α–β–γ sensitivities
 Appendix F: Predicting standard airborne fraction AF(t) from the generalized airborne fraction A(t)
 Appendix G: Transforming sensitivities defined for separate land and ocean temperatures to sensitivities defined for a single global temperature
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Generalized α–β–γ framework
 Predictive power of the generalized α–β–γ framework
 Timescale dependence of climate–carboncycle feedbacks and airborne fraction for weak perturbations in CMIP5 models
 Conclusions
 Discussion
 Outlook
 Appendix A: Calculation of generalized sensitivities for the MPIESM
 Appendix B: Calculation of the timescaledependent airborne fraction from emissiondriven simulations
 Appendix C: Calculation of χ_{ζ}(t) enforcing the constraint χ_{ζ}(0)=1
 Appendix D: Nonnegative and monotonically decreasing χ_{ζ}(t) implies $\mathrm{0}\le \stackrel{\mathrm{\u0303}}{A}\left(p\right)\le \mathrm{1}$
 Appendix E: Calculation of uncertainty range in α–β–γ sensitivities
 Appendix F: Predicting standard airborne fraction AF(t) from the generalized airborne fraction A(t)
 Appendix G: Transforming sensitivities defined for separate land and ocean temperatures to sensitivities defined for a single global temperature
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References