Quantifying memory and persistence in the atmosphere–land and ocean carbon system

Here we intend to further the understanding of the planetary burden (and its dynamics) caused by the effect of the continued increase in carbon dioxide (CO₂) emissions from fossil fuel burning and land use as well as by global warming from a new rheological (stress–strain) perspective. That is, we perceive the emission of anthropogenic CO₂ into the atmosphere as a stressor and survey the condition of Earth in stress–strain units (stress in units of Pa, strain in units of 1) – allowing access to and insight into previously unknown characteristics reflecting Earth's rheological status. We use the idea of a Maxwell body consisting of elastic and damping (viscous) elements to reflect the overall behavior of the atmosphere–land and ocean system in response to the continued increase in CO₂ emissions between 1850 and 2015. Thus, from the standpoint of a global observer, we see that the CO₂ concentration in the atmosphere is increasing (rather quickly). Concomitantly, the atmosphere is warming and expanding, while some of the carbon is being locked away (rather slowly) in land and oceans, likewise under the influence of global warming.

It is not known how reversible and how out of sync the latter process (uptake of carbon by sinks) is in relation to the former (expansion of the atmosphere). All we know is that the slower process remembers the influence of the faster one, which runs ahead. Important questions arise as to whether this global-scale memory – Earth's memory – can be identified and quantified, how it behaves dynamically, and, last but not least, how it interlinks with persistence by which we understand Earth's path dependency.

We go beyond textbook knowledge by introducing three parameters that characterize the system: delay time, memory, and persistence. The three parameters depend, ceteris paribus, solely on the system's characteristic viscoelastic behavior and allow deeper and novel insights into that system. The parameters come with their own limits which govern the behavior of the atmosphere–land and ocean carbon system, independently from any external target values (such as temperature targets justified by means of global change research). We find that since 1850, the atmosphere–land and ocean system has been trapped progressively in terms of persistence (i.e., it will become progressively more difficult to relax the system), while its ability to build up memory has been reduced. The ability of a system to build up memory effectively can be understood as its ability to respond still within its natural regime or, if the build-up of memory is limited, as a measure for system failures globally in the future. Approximately 60 % of Earth's memory had already been exploited by humankind prior to 1959. Based on these stress–strain insights we expect that the atmosphere–land and ocean carbon system will be forced outside its natural regime well before 2050 if the current trend in emissions is not reversed immediately and sustainably.

1 Motivation

Over the last century anthropogenic pressure on Earth became increasingly noticeable. Human activities turned out to be so pervasive and profound that the very life support system upon which humans depend is threatened (Steffen et al., 2004, 2015). The increase in emissions of greenhouse gases (GHGs) into the atmosphere is only one of several serious global threats, and their reduction is at the center of international agreements (Steffen et al., 2015; UN Climate Change, 2022; UN Sustainable Development Goals, 2022).

Here we intend to further the understanding of the planetary burden (and its dynamics) caused by the effect of the continued increase in GHG emissions and by global warming from a new rheological (stress–strain) perspective. That is, we perceive the emission of anthropogenic GHGs, notably carbon (CO₂), into the atmosphere as a stressor. This perspective goes beyond the global carbon mass balance perspective applied by the carbon community, which is widely referred to as the gold standard in assessing whether Earth will remain hospitable for life (Global Carbon Project, 2019). There, the condition of Earth is surveyed in units of petagrams of carbon per year (Pg C yr⁻¹), while we survey its condition in stress–strain units (stress in units of Pa, strain in units of 1) – allowing access to and insight into previously unknown characteristics reflecting Earth's rheological status.

We note that – although the focus is on the atmosphere–land and ocean carbon system – the stress–strain approach described herein should not be considered an appendix to a mass-balance-based carbon cycle model. Instead, it leads to a self-standing model belonging to the suite of reduced but still insightful models (such as radiation transfer, energy balance, or box-type carbon cycle models), which offer great benefits in safeguarding complex three-dimensional climate and global change models. A stress–strain model is missing in that suite of support models. Here we demonstrate the applicability and efficacy of such a model in an Earth systems context.

To develop a stress–strain systems perspective, we begin with the stress given by the CO₂ emissions from fossil fuel burning and land use between 1959 and 2015 (with the increase between 1850 and 1958 serving as antecedent or upstream emissions). Thus, from the standpoint of a global observer, we see that the CO₂ concentration in the atmosphere is increasing (rather quickly). Concomitantly, the atmosphere is warming (here combining the effect of tropospheric warming and stratospheric cooling) and expanding (by approximately 15–20 m in the troposphere per decade since 1990), while some of the carbon is being locked away (rather slowly) in land and oceans, likewise under the influence of global warming (Global Carbon Project, 2019; Lackner et al., 2011; Philipona et al., 2018; Steiner et al., 2011, 2020). We refer to these two processes together, the expansion of the atmosphere and the uptake of carbon by sinks, as the overall strain response of the atmosphere–land and ocean carbon system.

It is not known how reversible and how out of sync the latter process (uptake of carbon by sinks) is in relation to the former (expansion of the atmosphere) (Boucher et al., 2012; Dusza et al., 2020; Garbe et al., 2020; Schwinger and Tjiputra, 2018; Smith, 2012). All we know is that the slower process remembers the influence of the faster one, which runs ahead. Three (nontrivial) questions arise: (1) can this global-scale memory – Earth's memory – be quantified? (2) Can Earth's memory be compared with a buffer which is limited and negligently exploited; that is, what is the degree of depletion? And (3) does Earth's memory allow its persistence (path dependency) to be quantified, speculating that the two are not independent of each other? We answer these questions in the course of our paper.

This suggests, as the next step in developing a stress–strain systems perspective, getting a grip on Earth's memory. To this end, we focus on the slow-to-fast temporal offset inherent in the atmosphere–land and ocean system, while preferring an approach which is reduced to the highest possible extent; however, we do so without compromising complexity in principle. To this end, it is sufficient to resolve subsystems as a whole and to perceive their physical reaction in response to the increase in atmospheric CO₂ concentrations as a combined one (i.e., including effects such as that of global warming). From a temporal perspective, the subsystems' reactions, the expansion of the atmosphere by volume and the sequestration of carbon by sinks, can be considered sufficiently disjunct. Under optimal conditions (referring to the long-term stability of the temporal offset), the temporal offset view even suggests that we can refrain from disentangling the exchange of both thermal energy and carbon throughout the atmosphere–land and ocean system, as is done in climate–carbon models ranging from reduced to complex (Flato et al., 2013; Harman and Trudinger, 2014). The additional degree of reductionism, whilst preserving complexity, will prove an advantage in advancing our understanding of the temporal offset in terms of memory and persistence.

In view of the aforementioned questions, we chose a rheological stress–strain (σ–ε) model (Roylance, 2001; TU Delft, 2022): a Maxwell body (MB) consisting of an elastic element (its constant, traditionally denoted E for Young's modulus, is replaced by the compression modulus K) and a damping (viscous) element (the damping constant is denoted D) to capture the stress–strain behavior of the global atmosphere–land and ocean system (Fig. 1) and to simulate how humankind propelled that global-scale experiment historically. We note that the MB is a logical choice of model given the uninterrupted increase in atmospheric CO₂ concentrations since 1850 (Global Carbon Project, 2019).

Figure 1 Rheological model to capture the stress–strain behavior of the global atmosphere–land and ocean system as a Maxwell body consisting of elastic (atmosphere) and damping (land–ocean) elements. The stress (in units of Pa; known) is given by the carbon (CO₂) emissions from fossil fuel burning and land use, while the strain (in units of 1; assumed exponential, otherwise unknown) is given by the expansion of the atmosphere by volume and uptake of CO₂ by sinks. Independent estimates of K and D, the compression and damping characteristics of the MB, allow its stress–strain behavior to be captured and adjusted until consistency is achieved (see text).

In practice, rheology is principally concerned with extending continuum mechanics to characterize the flow of materials that exhibit a combination of elastic, viscous, and plastic behavior (that is, including hereditary behavior) by properly combining elasticity and (Newtonian) fluid mechanics. Limits (e.g., viscosity limits) exist beyond which basic rheological models are recommended to be refined. However, these limits are fluent, and basic rheological models also produce useful results beyond these limits (Malkin and Isayev, 2017; Mezger, 2006; TU Delft, 2022).

The mathematical treatment of an MB is standard. Depending on whether the strain (ε) or the stress (σ) is known (in addition to the compression and damping characteristics K and D), the stress–strain equation describing the MB between 0 and t can be applied in a stress-explicit form,

$$\begin{array}{}\text{(1a)}& \mathit{\sigma }\left(t\right)=\mathit{\sigma }\left(\mathrm{0}\right)\exp \left(-{\displaystyle \frac{K}{D}}t\right)+K\underset{\mathrm{0}}{\overset{t}{\int }}\dot{\mathit{\epsilon }}\left(\mathit{\tau }\right)\exp \left({\displaystyle \frac{K}{D}}(\mathit{\tau }-t)\right)\mathrm{d}\mathit{\tau },\end{array}$$

or in a strain-explicit form,

$$\begin{array}{}\text{(1b)}& \mathit{\epsilon }\left(t\right)=\mathit{\epsilon }\left(\mathrm{0}\right)+{\displaystyle \frac{\mathrm{1}}{K}}\left[\mathit{\sigma }\right(t)-\mathit{\sigma }(\mathrm{0}\left)\right]+{\displaystyle \frac{\mathrm{1}}{D}}\underset{\mathrm{0}}{\overset{t}{\int }}\mathit{\sigma }\left(\mathit{\tau }\right)\mathrm{d}\mathit{\tau },\end{array}$$

with σ(0) and ε(0) denoting initial conditions and a dot the derivative by time (Roylance, 2001; Bertram and Glüge, 2015).

Here, we focus on the application of these equations in an atmosphere–land and ocean carbon context. For an observer it is the overall strain response of that system (expansion of the atmosphere by volume and uptake of CO₂ by sinks) that is unknown. However, since atmospheric CO₂ concentrations have been observed to increase exponentially (quasi-continuously), the strain can be expected to be exponential or close to exponential. In addition, we provide independent estimates of the likewise unknown compression and damping characteristics of the MB. This a priori knowledge allows Eqs. (1a) and (1b) to be used stepwise in combination to narrow down our initial estimate of the $K/D$ ratio in particular. More accurate knowledge of this ratio is needed when we go beyond textbook knowledge by distilling three parameters – delay time (reflecting the temporal offset mentioned above), memory, and persistence – from the stress-explicit equation. The three parameters depend, ceteris paribus, solely on the system's characteristic $K/D$ ratio and allow deeper and novel insights into that system. We see the atmosphere–land and ocean system as being trapped progressively over time in terms of persistence. Given its reduced ability to build up memory, we expect system failures globally well before 2050 if the current trend in emissions is not reversed immediately and sustainably. Put differently, the stress–strain approach comes with its own internal limits which govern the behavior of the atmosphere–land and ocean carbon system independently from any external target values (such as temperature targets justified by means of global change research).

There is a wide range of other approaches which aim at exploring memory and persistence in Earth systems data, typically with the focus on individual Earth subsystems or processes (e.g., atmospheric temperature or carbon dioxide emissions). So far, applied approaches have mainly been based on classical time series and time–space analyses to uncover the memory or causal patterns contained in observational data (Barros et al., 2016; Belbute and Pereira, 2017; Caballero et al., 2002; Franzke, 2010; Lüdecke et al., 2013). However, these approaches come with well-known limitations which can all be attributed, directly or indirectly, to the issue of forecasting (more precisely, the conditions placed on the data to enable forecasting) or are not based on physics (Aghabozorgi et al., 2015; Darlington, 1996; Darlington and Hayes, 2016). By way of contrast, we do not forecast. We perpetuate long-term historical conditions which, in turn, allows the delay time in the atmosphere–land and ocean system to be expressed analytically in terms of memory and persistence. We are not aware of any scientific discipline or research area in which memory and persistence are defined other than statistically and are interlinked, if at all, other than via correlation.

Rheological approaches are common in Earth systems modeling as well. Typically, they are applied to mimic the long(er)-term behavior of Earth subsystems, e.g., its mantle viscosity, which is crucial for interpreting glacial uplift resulting from changes in planetary ice sheet loads (Müller, 1986; Whitehouse et al., 2019; Yuen et al., 1986). Yet, to the best of our knowledge, a rheological approach to unravel the memory and persistence behavior of the global atmosphere–land and ocean system in response to the long-lasting increase in atmospheric CO₂ emissions had not been applied before.

We describe our rheological model (MB) approach in detail in Sect. 2, while we provide an overview of the applied data and conversion factors in Sect. 3. In Sect. 4 we describe how we derive first-order estimates of the main characteristics of the atmosphere–land and ocean system (in terms of the MB's K and D characteristics) by using available knowledge. Although uncertain, these estimates become useful in Sect. 5 where we apply the aforementioned stress- and strain-explicit equations to quantify delay time, memory, and persistence of the atmosphere–land and ocean system. We conclude by taking account of our main findings in Sect. 6.

2 Method

This section provides an overview of how we process Eq. (1a) and how we distill delay time, memory, and persistence from this equation. To familiarize oneself with the details, the reader is referred to the Supplement.

To start with, we assume that we know the order of magnitude of both the $K/D$ ratio characteristic of the atmosphere–land and ocean system and the rate of change in the strain ε given by $\dot{\mathit{\epsilon }}\left(t\right)=\mathit{\alpha }\exp \left(\mathit{\alpha }t\right)$ with the exponential growth factor α>0. These first-order estimates permit Eqs. (1a) and (1b) to be used stepwise in combination.

Equation (1a): we vary both $K/D$ and α to reproduce the known stress σ given by the CO₂ emissions from fossil fuel burning (fairly well known) and land use (less known) (Global Carbon Project, 2019).

Equation (1b): we insert both the fine-tuned $K/D$ ratio and the known stress σ to compute the strain ε and check its derivative by time.

We consider this procedure a check of consistency, not a proof of concept.

Delay time, memory, and persistence are characteristic (functions) of the MB. They are contained in the integral on the right side of Eq. (1a) and are defined independently of initial conditions. These appear only in the lower boundary of that integral, which allows initial conditions other than zero to be considered by taking advantage of the integral's additivity. Thus, without loss of generality, we rewrite Eq. (1a) for σ(0)=0, which results in

$$\begin{array}{}\text{(2a)}& \mathit{\sigma }\left(t\right)={\displaystyle \frac{D}{\mathit{\beta }}}\dot{\mathit{\epsilon }}\left(t\right)\left(\mathrm{1}-q_{\mathit{\beta }}^t\right),\end{array}$$

(see Supplement Information 1) where $\mathit{\beta }=\mathrm{1}+\frac{D}{K}\mathit{\alpha }$ and $q_{\mathit{\beta }}^t=\exp \left(-\frac{K}{D}\mathit{\beta }t\right)$. The term $\frac{D}{K\mathit{\beta }}$ represents a time characteristic of the MB under (here) exponential strain (i.e., of the MB that responds to the stress acting upon it), whereas $\frac{D}{K}$ is the relaxation time of the MB (i.e., of the MB that relaxes unhindered after the stress causing that strain has vanished or that responds to strain held constant over time; also known as the relaxation test; Bertram and Glüge, 2015). However, to ensure that exponents still come in units of 1 after we split them up, we introduce the dimensionless time $n=\frac{t}{\mathrm{\Delta }t}$ globally (which will be discretized in the sequel when we refer to a temporal resolution of 1 year and set Δt=1 yr) such that, for example, $q^t=\exp {\left(-\frac{K}{D}\mathrm{\Delta }t\right)}^n$.

To understand the systemic nature of the MB, we explore its stress dependence on $q=\exp \left(-\frac{K}{D}\mathrm{\Delta }t\right)$, which contains the ratio of K and D, the two characteristic parameters of the MB, by way of derivation by q (while α is held constant). To this end, we transform Eq. (2a) further to

$$\begin{array}{}\text{(2b)}& {\mathit{\sigma }}_{\mathrm{D}}(q,t):={\displaystyle \frac{\mathrm{1}}{D}}\mathit{\sigma }\left(t\right)={\displaystyle \frac{\mathrm{1}}{D}}\mathit{\sigma }\left(n\right)=:{\mathit{\sigma }}_{\mathrm{D}}(q,n)\end{array}$$

and execute $\frac{\partial }{\partial q}{\mathit{\sigma }}_{\mathrm{D}}(q,n)$, the derivation by q of the system's rate of change σ_D (which is given in units of yr⁻¹). Doing so allows (what we call) delay time T to be distilled (see Supplement Information 2). It is defined as

$$\begin{array}{}\text{(3)}& T(q,n):={\displaystyle \frac{q_{\mathit{\beta }}}{S_n}}{\displaystyle \frac{\partial S_n}{\partial q_{\mathit{\beta }}}}=-{\displaystyle \frac{q_{\mathit{\beta }}^n}{\mathrm{1}-q_{\mathit{\beta }}^n}}n+{\displaystyle \frac{q_{\mathit{\beta }}}{\mathrm{1}-q_{\mathit{\beta }}}},\end{array}$$

where q_β=q_αq, $q_{\mathit{\alpha }}=\exp (-\mathit{\alpha }\mathrm{\Delta }t)$, and $S_n=S(q,n)=\frac{\mathrm{1}-q_{\mathit{\beta }}^n}{\mathrm{1}-q_{\mathit{\beta }}}$. The delay time behaves asymptotically for increasing n and approaches $T_{\mathrm{\infty }}=\underset{n\to \mathrm{\infty }}{lim}T=\frac{q_{\mathit{\beta }}}{\mathrm{1}-q_{\mathit{\beta }}}$. We further define

$$\begin{array}{}\text{(4)}& M:=S(q,n),\end{array}$$

with $M_{\mathrm{\infty }}=\frac{\mathrm{1}}{\mathrm{1}-q_{\mathit{\beta }}}$, and

$$\begin{array}{}\text{(5)}& P:=T(q,n{)}^{-\mathrm{1}},\end{array}$$

with $P_{\mathrm{\infty }}=\frac{\mathrm{1}}{T_{\mathrm{\infty }}}=\frac{\mathrm{1}-q_{\mathit{\beta }}}{q_{\mathit{\beta }}}$ as the MB's characteristic memory and persistence, respectively. As is commonly done, we keep the list of independent parameters minimal. (We only allow K and D – i.e., q – in addition to n; see Eqs. 2b and 3–5 in particular.)

T as given by Eq. (3) is not simply characteristic of the MB described by Eqs. (2a) and (2b); it can be shown to appear as delay time in the argument of any function dependent on current and previous times, with a weighting decreasing exponentially backward in time (see Supplement Information 3). Equation (4) reflects the history the MB was exposed to systemically prior to current time n (during which α was constant; see Supplement Information 4). Put simply, M can be understood as the depreciated (q-weighted) strain summed up backward in time. Equation (5) can be shortened to $T\cdot P=\mathrm{1}$. If we assume that q can be changed in retrospect at n=0, this equation tells us that if T – that is, ΔM per Δq (or, likewise, $\mathrm{\Delta }M/M$ per $\mathrm{\Delta }q/q$; see the first part of Eq. 3) – is small, P is great because the change in the system's characteristics (contained in q) hardly influences the MB's past, with the consequence that the past exhibits a great path dependency, and vice versa. We therefore perceive persistence and path dependency as synonymous.

An additional quantity to monitor is ln (M⋅P), which approaches ${\mathit{\lambda }}_{\mathit{\beta }}=\mathit{\lambda }\cdot \mathit{\beta }$ for increasing n, with $\mathit{\lambda }=\frac{K}{D}\mathrm{\Delta }t$ as the characteristic rate of change in the MB. The ratio $\mathit{\lambda }/\ln (M\cdot P)$ allows monitoring of how much the system's natural rate of change is exceeded as a consequence of the continued increase in stress (see Supplement Information 5).

3 Data and conversion factors

A detailed overview of the carbon data and conversion factors used in this paper (and also by the carbon community) is given in Supplement Information 6. The data pertain to atmosphere, land, and oceans:

• atmospheric CO₂ concentration (in ppm);

• CO₂ emissions from fossil fuel combustion and cement production (in Pg C yr⁻¹);

• land use change emissions (in Pg C yr⁻¹);

• net primary production (in Pg C yr⁻¹); and

• dissolved organic carbon (in µmol kg⁻¹).

They are given by source and time range and are also described briefly. The context within which they are used is revealed in each of the following sections. The conversion factors are standard; they are needed to convert C to CO₂ and parts per million by volume (ppmv) CO₂ to petagrams of carbon (Pg C) or pascal (Pa).

4 Independent estimates of D and K

In this section we provide independent estimates of the damping and compression characteristics of the atmosphere–land and ocean system, with D_L and D_O denoting the damping constants assigned to land and oceans, respectively, and K denoting the compression modulus assigned to the atmosphere. We capture the characteristics' right order of magnitude only – which can be done on physical grounds by evaluating the combined (net) strain response of each subsystem on the grounds of increasing CO₂ concentrations in the atmosphere. These first-order estimates are adequate as they allow sufficient flexibility for Sect. 5, where we narrow down our initial estimates by using Eqs. (1a) and (1b) stepwise in combination to achieve consistency.

4.1 Estimating the damping constant D_L

Increasing concentrations of CO₂ in the atmosphere trigger the uptake of carbon by the terrestrial biosphere. The intricacies of this process, including potential (positive and negative) feedback processes, are widely discussed (Dusza et al., 2020; Heimann and Reichstein, 2008; Smith, 2012). The crucial question is how we have observed the process of carbon uptake by the terrestrial biosphere taking place in the past. Compared to the reaction of the atmosphere to global warming (an expansion of the atmosphere by volume), we consider this process to be long(er) term in nature and perceive it as a Newton-like (damping) element.

Biospheric carbon uptake is described by the biotic growth factor:

$$\begin{array}{}\text{(6)}& {\mathit{\beta }}_{\mathrm{b}}={\displaystyle \frac{\mathrm{\Delta }\mathrm{NPP}/\mathrm{NPP}}{\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}}/{\mathrm{CO}}_{\mathrm{2}}}},\end{array}$$

which is used to approximate the fractional increase in net primary productivity (NPP) per unit increase in atmospheric CO₂ concentration (Wullschleger et al., 1995; Amthor and Koch, 1996; Luo and Mooney, 1996). Here we make use of the model-derived NPP time series (1900–2016) provided by O'Sullivan et al. (2019) to calculate β_b (O'Sullivan et al., 2019). To understand the uncertainty range underlying β_b for 1959–2018, we use the photosynthetic beta factor:

$$\begin{array}{}\text{(7)}& {\mathit{\beta }}_{\mathrm{Ph}}={\mathrm{CO}}_{\mathrm{2}}L=\left({\displaystyle \frac{\mathrm{dPh}}{\mathrm{Ph}}}\right)\left({\displaystyle \frac{{\mathrm{CO}}_{\mathrm{2}}}{\mathrm{d}{\mathrm{CO}}_{\mathrm{2}}}}\right),\end{array}$$

where L is the so-called leaf-level factor denoting the relative leaf photosynthetic response to a 1 ppmv change in the atmospheric concentration of CO₂, bounded by

$$\begin{array}{}\text{(8)}& L_{\mathrm{1}}\le L=f\left({\mathrm{CO}}_{\mathrm{2}}\right)\le L_{\mathrm{2}},\end{array}$$

(see below), and Ph is the global photosynthetic carbon influx (i.e., gross primary productivity). Equation (7) is similar to Eq. (6). In Eq. (6) β_b represents biomass production changes in response to CO₂ changes, whereas in Eq. (7) β_Ph describes photosynthesis changes in response to CO₂ changes (Luo and Mooney, 1996).

L can be shown to be independent of plant characteristics, light, and the nutrient environment and to vary little by geographic location or canopy position. Thus, L is virtually a constant across ecosystems and a function of time-associated changes in atmospheric CO₂ only (Luo and Mooney, 1996).

We use Eq. (7) to test whether β_b falls within the range of β_Ph given by the quantifiable photosynthetic limits L₁ (photosynthesis limited by electron transport) and L₂ (photosynthesis limited by RuBisCO activity). Figure 2 shows the biotic growth factors from O'Sullivan et al. (2019) that consider changes in NPP due to the combined effect of CO₂ fertilization, nitrogen deposition, climate change, and carbon–nitrogen synergy (β_{NPP_comb}) as well as due to CO₂ fertilization (${\mathit{\beta }}_{\mathrm{NPP}\mathrm{\_}{\mathrm{CO}}_{\mathrm{2}}}$) only. For 1960–2016, β_{NPP_comb} falls between β₁:=β_Ph(L₁) and β₂:=β_Ph(L₂), closer to β₁ than to β₂, whereas ${\mathit{\beta }}_{\mathrm{NPP}\mathrm{\_}{\mathrm{CO}}_{\mathrm{2}}}$ even falls below the lower β₁ limit.

Figure 2 Using the lower (β₁) and upper (β₂) limits of the photosynthetic beta factor to test the range of the biotic growth factor (β_b) for 1960–2016. The biotic growth factor is derived with the help of modeled net primary production (NPP) values accounting for CO₂ fertilization, nitrogen deposition, climate change, and carbon–nitrogen synergy. ${\mathit{\beta }}_{\mathrm{NPP}\mathrm{\_}{\mathrm{CO}}_{\mathrm{2}}}$ refers to O'Sullivan et al. (2019), who consider the change in NPP due to CO₂ fertilization only, and β_{NPP_comb} refers to the change in NPP due to the combined effect. All beta factors are in units of 1.

Rewriting Eq. (7) in the form

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{\mathrm{\Delta }{\mathrm{Ph}}_i}{\mathrm{Ph}}}=L_i\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}}(i=\mathrm{1},\mathrm{2})\end{array}$$

with Ph=120 Pg C yr⁻¹ indicates that the additional amount of annual relative photosynthetic carbon influx, stimulated by a yearly increase in atmospheric CO₂ concentration, can be estimated by L_i or the sequence of L_i if ΔCO₂ spans multiple years (see Supplement Information 7 and Supplement Data 1). Plotting $\mathrm{\Delta }{\mathrm{Ph}}_i/\mathrm{Ph}$ against time allows lower and upper slopes (rates of strain)

$$\begin{array}{}\text{(10a)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left({\displaystyle \frac{\mathrm{\Delta }{\mathrm{Ph}}_{\mathrm{1}}}{\mathrm{Ph}}}\right)\approx \mathrm{0.0019}{\mathrm{yr}}^{-\mathrm{1}}\end{array}$$

and

$$\begin{array}{}\text{(10b)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left({\displaystyle \frac{\mathrm{\Delta }{\mathrm{Ph}}_{\mathrm{2}}}{\mathrm{Ph}}}\right)\approx \mathrm{0.0041}{\mathrm{yr}}^{-\mathrm{1}}\end{array}$$

to be derived for 1959–2018. A linear fit works well in either case. The cumulative increase in atmospheric CO₂ concentration since 1959, $\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}}={\mathrm{CO}}_{\mathrm{2}}\left(t\right)-{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{1959}\right)$, exhibits a moderate exponential (close to linear) trend. Thus, plotting annual changes in CO₂, normalized on the aforementioned rates of strain, versus time allows the remaining (moderate) trends to be interpreted alternatively, namely as average photosynthetic damping constants with appropriate uncertainty given by half the maximal range (see Fig. 3 and Supplement Data 1).

$$\begin{array}{}\text{(11a)}& {\displaystyle }\begin{array}{rl}{D}_{\mathrm{1}}& \approx (\mathrm{815}\pm \mathrm{433})\mathrm{ppmv}\mathrm{yr}=(\mathrm{83}\pm \mathrm{44})\mathrm{Pa}\mathrm{yr}\\ & =(\mathrm{2606}\pm \mathrm{1383}){\mathrm{10}}^{\mathrm{6}}\mathrm{Pa}\mathrm{s}\end{array}\text{(11b)}& {\displaystyle }\begin{array}{rl}{D}_{\mathrm{2}}& \approx (\mathrm{378}\pm \mathrm{201})\mathrm{ppmv}\mathrm{yr}=(\mathrm{38}\pm \mathrm{20})\mathrm{Pa}\mathrm{yr}\\ & =(\mathrm{1207}\pm \mathrm{641}){\mathrm{10}}^{\mathrm{6}}\mathrm{Pa}\mathrm{s}.\end{array}\end{array}$$

Figure 3 Terrestrial carbon uptake perceived as damping (in ppmv yr) based on the limits of leaf photosynthesis (1960–2018: D₁ and D₂) and on model-derived changes in net primary production (NPP; 1960–2016) due to the combined effect of CO₂ fertilization, nitrogen deposition, climate change, and carbon–nitrogen synergy (D_{NPP_comb}) as well as CO₂ fertilization only ($D_{\mathrm{NPP}\mathrm{\_}{\mathrm{CO}}_{\mathrm{2}}}$). The linear trends of the four damping series are shown at the top. These are used to interpret damping as constants with appropriate uncertainty (given by half the maximal range).

Repeating the same procedure for 1959–2016 with O'Sullivan et al. (2019) model-derived NPP values considering the change in NPP due to CO₂ fertilization as well as the total change in NPP, we find

$$\begin{array}{}\text{(12a)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{\left({\displaystyle \frac{\mathrm{\Delta }\mathrm{NPP}}{\mathrm{NPP}}}\right)}_{{\mathrm{CO}}_{\mathrm{2}}}\approx \mathrm{0.0013}{\mathrm{yr}}^{-\mathrm{1}}\end{array}$$

and

$$\begin{array}{}\text{(12b)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{\left({\displaystyle \frac{\mathrm{\Delta }\mathrm{NPP}}{\mathrm{NPP}}}\right)}_{\mathrm{comb}}\approx \mathrm{0.0021}{\mathrm{yr}}^{-\mathrm{1}}\end{array}$$

(linear fits still work well), and consequently

$$\begin{array}{}\text{(13a)}& {\displaystyle }\begin{array}{rl}{D}_{{\mathrm{CO}}_{\mathrm{2}}}& \approx (\mathrm{1172}\pm \mathrm{617})\mathrm{ppmv}\mathrm{yr}=(\mathrm{119}\pm \mathrm{62})\mathrm{Pa}\mathrm{yr}\\ & =(\mathrm{3746}\pm \mathrm{1971}){\mathrm{10}}^{\mathrm{6}}\mathrm{Pa}\mathrm{s},\end{array}\text{(13b)}& {\displaystyle }\begin{array}{rl}{D}_{\mathrm{comb}}& \approx (\mathrm{726}\pm \mathrm{382})\mathrm{ppmv}\mathrm{yr}=(\mathrm{74}\pm \mathrm{39})\mathrm{Pa}\mathrm{yr}\\ & =(\mathrm{2319}\pm \mathrm{1220}){\mathrm{10}}^{\mathrm{6}}\mathrm{Pa}\mathrm{s}.\end{array}\end{array}$$

As before, these estimates are closer to the lower leaf-level factor (higher photosynthetic D) than to the higher leaf-level factor (lower photosynthetic D; Fig. 3).

Here we interpret the O'Sullivan et al. (2019) Earth systems model as a typical one, which means that the NPP changes it produces are common. We therefore (and sufficient for our purposes) choose the damping constant D₁ as a good estimator in light of the total change in NPP of the terrestrial biosphere since 1960. Hence,

$$\begin{array}{}\text{(14)}& \begin{array}{rl}{D}_{\mathrm{L}}& \approx (\mathrm{815}\pm \mathrm{433})\mathrm{ppmv}\mathrm{yr}=(\mathrm{83}\pm \mathrm{44})\mathrm{Pa}\mathrm{yr}\\ & =(\mathrm{2606}\pm \mathrm{1383}){\mathrm{10}}^{\mathrm{6}}\mathrm{Pa}\mathrm{s}.\end{array}\end{array}$$

D_L is on the order of viscosity indicated for bitumen and asphalt (Mezger, 2006).

4.2 Estimating the damping constant D_O

Increasing concentrations of CO₂ in the atmosphere trigger the uptake of carbon by the oceans (National Oceanic and Atmospheric Administration, 2015). Like the uptake of carbon by the terrestrial biosphere, we consider this process to behave like a Newton (damping) element in our MB because of the de facto irreversibility on the shorter timescale we are interested in (Schwinger and Tjiputra, 2018).

The Revelle (buffer) factor (R) quantifies how much atmospheric CO₂ can be absorbed by homogeneous reaction with seawater. R is defined as the fractional change in CO₂ relative to the fractional change in dissolved inorganic carbon (DIC):

$$\begin{array}{}\text{(15)}& R={\displaystyle \frac{\mathrm{\Delta }p{\mathrm{CO}}_{\mathrm{2}}/p{\mathrm{CO}}_{\mathrm{2}}}{\mathrm{\Delta }\mathrm{DIC}/\mathrm{DIC}}}.\end{array}$$

(Here, in contrast to before, atmospheric CO₂ is referred to in units of µatm and therefore indicated by pCO₂.) An R value of 10 indicates that a 10 % change in atmospheric CO₂ is required to produce a 1 % change in the total CO₂ content of seawater (Bates et al., 2014; Egleston et al., 2010; Emerson and Hedges, 2008).

DIC and R have been observed at seven ocean carbon time series sites for periods from 15 to 30 years (between 1983 and 2012) to change slowly and linearly with time (Bates et al., 2014):

$$\begin{array}{}\text{(16)}& {\displaystyle }& {\displaystyle \frac{\mathrm{\Delta }\mathrm{DIC}}{\mathrm{\Delta }t}}\approx [\mathrm{0.8};\mathrm{1.9}]\mathrm{\mu }\mathrm{mol}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}},\text{(17)}& {\displaystyle }& {\displaystyle \frac{\mathrm{\Delta }R}{\mathrm{\Delta }t}}\approx [\mathrm{0.01};\mathrm{0.03}]{\mathrm{yr}}^{-\mathrm{1}}\end{array}$$

(see also Supplement Data 2). Here it is sufficient to proceed with spatiotemporal averages. As before, the cumulative increase in atmospheric CO₂ concentration since 1983, $\mathrm{\Delta }p{\mathrm{CO}}_{\mathrm{2}}=p{\mathrm{CO}}_{\mathrm{2}}\left(t\right)-p{\mathrm{CO}}_{\mathrm{2}}\left(\mathrm{1983}\right)$, exhibits a moderate exponential (close to linear) trend. Thus, plotting annual changes in pCO₂, normalized on the rates of strain $\frac{\left(\mathrm{\Delta }\mathrm{DIC}/\mathrm{DIC}\right)}{\mathrm{\Delta }t}$, versus time allows the remaining (moderate) trend to be interpreted alternatively, namely as an average oceanic damping constant with appropriate uncertainty given by half the maximal range (see Fig. 4 and Supplement Data 2):

$$\begin{array}{}\text{(18)}& \begin{array}{rl}{D}_{\mathrm{O}}& \approx (\mathrm{3005}\pm \mathrm{588})\mathrm{ppmv}\mathrm{yr}=(\mathrm{304}\pm \mathrm{60})\mathrm{Pa}\mathrm{yr}\\ & =(\mathrm{9602}\pm \mathrm{1877}){\mathrm{10}}^{\mathrm{6}}\mathrm{Pa}\mathrm{s}.\end{array}\end{array}$$

D_O is on the order of viscosity indicated for bitumen and asphalt, yet approximately 3.7 times greater than D_L.

Figure 4 Oceanic carbon uptake perceived as damping (in ppmv yr) based on observations at seven ocean carbon time series sites for periods from 15 to 30 years (between 1983 and 2012). The linear trend in oceanic damping, shown at the bottom, is used to interpret damping as a constant with appropriate uncertainty (given by half the maximal range).

4.3 Estimating the compression modulus K

The long-lasting increase in GHG emissions has caused the CO₂ concentration in the atmosphere to increase and the atmosphere as a whole to warm (with tropospheric warming outstripping stratospheric cooling) and to expand (in the troposphere by approximately 15–20 m per decade since 1990) (Global Carbon Project, 2019; Lackner et al., 2011; Philipona et al., 2018; Steiner et al., 2011, 2020). Our whole-subsystem (net-warming) view does not invalidate the known facts that CO₂ in the atmosphere is well-mixed (except for very low altitudes at which deviations from uniform CO₂ concentrations are caused by the dynamics of carbon sources and sinks) and that the volume percentage of CO₂ in the atmosphere stays almost constant up to high altitudes (Abshire et al., 2010; Emmert et al., 2012).

Compared to the slow uptake of carbon by land and oceans, we assume the atmosphere to be represented well by a Hooke element in the MB and this to serve as a (sufficiently stable) surrogate physical descriptor for the reaction of the atmosphere as a whole (Sakazaki and Hamilton, 2020). However, in the case of a gas, Young's modulus E must be replaced by the compression modulus K, the reciprocal of which is compressibility κ. Both K and κ scale with altitude, which we get to grips with in the following. Compressibility is defined by

$$\begin{array}{}\text{(19)}& \mathit{\kappa }={\displaystyle \frac{\mathrm{1}}{K}}=-{\displaystyle \frac{\mathrm{1}}{V}}{\displaystyle \frac{\mathrm{d}V}{\mathrm{d}p}}\end{array}$$

(κ>0) (OpenStax, 2020). Depending on whether the compression happens under isothermal or adiabatic conditions, the compressibility is distinguished accordingly. It is defined by

$$\begin{array}{}\text{(20a)}& {\mathit{\kappa }}_{\mathrm{it}}={\displaystyle \frac{\mathrm{1}}{p}}\end{array}$$

in the isothermal case and

$$\begin{array}{}\text{(20b)}& {\mathit{\kappa }}_{\mathrm{ad}}={\displaystyle \frac{\mathrm{1}}{\mathit{\gamma }p}}\end{array}$$

in the dry adiabatic case, in which γ is the isentropic coefficient of expansion. Its value is 1.403 for dry air (1.310 for CO₂) under standard temperature (273.15 K) and pressure (1 atm; 101.325 kPa) (Wark, 1983). We consider a carbon-enriched atmosphere to also be air.

However, the observed expansion of the troposphere happens neither isothermally nor dry adiabatically but polytropically. Moreover, our ignorance of the exact value of κ is overshadowed by the uncertainty in altitude – or top of the atmosphere (TOA) – which we need as a reference for κ (thus K). As a matter of fact, there is considerable confusion as to which altitude the TOA refers to in climate models (CarbonBrief, 2018; NASA Earth Observatory, 2006).

To advance, we refer to the (dry adiabatic) standard atmosphere, which assigns a temperature gradient of −6.5 ^∘C per 1000 m up to the tropopause at 11 km, a constant value of −56.5 ^∘C (216.65 K) above 11 km and up to 20 km, and other gradients and constant values above 20 km (Cavcar, 2000; Mohanakumar, 2008). Guided by the distribution of atmospheric mass by altitude, we choose the stratopause as our TOA (at about 48 km of altitude and 1 hPa), with uncertainty ranging from middle to higher stratosphere (at about 43 km of altitude and 1.9 hPa) to mid-mesosphere (at about 65 km of altitude and 0.1 hPa) (Digital Dutch, 1999; International Organization for Standardization, 1975; Mohanakumar, 2008; Zellner, 2011). We assign the resulting uncertainty of 90 % in relative terms to

$$\begin{array}{}\text{(21)}& K=(\mathrm{1}\pm \mathrm{0.9})\mathrm{hPa}=(\mathrm{100}\pm \mathrm{90})\mathrm{Pa},\end{array}$$

which we consider sufficiently large to compensate for the unknown isentropic coefficient in the first place; that is,

$$\left[K_{\mathrm{ad},\min };K_{\mathrm{ad},\max }\right]\in \left[K_{\mathrm{it},\min };K_{\mathrm{ad},\max }\right]\in \left[K_{\min };K_{\max }\right].$$

For comparison, K_ad ranges from 400 to 412 hPa were the TOA allocated within the troposphere (exhibiting, for the reference used here, an expansion of 20 m; see Supplement Information 8).

5 Main findings

Equations (1a) (or 2a) and (1b) are used stepwise in combination to conduct three sets of stress–strain experiments including sensitivity experiments (SEs):
A. for the period 1959–2015 assuming zero stress and strain in 1959,
B. for the period 1959–2015 assuming zero stress and strain in 1900, and
C. for the period 1959–2015 assuming zero stress and strain in 1850
and, ultimately, also before 1850 (i.e., zero anthropogenic stress before that date).

The logic of the experiments is determined by both the availability of data (see Supplement Information 6) and the increasing complementarity from A to C (see below). The basic procedure is always the same. We insert into Eq. (1a) our first-order estimates of $D_{\mathrm{L}}\approx (\mathrm{83}\pm \mathrm{44})$ Pa yr and $D_{\mathrm{O}}\approx (\mathrm{304}\pm \mathrm{60})$ Pa yr: that is, $D=D_{\mathrm{L}}+D_{\mathrm{O}}\approx (\mathrm{387}\pm \mathrm{74})$ Pa yr, and $K\approx (\mathrm{100}\pm \mathrm{90})$ Pa. At the same time, we use the growth factor α_ppm=0.0043 yr⁻¹, which reflects the exponential increase in the CO₂ concentration in the atmosphere between 1959 and 2018 (see Supplement Data 1) as our first-order estimate for α in $\dot{\mathit{\epsilon }}=\mathit{\alpha }\exp \left(\mathit{\alpha }t\right)$, the rate of change in strain ε. We apply Eq. (1a) by varying both $K/D$ and α to reproduce the known stress σ on the left, given by the CO₂ emissions from fossil fuel burning and land use. To restrict the number of variation parameters to two, we let K and D deviate from their respective mean values equally in relative terms (i.e., we assume that our first-order estimates exhibit equal inaccuracy in relative terms) and express α as a multiple of α_ppm. This is easily possible with the introduction of suitable factors (see Supplement Data 3) that allow σ to be reproduced quickly and with sufficient accuracy. The main reason this works well is that the two factors pull the two exponential functions on the right side of Eq. (2a) – $\dot{\mathit{\epsilon }}\left(t\right)$ and ($\mathrm{1}-q_{\mathit{\beta }}^t$), which determine the quality of the fit – in different directions.

To A

This is our set of reference experiments, all for the period 1959–2015. This set comprises A.1 (a stress-explicit experiment), A.2 (three strain-explicit experiments), and A.3 (SEs expanding the strain-explicit experiments). The parameters α, λ, and λ_β are reported in units of yr⁻¹, as is commonly done.

To A.1

In this experiment we vary the ratio $K/D$ (λ in Table 1) and α to reproduce the monitored stress σ(t) on the left side of Eq. (2a) (see Supplement Data 3). This tuning process (hereafter referred to as Case 0) allows us to test whether K and D, in particular, stay within their estimated limits, namely, $K\in [\mathrm{10};\mathrm{190}]$ Pa, and $D\in [\mathrm{313};\mathrm{461}]$ Pa yr or, equivalently, $\mathit{\lambda }\in [\mathrm{0.0217};\mathrm{0.6078}]$ yr⁻¹. Column “Case 0” in Table 1 indicates that this case is practically identical to choosing $\mathit{\lambda }=(\mathrm{10}/\mathrm{461}){\mathrm{yr}}^{-\mathrm{1}}=\mathrm{0.0217}{\mathrm{yr}}^{-\mathrm{1}}$, the smallest ratio $K/D$ deemed possible. For Case 0 we find K=9.9 Pa and D=461.5 Pa yr (thus, $\mathit{\lambda }=K/D=\mathrm{0.0214}$ yr⁻¹) as well as, concomitantly, α=0.0247 yr⁻¹ (thus, ${\mathit{\lambda }}_{\mathit{\beta }}=(K/D)\mathit{\beta }=(K/D)+\mathit{\alpha }=\mathrm{0.0461}$ yr⁻¹).

Table 1 Overview of parameters in experiments A.1–A.3.

^a Given per year (yr⁻¹). ^b Derived for K and D deviating from their respective mean values equally in relative terms. The notation “n.a.” stands for not assessable.

Figure 5 Case 0: $K/D$ and α on the right side of Eq. (2a) are tuned to reproduce the stress σ(t) on the left side of that equation, given by the monitored (but cumulated) CO₂ emissions from fossil fuel burning and land use activities (in Pa). The value resulting for $K/D$ complies with its lower limit deemed possible based on the uncertainties derived for K and D in Sect. 4.

Figure 5 reflects the result of the tuning process graphically. It shows how well the monitored stress, given by the cumulated CO₂ emissions from fossil fuel burning and land use activities since 1959, can be reproduced by Eq. (2a). The quality of the tuning is observed by summing the squares of differences between monitored and reproduced stress from 1959 to 2015 using the SUMXMY2 command in Excel. (We stopped the tuning process with the sum at about 1.400 Pa² when changes in K and D became negligible, resulting in a correlation coefficient of 0.9998; see Supplement Data 3.)

Figure 5 also shows the parameters needed to describe the monitored stress by a second-order polynomial regression (see the grey box in the upper left corner of the figure). We have not yet used this regression but will do so in the strain-explicit experiments described next.

To A.2

We use Eq. (1b) with $\mathit{\sigma }\left(\mathrm{0}\right)=\mathit{\epsilon }\left(\mathrm{0}\right)=\mathrm{0}$ and $\mathit{\sigma }\left(t\right)=\mathrm{0.0028}{t}^{\mathrm{2}}+\mathrm{0.1811}t$, the second-order polynomial regression of the monitored stress (see Fig. 5), to conduct three experiments (hereafter referred to as Cases 1–3) to explore the spread in the strain ε. To this end, we let the ratio $K/D$ vary from minimum (Case 1) to mean (Case 2) to maximum (Case 3; see Table 1 and Supplement Data 4) irrespective of the outcome of the Case 0 experiment, which suggests that compared to Cases 2 and 3, Case 1 (K minimal: the atmosphere is rather compressible; D maximal: the uptake of carbon by land and oceans is rather viscous) appears to be more in conformity with reality than Cases 2 and 3.

Figure 6 Cases 1–3: the ratio $K/D$ is varied from minimum (Case 1: solid red) to mean (Case 2: solid black) to maximum (Case 3: solid blue) to explore the spread in the strain ε (in units of 1) on the left side of Eq. (1b), while the monitored stress is described by a second-order polynomial (see the text). These strain responses have to be shifted upward (so that they pass through 1 in 1959) to derive their rates of change if described by an exponential regression (here only demonstrated for Case 2). As is already illustrated in Case 0, the exponential regression in Case 1 is excellent (see the text), whereas second-order polynomial regressions provide better fits in Cases 2 and 3 (see the boxes in the figure; the polynomial regressions are not shown).

Figure 6 reflects these experiments graphically. It shows that the range of strain responses is encompassed by Case 1 ($K/D=(\mathrm{10}/\mathrm{461})$ yr⁻¹) and Case 2 ($K/D=(\mathrm{100}/\mathrm{387})$ yr⁻¹) but not by Case 1 and Case 3 ($K/D=(\mathrm{190}/\mathrm{313})$ yr⁻¹) – the solid blue line (Case 3) falls between the solid red (Case 1) and solid black (Case 2) lines – resulting from how K and D dominate the individual parts of Eq. (1b). These strain responses have to be shifted upward (so that they pass through 1 in 1959) to describe them by an exponential regression and to derive their rates of change. The exponential fit is excellent only in Case 1, as already illustrated in Case 0 (Case 0: λ=0.0214 yr⁻¹, Case 1: λ=0.0217 yr⁻¹), but inferior to the polynomial regressions, here of the second order, in Cases 2 and 3. However, a second-order polynomial approach to the strain has to be discarded because the stress derived with the help of Eq. (1a) would exhibit a linear behavior with increasing time and not be a polynomial of the second order as in Fig. 6 (see Supplement Information 9).

Figure 7 Cases 1–3: (a) delay time T and memory M (in units of 1), (b) persistence P (in units of 1), and (c) the ratio $\mathit{\lambda }/\ln (M\cdot P)$ (in %); all are versus time (in units of 1) as of n=0 (1959).

In this regard we note that a more targeted way forward would be to use a piecemeal approach. This approach requires the data series to be sliced into shorter time intervals, during which an exponential fit for the strain (which we assume to hold in principle in deriving Eq. 2a here) is sufficiently appropriate. Fortunately, as the SEs in A.3 indicate, we can hazard the consequences of using suboptimal growth factors resulting from suboptimal exponential regressions for the strain.

Equations (3) to (5) are used to determine delay time T, memory M, and persistence P (in units of 1) for Cases 1–3 as well as their characteristic limiting values T_∞, M_∞, and P_∞ (see Table 1 and Supplement Data 5 to 8). We recall that T, M, and P are characteristic functions of the MB and are defined independently of initial conditions; these only specify the reference time for n=0 (here 1959). Figure 7a and b reflect the behavior of T, M, and P over time (in units of 1). For a better overview, Table 1 lists the times when these parameters exceed 50 % or 95 % of their limiting values (without indicating whether these levels go hand in hand with, e.g., global-scale ecosystem changes of equal magnitude). In the table we also specify the ratio $\mathit{\lambda }/\ln (M\cdot P)$ for each of these times (see also Fig. 7c). The ratio approaches $\mathit{\lambda }/{\mathit{\lambda }}_{\mathit{\beta }}$ for n→∞ and indicates (as a percentage) how much smaller the system's natural rate of change in the numerator turns out to be compared to the system's rate of change in the denominator under the continued increase in stress. As is illustrated, in particular by Case 1 in the figure, the ratio does not increase at a constant pace as n increases, which shows the nonlinear strain response of the atmosphere–land and ocean system.

To A.3

Three sets of SEs serve to assess the influence of the exponential growth factor on the strain-explicit experiments described above.
SE1: α₁=0.0248 yr⁻¹ as in Case 1 (see Fig. 6) is also used in Cases 2 and 3 (hereafter referred to as Cases 21 and 31).
SE2: α₂=0.0158 yr⁻¹ as in Case 2 (see Fig. 6) is also used in Cases 1 and 3 (hereafter referred to as Cases 12 and 32).
SE3: α₃=0.0174 yr⁻¹ as in Case 3 (see Fig. 6) is also used in Cases 1 and 2 (hereafter referred to as Cases 13 and 23).

Table 1 shows that the influence of a change in the exponential growth factor is small vis à vis the dominating influence of K and D and the quality in the estimates of T, M, and P. For instance, the dimensionless time n at $M/M_{\mathrm{\infty }}=\mathrm{0.5}$ ranges from 15 to 19 in Case 1 and experiments related to Case 1 (small persistency) and from 2 to 3 in Case 2 and experiments related to Case 2 (great persistency); in Case 3 and experiments related to Case 3, it does not exhibit a range at all (n≈1; very great persistency). These ranges for n tell us how long it takes to build up 50 % of the memory with time running as of n=0 (1959).

Table 2 Cases 1–3 and related experiments: build-up of memory (%) as of n=0 (1959).

^* Start year: ${\mathit{\sigma }}_{\mathrm{0}}={\mathit{\epsilon }}_{\mathrm{0}}=\mathrm{0}$.

Alternatively, we can ask how much memory has been built up until a given year. Table 2 tells us that after 56 years (i.e., in 2015) memory is still building up only in Case 1 and experiments related to Case 1, which means that the system still responds in its own characteristic way (as a result of a small K and a great D) to the continuously increasing stress; this is not so in Cases 2 and 3 (and related experiments). In the latter two cases today's uptake of carbon by land and oceans happens de facto outside the system's natural regime and solely in response to the sheer, continuously increasing stress imposed on it, whereas in Case 1 and experiments related to Case 1 the limits of the natural regime are not yet reached. This interpretation of Cases 1–3 (and related experiments) does not depend on how much carbon the system already took up before 1959. M is additive and defined independently of initial conditions; these only specify 1959 as reference time for n=0. This means by implication that the current M value (or its perpetuation) is contained in the M value (or is part of that value's perpetuation) which starts accruing from an earlier point in time (see also experiments B and C below).

Finally, it is important to note that it is prudent to expect that natural elements (like land and oceans) will not continue to maintain their damping (i.e., carbon uptake) capacity – or their capacity to embark on a, most likely, hysteretic downward path in the case of a sustained decrease in emissions – even well before they reach the limits of their natural regimes. They may simply collapse globally when reaching a critical threshold. We note that our choice of model binds us to the global scale and also does not allow “failure” to be specified further, e.g., with respect to when exactly a critical threshold will occur and in terms of whether carbon uptake decreases only or even ceases upon reaching the threshold.

To B and C

We report on the sets of stress–strain experiments B and C in combination. They can be understood as a repetition of the 1959–2015 Case 0 experiment (see A.1) but with the difference that now upstream emissions as of 1900 (B) or 1850 (C) are considered. This allows initial conditions for 1959 other than zero, as in the Case 0 experiment, to be considered (see Supplement Information 10 and Supplement Data 9 to 16).
Case 0: 1959–2015.
B: 1900–1958 (upstream emissions), 1959–2015.
C: 1850–1958 (upstream emissions), 1959–2015.

The experiments can be ordered consecutively in terms of time with the three 1959–2015 periods comprising a min–max interval to facilitate the drawing of a number of robust results in spite of the uncertainty underlying these stress–strain experiments (see Supplement Information 10). Between 1850 and 1959–2015 (i) the compression modulus K increased from ∼2 to 10–13 Pa (the atmosphere became less compressible), while (ii) the damping constant D decreased from ∼468 to 459–462 Pa yr (the uptake of carbon by land and oceans became less viscous), with the consequence that (iii) the ratio $\mathit{\lambda }=K/D$ increased from ∼0.004–0.005 yr⁻¹ to 0.021–0.028 yr⁻¹ (i.e., by a factor of 4–6). Likewise, (iv) delay time T_∞ decreased (hence, persistence P_∞ increased) from ∼51 (∼0.02) to 18–21 (0.047–0.055), while (v) memory M_∞ decreased from ∼52 to 19–22 on the dimensionless timescale.

6 Account of the findings

Here we discuss our main findings in greater depth, recollect the assumptions underlying our global stress–strain approach, and conclude by returning to the three questions posed at the beginning.

We make use of an MB to model the stress–strain behavior of the global atmosphere–land and ocean carbon system and to simulate how humankind propelled that global-scale experiment historically, here as of 1850. The stress is given by the CO₂ emissions from fossil fuel burning and land use, while the strain is given by the expansion of the atmosphere by volume and uptake of CO₂ by sinks. The MB is a logical choice of stress–strain model given the uninterrupted increase in atmospheric CO₂ concentrations since 1850.

The stress–strain model is unique and a valuable addendum to the suite of models (such as radiation transfer, energy balance, or box-type carbon cycle models), which are highly reduced but do not compromise complexity in principle. These models offer great benefits in safeguarding complex three-dimensional global change models. Here, too, the proposed stress–strain approach allows three system-characteristic parameters to be distilled from the stress-explicit equation – delay time, memory, and persistence – and new insights to be gained. What we consider most important is that these parameters come with their own internal limits, which govern the behavior of the atmosphere–land and ocean carbon system. These limits are independent from any external target values (such as temperature targets justified by means of global change research).

Knowing these limits is precisely the reason why we can advance the discussion and draw some preliminary conclusions. To start with, we look at the Case 0 experiment and the stress–strain experiments B and C in combination. The values of the Case 0 parameters T_∞ and M_∞, in particular, are at the upper end of the respective 1959–2015 min–max intervals (see Supplement Information 10). That is, the respective characteristic ratios $T/T_{\mathrm{\infty }}$ and $M/M_{\mathrm{\infty }}$ reach specified levels (e.g., 0.5 or 0.95; see Fig. 7a) slightly sooner than when T_∞ and M_∞ take on values at the lower end of the 1959–2015 min–max intervals. Given that Case 0 is well represented by Case 1, we can use the parameter values of the latter. According to column “Case 1” in Table 1, $M/M_{\mathrm{\infty }}$ and $T/T_{\mathrm{\infty }}$ reached their 0.5 levels after about 15- and 28-year-equivalent units on the dimensionless timescale (which was in 1974 and 1987), whereas they will reach their 0.95 levels after about 64- and 98-year-equivalent units (which will be in 2023 and 2057) – if the exponential growth factor α remains unchanged in the future.

This not unthinkable worst case provides a reference, as follows: we understand, in particular, the ability of a system to build up memory effectively as its ability to still respond to stress in its own characteristic way (i.e., within its natural regime). Therefore, it appears precautionary to prefer memory over delay time in avoiding potential system failures globally in the future. These we expect to happen well before 2050 if the current trend in emissions is not reversed immediately and sustainably. However, we reiterate that our choice of model binds us to the global scale and also does not allow failure to be specified further.

We consider our precautionary statement robust given both the uncertainties we dealt with in the course of our evaluation and the restriction of our variation parameters to two. One of the two variation parameters (λ) presupposes knowing K and D with equal inaccuracy in relative terms. This procedural measure in treating λ, in particular, offers a great applicational benefit but no serious restriction given that (while, ideally, α is constant) it is the $K/D$ ratio that matters and whose ultimate value is controlled by consistency – which comes in as a powerful rectifier. As a matter of fact, fulfilling consistency results in a $K/D$ ratio that ranges close to the lower uncertainty boundary which we deem adequate based on our preceding assessment. That means a smaller K, i.e., the atmosphere is more compressible than previously thought, and a greater D, i.e., the uptake of carbon by land and oceans is more viscous than previously thought (see Cases 1–3 in Table 1). However, the overall effect of the continued release of CO₂ emissions since 1850 on the $K/D$ ratio is unambiguous – the ratio increased (see λ in Table SI10-2) by a factor of 4–6 (K increased: the atmosphere became less compressible; D decreased: the uptake of carbon by land and oceans became less viscous).

By way of contrast, persistence is less intelligible. Equation (5) allows persistence (as well as its systemic limit) to be followed quantitatively. However, it is conducive to understand persistence as path dependency and in qualitative terms, i.e., whether it increased or decreased. Thus, we see that P_∞ increased since 1850 by a factor of 2–3 (see P_∞ in Table SI10-2), which indicates that the atmosphere–land and ocean system is progressively trapped from a path dependency perspective. This, in turn, means that it will become progressively more difficult to (strain) relax the entire system (i.e., the atmosphere including land and oceans) – a mere 1-year decrease of a few percentage points in CO₂ emissions, as reported recently for 2020, will have virtually no impact (Global Carbon Project, 2020).

To conclude, we return to the three questions posed in the beginning. These can be answered unambiguously.

Memory, just as persistence, is a characteristic (function) of the MB. Mathematically spoken, it is contained in the integral on the right side of Eq. (1a) and is defined independently of initial conditions. These appear only in the lower boundary of that integral, which allows initial conditions other than zero to be considered by taking advantage of the integral's additivity.

The memory of the atmosphere–land and ocean carbon system – Earth's memory – can be quantified. It can be understood as the depreciated strain summed up backward in time. We let memory extend backward in time to 1850, assuming zero anthropogenic stress before that date. Memory is measured in units of 1 and accrues continually over time (here as the result of the uninterrupted increase in stress).

Memory is constrained. It can be compared with a limited buffer, approximately 60 % of which humankind had already exploited prior to 1959 (see M_∞ in Table SI10-2). We understand the effective build-up of memory as Earth's ability to still respond within its own natural stress–strain regime. However, this ability declines considerably with memory reaching high levels of exploitation (see $M/M_{\mathrm{\infty }}\ge \mathrm{0.95}$ in Table 1) – which we anticipate happening in the foreseeable future if CO₂ emissions continue to increase globally as before.

Finally, we can also quantify the persistence of the atmosphere–land and ocean carbon system. It is also measured in units of 1. Persistence can be understood intuitively as path dependency and in qualitative terms. Concomitantly with the exploitation of memory, we see that P_∞ increased since 1850 by approximately a factor of 2–3 – and can be expected to increase further if the release of CO₂ emissions globally continues as before.

Based on these stress–strain insights we expect that the atmosphere–land and ocean carbon system will be forced outside its natural regime well before 2050 if the current trend in emissions is not reversed immediately and sustainably.

Appendix A: Acronyms and nomenclature

If terms or symbols are used in more than one way, we make them unambiguous by specifying (in parentheses) how they are used in the paper (e.g., CO₂ as a chemical formula in the text or as a physical parameter in units of ppmv in mathematical equations). As a basic rule, physical parameters are always specified by their units.

ad	adiabatic
C	carbon
comb	combined
CO2	carbon dioxide (chemical formula)
CO2	atmospheric CO2 concentration (in ppmv;
	parameter)
D	damping constant (in Pa yr)
DIC	dissolved inorganic carbon (in µmol kg−1)
E	Young's modulus (in Pa)
GHG	greenhouse gas
h	altitude (in m)
it	isothermal
K	compression modulus (in Pa)
L	land (index)
L	leaf-level factor (in ppmv−1; parameter)
M	memory (in units of 1)
MB	Maxwell body
n.a.	not assessable

NPP	net primary productivity and production
	(in Pg C yr−1)
O	oceans
p	atmospheric pressure (in hPa)
pCO2	partial pressure of atmospheric
	CO2 (in µatm)
P	persistence (in units of 1)
Ph	global photosynthetic carbon influx
	(in Pg C yr−1)
q	auxiliary quantity (in units of 1)
R	Revelle (buffer) factor (in units of 1)
SD	supplementary data
SE	sensitivity experiment
SI	supplementary information
t	time (in yr)
T	delay time (in units of 1)
TOA	top of the atmosphere
w	weight(ed)
α	exponential growth factor of the strain
	(in yr−1)
αppm	exponential growth factor of the atmospheric
	CO2 concentration (in yr−1)
β	auxiliary quantity (in units of 1)
βb	biotic growth factor (in units of 1)
βPh	photosynthetic beta factor (in units of 1)
ε	strain (referring to atmospheric expansion by
	volume and CO2 uptake by sinks; in units
	of 1)
γ	isentropic coefficient of expansion (in units
	of 1)
κ	(in Pa−1)
σ	stress (atmospheric CO2 emissions from fossil
	fuel burning and land use; in Pa)