The climate benefit of carbon sequestration

Ecosystems play a fundamental role in climate change mitigation by photosynthetically fixing carbon from the atmosphere and storing it for a period of time in organic matter. Although climate impacts of carbon emissions by sources can be quantified by global warming potentials, the appropriate formal metrics to assess climate benefits of carbon removals by sinks are unclear. We introduce here the climate benefit of sequestration (CBS), a metric that quantifies the radiative effect of fixing carbon dioxide from the atmosphere and retaining it for a period of time in an ecosystem before releasing it back as the result of respiratory processes and disturbances. In order to quantify CBS, we present a formal definition of carbon sequestration (CS) as the integral of an amount of carbon removed from the atmosphere stored over the time horizon it remains within an ecosystem. Both metrics incorporate the separate effects of (i) inputs (amount of atmospheric carbon removal) and (ii) transit time (time of carbon retention) on carbon sinks, which can vary largely for different ecosystems or forms of management. These metrics can be useful for comparing the climate impacts of carbon removals by different sinks over specific time horizons, to assess the climate impacts of ecosystem management, and to obtain direct quantifications of climate impacts as the net effect of carbon emissions by sources versus removals by sinks.

1 Introduction

Terrestrial ecosystems exchange carbon with the atmosphere at globally significant quantities, thereby influencing Earth's climate and potentially mitigating warming caused by increasing concentrations of CO₂ in the atmosphere. Carbon fixed during the process of photosynthesis remains stored in the terrestrial biosphere over a range of timescales, from days to millennia – timescales of relevance for affecting the concentration of greenhouse gases in the atmosphere (Archer et al., 2009; IPCC, 2014; Joos et al., 2013). During the time carbon is stored in the terrestrial biosphere, it is removed from the radiative forcing effect that occurs in the atmosphere; thus, it is of scientific and policy relevance to understand the timescale of carbon storage in ecosystems, i.e., for how long newly fixed carbon is retained in an ecosystem before it is released back to the atmosphere.

Timescales of element cycling and storage are unambiguously characterized by the concepts of system age and transit time (Bolin and Rodhe, 1973; Rodhe, 2000; Rasmussen et al., 2016; Sierra et al., 2017; Lu et al., 2018). In a system of multiple interconnected compartments, system age characterizes the time that the mass of an element observed in the system has remained there since its entry. Transit time characterizes the time that it takes element masses to traverse the entire system, from the time of entry until they are released back to the external environment (Sierra et al., 2017). Both metrics are excellent system-level diagnostics of the dynamics and timescales of ecosystem processes. Because system age and transit time both can be reported as mass or probability distributions, they provide different information about an ecosystem over a wide range in the time domain.

System age and transit time are closely related to the complexity of the ecosystem and its process rates, which are affected by the environment (Luo et al., 2017; Rasmussen et al., 2016; Sierra et al., 2017; Lu et al., 2018). Mean system ages of carbon are consistently greater than mean transit time (Lu et al., 2018; Sierra et al., 2018b), suggesting that once a mass of carbon enters an ecosystem a large proportion gets quickly released back to the atmosphere, but a small proportion remains for very long times. Furthermore, differences in transit times across ecosystems suggest that not all carbon sequestered in the terrestrial biosphere spends the same amount of time stored; e.g., one unit of photosynthesized carbon is returned back to the atmosphere faster in a tropical than in a boreal forest (Lu et al., 2018). Therefore, not all carbon drawn down from the atmosphere should be treated equally for the purpose of quantifying the climate mitigation potential of sequestering carbon in ecosystems as it is currently recommended in accounting methodologies (IPCC, 2006).

Global warming potentials (GWPs; see definition in Sect. 2) quantify the radiative effects of greenhouse gases emitted to the atmosphere (Fig. 1) but do not consider the avoided radiative effect of storing carbon in ecosystems (Neubauer and Megonigal, 2015). GWPs are computed using the age distribution of CO₂ and other greenhouse gases in the atmosphere (Rodhe, 1990; Joos et al., 2013) but do not consider age or transit times of carbon in ecosystems in the case of sequestration. Transit time distributions, in particular, can better inform us about the time newly sequestered carbon will be removed from radiative effects in the atmosphere.

For more comprehensive accounting of the contribution of carbon sequestration to climate change mitigation, it is necessary to quantify the avoided warming effects of sequestered carbon in ecosystems over the timescale the carbon is stored. The GWP metric is inappropriate to quantify avoided warming potential as a result of sequestration. A metric that can capture this avoided warming effect could have applications for (1) comparing different carbon sequestration activities considering the time carbon is stored in ecosystems and (2) providing better accounting methods for the effect of removals by sinks in climate policy. Currently, the Intergovernmental Panel on Climate Change (IPCC) recommends that countries and project developers report only emissions by sources and removals by sinks of greenhouse gases (GHGs), treating all removals equally in terms of their fate (IPCC, 2006).

Problems with applying GWPs to compute climate benefits of sequestering carbon in ecosystems are well documented (Moura Costa and Wilson, 2000; Fearnside et al., 2000; Brandão et al., 2013; Neubauer and Megonigal, 2015). Several approaches have been proposed to deal with the issue of timescales (Brandão et al., 2013), many of which deal with time as some form of delay in emissions. However, to our knowledge, no solution proposed thus far explicitly accounts for the time carbon is sequestered in ecosystems, from the time of photosynthetic carbon fixation until it is returned back to the atmosphere by autotrophic and heterotrophic respiration, and fires.

Therefore, the main objective of this paper is to introduce a metric to assess the climate benefits of carbon sequestration while accounting for the time carbon is stored in ecosystems. We first present the theoretical framework for the development of the metric, then provide simple examples for its computation, and discuss potential applications for ecosystem management and for climate change mitigation.

2 Theoretical framework

2.1 Absolute global warming potential (AGWP)

The direction of carbon flow, into or out of ecosystems, is of fundamental importance to understand and quantify their contribution to climate change mitigation. The absolute global warming potential (AGWP) of carbon dioxide quantifies the radiative effects of a unit of CO₂ emitted to the atmosphere during its life time – in the direction land → atmosphere. It is expressed as (Lashof and Ahuja, 1990; Rodhe, 1990)

$$\begin{array}{}\text{(1)}& \mathrm{AGWP}(T,t_{\mathrm{0}})=\underset{t_{\mathrm{0}}}{\overset{t_{\mathrm{0}}+T}{\int }}{k}_{{\mathrm{CO}}_{\mathrm{2}}}\left(t\right)M_{\mathrm{a}}\left(t\right)\mathrm{d}t,\end{array}$$

where $k_{{\mathrm{CO}}_{\mathrm{2}}}\left(t\right)$ is the radiative efficiency or greenhouse effect of one unit of CO₂ (in mole or mass) in the atmosphere at time t, and M_a(t) is the amount of gas present in the atmosphere at time t (Rodhe, 1990; Joos et al., 2013). The AGWP quantifies the amount of warming produced by CO₂, while it stays in the atmosphere since the time the gas is emitted at time t₀ over a time horizon T. The function M_a(t) quantifies the fate of the emitted carbon in the atmosphere and can be written in general form as

$$\begin{array}{}\text{(2)}& M_{\mathrm{a}}\left(t\right)=h_{\mathrm{a}}(t-t_{\mathrm{0}})M_{\mathrm{a}}\left(t_{\mathrm{0}}\right)+\underset{t_{\mathrm{0}}}{\overset{t}{\int }}{h}_{\mathrm{a}}(t-\mathit{\tau })Q\left(\mathit{\tau }\right)\mathrm{d}\mathit{\tau },\end{array}$$

where h_a(t−t₀) is the impulse response function of atmospheric CO₂ released into the atmosphere, M_a(t₀) is the content of atmospheric CO₂ at time t₀, and Q(τ) is the perturbation of new incoming carbon to the atmosphere between t₀ and t.

For a pulse, or instantaneous emission of CO₂, M_a(t₀)=E₀, and

$$\begin{array}{}\text{(3)}& M_{\mathrm{a}}\left(t\right)=h_{\mathrm{a}}(t-t_{\mathrm{0}})E_{\mathrm{0}},\end{array}$$

assuming no additional carbon enters the atmosphere after the pulse. If the pulse is equivalent to 1 kg or mole of CO₂, then E₀=1 and $M_{\mathrm{a}}\left(t\right)=h_{\mathrm{a}}(t-t_{\mathrm{0}})$. For a pulse emission of any arbitrary size, and assuming constant radiative efficiency (see details about this assumption in Sect. 2.2),

$$\begin{array}{}\text{(4)}& \mathrm{AGWP}(T,E_{\mathrm{0}},t_{\mathrm{0}})=k_{{\mathrm{CO}}_{\mathrm{2}}}{E}_{\mathrm{0}}\underset{t_{\mathrm{0}}}{\overset{t_{\mathrm{0}}+T}{\int }}{h}_{\mathrm{a}}(t-t_{\mathrm{0}})\mathrm{d}t.\end{array}$$

The AGWP can be computed for any other greenhouse gas using their respective radiative efficiencies and fate in the atmosphere (impulse response function). To compare different gases, the global warming potential (GWP) is defined as the AGWP of a particular gas divided by the AGWP of CO₂ (Shine et al., 1990; Lashof and Ahuja, 1990). Our interest in this paper is on carbon fixation and respiration in the form CO₂; therefore, we primarily concentrate here on AGWP.

The impulse response function h_a(t−t₀) plays a central role within the AGWP framework. The function encodes information about the fate of a gas once it enters the atmosphere and determines for how long the gas will remain. Therefore, it can be interpreted as a density distribution for the transit time of a gas, since the time of emission until it is removed by natural sinks (e.g., CO₂) or by chemical reactions (e.g., CH₄).

The function typically is assumed to be static; i.e., the time at which the gas enters the atmosphere is not relevant, only the time it remains there (t−t₀). However, this function can be time-dependent, expressing different shapes depending on the time the gas enters the atmosphere, i.e., $h_{\mathrm{a}}(t_{\mathrm{0}},t-t_{\mathrm{0}})$. For example, when natural sinks saturate, faster accumulation of CO₂ and longer transit times of carbon in the atmosphere are observed (Metzler et al., 2018). In this situation, the specific time of an emission would lead to different response functions in the atmosphere. Because current research on impulse response functions primarily considers the static time-independent case (see Millar et al., 2017, for an exception), we will consider only the static case for the remainder of this paper.

2.2 The radiative efficiency of CO₂ and its impulse response function

The radiative efficiency of CO₂ is a function of the concentration of this gas and the concentration of other gases in the atmosphere with overlapping absorption bands (Lashof and Ahuja, 1990; Shine et al., 1990). Therefore, $k_{{\mathrm{CO}}_{\mathrm{2}}}$ changes as the concentration of GHGs change in the atmosphere. For most applications however, the radiative efficiency of CO₂ has been assumed constant in the limit of a small perturbation at a specific background concentration (Lashof and Ahuja, 1990; Shine et al., 1990; Joos et al., 2013; Myhre et al., 2013).

Here, we use a constant value of $k_{{\mathrm{CO}}_{\mathrm{2}}}=\mathrm{6.48}\times {\mathrm{10}}^{-\mathrm{12}}$ W m⁻² per megagram of carbon based on results reported by Joos et al. (2013) for an atmospheric background of 389 ppm (∼ present day). This radiative efficiency represents the change in radiative forcing caused by a change of 1 Mg of carbon in the atmosphere in the form of CO₂ in units of rate of energy transfer (watt) per square meter of surface.

Joos et al. (2013) have also derived impulse response functions (IRFs) of CO₂ in the atmosphere using coupled carbon–climate models that include multiple feedbacks among Earth system processes. One function was obtained by emitting a pulse of 100 Gt of carbon to a pre-industrial atmosphere with a background concentration of 280 ppm (PI100 function from here on), and another function was obtained by emitting 100 Gt of carbon to a present-day atmosphere with a background of 389 ppm (PD100 from here on). The functions they report are averages from the numerical output of multiple models fitted to a sum of exponential functions that include an intercept term. This intercept implies that a proportion of the added CO₂ never leaves from the atmosphere–ocean–terrestrial system to long-term geological reservoirs. Following Millar et al. (2017), we added a timescale of 1 million years that corresponds to the intercept term in the IRFs. The addition of this timescale has no effect on the results presented here, which are focused on much shorter timescales, but they avoid the mathematical problem that the integrals of the original functions go to infinity with time (Lashof and Ahuja, 1990; Millar et al., 2017).

2.3 Carbon sequestration CS and the climate benefit of carbon sequestration (CBS)

GWPs are useful to quantify the climate impacts of increasing or reducing emissions of GHGs to the atmosphere. However, it is also necessary to quantify the climate benefits of carbon flows in the opposite direction, atmosphere → land. Furthermore, it is also important to quantify not only how much and how fast carbon enters ecosystems, but also for how long the carbon stays (Körner, 2017).

Carbon taken up from the atmosphere through the process of photosynthesis is stored in multiple ecosystem reservoirs for a particular amount of time. Carbon sequestration can be defined as the process of capture and long-term storage of CO₂ (Sedjo and Sohngen, 2012). We define here carbon sequestration CS over a time horizon T as

$$\begin{array}{}\text{(5)}& \mathrm{CS}(T,S_{\mathrm{0}},t_{\mathrm{0}}):=\underset{t_{\mathrm{0}}}{\overset{t_{\mathrm{0}}+T}{\int }}{M}_{\mathrm{s}}(t-t_{\mathrm{0}})\mathrm{d}t,\end{array}$$

where M_s(t−t₀) represents the fate of a certain amount of carbon S₀ taken up by the sequestering system at a time t₀. Notice that this definition of carbon sequestration is very similar to that of AGWP for an emission, with the exception that the radiative efficiency term is omitted.

To obtain the fate of sequestered carbon over time, we represent carbon cycling and storage in ecosystems using the theory of compartmental dynamical systems (Luo et al., 2017; Sierra et al., 2018a). In their most general form, we can write carbon cycle models as

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\mathrm{d}\mathit{x}\left(t\right)}{\mathrm{d}t}}=\dot{\mathit{x}}\left(t\right)=\mathit{u}(\mathit{x},t)+\mathbf{B}(\mathit{x},t)\mathit{x},\end{array}$$

where x(t)∈ℝⁿ is a vector of n ecosystem carbon pools, $\mathit{u}(\mathit{x},t)\in {\mathbb{R}}^n$ is a time-dependent vector-valued function of carbon inputs to the system, and $\mathbf{B}(\mathit{x},t)\in {\mathbb{R}}^{n\times n}$ is a time-dependent compartmental matrix. The latter two terms can depend on the vector of states, in which case the compartmental system is considered nonlinear. In case the input vector and the compartmental matrix have fixed coefficients (no time dependencies), the system is considered autonomous, and it is considered non-autonomous otherwise (Sierra et al., 2018a). This distinction of models with respect to linearity and time dependencies (autonomy) is fundamental to distinguish important properties of models. For instance, models expressed as autonomous linear systems have a steady-state solution given by ${\mathit{x}}^{\ast }=-{\mathbf{B}}^{-\mathrm{1}}\mathit{u}$, where x^∗ is a vector of steady-state contents for all ecosystem pools. Non-autonomous models have no steady-state solution.

The fate of the fixed carbon for the general nonlinear non-autonomous case can be obtained as

$$\begin{array}{}\text{(7)}& M_{\mathrm{s}}(t-t_{\mathrm{0}})=\parallel \mathbf{\Phi }(t,t_{\mathrm{0}})\mathit{\beta }\left(t_{\mathrm{0}}\right)S_{\mathrm{0}}\parallel ,\end{array}$$

where β(t₀)S₀=u(t₀), and β(t₀) is an n-dimension vector representing the partitioning of the total sequestered carbon among n ecosystem carbon pools (Ceballos-Núñez et al., 2020). The n×n matrix Φ(t,t₀) is the state-transition operator, which represents the dynamics of how carbon moves in a system of multiple interconnected compartments (see details in Appendix B). Throughout this document, we use the symbol $\Vert \Vert$ to denote the 1-norm of a vector, i.e., the sum of the absolute values of all elements in a vector.

Because ecosystems and most reservoirs are open systems, the sequestered carbon S₀ returns back to the atmosphere, mostly as CO₂ due to ecosystem respiration and fires. Carbon release r(t) from ecosystems can be obtained according to

$$\begin{array}{}\text{(8)}& r\left(t\right)=-{\mathbf{1}}^{\mathit{⊺}}\mathbf{B}\left(t\right)\mathbf{\Phi }(t,t_{\mathrm{0}})\mathit{\beta }\left(t_{\mathrm{0}}\right)S_{\mathrm{0}},\end{array}$$

where 1^⊺ is the transpose of the n-dimensional vector containing only 1s. The state-transition matrix captures the entire fate and dynamics of the sequestered carbon, from the time it enters t₀ until release at any t.

The link between the time it takes sequestered carbon S₀ to appear in the release flux r(t) is established by the concept of transit time (Metzler et al., 2018). In particular, we define the forward transit time (FTT) as the age that fixed carbon will have at the time it is released back to the atmosphere, or how long a mass fixed now will stay in the system. The backward transit time (BTT) is defined as the age of the carbon in the output flux since the time it was fixed, or how long the mass leaving the system now had stayed. This implies that

$$\begin{array}{}\text{(9)}& r\left(t\right)=p_{\mathrm{BTT}}(t-t_{\mathrm{0}},t)=p_{\mathrm{FTT}}(t-t_{\mathrm{0}},t_{\mathrm{0}}),\end{array}$$

where $p_{\mathrm{BTT}}(t-t_{\mathrm{0}},t)$ is the backward transit time distribution of carbon leaving the system at time t with an age t−t₀, while $p_{\mathrm{FTT}}(t-t_{\mathrm{0}},t_{\mathrm{0}})$ is the forward transit time distribution of carbon entering the system at time t₀ and leaving with an age t−t₀. For systems in equilibrium, both quantities are equal (Metzler et al., 2018). For systems not in equilibrium, semi-explicit formulas for their distributions are given in Appendix B.

For the atmosphere, carbon sequestration is a form of negative emission, and we can represent its fate in the atmosphere as

$$\begin{array}{}\text{(10)}& M_{\mathrm{a}}^{\prime }\left(t\right)=-h_{\mathrm{a}}(t-t_{\mathrm{0}})S_{\mathrm{0}}+\underset{t_{\mathrm{0}}}{\overset{t}{\int }}{h}_{\mathrm{a}}(t-\mathit{\tau })r\left(\mathit{\tau }\right)\mathrm{d}\mathit{\tau },\end{array}$$

where the prime symbol represents a perturbed atmosphere as an effect of sequestration. The first term in the right-hand side represents the response of the atmosphere to an instantaneous sequestration S₀ at t₀, and the second term represents the perturbation in the atmosphere of the carbon returning back from the terrestrial biosphere. Notice that the integral in this equation can be written as a convolution (h_a⋆r)(t) between the impulse response function of atmospheric CO₂ and the carbon returning from ecosystems to the atmosphere.

We define now the climate benefit of sequestration for a pulse of CO₂ into an ecosystem as

$$\begin{array}{}\text{(11)}& \begin{array}{rl}& \mathrm{CBS}(T,S_{\mathrm{0}},t_{\mathrm{0}}):=\underset{t_{\mathrm{0}}}{\overset{t_{\mathrm{0}}+T}{\int }}{k}_{{\mathrm{CO}}_{\mathrm{2}}}{M}_{\mathrm{a}}^{\prime }\left(t\right)\mathrm{d}t,\\ & =-k_{{\mathrm{CO}}_{\mathrm{2}}}\underset{t_{\mathrm{0}}}{\overset{t_{\mathrm{0}}+T}{\int }}\left(h_{\mathrm{a}}(t-t_{\mathrm{0}})S_{\mathrm{0}}-(h_{\mathrm{a}}\star r)\left(t\right)\right)\mathrm{d}t.\end{array}\end{array}$$

This metric integrates over a time horizon T the radiative effect avoided by sequestration of an amount of carbon S₀ taken up at time t₀ by an ecosystem. It captures the timescale at which the carbon is stored and gradually returns back to the atmosphere. It can also be interpreted as the atmospheric response to carbon sequestration in the form of a negative emission of CO₂ during a time horizon of interest. It relies on knowledge of the atmospheric response to perturbations in the form of an impulse response function and the transit time of carbon in an ecosystem.

Figure 1 Contrast between current approach to quantification of climate effects of emissions and sequestration (a), and the proposed approach for sequestration (b). Plots and equations represent the concepts of absolute global warming potential (AGWP) of an emission of CO₂, carbon sequestration (CS), and climate benefits of sequestration (CBS). AGWP integrates over a time horizon T the fate of an instant emission at time t₀ of a gas (M_a(t)) and multiplies by the radiative efficiency k of the gas. A similar idea can be used to define CS as the integral of the fate M_s(t) of an instant amount of carbon uptake S₀ over T. The CBS captures the atmospheric disturbance caused by CO₂ uptake and subsequent release by respiration as the integral over T of the fate of sequestered carbon $M_{\mathrm{a}}^{\prime }\left(t\right)$ multiplied by the radiative efficiency of CO₂.

2.4 Ecosystems in equilibrium: the linear, steady-state case

The computation of CS and CBS is simplified for systems in equilibrium. For linear systems at a steady state, the time at which the carbon enters the ecosystem is irrelevant (Kloeden and Rasmussen, 2011; Rasmussen et al., 2016); one only needs to know for how long the carbon has been in the system to predict how much of it remains. Mathematically, this implies

$$\begin{array}{}\text{(12)}& \mathbf{\Phi }(t,t_{\mathrm{0}})=e^{a\cdot \mathbf{B}}\text{for all}{t}_{\mathrm{0}}\le t\text{and}a=t-t_{\mathrm{0}}.\end{array}$$

Therefore, for linear systems at a steady state, we have the special cases

$$\begin{array}{}\text{(13)}& M_{\mathrm{s}}\left(a\right)=\Vert e^{a\cdot \mathbf{B}}\mathit{u}\Vert ,\end{array}$$

and

$$\begin{array}{}\text{(14)}& M_{\mathrm{s}\mathrm{1}}\left(a\right)=\left|\left|e^{a\cdot \mathbf{B}}{\displaystyle \frac{\mathit{u}}{\Vert \mathit{u}\Vert }}\right|\right|,\end{array}$$

where M_s1 represents the fate of one unit of fixed carbon, which can also be interpreted as the proportion of carbon remaining after the time of fixation.

The amount of released carbon returning to the atmosphere is therefore

$$\begin{array}{}\text{(15)}& r\left(a\right)=-{\mathbf{1}}^{\mathit{⊺}}\mathbf{B}{e}^{a\cdot \mathbf{B}}\mathit{u},\end{array}$$

which for one unit of fixed carbon is equal to the transit time density distribution f(τ) of a linear system (Metzler and Sierra, 2018, see also Appendix B)

$$\begin{array}{}\text{(16)}& r_{\mathrm{1}}\left(a\right)=-{\mathbf{1}}^{\mathit{⊺}}\mathbf{B}{e}^{a\cdot \mathbf{B}}{\displaystyle \frac{\mathit{u}}{\Vert \mathit{u}\Vert }},\end{array}$$

where r₁(a)=f(τ), with mean (expected value) transit time given by

$$\begin{array}{}\text{(17)}& \mathbb{E}\left(\mathit{\tau }\right)=-{\mathbf{1}}^{\mathit{⊺}}{\mathbf{B}}^{-\mathrm{1}}{\displaystyle \frac{\mathit{u}}{\Vert \mathit{u}\Vert }}={\displaystyle \frac{\Vert {\mathit{x}}^{\ast }\Vert }{\Vert \mathit{u}\Vert }}.\end{array}$$

We can now derive the steady-state expression of CS as

$$\begin{array}{}\text{(18)}& \mathrm{CS}\left(T\right)=\underset{\mathrm{0}}{\overset{T}{\int }}\Vert e^{a\cdot \mathbf{B}}\mathit{u}\Vert \mathrm{d}a.\end{array}$$

Furthermore, it is possible to find a closed-form expression for this integral:

$$\begin{array}{}\text{(19)}& \mathrm{CS}\left(T\right)=\Vert {\mathbf{B}}^{-\mathrm{1}}\left(e^{T\cdot \mathbf{B}}-\mathbf{I}\right)\mathit{u}\Vert ,\end{array}$$

where $\mathbf{I}\in {\mathbb{R}}^{n\times n}$ is the identity matrix. Similarly, for one unit of carbon entering a steady-state system at any time, we define CS₁ as

$$\begin{array}{}\text{(20)}& {\mathrm{CS}}_{\mathrm{1}}\left(T\right)=\underset{\mathrm{0}}{\overset{T}{\int }}∥e^{a\cdot \mathbf{B}}{\displaystyle \frac{\mathit{u}}{\Vert \mathit{u}\Vert }}∥\mathrm{d}a,\end{array}$$

which by integration gives

$$\begin{array}{}\text{(21)}& {\mathrm{CS}}_{\mathrm{1}}\left(T\right)=∥{\mathbf{B}}^{-\mathrm{1}}\left(e^{T\cdot \mathbf{B}}-\mathbf{I}\right){\displaystyle \frac{\mathit{u}}{\Vert \mathit{u}\Vert }}∥.\end{array}$$

These steady-state expressions can be very useful to compare different systems or changes to a particular system if the steady-state assumption is justified. Furthermore, it can be shown that in the long term, as the time horizon T goes to infinity (∞), the term $(e^{T\cdot \mathbf{B}}-\mathbf{I})$ converges to −I, and therefore Eq. (19) converges to the expression

$$\begin{array}{}\text{(22)}& \underset{T\to \mathrm{\infty }}{lim}\mathrm{CS}\left(T\right)=\Vert {\mathit{x}}^{\ast }\Vert ,\end{array}$$

which means that the total amount of carbon at a steady state is equal to the long-term carbon sequestration of an instantaneous amount of fixed carbon at an arbitrary time.

Similarly, for one unit of carbon entering a system at a steady state, the long-term CS₁ from Eq. (21) can be obtained simply as

$$\begin{array}{}\text{(23)}& \underset{T\to \mathrm{\infty }}{lim}{\mathrm{CS}}_{\mathrm{1}}\left(T\right)=\mathbb{E}\left(\mathit{\tau }\right)\end{array}$$

by using the definition of mean transit time of Eq. (17). This means that long-term sequestration of one unit of CO₂ converges to the mean transit time of carbon in an ecosystem.

2.5 Dynamic ecosystems out of equilibrium: the continuous sequestration and emissions case

In addition of considering isolated pulses of emissions E₀ or sequestrations S₀, we can also consider permanently ongoing emissions $e:t\mapsto E\left(t\right)$ and sequestration $s:t\mapsto S\left(t\right)$, respectively. Hence,

$$\begin{array}{}\text{(24)}& \mathrm{CS}(T,s,t_{\mathrm{0}}):=\underset{t_{\mathrm{0}}}{\overset{t_{\mathrm{0}}+T}{\int }}{M}_{\mathrm{s}}\left(t\right)\mathrm{d}t,\end{array}$$

where

$$\begin{array}{}\text{(25)}& M_{\mathrm{s}}\left(t\right)=\underset{t_{\mathrm{0}}}{\overset{t}{\int }}\Vert \mathbf{\Phi }(t,\mathit{\tau })\mathit{\beta }\left(\mathit{\tau }\right)s\left(\mathit{\tau }\right)\Vert \mathrm{d}\mathit{\tau }.\end{array}$$

Here s(τ) is a scalar flux of sequestration at time τ. This leads to

$$\begin{array}{}\text{(26)}& r\left(t\right)=-{\mathbf{1}}^{\mathit{⊺}}\mathbf{B}\left(t\right)\underset{t_{\mathrm{0}}}{\overset{t}{\int }}\mathbf{\Phi }(t,\mathit{\tau })\mathit{\beta }\left(\mathit{\tau }\right)s\left(\mathit{\tau }\right)\mathrm{d}\mathit{\tau }.\end{array}$$

The fate of sequestered carbon, for the atmosphere in the form of a balance between simultaneous sequestration and return of carbon, can now be obtained as

$$\begin{array}{}\text{(27)}& \begin{array}{rl}{M}_{\mathrm{a}}^{\prime }\left(t\right)& =-\underset{t_{\mathrm{0}}}{\overset{t}{\int }}{h}_{\mathrm{a}}(t-\mathit{\tau })s\left(\mathit{\tau }\right)\mathrm{d}\mathit{\tau }+\underset{t_{\mathrm{0}}}{\overset{t}{\int }}{h}_{\mathrm{a}}(t-\mathit{\tau })r\left(\mathit{\tau }\right)\mathrm{d}\mathit{\tau }\\ & =-\underset{t_{\mathrm{0}}}{\overset{t}{\int }}{h}_{\mathrm{a}}(t-\mathit{\tau })\left[s\left(\mathit{\tau }\right)-r\left(\mathit{\tau }\right)\right]\mathrm{d}\mathit{\tau }\\ & =-(h_{\mathrm{a}}\star (s-r\left)\right)\left(t\right).\end{array}\end{array}$$

We can now define the climate benefit of sequestration for a dynamic ecosystem with continuous sequestration and respiration as

$$\begin{array}{}\text{(28)}& \begin{array}{rl}\mathrm{CBS}(T,s,t_{\mathrm{0}})& :=\underset{t_{\mathrm{0}}}{\overset{t_{\mathrm{0}}+T}{\int }}{k}_{{\mathrm{CO}}_{\mathrm{2}}}{M}_{\mathrm{a}}^{\prime }\left(t\right)\mathrm{d}t,\\ & =-k_{{\mathrm{CO}}_{\mathrm{2}}}\underset{t_{\mathrm{0}}}{\overset{t_{\mathrm{0}}+T}{\int }}(h_{\mathrm{a}}\star (s-r\left)\right)\left(t\right)\mathrm{d}t.\end{array}\end{array}$$

This expression of CBS accounts for the dynamic behavior of inputs and outputs of carbon in ecosystems, and it can be used to represent time dependencies resulting from environmental changes and disturbances or produced by emission scenarios or scheduled management activities. This time-dependent CBS is computed for a time horizon T starting at any initial time t₀. In other words, it can be used to analyze specific time windows of interest, accounting for the fate of all carbon sequestered during specific time intervals.

3 Example 1: CS and CBS for linear systems in equilibrium

3.1 The fate of a pulse of inputs through the system

A simple ecosystem carbon model, the terrestrial ecosystem model (TECO), will now demonstrate an application of the theory to compute CS and CBS assuming a linear system at a steady state (i.e., in equilibrium). We used a modified version of the TECO model, originally described by Weng and Luo (2011) with parameter values obtained through data assimilation using observations from the Duke Forest in North Carolina, USA. It contains eight main compartments: foliage x₁, woody biomass x₂, fine roots x₃, metabolic litter x₄, structural litter x₅, fast soil organic matter (SOM) x₆, slow SOM x₇, and passive SOM x₈ (Fig. 2). The model represents the dynamics of carbon at a temperate forest dominated by loblolly pine. We chose this model due to its simplicity and tractability, but the framework presented in Sect. 2 can be applied to more complex models and for other ecosystems (see reference in Sect. “Executable research compendium (ERC)” for an example with a nonlinear model). In addition to its simplicity and tractability, there are two advantages of using this model over others: (1) it provides reasonable predictions of net ecosystem carbon fluxes and biometric pool data (Weng and Luo, 2011); (2) it is commonly used to express complex ecosystem-level concepts such as the matrix generalization of carbon cycle models, their traceability, and transient behavior (e.g., Luo and Weng, 2011; Luo et al., 2012; Xia et al., 2013; Luo et al., 2017; Sierra, 2019).

Figure 2 Graphical representation of the terrestrial ecosystem model (TECO) described in Weng and Luo (2011) and Luo et al. (2012). Carbon enters the ecosystem through canopy photosynthesis and is allocated to three biomass pools: foliage, woody biomass, and fine roots. From these pools, carbon is transferred to metabolic and structural litter pools, from where it can be respired as CO₂ or transferred to the soil organic matter (SOM) pools. Blue arrows represent transfers among compartments, and red arrows release to the atmosphere in the form of CO₂.

The model is commonly expressed as

$$\begin{array}{}\text{(29)}& {\displaystyle \frac{\mathrm{d}\mathit{X}\left(t\right)}{\mathrm{d}t}}=\mathit{b}U\left(t\right)+\mathit{\xi }\left(t\right)\mathbf{AC}\mathit{X}\left(t\right),\end{array}$$

where X is a vector of ecosystem carbon pools, C is a diagonal matrix with cycling rates for each pool, A is a matrix of transfer coefficients among pools, and b is a vector of allocation coefficients to plant parts. We modified the entries of matrix A to allow autotrophic respiration to be computed from the vegetation pools and not from the GPP flux as in the original model (see details in Appendix C). The function U(t) determines the carbon inputs to the system as gross primary production (GPP), and ξ(t) is a time-dependent function that modifies ecosystem cycling rates according to changes in the environment.

For this steady-state example, we assume constant inputs (U(t)=U) and constant rates (ξ(t)=1). Furthermore, defining $\mathbf{B}:=\mathbf{AC}$, and u:=bU, we can write this model as a linear, autonomous compartmental system of the form

$$\begin{array}{}\text{(30)}& \dot{\mathit{x}}=\mathit{u}+\mathbf{B}\mathit{x},\end{array}$$

with values for B and u as described in Appendix C.

The fate of a pulse of carbon input entering the ecosystem at an arbitrary time when the system is in equilibrium can be obtained by applying Eqs. (13) and (14) (Fig. 3). Carbon enters the ecosystem through foliage, wood, and fine-root pools. A large proportion of this carbon is quickly transferred from these pools to the fine and metabolic litter pools. Subsequently, the carbon moves to the SOM pools with important respiration losses during these transfers. Most carbon is returned back to the atmosphere with a mean transit time of 30.4 years for the whole system. Half of the sequestered carbon is returned back to the atmosphere in 7.6 years and 95 % in 124 years.

Figure 3 Fate of carbon (M_s(t), left axis; and M_s1(t), right axis) entering the ecosystem according to the TECO model parameterized for the Duke Forest and calculated using Eq. (13) for the upper panel, and respired carbon (r(t)) returning back to the atmosphere calculated using Eq. (15).

Ecosystem-level CS, i.e., the area under the curve of the amount of remaining carbon over time (area under dotted line in Fig. 3, upper panel), increases towards an asymptote as the time horizon of integration increases (Fig. 4a). Here, CS is reported in units of Mg C ha⁻¹ yr, because this is the amount of carbon retained in organic matter over a fixed time horizon. For relevant time horizons of 50, 100, 500, and 1000 years, CS was 233.51, 317.68, 371.64, and 373.42 Mg C ha⁻¹ yr, respectively. In the long term (i.e., as the time horizon goes to infinity), CS converges to the steady-state carbon stock predicted by the model of 373.67 Mg C ha⁻¹.

A similar computation can be made for one unit of fixed carbon (CS₁). In this case CS₁ was 18.98, 25.83, 30.21, and 30.36 years for time horizons of 50, 100, 500, and 1000 years, respectively. In the long term, CS₁ converges to the mean transit time of carbon: 30.4 years (Fig. 4b).

Figure 4 Carbon sequestration (CS) and climate benefit of sequestration (CBS) for instantaneous carbon uptake at any given time. (a) CS due to the uptake of 12.3 Mg C ha⁻¹, which corresponds to GPP of 1 year. (b) CS due to the uptake of one unit of carbon (CS₁). (c) CBS due to the uptake of 12.3 Mg C ha⁻¹ for two different impulse response functions (pre-industrial atmosphere with a pulse of 100 Gt of carbon: PI100, and present-day atmosphere with a pulse of 100 Gt of carbon: PD100). (d) CBS due to the uptake of one unit of carbon (CBS₁) for two different impulse response function. Dotted lines in (a) and (b) represent steady-state carbon storage and mean transit time, respectively.

Figure 5 Absolute global warming potential (AGWP) due to the emission of 1 Mg of CO₂–C to the atmosphere for the two different IRFs (pre-industrial atmosphere with a pulse of 100 Gt of carbon: PI100, and present-day atmosphere with a pulse of 100 Gt of carbon: PD100) reported by Joos et al. (2013).

Due to sequestration at t₀, the CBS shows a rapid negative increase in radiative forcing, which decreases as the time horizon increases due to the return of carbon to the atmosphere as an effect of respiration (Fig. 4c). The shape of the curve, however, depends strongly on the IRF for atmospheric CO₂. CBS is larger over the long term (>200 years) for the present-day (PD100) curve proposed by Joos et al. (2013) than for the pre-industrial curve (PI100). During the pre-industrial period, perturbations of CO₂ in the atmosphere are lower than in the present-day period due to higher uptake of carbon from the oceans and the land biosphere (Joos et al., 2013). Therefore, the benefits of carbon sequestration are larger under present-day conditions based on these IRF curves. Impulse response functions depend strongly on the magnitude and timing of the pulse (Joos et al., 2013; Millar et al., 2017). Therefore, estimates of climate impacts of emissions (AGWP, Fig. 5) and climate benefits of sequestration (CBS, Fig. 4c, d) depend strongly on the choice of the IRF. For the purpose of this paper, we will use the present-day curve (PD100) from here on.

Because AGWP and CBS are based on similar concepts and share similar units, it becomes possible to directly compare one another (Fig. 6) and obtain an estimate of the climate impact of emissions versus sequestration. This can be done either as the ratio of the absolute value of CBS to AGWP, i.e., $\left|\mathrm{CBS}\right|/\mathrm{AGWP}$ (unitless), or as the net radiative balance CBS + AGWP (W m⁻² yr). It is possible to compute these relations using the CBS for one unit of sequestered carbon, which provides a direct estimate of the impact of one unit of sequestration versus one unit of emission, or corresponding to the amount of GPP sequestered in 1 year (12.3 Mg C ha⁻¹ yr⁻¹ for Duke Forest).

In our example, the emission of 1 Mg of carbon to the atmosphere has a predominant warming effect that cannot be compensated for by the sequestration of 1 Mg of carbon at the Duke Forest (Fig. 6). However, the sequestration of the equivalent of GPP in 1 year can have a significant climate benefit compared to the emission of 1 Mg of carbon, depending on the time horizon of analysis. When one integrates in time horizons lower than 200 years, CBS outweighs AGWP in this example. However, because the lifetime of an emission of CO₂ is much longer in the atmosphere than the transit time of carbon through a forest ecosystem, AGWP outweighs CBS on longer timescales.

The time of integration in the computation of GWP has been a heavily debated topic in the past, and this is related to the topic of “permanence” of sequestration in carbon accounting and climate policy (Moura Costa and Wilson, 2000; Noble et al., 2000; Sedjo and Sohngen, 2012). One problem in these previous debates is that the timescale of carbon in ecosystems was not considered explicitly while the timescale of carbon in the atmosphere was. With the approach proposed here, both are explicitly taken into account and can better inform management and policy debates about sequestration of carbon in natural and man-made sinks.

Figure 6 Relations between CBS and AGWP for the IRF PD100 as a function of time horizon T. (a) Ratio between the absolute value of CBS and AGWP, based on a total sequestration of 12.3 Mg of carbon (back line, GPP equivalent for 1 ha and 1 year at Duke Forest) versus a sequestration of 1 Mg of carbon (dashed green line). (b) Radiative balance (net difference) between CBS and AGWP for the sequestration of 12.3 Mg of carbon (black line) and 1 Mg of carbon (dashed green line).

3.2 Carbon management to maximize the climate benefit of carbon sequestration

In the context of climate change mitigation, management of ecosystems may be oriented to increase carbon sequestration and its climate benefit. In the recent past, scientists and policy makers have advocated increasing the amount of inputs to ecosystems as an effective form of carbon management (e.g., Silver et al., 2000; Grace, 2004; Lal, 2004; Chabbi et al., 2017; Minasny et al., 2017). Although increases in carbon inputs can increase the amount of stored carbon in an ecosystem with related climate benefits, it does not necessarily increase the amount of time the sequestered carbon will stay in the system. Therefore, strategies that focus on increasing carbon inputs alone do not take full advantage of the potential of ecosystems to mitigate climate change.

We can conceptualize any management activity that increases or reduces carbon inputs to an ecosystem by a factor γ, so the new inputs are given by the product γ u. For example, if we increase carbon inputs to an ecosystem by 10 %, γ=1.1. Increasing carbon inputs by a proportion γ>1 increases carbon storage at a steady state by an equal proportion since

$$\begin{array}{}\text{(31)}& \begin{array}{rl}-{\mathbf{B}}^{-\mathrm{1}}\left(\mathit{\gamma }\mathit{u}\right)& =\mathit{\gamma }(-{\mathbf{B}}^{-\mathrm{1}}\mathit{u}),\\ & =\mathit{\gamma }{\mathit{x}}^{\ast }.\end{array}\end{array}$$

Similarly, a decrease in carbon inputs by a proportion γ<1 decreases steady-state carbon storage by an equal proportion. However, the time carbon requires to travel through the ecosystem is still the same since the transit time does not change, as we can see from the mean transit time expression

$$\begin{array}{}\text{(32)}& -{\mathbf{1}}^{\mathit{⊺}}{\mathbf{B}}^{-\mathrm{1}}{\displaystyle \frac{\mathit{\gamma }\mathit{u}}{\Vert \mathit{\gamma }\mathit{u}\Vert }}=\mathbb{E}\left(\mathit{\tau }\right).\end{array}$$

Both the transit time distribution (Eqs. B4 and 16) and the mean transit time (Eq. 17) only take into account the proportional distribution of the carbon inputs to the different pools $(\mathit{u}/\Vert \mathit{u}\Vert )$ but not the total amount of inputs. Therefore, a unit of carbon that enters an ecosystem stays there for the same amount of time independent of how much carbon is entering the system. Although these results only apply to linear systems at a steady state, they provide some intuition about what might be the case in systems out of equilibrium.

Carbon management can also be oriented to modify process rates in ecosystems as encoded in the matrix B. A proportional decrease in process rates by a factor ξ<1 not only increases carbon storage as

$$\begin{array}{}\text{(33)}& \begin{array}{rl}-(\mathit{\xi }\mathbf{B}{)}^{-\mathrm{1}}\mathit{u}=& {\displaystyle \frac{\mathrm{1}}{\mathit{\xi }}}(-{\mathbf{B}}^{-\mathrm{1}}\mathit{u}),\\ =& {\displaystyle \frac{{\mathit{x}}^{\ast }}{\mathit{\xi }}},\end{array}\end{array}$$

but also increases the mean transit time as

$$\begin{array}{}\text{(34)}& -{\mathbf{1}}^{\mathit{⊺}}(\mathit{\xi }\mathbf{B}{)}^{-\mathrm{1}}{\displaystyle \frac{\mathit{u}}{\Vert \mathit{u}\Vert }}={\displaystyle \frac{\mathbb{E}\left(\mathit{\tau }\right)}{\mathit{\xi }}}.\end{array}$$

A proportional change in the opposite direction (ξ>1) causes the opposite effect; a proportional increase in process rates decreases carbon storage and decreases mean transit time.

Based on these results, it is now clear that carbon management to increase carbon inputs alone can only increase CS but not CS₁; i.e., the new carbon inputs have a sequestration benefit only through increase of carbon storage but not through a longer transit time in ecosystems. Management to decrease process rates, on the contrary, can increase both CS and CS₁ because the new carbon entering the system stays there for longer.

We can see these effects of carbon management on CS by running simulations using the TECO model at a steady state (Fig. 7). Now, we modified carbon inputs and process rates by either increasing them by 10 % and 50 % (γ, ξ=1.1, 1.5) or decreasing them by 10 % and 50 % (γ, ξ=0.9, 0.5). The simulations showed that increasing or decreasing carbon inputs increase or decrease CS for any time horizon (Fig. 7a), but it does not modify the behavior of one unit of sequestered carbon (CS₁) (Fig. 7b). On the contrary, decreasing or increasing process rates increase or decrease both CS (Fig. 7c) and CS₁ (Fig. 7d).

Figure 7 Different carbon management strategies and their effect on the CS and CS₁. Management to increase or decrease carbon inputs in the vector u by specific proportions γ is shown in panel (a) and (b). Management to increase or decrease process rates in the matrix B by a proportion ξ is shown in panels (c) and (d). Since CS₁ quantifies carbon sequestration of one unit of carbon, management of the amount of carbon inputs does not modify CS₁ in panel (b), and all lines overlap.

The resultant effects of changes in management of inputs or process rates on CBS can differ substantially. Increases or decreases of carbon inputs have similar proportional effects on CBS, but differences in processes rates are not equally proportional. While an increase in inputs by 50 % would increase CBS by 50 %, a decrease in process rates by 50 % would have an increase in CBS by more than 100 % for time horizons longer than 300 years (Fig. 8). Similarly, while a decrease in inputs by 50 % would reduce CBS by 50 %, an increase in process rates by 50 % would decrease CBS by only ∼ 40 %.

Figure 8 Effects of different management strategies on CBS. (a) Effect of increasing or decreasing carbon inputs by a proportion γ on CBS; (b) same effect of γ expressed as a ratio with respect to the reference case of γ=1. (c) Effects of decreasing or increasing process rates in the matrix B by a proportion ξ on CBS; (d) same effect of ξ expressed as a ratio with respect to the reference case ξ=1.

These results show that management of transit time, e.g., by decreasing process rates, may lead to stronger climate benefits than managing carbon inputs alone. Furthermore, one could think about optimization scenarios in which both inputs and transit times are managed to achieve larger climate benefits given certain constraints. The concept of CBS is thus a useful mathematical framework to formally pose such an optimization problem.

We can also use these results to infer differences in CS and CBS for different ecosystem types. Without management, we would expect large variability of CS and CBS in the terrestrial biosphere. Inputs and process rates vary considerably for terrestrial ecosystems as previously reported in other studies. For instance, gross primary productivity can range from about 1 to >30 Mg C ha⁻¹ yr⁻¹ from high- to low-latitude ecosystems (Jung et al., 2020). Based on simulations from the CABLE model, Lu et al. (2018) found a range of mean transit times between 13 and 341 years from low- to high-latitude ecosystems. These large ranges of variability for GPP and mean transit time suggest that CS and CBS may vary among ecosystems by large proportions (>20 times larger or smaller depending on the ecosystems being compared).

4 Example 2: CS and CBS for dynamic systems out of equilibrium

4.1 Pulses entering at different times and experiencing different environments

The steady-state examples above are useful to gain some intuition about potential long-term patterns in CS and CBS, but for real-world applications it is necessary to consider systems out of equilibrium and driven by specific time-dependent signals. We will consider now the case of the temperate ecosystem of our previous example driven by increases in atmospheric CO₂ concentrations that lead to higher photosynthetic uptake and increasing temperatures that lead to faster cycling rates. We will thus consider a non-autonomous version of the TECO model that follows the general form

$$\begin{array}{}\text{(35)}& \dot{\mathit{x}}\left(t\right)=\mathit{\gamma }\left(t\right)\cdot \mathit{u}+\mathit{\xi }\left(t\right)\cdot \mathbf{B}\cdot \mathit{x}\left(t\right),\end{array}$$

where the time-dependent function γ(t) incorporates the effects of temperature and atmospheric CO₂ on primary production, and the function ξ(t) incorporates the effects of temperature on respiration rates. Specific shapes for these functions were taken from Rasmussen et al. (2016) and are described in detail in Appendix C. When applied to the CASA model in Rasmussen et al. (2016), these functions predicted an increase in primary production and an increase in process rates, which resulted in a decrease in transit times over a simulation of 600 years.

We used the same simulation setup here starting from an empty system (x(0)=0) and obtained similar results in terms of primary production and transit times as in Rasmussen et al. (2016). We used these simulation results to compute CS and CBS for carbon entering the ecosystem at different times during the simulation window. In particular, we considered the case of the amount of carbon sequestered at years 100 and 300 after the start of the simulation; i.e., we considered the cases t₀=100 and t₀=300 (Fig. 9a) and computed the fate of this carbon ($M_{\mathrm{s}}(t,t_{\mathrm{0}},{\mathit{u}}_{\mathrm{0}})$), its carbon sequestration (CS$(T,{\mathit{u}}_{\mathrm{0}},t_{\mathrm{0}})$) and the climate benefit of sequestration (CBS$(T,{\mathit{u}}_{\mathrm{0}},t_{\mathrm{0}})$) for different time horizons T.

Figure 9 Prediction of CS and CBS for a non-steady-state case with time-dependent inputs u(t) controlled by CO₂ fertilization and temperature and process rates controlled by temperature modified by a time-dependent factor ξ(t). (a) Predicted time-dependent inputs u(t), and the fate of carbon entering the ecosystem at simulation year 100 ($M_{\mathrm{s}}(t,t_{\mathrm{0}}=\mathrm{100})$) and simulation year 300 ($M_{\mathrm{s}}(t,t_{\mathrm{0}}=\mathrm{300})$). (b) Predicted carbon accumulation in the ecosystem ($\parallel \mathit{x}\left(t\right)\parallel$) for the entire simulation period. (c) Carbon sequestration for the amount of inputs entering at simulation years 100 and 300 calculated for different time horizons T. (d) Climate benefit of sequestration for carbon entering the ecosystem at simulation years 100 and 300 integrated for different time horizons T.

Although more carbon enters the ecosystem at simulation year 300 than at year 100 due to the CO₂ fertilization effect, it is lost much faster because of higher temperatures that result in faster transit times for simulation times above 300 years (Fig. 9a). The slower transit times experienced by the carbon that enters at year 100 due to lower temperature result then in much higher values of CS for time horizons T>100 years (Fig. 9c). Similarly for CBS, where differences are evident much earlier, lower temperatures lead to higher values of CBS for time horizons T>50 years (Fig. 9d).

This simple example highlights the importance of time-dependent transit times in determining CS and CBS. If changes in climate lead to faster carbon processing rates, we would thus expect carbon to transit faster through the ecosystem, returning faster to the atmosphere, and therefore with lower values for carbon sequestration and its climate benefit.

4.2 Continuous inputs into a changing environment

In the previous example, we considered the case of two single pulses entering the ecosystem at different times under changing environmental conditions during a simulation. A consolidated view can be obtained by taking all single pulses and integrating them continuously in time to compute CS and CBS using Eqs. (24) and (28), respectively. In this case, CS increases monotonically, and CBS decreases monotonically with time horizon (Fig. 10, continuous black lines), which is somewhat obvious because as the ecosystem accumulates carbon, more of it is retained in the ecosystem and is isolated from atmospheric radiative effects. However, this simulation only considers carbon that enters the ecosystem from the beginning of the simulation until the end of the time horizon, from t₀ to t₀+T. An important aspect to consider is the role of carbon already present in the ecosystem at t₀.

Figure 10 Computation of CS and CBS for continuous inputs and release of carbon in simulations with different initial conditions x₀: in one simulation the ecosystem develops from empty pools (x(0)=0, i.e., bare ground, black lines), and in the second simulation the ecosystem develops from existing litter and SOM pools but empty biomass pools ($\Vert \mathit{x}\left(\mathrm{0}\right)\Vert =\mathrm{149.04}$ Mg C ha⁻¹, dashed magenta color lines). (a) Inputs u(t) and release fluxes r(t) along the simulation time. (b) Carbon stocks predicted by the model along the simulation time. (c) Carbon sequestration CS for a sequence of time horizons. (d) Climate benefit of sequestration CBS for a sequence of time horizons.

We will consider now the case of continuous sequestration and release of carbon with differences in the initial conditions in the simulation, which can vary according to land use changes. For example, when changing land use from agriculture to forest, or from natural forest to plantation, there are carbon legacies that have an influence on future carbon trajectories (Harmon et al., 1990; Janisch and Harmon, 2002; Sierra et al., 2012). These carbon legacies are usually dead biomass and detritus, which cause ecosystems to lose carbon via decomposition before photosynthesis from new biomass compensates for the losses. In these initial stages of recovery, ecosystems are usually net carbon sources, but they still may store more carbon than an ecosystem developing from bare ground.

The CS and CBS concepts can be very useful to compare contrasting trajectories of ecosystem development and assess their role in terms of carbon sequestration alone and their climate impact. For this purpose, we performed an additional simulation in which at the starting time there is no living biomass, but the detritus pools and the SOM pools are 1.5 and 1.0 times as large as in the equilibrium case, respectively ($\Vert \mathit{x}\left(\mathrm{0}\right)\Vert =\mathrm{149.04}$ Mg C ha⁻¹). In this simulation, the ecosystem losses a significant amount of carbon in the early stages of development, and respiration is much larger than primary production ($r\left(t\right)>\Vert \mathit{u}\left(t\right)\Vert$) (Fig. 10a, dashed magenta line). Because soils are already close to an equilibrium value, the ecosystem already has a large amount of carbon stored; therefore in the computation of the fate of carbon M_s(t,t₀) there is already a larger amount of carbon to consider, which causes CS to be larger for the land-use-change case than for the bare ground case (Fig. 10c). On the contrary, because there are more emissions from the ecosystem in early development stages, CBS is lower for the land-use-change case than for the bare ground case (Fig. 10d).

These contrasting results between CS and CBS for the continuous case with contrasting initial conditions can be very useful to address debates and controversies about the role of land use change and baselines in carbon accounting. The results show that carbon sequestration can still be high in ecosystems where emission fluxes are large, but climate impacts can differ significantly. By using two different metrics, these two different aspects of carbon sequestration can be discussed separately.

5 Discussion

The metrics introduced here, carbon sequestration (CS) and the climate benefit of sequestration (CBS), integrate both the amount of carbon entering an ecosystem and the time it is stored there, thus avoiding radiative effects in the atmosphere. Disproportionate attention is given to quantifying sources and sinks of carbon in ongoing debates about the role of ecosystems in climate change mitigation, with much less attention paid to the fate of carbon once it enters an ecosystem. The time carbon remains in an ecosystem, encapsulated in the concept of transit time, is critical for climate change mitigation because during this time carbon is removed from radiative effects in the atmosphere.

The CS and CBS concepts unify atmospheric and ecosystem approaches to quantifying the greenhouse effect. The CBS concept builds on that of the absolute global warming potential (AGWP) of a greenhouse gas. The main difference is that CBS quantifies avoided warming during the time carbon is stored in an ecosystem, while AGWP quantifies potential warming when the carbon enters the atmosphere. Both metrics rely on the quantification of the fate of carbon (or other GHGs for AGWP) once it enters the particular system. For atmospheric systems, a significant amount of work has been done in determining the fate of GHGs once they enter the atmosphere after emissions (e.g., Rodhe, 1990; O'Neill et al., 1994; Prather, 1996; Archer et al., 2009; Joos et al., 2013). For terrestrial ecosystems; however, robust methods to quantify the fate of carbon as it flows through terrestrial system components have been developed only recently (Rasmussen et al., 2016; Metzler and Sierra, 2018; Metzler et al., 2018).

Global warming potential (GWP), or the climate impact of an emission of a certain gas in relation to the impact of an emission of CO₂, is often used to assess climate impacts of actions, e.g., avoided deforestation, land use change, and even enhanced carbon sequestration. However, this metric has two limitations when applied to carbon sequestration and in comparison to the combined use of CBS and AGWP we advocate here: (1) it only quantifies the climate effects of emissions but not of sequestration and treats all fixed carbon equally independent of its transit time in the ecosystem and (2) it is a relative measure with respect to the emission of CO₂. GWPs are commonly reported in units of CO₂ equivalents, which only address indirectly the effect of a gas in producing warming. In contrast, CBS quantifies the effects of avoided warming in units of W m⁻² over the period of time carbon is retained.

Other concepts have been proposed in the past to account for the temporary nature of carbon sequestration (see review by Brandão et al., 2013, and references therein), with special interest in accounting for credits in carbon markets. In fact, “ton-year” accounting methods (Noble et al., 2000) resemble our definition of carbon sequestration; however, none of these previous concepts explicitly considers the time carbon is retained in the ecosystem. Instead, these approaches relate carbon sequestration to delay in fossil fuel emissions (Fearnside et al., 2000), or as the equivalence of the amount of carbon storage to AGWP (Moura Costa and Wilson, 2000). The concepts of sustained global warming potential (SGWP) and sustained global cooling potential (SGCP) proposed by Neubauer and Megonigal (2015) are notable exceptions. The CBS concept captures some of the ideas of the SGCP concept but differs in some fundamental assumptions related to the interpretation of the impulse response functions, the treatment of time-dependent fluxes and rates, and reporting. While SGCP reports values in reference to CO₂ as is commonly done for GWP, we report CBS for individual gases as it is done for AGWP. Appendix A elaborates on other aspects of the SGWP and SGCP concepts.

The concept of CBS improves our ability to address some of the existing debates about the role of ecosystems in mitigating climate change and enhances our potential to provide decision support. In combination with quantifications of AGWP, CBS provides the net climate effect of an ecosystem or some management. For example, CBS can be used to better understand the climate impacts of storing carbon in long-term reservoirs such as soils and wood products, as well as the climate benefits of increasing the transit time in these systems. CBS can be used to better quantify the climate benefits of using biofuels as fossil fuel substitution by computing the CBS of the whole bioenergy production system and adding the negative AGWP attributed to the avoided emission. Similarly, it can be incorporated in assessments of sequestration in industrial systems with associated carbon capture and storage.

Carbon management of ecosystems can maximize CS and/or CBS by not only increasing carbon inputs, but also by increasing the transit time of carbon. There are many ways in which the transit time of carbon can be increased – for instance, by increasing transfers of carbon to slow cycling pools such as the case of increasing wood harvest allocation to long-duration products (Schulze et al., 2019), or addition of biochar to soils, or by reducing cycling rates of organic matter such as the case of soil flipping (Schiedung et al., 2019). Independently of the management activity, CS and CBS can be powerful metrics to quantify their climate benefits, make comparisons among them, and compare against baselines or no-management scenarios.

The examples we provided in this paper illustrate the use and interpretation of CS and CBS metrics under the assumptions of linearity, steady state, or time dependencies in carbon cycle dynamics with subsequent consequences for carbon sequestration and its climate benefits. The computation of the CBS relies on a model, which can be as simple as a one-pool model or a state-of-the-science land surface model. The TECO model is an excellent tool to illustrate ecosystem-level concepts because of its simplicity and tractability, but other models with more accurate parameterizations and including more processes should be considered for practical applications. The formulas and formal theory developed in Sect. 2 are general enough to deal with the non-steady-state case as well as with models with nonlinear interactions among state variables. In Sierra (2020), we provide an example in the form of a Jupyter Notebook to compute CS and CBS for a nonlinear model (see Sect. “Executable research compendium (ERC)” for details).

The concepts of CS and CBS present improvements to the current guidelines for carbon inventories that treat all carbon removals by sinks equally (IPCC, 2006) by explicitly considering the transit time of carbon in ecosystems. Therefore, these new concepts have potential for being incorporated in revised policies for carbon accounting in the context of international climate agreements and carbon markets. CS and CBS can aid in the economic valuation of carbon by adding economic incentives to sequestration activities that retain carbon in ecosystems for longer times. In addition, the concepts can help in dealing with the issue of permanence of carbon by explicitly quantifying climate benefits of sequestration that can be compared directly with the climate impacts of emissions on a similar time horizon.

Two potential limitations to apply the concepts of CS and CBS are that they rely (1) on a model that tracks the fate of the fixed carbon and (2) on an impulse response function of CO₂ in the atmosphere. Reliable models may not be available for certain types of ecosystems or may include large uncertainties that propagate to CS and CBS estimates. Also, estimates of impulse response functions for atmospheric CO₂ seem to also have uncertainties, particularly related to the size of the emission pulse, the atmospheric background at which the pulse is applied, and the long-term behavior of the curve for timescales longer than 1000 years (Archer et al., 2009; Lashof and Ahuja, 1990; Joos et al., 2013; Millar et al., 2017). However, one advantage of the functions proposed by Joos et al. (2013) is that they are derived from coupled climate–carbon models that include multiple feedbacks. Therefore, when computing CS and CBS for small perturbations of the carbon cycle, it is not necessary to explicitly compute carbon–climate feedbacks. Also, when comparing two different systems with a CBS ratio as in Fig. (8) or a ratio CBS to AGWP (Fig. 6), uncertainties in the IRFs would tend to cancel each other out. Nevertheless, advances in our understanding of the fate of emitted CO₂ to the atmosphere will consequently derive better estimates of the climate benefits of carbon sequestration.

6 Conclusions

Analyses of carbon sequestration for climate change mitigation purposes must consider both the amount of carbon inputs and the transit time of carbon. Both concepts are encapsulated in the unifying concepts of carbon sequestration (CS) and climate benefit of sequestration (CBS) that we propose. Carbon management can be oriented to maximize CS and CBS, which can be achieved by managing both rates of carbon input and process rates in ecosystems. We believe the use of these metrics can help to better deal with current discussions about the role of ecosystems in mitigating climate change, provide better estimates of avoided or human-induced warming, and have the potential to be included in accounting methods for climate policy.

Appendix A: Comment on Neubauer and Megonigal (2015)

Neubauer and Megonigal (2015) proposed two metrics, the sustained global warming potential (SGWP) and the sustained global cooling potential (SGCP), to overcome issues with GWP. However, there is an important misconception in their study that we would like to address here. In particular, these authors state “… GWPs requires the implicit assumption that greenhouse gas emissions occur as a single pulse; this assumption is rarely justified in ecosystem studies”. The use of pulse emissions in computing AGWP, as shown in Eq. (3), is done with the purpose of obtaining a representation of the fate of a unit of emissions under the assumption that the system is in equilibrium. This is a mathematical property of linear time-invariant dynamical systems by which an impulse response function can provide a full characterization of the dynamics of the system (Hespanha, 2009). In other words, the emission pulse is a mathematical method to obtain a description of the fate of incoming mass into the system, but it is not an assumption imposed on the system.

To use impulse response functions, it is necessary to assume that a system is in equilibrium and that all rates remain constant for all times. It is this assumption that is problematic and difficult to impose on ecosystems and not the pulse emission because it is simply a method. Therefore, we are of the opinion that the sustained-flux global warming potential metric proposed by these authors is unjustified on the argument that it removes the assumption of pulse emissions.

One interesting characteristic of the study of Neubauer and Megonigal (2015) is that it uses a model that couples an ecosystem compartment with the atmosphere, and their computation of SGWP and SGCP captures the interactions between these two reservoirs similarly as in the framework described here in Sect. 2. The SGCP is very similar in spirit to the CBS. However, their approach differs from the approach we present here in that our mathematical framework is general enough to deal with ecosystem models of any level of complexity and not restricted to a one-pool model and constant parameters and sequestration rates. Furthermore, we abstain from proposing a metric that is relative to CO₂. We are rather interested in an absolute metric that quantifies the effect of CO₂ sequestration on radiative forcing and not in equivalents to sequestration or emissions of other gases.

Appendix B: Fate and timescales of carbon in compartmental systems

Carbon cycling in the terrestrial biosphere is well characterized by a particular type of dynamical systems called compartmental systems (Anderson, 1983; Jacquez and Simon, 1993). These systems of differential equations generalize mass-balanced models and therefore generalize element and carbon cycling models in ecosystems (Rasmussen et al., 2016; Luo et al., 2017; Sierra et al., 2018a). In their most general form, we can write carbon cycle models as

$$\begin{array}{}\text{(B1)}& {\displaystyle \frac{\mathrm{d}\mathit{x}\left(t\right)}{\mathrm{d}t}}=\dot{\mathit{x}}\left(t\right)=\mathit{u}(\mathit{x},t)+\mathbf{B}(\mathit{x},t)\mathit{x},\end{array}$$

where x(t)∈ℝⁿ is a vector of ecosystem carbon pools, $\mathit{u}(\mathit{x},t)\in {\mathbb{R}}^n$ is a time-dependent vector-valued function of carbon inputs to the system, and $\mathbf{B}(\mathit{x},t)\in {\mathbb{R}}^{n\times n}$ is a time-dependent compartmental matrix. The latter two terms can depend on the vector of states, in which case the compartmental system is considered nonlinear. In case the input vector and the compartmental matrix have fixed coefficients (no time dependencies), the system is considered autonomous, and it is considered non-autonomous otherwise (Sierra et al., 2018a). At a steady state, the autonomous linear system has the general solution ${\mathit{x}}^{\ast }=-{\mathbf{B}}^{-\mathrm{1}}\mathit{u}$.

The probability density function (pdf) for system age of linear autonomous models at a steady state can be computed by the following expression (Metzler and Sierra, 2018)

$$\begin{array}{}\text{(B2)}& f\left(a\right)=-{\mathbf{1}}^{\mathit{⊺}}\mathbf{B}{e}^{a\cdot \mathbf{B}}{\displaystyle \frac{{\mathit{x}}^{\ast }}{\Vert {\mathit{x}}^{\ast }\Vert }},a\ge \mathrm{0},\end{array}$$

where a is the random variable age, 1^⊺ is the transpose of the n-dimensional vector containing ones, e^a⋅B is the matrix exponential computed for each value of a, and $\Vert {\mathit{x}}^{\ast }\Vert$ is the sum of the stocks of all pools at a steady state.

The mean, i.e., the expected value, of the age pdf can be computed by the expression

$$\begin{array}{}\text{(B3)}& \mathbb{E}\left(a\right)=-{\mathbf{1}}^{\mathit{⊺}}{\mathbf{B}}^{-\mathrm{1}}{\displaystyle \frac{{\mathit{x}}^{\ast }}{\Vert {\mathit{x}}^{\ast }\Vert }}={\displaystyle \frac{\Vert {\mathbf{B}}^{-\mathrm{1}}{\mathit{x}}^{\ast }\Vert }{\Vert {\mathit{x}}^{\ast }\Vert }}.\end{array}$$

The pdf of the transit time variable τ for linear autonomous systems in equilibrium is given by (Metzler and Sierra, 2018)

$$\begin{array}{}\text{(B4)}& f\left(\mathit{\tau }\right)=-{\mathbf{1}}^{\mathit{⊺}}\mathbf{B}{e}^{\mathit{\tau }\cdot \mathbf{B}}{\displaystyle \frac{\mathit{u}}{\Vert \mathit{u}\Vert }},\mathit{\tau }\ge \mathrm{0}\end{array}$$

and the mean transit time by

$$\begin{array}{}\text{(B5)}& \mathbb{E}\left(\mathit{\tau }\right)=-{\mathbf{1}}^{\mathit{⊺}}{\mathbf{B}}^{-\mathrm{1}}{\displaystyle \frac{\mathit{u}}{\Vert \mathit{u}\Vert }}={\displaystyle \frac{\Vert {\mathit{x}}^{\ast }\Vert }{\Vert \mathit{u}\Vert }}.\end{array}$$

For the most general case of nonlinear non-autonomous systems, we follow the approach described in Metzler et al. (2018). For these systems, the age distribution of mass is given by

$$\begin{array}{l}\text{Mass in the}\\ \text{system at}\\ \text{time}t\text{with}\\ \text{age}a\end{array}=\left\{\begin{array}{ll}\mathbf{\Phi }(t,t-a)\mathit{u}(t-a),& a<t-t_{\mathrm{0}},\\ \mathbf{\Phi }(t,t_{\mathrm{0}}){\mathit{f}}^{\mathrm{0}}(a-(t-t_{\mathrm{0}}\left)\right),& a\ge t-t_{\mathrm{0}},\end{array}\right.$$

where Φ is a state-transition matrix, and f⁰ is an initial age density distribution at initial time t₀. We obtain Φ by taking advantage of an existing numerical solution x(t), which we plug into the original system, obtaining a new compartmental matrix $\stackrel{\mathrm{̃}}{\mathbf{B}}\left(t\right):=\mathbf{B}\left(\mathit{x}\right(t),t)$ and a new input vector $\stackrel{\mathrm{̃}}{\mathit{u}}:=\mathit{u}\left(\mathit{x}\right(t),t)$. Then, the new linear non-autonomous compartmental system,

$$\begin{array}{}\text{(B6)}& \dot{\mathit{y}}\left(t\right)=\stackrel{\mathrm{̃}}{\mathbf{B}}\left(t\right)\mathit{y}\left(t\right)+\stackrel{\mathrm{̃}}{\mathit{u}}\left(t\right),t>t_{\mathrm{0}},\end{array}$$

has the unique solution y(t)=x(t), which emerges from the fact that both systems are identical. The solution of the system is then given by

$$\begin{array}{}\text{(B7)}& \mathit{x}\left(t\right)=\mathbf{\Phi }(t,t_{\mathrm{0}}){\mathit{x}}^{\mathrm{0}}+\underset{t_{\mathrm{0}}}{\overset{t}{\int }}\mathbf{\Phi }(t,s)\mathit{u}\left(s\right)\mathrm{d}s,\end{array}$$

where ${\mathit{x}}^{\mathrm{0}}={\int }_{\mathrm{0}}^{\mathrm{\infty }}{\mathit{f}}^{\mathrm{0}}\left(a\right)\mathrm{d}a$ is the initial vector of carbon stocks. We obtain the state-transition matrix as the solution of the following matrix differential equation:

$$\begin{array}{}\text{(B8)}& {\displaystyle \frac{\mathbf{\Phi }(t,t_{\mathrm{0}})}{\mathrm{d}t}}=\mathbf{B}\left(t\right)\mathbf{\Phi }(t,t_{\mathrm{0}}),t>t_{\mathrm{0}},\end{array}$$

with initial condition

$$\begin{array}{}\text{(B9)}& \mathbf{\Phi }(t_{\mathrm{0}},t_{\mathrm{0}})=\mathbf{I},\end{array}$$

where $\mathbf{I}\in {\mathbb{R}}^{n\times n}$ is the identity matrix. For the special case in which the time-dependent matrix can be expressed as a product between a time-dependent scalar factor ξ(t) and a constant value matrix B, i.e., B(t)=ξ(t)B, we obtain the state-transition matrix as

$$\begin{array}{}\text{(B10)}& \mathbf{\Phi }(t,t_{\mathrm{0}})=\exp \left(\underset{t_{\mathrm{0}}}{\overset{t}{\int }}\mathit{\xi }\left(\mathit{\tau }\right)\mathrm{d}\mathit{\tau }\cdot \mathbf{B}\right).\end{array}$$

These formulas can be applied to any carbon cycle model represented as a compartmental system to obtain the fate of carbon once it enters the ecosystem as well as timescale metrics such as age and transit time distributions.

Computation of the mass remaining in the system

From Eq. (B7), we can see from the first term that the initial amount of carbon in the system x⁰ changes over time according to the term Φ(t,t₀) x⁰. Rasmussen et al. (2016) showed that, under certain circumstances, Eq. (B7) is exponentially stable as long as B is invertible, and the state-transition operator acts as a term that exponentially “decomposes” the initial amount of carbon. Furthermore, the state-transition operator tracks the dynamics of the incoming carbon and how it is transferred among the different pools before it is respired. Therefore, this operator can be used to compute the fate of an amount of carbon sequestered at time t_s as

$$\begin{array}{}\text{(B11)}& M_{\mathrm{s}}(t-t_{\mathrm{s}})=M_{\mathrm{s}}\left(a\right)=\Vert \mathbf{\Phi }(t,t_{\mathrm{s}})\mathit{u}\left(t_{\mathrm{s}}\right)\Vert ,a=t-t_{\mathrm{s}}.\end{array}$$

Similarly, the fate of one unit of sequestered carbon at time t_s can be computed as

$$\begin{array}{}\text{(B12)}& M_{\mathrm{s}\mathrm{1}}\left(a\right)=∥\mathbf{\Phi }(t,t_{\mathrm{s}})\cdot {\displaystyle \frac{\mathit{u}\left(t_{\mathrm{s}}\right)}{\Vert \mathit{u}\left(t_{\mathrm{s}}\right)\Vert }}∥,\end{array}$$

where the subscript 1 denotes that the function predicts the fate of one unit of carbon.

Appendix C: Detailed representation of the TECO model and the transient simulations used in examples

The terrestrial ecosystem model (TECO) described in Weng and Luo (2011) and Luo et al. (2012) has eight pools to simulate ecosystem-level carbon dynamics, with a parameterization for the Duke Forest, a temperate forest in North Carolina, USA. The annual amount of photosynthetically fixed carbon predicted by the model in this forest (GPP) is U=12.3 Mg C ha⁻¹ yr⁻¹. The vector of carbon allocation is given by

$$\mathit{b}=\left(\begin{array}{c}\mathrm{0.14}\\ \mathrm{0.26}\\ \mathrm{0.14}\\ \mathrm{0}\\ \mathrm{0}\\ \mathrm{0}\\ \mathrm{0}\\ \mathrm{0}\end{array}\right),$$

which shows that from all photosynthetically fixed carbon, 14 % is allocated to foliage, 26 % to woody biomass, and 14 % to roots.

Each pool in the model cycles at annual rates given by the diagonal elements of the matrix

$$\mathbf{C}=\left(\begin{array}{cccccccc}\mathrm{0.942}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}\\ \mathrm{0.000}& \mathrm{0.021}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}\\ \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.872}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}\\ \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{3.978}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}\\ \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.347}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}\\ \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{3.833}& \mathrm{0.000}& \mathrm{0.000}\\ \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.036}& \mathrm{0.000}\\ \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.000}& \mathrm{0.004}\end{array}\right),$$

with a matrix of transfer coefficients as

$$\mathbf{A}=\left(\begin{array}{cccccccc}-\mathrm{1.00}& \mathrm{0.00}& \mathrm{0.00}& \mathrm{0.00}& \mathrm{0.00}& \mathrm{0.00}& \mathrm{0.00}& \mathrm{0.00}\\ \mathrm{0.00}& -\mathrm{1.00}& \mathrm{0.00}& \mathrm{0.00}& \mathrm{0.00}& \mathrm{0.00}& \mathrm{0.00}& \mathrm{0.00}\\ \mathrm{0.00}& \mathrm{0.00}& -\mathrm{1.00}& \mathrm{0.00}& \mathrm{0.00}& \mathrm{0.00}& \mathrm{0.00}& \mathrm{0.00}\\ \mathrm{0.82}& \mathrm{0.00}& \mathrm{0.12}& -\mathrm{1.00}& \mathrm{0.00}& \mathrm{0.00}& \mathrm{0.00}& \mathrm{0.00}\\ \mathrm{0.02}& \mathrm{0.85}& \mathrm{0.72}& \mathrm{0.00}& -\mathrm{1.00}& \mathrm{0.00}& \mathrm{0.00}& \mathrm{0.00}\\ \mathrm{0.00}& \mathrm{0.00}& \mathrm{0.00}& \mathrm{0.45}& \mathrm{0.28}& -\mathrm{1.00}& \mathrm{0.42}& \mathrm{0.45}\\ \mathrm{0.00}& \mathrm{0.00}& \mathrm{0.00}& \mathrm{0.00}& \mathrm{0.28}& \mathrm{0.30}& -\mathrm{1.00}& \mathrm{0.00}\\ \mathrm{0.00}& \mathrm{0.00}& \mathrm{0.00}& \mathrm{0.00}& \mathrm{0.00}& \mathrm{0.00}& \mathrm{0.01}& -\mathrm{1.00}\end{array}\right).$$

Matrix A was modified from the original publication (Luo et al., 2012) by decreasing the proportion that is transferred from vegetation to litter pools in a proportion consistent with the proportions of carbon that are respired by autotrophic respiration. In other words, autotrophic respiration is not computed here as in the original publication where the inputs enter the ecosystem in the form of net primary production, i.e., U = GPP − R_a. We compute autotrophic respiration as the proportion that leaves the autotrophic pools and that is not transferred to the litter pools. In this way, U = GPP, and all carbon that is fixed enters the vegetation pools from where it is subsequently respired or added to the litter pools.

Defining B:=AC and u=bU, we obtained the steady-state solution as

$$\begin{array}{}\text{(C1)}& \begin{array}{rl}{\mathit{x}}^{\ast }& =-{\mathbf{B}}^{-\mathrm{1}}\cdot \mathit{u}\\ & =\left(\begin{array}{c}\mathrm{3.83}\\ \mathrm{237.70}\\ \mathrm{4.14}\\ \mathrm{0.86}\\ \mathrm{20.18}\\ \mathrm{1.28}\\ \mathrm{92.96}\\ \mathrm{12.72}\end{array}\right).\end{array}\end{array}$$

For the simulation with initial conditions as in a land-use-change case, the initial conditions x₀ of the simulation were set as

$$\begin{array}{}\text{(C2)}& {\mathit{x}}_{\mathrm{0}}={\left(\begin{array}{cccccccc}\mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1.5}& \mathrm{1.5}& \mathrm{1.0}& \mathrm{1.0}& \mathrm{1.0}\end{array}\right)}^{\mathit{⊺}}\circ {\mathit{x}}^{\ast },\end{array}$$

where the symbol ∘ represents entry-wise multiplication.

For the transient simulations, we derived time-dependent modifiers for inputs γ(t) and for process rates ξ(t) following the approach described in Rasmussen et al. (2016). Atmospheric CO₂ concentrations increase following a sigmoid curve given by

$$\begin{array}{}\text{(C3)}& x_{\mathrm{a}}\left(t\right)=\mathrm{284}+\mathrm{1715}\exp \left({\displaystyle \frac{\mathrm{0.0305}t}{(\mathrm{1715}+\exp (\mathrm{0.0305}t)-\mathrm{1})}}\right),\end{array}$$

and surface air temperature increases with CO₂ concentrations according to

$$\begin{array}{}\text{(C4)}& T_{\mathrm{s}}\left(t\right)=T_{\mathrm{s}\mathrm{0}}+{\displaystyle \frac{\mathit{\sigma }}{\ln \left(\mathrm{2}\right)}}\ln \left(x_{\mathrm{a}}\right(t)/\mathrm{285}).\end{array}$$

The combined effect of CO₂ concentrations and air surface temperature on primary production is then computed as

$$\begin{array}{}\text{(C5)}& \mathit{\gamma }\left(t\right)=(\mathrm{1}+\mathit{\beta }(x_{\mathrm{a}}\left(t\right),T_{\mathrm{s}}\left(t\right)\left)\ln \right(x_{\mathrm{a}}\left(t\right)/\mathrm{285}\left)\right),\end{array}$$

with

$$\begin{array}{}\text{(C6)}& \begin{array}{rl}& \mathit{\beta }\left(x_{\mathrm{a}}\right(t),T_{\mathrm{s}}(t\left)\right)=\\ & {\displaystyle \frac{\mathrm{3}\mathit{\rho }{x}_{\mathrm{a}}\left(t\right)\mathrm{\Gamma }\left(T_{\mathrm{s}}\right(t\left)\right)}{\left(\mathit{\rho }{x}_{\mathrm{a}}\right(t)-\mathrm{\Gamma }(T_{\mathrm{s}}\left(t\right)\left)\right)\left(\mathit{\rho }{x}_{\mathrm{a}}\right(t)+\mathrm{2}\mathrm{\Gamma }(T_{\mathrm{s}}\left(t\right)\left)\right)}},\end{array}\end{array}$$

where β(x_a(t),T_s(t)) is the sensitivity of primary production with respect to atmospheric CO₂ and air surface temperature, and ρ=0.65 is the ratio of intracellular CO₂ to x_a(t). The response function with respect to temperature Γ(T_s(t)) is given by

$$\begin{array}{}\text{(C7)}& \mathrm{\Gamma }\left(T_{\mathrm{s}}\right(t\left)\right)=\mathrm{42.7}+\mathrm{1.68}\left(T_{\mathrm{s}}\right(t)-\mathrm{25})+\mathrm{0.012}\left(T_{\mathrm{s}}\right(t)-\mathrm{25}{)}^{\mathrm{2}}.\end{array}$$

The separate effect of air surface temperatures on process rates is computed with a power function of the form

$$\begin{array}{}\text{(C8)}& \mathit{\xi }\left(T_{\mathrm{s}}\right(t\left)\right)={\mathit{\xi }}_b^{\mathrm{0.1}{T}_{\mathrm{s}}\left(t\right)-\mathrm{1.5}},\end{array}$$

with ξ_b=2.