Reservoir Routing and Its Application to Atmospheric Carbon Dioxide Balance

1. Introduction

Thanks to their links with real-world problems and their engineering approach in seeking rational solutions to them, hydrology and hydraulics have ever been in close contact with reality. This prevented them from practicing hypocrisy and vagueness, similar to the mathematics (cf. the motto above). Historically, the solutions to real-world problems preceded the development of hydrology and hydraulics as sciences, and, thus, could not be based on scientific knowledge. Rather, they were based on common sense, which historically has provided the foundation of philosophical and scientific knowledge, while in recent decades is being abandoned if favour of its opposite, which euphemistically has been called “political correctness”.

Among the solutions to real-world problems, the technology of storage, implemented by cisterns at small scale and reservoirs at dammed rivers at large scale, has been most determinant. Perhaps the first dam with a reservoir storing water is that in Jawa, in the desert of Northeastern Jordan, a 5 m high and 80 m long dam dated between 3500 and 3400 BC [2]. Another prehistoric example is the so-called Sadd el Kafara dam in Wadi Garawi, about 50 km southeast of Cairo, dated, into the old kingdom of Egypt, around 2650 BC [3].

Apparently, storage facilities for surpluses of goods, in particular grains, to be consumed in periods of deficit, preceded that of water. Granaries date as far back as 9500 BC in the Jordan Valley, coinciding with the dawn of agriculture, and they were initially at the household scale, later expanding to larger scales for collective storage of cereal [4].

We may speculate that this type of storage, whether for water or for cereal, was for seasonal (intra-annual) regulation and is far different from modern facilities for overyear regulation. Yet we have written information of the latter type of regulation, which is necessary to deal with famines, in the biblical story of Joseph and the pharaoh’s dream of the seven fat and the seven skinny cows. The famine, which was predicted by Joseph by interpreting pharaoh’s dream, indeed occurred, but Egypt was prepared for it:

The famine became bad everywhere in Egypt, so Joseph opened the storehouses and sold the grain to the Egyptians. People from all over the world came to Egypt to buy grain, because the famine was so severe in their countries. (Καὶ ὁ λιμὸς ἦν ἐπὶ προσώπου πάσης τῆς γῆς· ἀνέῳξε δὲ Ἰωσὴφ πάντας τοὺς σιτοβολῶνας, καὶ ἐπώλει πᾶσι τοῖς Aἰγυπτίοις. Καὶ πᾶσαι αἱ χῶραι ἦλθον εἰς Aἴγυπτον, ἀγοράζειν πρὸς Ἰωσήφ· ἐπεκράτησε γὰρ ὁ λιμὸς ἐν πάσῃ τῇ γῇ).[5]

Bell [6] and Said [7] confirm this biblical story and date it to an uncertain number of years preceding 1740 BC. Specifically, Bell [6] links this story with an inscription in the tomb of Sobek-nakht, which was decorated in the three-year reign of Sobekhotep III, just prior to ca. 1740, and reads:

I was a man who protected the afflicted against the powerful […] who supplied the granaries of the god […] who summoned his entire energy every time he saw an insufficient flood.

She also refers to another (undated) tomb inscription which reads:

I gave grain to the entire country, I saved my town from famine […] no one has done what I did.

This biblical story is of great hydrological and climatic importance as it is perhaps the oldest text referring to a long-lasting drought, as well as to the natural behavior of clustering in time of similar events (such as dry or wet years), which in recent times was quantified by Hurst [8] in 1951. Mandelbrot and Wallis [9] used the term Joseph effect for this behavior, which today is more often called Hurst-Kolmogorov dynamics [10]. The importance of the story extends to good management practices, in which the storage of goods (in this case in granaries) during periods of abundance can mitigate the adverse consequences in periods of shortage.

Yet storage is not a human invention, because Nature uses it extensively. Numerous natural processes include storage and their modelling needs to properly represent it. On the other hand, the need for design and management of human storage facilities, particularly water reservoirs, enabled deeper insights and more concrete knowledge on modelling storage and the related processes. This modelling typically follows a systems approach, with an input (inflow), an output (outflow) and a state variable (stored quantity). These three quantities are interlinked with a first principle, i.e. mass conservation, implemented as a linear differential or a difference equation. This, however, does not suffice for complete modelling. One more equation is needed. In detailed physical modelling this can be offered by another physical principle, such as conservation of momentum, but in macroscopic modelling this is not possible and the second equation is macroscopic with some conceptual or statistical meaning, depending on the nature of the system. In particular, if this second equation is a linear relationship between inflow and storage, then there exists an analytical solution of the resulting differential equation, which can directly be utilized to find the outflow from known inflow, and to obtain macroscopic characteristics of the process, such as response and residence times, and their distribution functions.

Here we revisit and streamline the reservoir routing framework, and extend it to find approximate solutions for nonlinear cases (Section 2). Furthermore, we apply this framework to a topical issue of climatic interest, the mass balance of atmospheric carbon dioxide (Section 2). The framework allows determining characteristic residence times, which have been an issue of controversy. Application of the theoretical framework results in excellent agreement with real world data. In this manner, we easily quantify the atmospheric carbon exchanges, and obtain reliable and intuitive results, without the need to resort to complex climate models.

2. Theoretical Analysis

2.1 System Components and Determination of Their Temporal Evolution

We consider a system that can be represented as a reservoir, with input $I(t)$, output $Q(t)$, and storage S$(t)$. These three quantities are connected by the continuity equation, which expresses the conservation of mass and is written in differential form as:

$$\frac{dS(t)}{dt}+Q(t)=I(t)$$

(1)

In hydrosystems, the inflow is typically the upstream flow discharge and the outflow is the downstream flow discharge or a reservoir spill, and are usually expressed in terms of volumes instead of masses (but unless we deal with uncompressible fluids, we should stick with mass units). The reservoir could be a real one, upstream of a dam, but the use of imaginary reservoirs is not uncommon. For example, an entire catchment is often represented as a single imaginary reservoir with the inflow being the precipitation [10,11] or as a cascade of reservoirs [12], in which outputs from upstream reservoirs are inputs to the downstream reservoirs in the cascade. Similar representations may be appropriate for groundwater stored in aquifers [13].

Storage, however, is not limited to locally resolved processes of surface- or ground-water. The global hydrological cycle can be viewed in terms of mass exchange with the main storage being the oceans, but with the atmosphere having a crucial role when evaporation and precipitation are studied [14]. Modelling of the carbon cycle can also be made in a similar setting, as will be shown in section 3.

The typical problem is to determine the outflow $Q(t)$ for known inflow $I(t)$. To fully describe the transformation of inflow to outflow, we need one more equation. The solution of the two equations is commonly known as reservoir routing. In hydrology and hydraulics we usually construct the second equation by means of a stage-discharge and a stage-storage relationship. Different equations can be formulated, depending on the system dynamics, but the following power type (combined) equation is representative for most problems:

$$Q\left(S\right)=Q_0{\left(\frac{S}{S_0}\right)}^b$$

(2)

where $b$ is a dimensionless parameter (exponent) and $S_0$ and $Q_0$ are parameters with units of mass and mass flow, respectively. These are necessary to make the equation physically (dimensionally) consistent and here will be regarded to be the initial conditions (at time zero, i.e. $I(0)$, and $Q(0)$, respectively). Notice that these parameters are not unique; for example we can multiply $S_0$ and $Q_0$ by $c$ and $c^{-b}$, respectively, and get another pair of valid parameters.

It is well known that when $b=1$, we get a first-order linear differential equation with constant coefficients, i.e.,

$$\frac{dS(t)}{dt}+\frac{1}{W_0}S\left(t\right)=I\left(t\right),W_0:=\frac{S_0}{Q_0}$$

(3)

where $W_0$ is a characteristic time, whose meaning will be discussed below. This admits a closed solution:

$$S\left(t\right)=S_0{e}^{-t/W_0}+\frac{1}{W_0}{\int }_0^t{e}^{-(t-s)/W_0}I\left(s\right)ds, Q\left(t\right)=Q_0{e}^{-t/W_0}+\frac{1}{W_0}{\int }_0^t{e}^{-(t-s)/W_0}I\left(s\right)ds$$

(4)

This case ($b=1$) is known as the linear reservoir [15,16,17,18,19,20] and the fact that it admits a general analytical solution makes it a useful tool. We will call a reservoir with $b<1$ a sublinear reservoir (as the discharge function $Q\left(S\right)$ is sublinear, meaning that $\underset{S\to \infty }{\lim }Q\left(S\right)/S=0$). Conversely, we will call a reservoir with $b>1$ a superlinear reservoir. For both sublinear and superlinear reservoirs no general analytical solution can be found, except for very few special cases (e.g. [21,22,23]), some of which will be discussed below and in Appendix A. Here we will seek general approximate solutions for $b\ne 1$.

As a first step, we standardize the equation in the following form, in which all variables are dimensionless:

$$\frac{\mathrm{d}s\left(\tau \right)}{\mathrm{d}\tau }+{\left(s\left(\tau \right)\right)}^b=\frac{i(\tau )}{q_0}$$

(5)

with initial conditions

$$s\left(0\right)=q\left(0\right)=1$$

(6)

where

$$\tau =\frac{t}{W_0}, s\left(\tau \right)=\frac{S\left(\tau W_0\right)}{S_0}, q\left(\tau \right)=\frac{Q\left(\tau W_0\right)}{Q_0}, i\left(\tau \right)=\frac{I\left(\tau W_0\right)}{I_0}, q_0=\frac{Q_0}{I_0}$$

(7)

with $S_0, I_0$, $Q_0$ denoting the values of the respective variable at time $t=\tau =0$. Notice the use of lower- and upper-case symbols for dimensionless and dimensional quantities, respectively, with the exception of dimensional time, where we use the common symbol t, while we use the Greek symbol τ for dimensionless time.

Now we make a first order linear approximation of ${\left(s\left(\tau \right)\right)}^b$ at $s=1$ (similar to Basha [24]), as

$${\left(s\left(\tau \right)\right)}^b=bs\left(\tau \right)+1-b+O\left({\left(s-1\right)}^2\right)$$

(8)

and neglecting the high order terms we get

$$\frac{\mathrm{d}s\left(\tau \right)}{\mathrm{d}\tau }+bs\left(\tau \right)=\frac{i(\tau )}{q_0}+b-1$$

(9)

whose solution is

$$s\left(\tau \right)=e^{-\tau }+{\int }_0^{\tau }{e}^{-\left(\tau -s\right)}\left(\frac{i(\tau )}{q_0}+b-1\right)ds, q\left(\tau \right)=b{e}^{-\tau }-b+1+b{\int }_0^{\tau }{e}^{-\left(\tau -s\right)}\left(\frac{i(\tau )}{q_0}+b-1\right)ds$$

(10)

If the dimensionless inflow is constant (i.e. $i\left(\tau \right)=1$, or dimensional inflow $I\left(t\right)=I_0=Q_0/q_0$), this simplifies to

$$s\left(\tau \right)=1-\frac{1}{b}\left(1-\frac{1}{q_0}\right)\left(1-e^{-b\tau }\right), q\left(\tau \right)=1-\left(1-\frac{1}{q_0}\right)(1-e^{-b\tau })=\frac{1}{q_0}+\left(1-\frac{1}{q_0}\right)e^{-b\tau }$$

(11)

where, as $t\to \infty$, the dimensionless outflow becomes $1/q_0$ and the dimensional output becomes $Q_0 (1/q_0)=I_0$.

It is important to note that, despite being an approximation, this preserves the mass balance precisely. Indeed, the time integral of outflow minus inflow is

$$\underset{0}{\overset{t}{\int }}{Q}_{0}q\left(\frac{u}{W_0}\right)du-I_{0}t=\frac{\left(q_0-1\right)Q_0{W}_0}{b{q}_0}\left(1-e^{-\frac{bt}{W_0}}\right)$$

(12)

while the change in storage is

$$S_0-S_{0}s\left(\frac{t}{W_0}\right)=S_0\frac{1}{b}\left(1-\frac{1}{q_0}\right)\left(1-e^{-\frac{bt}{W_0}}\right)$$

(13)

and upon substituting $W_0{Q}_0\mathrm{for}{S}_0$ we confirm that the right-hand sides in the above two equations are identical.

It is readily seen that if $q_0=1$, then $s\left(\tau \right)=q\left(\tau \right)=1$, which describes a steady state flow. In this particular case, the solution is exact. When the inlow is zero, a case that is represented as $q_0=\infty$, the solution becomes

$$s\left(\tau \right)=1-\frac{1}{b}\left(1-e^{-b\tau }\right), q\left(\tau \right)=e^{-b\tau }$$

(14)

but it is not exact; the exact one will be discussed below. For $b>1$ and $\tau \to \infty ,$this results in $s(\infty )=1-1/b>0$, meaning that the reservoir never empties. This looks absurd, but it is consistent with the linear approximation made, as seen in Equation (8), according to which no outflow occurs for $s(t)\le 1-1/b$. For $b<1$, the resevoir empties at $\tau =-\ln (1-b)/b$ and beyond that time $q\left(\tau \right)=0$.

For inflow linearly varying with time, i.e. $i\left(\tau \right)=1+a\tau$, the solution becomes

$$s\left(\tau \right)=1+\frac{a\tau }{b{q}_0}-\frac{1}{b}\left(1-\frac{1}{q_0}+\frac{a}{b{q}_0}\right)\left(1-{\mathrm{e}}^{-b\tau }\right), q\left(\tau \right)=\frac{1}{q_0}-\frac{a}{b{q}_0}+\frac{a\tau }{q_0}+\left(1-\frac{1}{q_0}+\frac{a}{b{q}_0}\right){\mathrm{e}}^{-b\tau }$$

(15)

where the rightmost equation can also be written as

$$q\left(\tau \right)=\frac{i\left(\tau -1/b\right)}{q_0}+\left(1-\frac{i\left(-1/b\right)}{q_0}\right){\mathrm{e}}^{-b\tau }$$

(16)

If $a\ge 0$, the solution is valid for any τ; if $a<0$, it is valid for $\tau \le -1/a$.

Some special cases of $b$ and $i\left(\tau \right)$ also admit exact analytical solutions, like in the linear reservoir problem. A first case is when the inflow is zeroed at time $\tau =0$, and continues to be zero ($i(\tau )=0$) at any $\tau >0$. In this case the solution of the differential equation (5), which is now homogenous, with initial conditions (6) is easily found to be

$${s\left(\tau \right)=\left(\left(b-1\right)\tau +1\right)}^{\frac{1}{1-b}}, {q\left(\tau \right)=\left(\left(b-1\right)\tau +1\right)}^{\frac{b}{1-b}}$$

(17)

and it apparently holds for $\left(1-b\right)\tau \le 1$, that is, for any $\tau$ if $b\ge 1$, but for $\tau \le 1/\left(1-b\right)$ if $b<1$; beyond that time $s\left(\tau \right)=q\left(\tau \right)=0$. As$b\to 1$, we have the limiting case:

$$s\left(\tau \right)=q\left(\tau \right)=e^{-\tau }$$

(18)

We also examine the following two special cases, which are useful as benchmarks. The first benchmark case is superlinear, $q\left(\tau \right)={\left(s\left(\tau \right)\right)}^2,i\left(\tau \right)=1$, for which:

$$q\left(\tau \right)=\frac{1}{q_0}{\left(1-\frac{2\left(1-\sqrt{q_0}\right)}{\left(1+\sqrt{q_0}\right)e^{\frac{2\tau }{\sqrt{q_0}}}+(1-\sqrt{q_0})}\right)}^2, s\left(\tau \right)=\sqrt{q(\tau )}$$

(19)

The second is sublinear, $q\left(\tau \right)={\left(s\left(\tau \right)\right)}^{1/2},i\left(\tau \right)=1$, with

$$q\left(t\right)=\frac{1}{q_0} \left(\mathrm{W}\left(\left(q_0-1\right){\mathrm{e}}^{-\frac{q_0\tau }{2}+q_0-1}\right)+1\right), s\left(\tau \right)={\left(q\left(\tau \right)\right)}^2$$

(20)

where $\mathrm{W}\left(z\right)$ denotes the Lambert W function, defined as the inverse of the function $z=f\left(w\right)≔w{\mathrm{e}}^w$. The case $q\left(\tau \right)={\left(s\left(\tau \right)\right)}^2$ also has an analytical solution for changing inflow as a linear function but this is too complicated to write down. It is thus preferable to use numerical integration, which can work in every case, yet the analytical approximations are quite useful, as will be seen in the next sections. Additional cases of analytical solutions and different approximations are studied in Appendix A.1.

The good performance of the approximate solution is illustrated in Figure 1 for constant inflow and in Figure 2 for input linearly changing with time. From these figures we infer that, as long as the exponent $b$ is between 1/2 and 2, and the initial dimensionless output $q_0$ between 0.8 and 1.25, the linear approximation is satisfactory. In practical tasks these ranges cover the most typical applications. An exception is the determination of the response time (see section 2.3), in which the input is zero (and hence $q_0=\infty$) but for this case we have an exact solution (Equation (17)).

2.2 Residence Time

The above mathematical framework allows us to determine characteristic times, such as residence and response times and probabilities thereof. In this subsection we deal with the former, starting with its definition, and in the next one with the latter.

Definition 1: 

The residence time, W, is defined to be the time duration that a particle (molecule) spends in the reservoir from since it entered it up to the time it left it.

As a starting point, we consider a simple reservoir in which the flow is steady, laminar and one-dimensional, i.e., the flow properties change only along the flow direction, described with a variable $x\in \left[0,L\right]$, while the cross-sectional area may change with x and is $A(x)$. As the flow is steady, the discharge$Q_0$ is constant in time and the same for all x (inflow equal to outflow), while the velocity changes with $x$ as $V\left(x\right)=Q_0/A\left(x\right)$. Due to the laminar flow a molecule flowing into the reservoir at point $x=0$ will follow a smooth trajectory and will travel a distance $\mathrm{d}x$ in time $\mathrm{d}t=\mathrm{d}x/V\left(x\right)=A\left(x\right)\mathrm{d}x/Q_0$. Therefore, by integration, the total time that the molecule spends inside the reservoir, i.e., the residence time, denoted as $W_0$, is ${\int }_0^{L}A\left(x\right)\mathrm{d}x/Q_0$ and in this case is the same for all particles. The integral is the reservoir volume $S_0$ and hence

$$W_0=\frac{S_0}{Q_0}$$

(21)

As we will see below, this residence time, $W_0$, despite being determined above by assuming a simplistic and unrealistically regular system, is characteristic for any reservoir, however complex and chaotic.

Next, we examine a fully random system, in which particles are well mixed and the residence time can only be modelled by a stochastic approach. Let $\underset{\_}{W}$ be a stochastic variable (random variable) representing the time a molecule left the reservoir, once it entered at time $t=0$. (Notice that we use the Dutch notational convention to underline the stochastic variables.) At time $t$ the mass stored in the reservoir is $S(t)$. At time $t+\mathrm{d}t$ the mass is $S\left(t\right)+I\left(t\right)\mathrm{d}t-Q(t)\mathrm{d}t$. Of the molecules contained in $S(t)$ at time $t$, the proportion of those removed at time $t+\mathrm{d}t$ is found, by neglecting second order differentials, to be $Q(t)\mathrm{d}t/S(t)$.

The event that a molecule, which entered the reservoir at $t=0$, has remained at time $t=W$, can be denoted as $\{\underset{\_}{W}>W$}. The event that it leaves the reservoir in the next elementary interval $\mathrm{d}W$ is $\{\underset{\_}{W}\le W+\mathrm{d}W\}$ and its probability, conditional on $\{\underset{\_}{W}>W\}$, is equal to the ratio of the mass leaving during this interval, $Q\left(W\right)\mathrm{d}W$, to the total mass:

$$P\left\{\underset{\_}{W}\le W+\mathrm{d}W|\underset{\_}{W}>W\right\}=\frac{Q\left(W\right)\mathrm{d}W}{S(W)}$$

(22)

Applying the definition of conditional probability, we write

$$\frac{P\left\{W<\underset{\_}{W}<W+dW\right\}}{P\left\{\underset{\_}{W}>W\right\}}=\frac{Q\left(W\right)\mathrm{d}W}{S(W)}$$

(23)

The numerator is the probability density multiplied by $\mathrm{d}W$, i.e. the derivative of the distribution function $F_{\underset{\_}{W}}\left(W\right)$. Hence

$$\frac{F_{\underset{\_}{W}}^{\prime }\left(W\right)}{1-F_{\underset{\_}{W}}\left(W\right)}=\frac{Q\left(W\right)}{S(W)}$$

(24)

The solution of this differential equation is

$$F_{\underset{\_}{W}}\left(W\right)=1-\exp \left(-{\int }_0^W\frac{Q\left(t\right)}{S(t)}dt\right)$$

(25)

and in dimensionless form, with

$$\underset{\_}{w} :=\frac{\underset{\_}{W}}{W_0}$$

(26)

the distribution function of the dimensionless residence time is

$$F_{\underset{\_}{w}}\left(w\right)=1-\exp \left(-{\int }_0^w\frac{q\left(\tau \right)}{s\left(\tau \right)}d\tau \right)$$

(27)

or

$$F_{\underset{\_}{w}}\left(w\right)=1-\exp \left(-{\int }_0^w{\left(q\left(\tau \right)\right)}^{1-\frac{1}{b}}d\tau \right)$$

(28)

This is simplified to the following expressions for the indicated special cases:

• Linear reservoir (in which $q\left(\tau \right)/s\left(\tau \right)=1$), any inflow:

  $$F_{\underset{\_}{w}}\left(w\right)=1-{\mathrm{e}}^{-w}$$

  (29)
• Superlinear benchmark reservoir, $q\left(t\right)={\left(s\left(t\right)\right)}^2$, constant inflow:

  $$F_{\underset{\_}{w}}\left(w\right)=1-\frac{2e^{\frac{w}{\sqrt{q_0}}}}{\left(2+\left(e^{\frac{2w}{\sqrt{q_0}}}-1\right)\left(1+\sqrt{q_0}\right)\right)}$$

  (30)
• Sublinear benchmark reservoir, $q\left(t\right)=\sqrt{s\left(t\right)}$, constant inflow:

  $$F_{\underset{\_}{w}}\left(w\right)=1-\frac{{\mathrm{W}\left(\left(q_0-1\right){\mathrm{e}}^{-\frac{q_{0}w}{2}+q_0-1}\right)}^2}{{(q_0-1)}^2}$$

  (31)

Equation (29) is the standard exponential distribution. The equations for the benchmark reservoirs are non-standard and look complicated, yet if we plot them we see that they do not differ substantially from the exponential distribution. To make the differences more visible we choose to plot the probability ratio $F_{\underset{\_}{w}}\left(w\right)/\left(1-F_{\underset{\_}{w}}\left(w\right)\right)$ on a logarithmic axis, instead of the distribution function $F_{\underset{\_}{w}}\left(w\right)$, which would not show any visible difference. This is seen in Figure 3, where indeed the body of each distribution is almost identical to each other, and differences appear in the tails ($w>2$, noting that both the mean and standard deviation of the standard exponential distribution equal 1).

Next we proceed to forming an approximation of the distribution function for the general case, i.e. not limited to the benchmark cases. Specifically, for the ratio $q\left(\tau \right)/s\left(t\right)$ appearing in Equation (27) we propose the approximation:

$$r\left(\tau \right)≔\frac{q\left(\tau \right)}{s\left(\tau \right)}={\left(q\left(\tau \right)\right)}^{1-\frac{1}{b}}\approx d-\frac{Ac}{1+c\tau }=:r_{\mathrm{A}}(\rho )$$

(32)

where $A,c,d$ are parameters. Then the integral in Equation (27) becomes

$$\underset{0}{\overset{w}{\int }}\frac{q\left(\tau \right)}{s\left(\tau \right)}d\tau \approx \underset{0}{\overset{w}{\int }}\left(d-\frac{Ac}{1+c\tau }\right)d\tau =dw-\ln {\left(1+cw\right)}^A$$

(33)

so that Equation (27) becomes

$$F_{\underset{\_}{w}}\left(w\right)=1-e^{-dw}{\left(1+cw\right)}^A$$

(34)

The distribution is interesting, and its domain is $w\ge 0$ if $c\ge 0$ and $0\le w<-1/c$ if $c\ge 0$. In both cases it has all its moments finite if $d\ne 0$. For the special case $d=0$, the distribution becomes

$$F_{\underset{\_}{w}}\left(w\right)=1-{\left(1+cw\right)}^A$$

(35)

To find the parameters, we match the true values of${\left(q\left(\tau \right)\right)}^{1-1/b}$ at $\tau =0$ and $\infty$ and its derivative at $\tau =0$. The latter is found by approximating $q\left(\tau \right)$ from Equation (11). The obtained equations are

$$1=d-Ac, {\left(\frac{1}{q_0}\right)}^{1-\frac{1}{b}}=d, \left(b-1\right)\left(\frac{1}{q_0}-1\right)=A{c}^2$$

(36)

and their solution is

$$d=q_0^{\frac{1}{b}-1} , c=\frac{b-1}{d-1}\left(\frac{1}{q_0}-1\right), A=\frac{d-1}{c}$$

(37)

As seen in Figure 3, the approximation compares very well to the exact solutions for the benchmark cases (Equations (29)-(31)).

By integration of $w{F}_{\underset{\_}{w}}^{\prime }\left(w\right)$ we find that the mean is

$${\mu }_w=\left\{\begin{array}{ll}\frac{1}{d},& A=0\\ \frac{{\mathrm{e}}^{d/c}\mathsf{\Gamma }\left(A+1,d/c\right)}{c{\left(d/c\right)}^{A+1}},& c>0,d>0\\ {d\mathrm{e}}^{-d}& c=0,d>0\\ \frac{{\mathrm{e}}^{d/c}\left(\mathsf{\Gamma }\left(A+1,d/c\right)-\mathsf{\Gamma }\left(A+1\right)\right)}{c{\left(d/c\right)}^{A+1}},& c<0,d>0,A>0\\ \frac{-1}{\left(A+1\right)c},& d=0 \mathrm{a}\mathrm{n}\mathrm{d}\left\{\begin{array}{l}c>0,A<-1\\ \mathrm{o}\mathrm{r}c<0,A>0 \end{array}\right.\end{array}\right.$$

(38)

while in subcases not included above (e.g. $c>0,A>-1$), the mean diverges to infinity. From the equation $F_{\underset{\_}{w}}\left(w_{1/2}\right)=1/2$ the median is found to be

$$w_{1/2}=\left\{\begin{array}{ll}\frac{\ln 2}{d},& c=0 \mathrm{o}\mathrm{r} A=0\\ \frac{2^{-\frac{1}{A} }-1}{c},& d=0\\ -\frac{A}{d}W\left(-\frac{2^{-\frac{1}{A} }d{e}^{-\frac{d}{Ac}}}{Ac}\right)-\frac{1}{c},& \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}\right.$$

(39)

It is useful to note that for a linear reservoir ($b=1$) the above approximations result in $d=1,A=0$ and, hence, recover the exact case of an exponential distribution, in which ${\mu }_w=1, w_{1/2}=\ln 2$. In the exponential distribution we have $F_{\underset{\_}{w}}\left({\mu }_w\right)=1-{\mathrm{e}}^{-1}$ and this property has been used to define other variants of characteristic times such as the so-called e-time (used by Berry [25]), relaxation time, e-folding time, half-life, or decay constant (used by IPCC [26,27], but without providing definitions). Here we prefer not to use such terms, but stick to the probabilistically meaningful terms mean and median (dimensionless) residence times. The dimensional mean and median are found after multiplication by $W_0,\mathrm{i}.\mathrm{e}.{\mu }_W=W_0{\mu }_w, W_{1/2}=W_0{w}_{1/2}$.

We stress that, excepting some special cases, the above formulae provide approximations rather than exact values. To find a second approximation, simpler to apply, we conducted a numerical investigation and concluded with the following simple formulae, designated as approximation 2:

$${\mu }_w=b^{0.45\ln \left(q_0-a\right)}, w_{1/2}=\ln 2b^{\gamma \ln q_0}, \gamma =\left\{\begin{array}{cc}0.2,& b\le 1\\ 0.25,& b\ge 1\end{array}\right.$$

(40)

where the approximation of ${\mu }_w$ also includes the case of a linear change of inflow with rate $a$ (dimensionless).

Comparisons of the approximations of the mean residence time with the exact values for the benchmark cases and constant inflow are given in Figure 4. It is observed that approximation 2 is somewhat better than approximation 1—and also simpler. In addition, Figure 5 shows comparisons for changing inflows. Generally, the approximations perform well. Most importantly, both figures show that for the common ranges of $b$ and $q_0$, the dimensionless residence time is close to 1, not differing more than ±10%. Therefore, this thorough theoretical investigation results in the simple conclusion that in practical problems it is sufficient to take ${\mu }_w=1$ and $w_{1/2}=\ln 2=0.693$. Consequently, considering the dimensional time, in practice we may take the mean residence time as $W_0=S_0/Q_0$ and the median one as 70% of $W_0$.

2.3 Response Time

The response time is conceptually different from the residence time as it is related to the impulse response function (IRF), defined as follows:

Definition 2: 

The impulse response function (IRF), g(h), of a system is defined to be the system output at time lag h after the system is stimulated by an input that is an (instantaneous) impulse of unit mass (a Dirac delta function).

A thorough and general presentation of the IRF concept in a causality framework has been presented by Koutsoyiannis et al. [28,29,30]. In our case, the following particular considerations are taken into account for a reservoir: (a) the system can be studied in terms of the dimensionless quantities and the dimensional ones can then be obtained through Equation (7); (b) the response function $g\left(\eta \right)$ for dimensionless time lag $\eta =h/W_0$is identical to the output function $q\left(\tau \right)$ for $\eta =\tau$, which results from the impulse of an otherwise empty reservoir; and (c) the system is causal, which means that $g\left(\eta \right)=0$ for $\eta <0$.

Based on these considerations and the framework in [28], we proceed to the following:

Definition 3: 

Based on the dimensionless IRF of a reservoir, g(η), we define the mean and the median dimensionless response time as the mean ${\mu }_{\eta }$ and the median ${\eta }_{1/2}$ of the function g(η), namely:

$${\mu }_{\eta }≔\underset{0}{\overset{\infty }{\int }}\eta g\left(\eta \right)\mathrm{d}\eta , \underset{0}{\overset{{\eta }_{1/2}}{\int }}g\left(\eta \right)\mathrm{d}\eta =\frac{1}{2}$$

(41)

Notice that the median ${\eta }_{1/2}$ is defined by an implicit equation. To find the dimensional time lags we need to multiply those by the characteristic time $W_0=S_0/Q_0$.

To conceptualize the impulse, we imagine an empty reservoir (in dimensionless setting) with zero inflow, which at time zero receives instantaneously a unit input, increasing the storage to $s(0)=1$ and causing an outflow $q(0)=1$. Subsequently, the input becomes zero and the outflow is decreased. Hence, the outflow is given by Equation (17) and, consequently:

$${g\left(\eta \right)=\left(\left(b-1\right)\eta +1\right)}^{\frac{b}{1-b}}$$

(42)

It is then easy to find from Equation (41) the characteristic times, which are:

$${\mu }_{\eta }=\frac{1}{2-b}, {\eta }_{1/2}=\frac{2^{b-1}-1}{b-1}$$

(43)

Notice that, for $b\ge 2$, the mean ${\mu }_{\eta }$ diverges to infinity. It is stressed that the quantity ${\eta }_{1/2}$ is the dimensionless time required to empty half of the reservoir storage, and not that required for the outflow to be reduced to half its initial value (denoted as ${\zeta }_{1/2}$ and implicitly defined as $q\left({\zeta }_{1/2}\right)=1/2$). For completeness, the resulting expression for the latter is

$${\zeta }_{1/2}=\frac{2^{\frac{b-1}{b}}-1}{b-1}$$

(44)

While the median ${\eta }_{1/2}$ (similar to the mean ${\mu }_{\eta }$) is an increasing function of $b$, ${\zeta }_{1/2}$ is a decreasing one (for $b>0.21$), taking identical values only for $b=1$ (linear reservoir). The particular values for $b=1$ (linear reservoir) are identical to the mean and median residence times:

$${\mu }_{\eta }=1, {\eta }_{1/2}={\zeta }_{1/2}=\ln 2$$

(45)

It is interesting to note that the original definitions of the characteristic lag times in causality framework of Koutsoyiannis et al. [28] were based on a linearity assumption. On the other hand, Equations (27) and (17) did not assume linearity and thus the results found do hold for a nonlinear system (reservoir). This offers a basis for an extension of that causality framework. We note, however, that, while in the linear case the function $g\left(\eta \right)$ suffices to determine the output for any input through a convolution equation, this is not the case if the dynamics is nonlinear.

It is also easy to find the residence time in the case that the input is an impulse at time $\tau =w=0$. Using Equations (27) and (17) we find that its exact distribution function is

$$F_{\underset{\_}{w}}\left(w\right)=1-{\left(1+\left(b-1\right)w\right)}^{\frac{1}{1-b}}$$

(46)

It is readily seen that Equation (46) fully corresponds to Equation (34), if the parameters of the latter are $d=0,c=b-1,A=1/(1-b)$. We note though that these parameter values are not determined from Equation (37), as this does not work for zero inflow.

For $b=1$, this corresponds the exponential distribution; indeed, by taking the limit as $b\to 1$, we recover equation (29). For$b<1$ the variable $w$ becomes bounded from above by $1/(1-b)$, which means that the reservoir empties (and the outflow ceases) at a finite time. For $b>1$ the time to full emptying is infinite and for $b>2$ the mean of $\underset{\_}{w}$ is also infinite.

Furthermore, it can be observed that

$$F_{\underset{\_}{w}}\left(w\right)=1-s\left(w\right)=\underset{0}{\overset{w}{\int }}q\left(\tau \right)\mathrm{d}t$$

(47)

and if we take the derivatives, we find that the probability density function is

$$f_{\underset{\_}{w}}\left(w\right)=q\left(w\right)$$

(48)

This is not a coincidence, but it can be easily shown that it holds for any function $q\left(\tau \right)=f\left(s\left(\tau \right)\right)$. We can thus state this result as follows:

Theorem 1: 

The mean and median response time equal the mean and median residence time for the case that the input is an impulse function.

For different input, though, the residence time is different from the response time. The former depends on the input, while the latter depends only on the exponent $b$ (or more generally on the function $q\left(\tau \right)=f\left(s\left(\tau \right)\right)$). Illustrations of the different residence and response times for characteristic values of $b$ and $q_0$ are shown in Figure 6, which allows formulating the following:

• Observation 1: For a linear reservoir, the mean residence time and the mean response time are equal to each other, with a value $W_0=S_0/Q_0$. The median residence time and the median response are also equal to each other and smaller than $W_0$ by a factor ln 2 = 0.69.
• Observation 2: For a sublinear reservoir, the mean and median response times are generally smaller than the mean and median residence time, respectively, and can only become equal if the input is zero.
• Observation 3: From a practical point of view and for a reservoir that is not superlinear, the response times are smaller than the characteristic value $W_0=S_0/Q_0$, and the residence times can only slightly (by < 10%) exceed this value (for highly sublinear reservoirs and initial inflow higher than outflow).

It is useful to notice that the linear reservoir approximation, which in this case takes the form of Equation (14), is not satisfactory when inflow is zero. This is illustrated in Figure 7, where the approximation is compared to the exact one. The same is the case for the distribution function of residence time when inflow is zero, as seen in Figure 8. However, this is not a problem in our framework, since the exact solutions in this case are explicit (Equation (17) for outflow and (46) for residence time).

2.4 Parameters and their Estimation

The reservoir model, as a system with a power law connecting the outflow and storage, is quite simple and involves two parameters, the dimensionless exponent $b$ and a dimensional parameter. For a linear reservoir, the characteristic dimensional parameter is the mean residence time, $W_0=S_0/Q_0$, which is invariant (does not change if we change the time origin). However, if the reservoir is nonlinear this is not invariant. Instead, the invariant dimensional parameter is:

$$\Omega :=\frac{S^b}{Q}=\frac{S_0^b}{Q_0}=W_0{S}_0^{b-1}$$

(49)

which for $b=1$ becomes identical to $W_0$. The dimensions of $\Omega$ are [${\mathrm{M}}^{b-1}\mathrm{T}]$.

To estimate the parameters, we need at least a couple of simultaneous observations of $Q$ and $S$. However, these would not give any information about the appropriateness of the model. The ideal case is to have systematic measurements (time series) of all three processes $I,Q$ and $S$, and fit the model by minimizing the error on processes $Q$ and $S$. If we have observations of the storage process only, which is the easiest to measure, the model cannot be fitted, unless we have some additional estimates of the balance at some time scale. An example will be discussed in section 3.

3. Application to the Carbon Cycle

3.1 A Summary of the Established Approach

Typically, the carbon cycle is modelled in complex approaches, such as in climate models, which do not allow transparency and easy understanding. The established approach is reflected in the Assessment Reports of the Intergovernmental Panel on Climate Change (IPCC), among which here we refer to the Fifth Assessment Report (AR5) [26] and the Sixth Assessment Report (AR6) [27]. Comprehensive critiques on the established approach, accompanied with alternative approaches and quantifications, have been provided by Salby [31], Humlum et al. [32], Harde [33,34], Berry [25], Poyet [35] and Stallinga [36].

The approach presented here, with its simplicity and transparency, can shed light on the unclear issues and highlight the problems in the established approach. Among these problems, here we emphasize three, discussed in the following subsections.

3.1.1 Improper Terminology, Obscureness, Ambiguity and Vagueness

The terminology within the IPCC reports is different than used here, as seen in the following extract from the latest IPCC report [27 (Glossary; pp. 2237, 2246)]:

• Lifetime is a general term used for various time scales characterizing the rate of processes affecting the concentration of trace gases. The following lifetimes may be distinguished:
• […] Response time or adjustment time (T_a) is the time scale characterizing the decay of an instantaneous pulse input into the reservoir. The term adjustment time is also used to characterize the adjustment of the mass of a reservoir following a step change in the source strength. Half-life or decay constant is used to quantify a first-order exponential decay process. […]
• The term lifetime is sometimes used, for simplicity, as a surrogate for adjustment time.
• In simple cases, where the global removal of the compound is directly proportional to the total mass of the reservoir, the adjustment time equals the turnover time: T = T_a.
• […]
• Turnover time (T) (also called global atmospheric lifetime) is the ratio of the mass M of a reservoir (e.g., a gaseous compound in the atmosphere) and the total rate of removal S from the reservoir: T = M/S.
• […]
• Response time or adjustment time In the context of climate variations, the response time or adjustment time is the time needed for the climate system or its components to re-equilibrate to a new state, following a forcing resulting from external processes. It is very different for various components of the climate system. The response time of the troposphere is relatively short, from days to weeks, whereas the stratosphere reaches equilibrium on a time scale of typically a few months. […] In the context of lifetimes, response time or adjustment time (T_a) is the time scale characterizing the decay of an instantaneous pulse input into the reservoir.

We notice in the above definitions the terms lifetime, turnover time, global atmospheric lifetime, response time, adjustment time, half-life or decay constant, none of which is clear enough to allow quantification and even to allow distinguishing which one is referred to each time. The most general term appears to be “lifetime”, which is typically used to describe decaying entities such as unstable atoms (radioactive) or particles, but for the processes in question may be inappropriate. For example, the atmospheric carbon dioxide may not die or decay (e.g. in the ocean-atmosphere exchange), and therefore it is meaningless to speak about its lifetime. In contrast, it is meaningful to speak about its residence time, the time between entering and leaving the atmosphere. As explained in Section 2, the residence time can be represented as a stochastic variable ranging from zero to infinity, yet we can quantify and specify it through, either (a) its distribution function, (b) its average, (c) its median, (d) any other statistic (e.g. standard deviation, skewness, etc.) that would be relevant in each case. No such information is given or hinted in the above definitions by IPCC.

The ambiguity in the terminology and definitions is manifest also in the assessments and results provided. Thus, the IPCC AR5 contains the statements:

the concept of a single, characteristic atmospheric lifetime is not applicable to CO2 [26 (p. 473)].

No single lifetime can be given [for CO₂]. The impulse response function for CO₂ from Joos et al. (2013) [37] has been used. [26 (p. 737)].

Likewise, the IPCC AR6 refers to “multiple lifetimes for CO₂” without specifying which ones [27 (p. 302, Table 2.2; p. 1017, Table 7.15)].

Apparently, the residence time (and IPCC’s “lifetime”) may take any positive real value, if modelled as a stochastic variable, yet it has certain statistics, such as a mean, which IPCC avoids specifying, preferring to report that the values are multiple. It is interesting that the same reports give specific values for other substances. The reasons for this special treatment of CO₂ by IPCC, may be inferred from what follows.

3.1.2 Separate Treatment of CO₂ Depending on its Origin

The ambiguity is accompanied with inappropriate assumptions and speculations, the weirdest of which is that the behaviour of the CO₂ in the atmosphere depends on its origin and that CO₂ emitted by anthropogenic fossil fuel combustion has higher residence time than naturally emitted. This is clear in the IPCC AR5:

Simulations with climate–carbon cycle models show multi-millennial lifetime of the anthropogenic CO₂ in the atmosphere [26 (p. 435)].

It is also repeated in IPCC AR6:

This delay between a peak in emissions and a decrease in concentration is a manifestation of the very long lifetime of CO₂ in the atmosphere; part of the CO₂ emitted by humans remains in the atmosphere for centuries to millennia [27 (p. 642 FAQ 4.2)].

This weird idea has a long history, as it was thought from the beginning of climate modelling that the fate of anthropogenic CO₂, is different from that of the natural CO₂. For example, Joos et al. [38] stated:

When considering the fate of anthropogenic CO₂, the emission into the atmosphere can be considered as a series of consecutive pulse inputs.

More recently, in their study entitled “The millennial atmospheric lifetime of anthropogenic CO₂”, Archer and Brovkin [39] stated:

The largest fraction of the CO₂ recovery will take place on time scales of centuries, as CO₂ invades the ocean, but a significant fraction of the fossil fuel CO₂, ranging in published models in the literature from 20–60%, remains airborne for a thousand years or longer.

Also, Archer et al. [40] stated:

The models agree that 20-35% of the CO₂ remains in the atmosphere after equilibration with the ocean (2-20 centuries).

The idea is also redundantly repeated in grey literature (and more recently promoted by artificial intelligence chatbots), including in publications by universities and research organizations, such as the following by the Massachusetts Institute of Technology (MIT) and the National Aeronautics and Space Administration (NASA), respectively:

estimates for how long carbon dioxide (CO₂) lasts in the atmosphere […] are often intentionally vague, ranging anywhere from hundreds to thousands of years. […] As it stands, says [Ed] Boyle, human-generated carbon dioxide is expected to continue warming the planet for tens of thousands of years [41].

Once [carbon dioxide is] added to the atmosphere, it hangs around, for a long time: between 300 to 1,000 years. Thus, as humans change the atmosphere by emitting carbon dioxide, those changes will endure on the timescale of many human lives [42].

We may highlight in the former quotation the phrase “intentionally vague”, which faithfully conveys the fact that, behind all this vagueness, there are intensions. The reader interested in some amusement, may also see a summarizing depiction of the idea in a cartoon by Berry [43] depicting a demon that separates and delays the human CO₂ molecules in the atmosphere.

3.1.2 Inappropriate Modelling and Inconsistent Results

IPCC’s methodology in modelling the atmospheric CO₂ exchange is based on the so-called Bern modelling approach (Joos et al. [37,38]; Myhre et al. [44]; Strassmann and Joos [45]; Luo et al. [46]). It is reflected in the following expression of the IRF as the sum of three exponential functions and a constant term:

$$g\left(h\right)=a_0+\sum _{i=1}^3{a}_i{\mathrm{e}}^{-h/W_i}$$

(50)

The fitting of the parameters of this expression has been based on climate model results. This is made clear by Joos et al. [37] who stated (below their equation (11) and in the caption of their table 5) that they fitted on the mean of multimodel mean in future studies. In other words, the parameters were not obtained from observed data.

There are several problems with this methodology, in addition to the fact that it is based on imaginary data. These are discussed in general mathematical terms in Appendix A.2, as well as in numerical terms, with the specified values of the parameters also given in Appendix A.2, which were used in IPCC AR5 and IR6. In particular, the form of the equation is arbitrary and does not correspond to a reservoir’s dynamics. The inclusion of the constant term ($a_0$) results in theoretically infinite mean response time. Even if the constant term is excluded, the resulting mean response time is 353 years. With the inclusion of this term, even if we replace the nominal upper limit of integration, which is infinity, with 1000 years (the duration considered by Joos et al. [37] for their model fitting), the mean response time is no less than 432 years. These values can hardly be reconciled with the fact that the residence time of CO₂ is no more than 4 years, as admitted even by IPCC ([27 (p. 2237)]:

Carbon dioxide (CO₂) is an extreme example. Its turnover time is only about 4 years because of the rapid exchange between the atmosphere and the ocean and terrestrial biota. However, a large part of that CO₂ is returned to the atmosphere within a few years. The adjustment time of CO₂ in the atmosphere is determined from the rates of removal of carbon by a range of processes with time scales from months to hundreds of thousands of years. As a result, 15 to 40% of an emitted CO₂ pulse will remain in the atmosphere longer than 1,000 years, 10 to 25% will remain about ten thousand years, and the rest will be removed over several hundred thousand years.

In the next sections we will show that the first part of this quotation (referring to a 4 years’ time) is correct, while the last part is blatantly incorrect as not even one CO₂ molecule remains in the atmosphere for such long times.

3.2 Data

Systematic measurements of the atmospheric CO₂ have been made since 1958 [47] by the Scripps CO₂ Program of the Scripps Institution of Oceanography, University of California, and are available online [48,49,50]. The data include observations of CO₂ concentration (in micro-moles CO₂ per mole, or parts per million—ppm), and are processed to extract monthly values, filled in in case of missing data. Here the monthly time series have been retrieved and processed for two stations namely Mauna Loa Observatory, Hawaii (19.5° N, 155.6° W, 3397 m a.s.l.,1958 – present) and Barrow (recently renamed to Utqiagvik), Alaska (71.3° N, 156.6° W, 11 m a.s.l., 1961 – present).

Data of global human carbon emissions are also available online for years 1850 – 2022 and have been retrieved from [51,52]. The value of 2023 was taken as that of 2022 increased by 1.01, according to the International Energy Agency’s report [53 (p.3)].

For conversion of different units, we use the following coefficients:

• From mass of C to mass of CO₂, we multiply by 44/12 = 3.67 kg CO₂ / kg C (where 44 and 12 are the molecular masses of CO₂ and C);
• From atmospheric CO₂ concentration in ppm to total atmospheric mass in Gt CO₂ we multiply by 7.8 Gt CO₂ / ppm CO₂.

3.3 Premises of the Application

Based on the IPCC AR6 estimates of the global carbon balance [27 (Figure 5.12)], Koutsoyiannis et al. [30] compiled a summary graph of total carbon emissions and sinks, distinguishing the preindustrial quantities (before 1750) and modern additions. This graph is reproduced here as Figure 9, after conversion from Gt C to ppm CO₂.

Based on this graph we make the following observations, which are important for the modelling of the CO₂ exchanges that follow:

• Human activities are responsible for only 4% of carbon emissions.
• The vast majority of changes in the atmosphere since 1750 (red bars in the graph) are due to natural processes, respiration and photosynthesis.
• The increases of both CO₂ emissions and sinks are due to the temperature increase, which expands the biosphere and makes it more productive.
• The terrestrial biosphere processes are much stronger than the maritime ones in terms of both production and absorption of CO₂.
• The CO₂ emissions by merely the ocean biosphere are much larger than human emissions.
• The modern (post 1750) CO₂ additions to pre-industrial quantities (red bars in the right half of the graph, corresponding to positive values) exceed the human emissions by a factor of ~4.5. In the most recent 65 years, covered by measurements, the rate of natural emissions is ~3.5 times greater than the CO₂ emissions from fossil fuels.

Point 3 above implies a causality direction between temperature and CO₂ concentration that is opposite to the popular one, which is the one also assumed and embedded in climate models. Specifically, the conventional wisdom that it is the increased atmospheric carbon dioxide concentration ([CO₂]) that caused the increase of temperature (T) was questioned by Koutsoyiannis and Kundzewicz [54], while later Koutsoyiannis et al. [28,29,30] provided evidence, based on analyses of instrumental measurements of the last seven decades, for a unidirectional, potentially causal link between T as the cause and [CO₂] as the effect.

The effect of [CO₂] to the CO₂ inflow to the atmosphere is depicted in Figure 10. The two end points, corresponding to 1750 and 2017, were determined from the quantities given in Figure 9, with the additional information that [CO₂] was 280.0 and 404.6 in these two years, respectively.

The remaining points are obtained from the so-called $Q_{10}$ model [55], as in Koutsoyiannis et al. [30 (Appendix A.1)]. This is based on the equation:

$$\frac{R_t}{R_0}=Q_{10}^{\left(T_t-T_0\right)/(10\mathrm{K})}$$

(51)

where $R_t$ and $R_0$ denote the respiration rate at times $t$ and $t_0$, respectively, $T_t$ and $T_0$ are the temperatures at these times, and $Q_{10}$ is a dimensionless parameter, characteristic of the species. Assuming a linear trend $c_T$ in temperature, we get

$$\frac{R_t}{R_0}=Q_{10}^{c_T(t-t_0)/(10\mathrm{K})}$$

(52)

The linear trends for the last 65-year period, were calculated in [30] from the NCEP/NCAR Reanalysis data at 0.26 °C/decade for the terrestrial and 0.12 °C for the maritime part. The literature gives representative average $Q_{10}$ values of 3.05 for the terrestrial respiration [55] and 4.07 for the maritime respiration [56]. For these values and after anchoring the calculations at year 2017, the resulting total inflow to the atmosphere (the sum of respiration from the terrestrial and the maritime part) are shown in Figure 10 against the [CO₂] observed values at a decadal time step. The figure also shows a power law between [CO₂] and input (emitted) CO₂, with an exponent of 0.7, determined from the two end points. The intermediate points do not perfectly agree with this power law, although they show a rising trend. If we increased the $Q_{10}$ values by 30%, the agreement would be better, as also shown in the figure. In addition, the figure shows two points of CO₂ outflow from the atmosphere corresponding to 1750 and 2017, which were again determined from the quantities given in Figure 9. These points show a similarity with those of input, with a greater rise, yielding an exponent of a power law equal to 0.77.

3.4 Model and its Fitting Methodology

To apply the reservoir routing methodology to the atmospheric CO₂, we represent the storage $S(t)$ as the atmospheric [CO₂] (in ppm) and the inflow $I(t)$ and outflow $Q(t)$ as the emissions and sinks, respectively (ppm/year). An eminent characteristic of the atmospheric CO₂ exchange is its seasonality, implying seasonal variation of the characteristic residence time. To take seasonality into account in a parsimonious manner, we modify equation (2) by substituting $S_0/W_0$ for $Q_0$ and then replacing $W_0$ with a periodic function of time, $W\left(t\right)$:

$$Q\left(S\right)=\frac{S_0}{W\left(t\right)}{\left(\frac{S}{S_0}\right)}^b$$

(53)

In a first phase (Phase 1), we made several preliminary model runs with different parameterizations, including mathematical expressions of $Q\left(S\right)$ different from the power law, but here we present final runs. We concluded that the best results are obtained by the power-law relationship and with the following simple and parsimonious mathematical form for $W\left(t\right)$:

$$W\left(t\right)=A{\left(\cos \left(2\pi t+\phi \right)+\psi \right)}^b$$

(54)

where the parameter $A$ has dimensions of time (as does $W\left(t\right)$) and the parameters $\phi$ (phase) and $\psi$ are dimensionless. Hence:

$$Q\left(t\right)=\frac{S_0}{A}{\left(\frac{S(t)}{S_0\left(\cos \left(2\mathsf{\pi }t+\phi \right)+\psi \right)}\right)}^b$$

(55)

For two times $t_1$ and $t_2$ differing by $2k\mathsf{\pi }t$, where $k$ is an integer (meaning: referring to the same date in different years), we have $Q\left(t_1\right)/S{\left(t_1\right)}^b=Q\left(t_2\right)/S{\left(t_2\right)}^b$, which is a periodic extension of the invariance condition in Equation (49).

Given the similar behaviours in the input and outflow, as seen in Figure 10, we choose a similar expression for the natural inflow, $I_{\mathrm{N}}\left(t\right)$, which is not measured:

$$I_{\mathrm{N}}\left(t\right)=\frac{S_0}{A_I}{\left(\frac{S(t)}{S_0\left(\cos \left(2\mathsf{\pi }t+{\phi }_I\right)+{\psi }_I\right)}\right)}^{b_I}$$

(56)

with parameters $b_I,A_I,{\phi }_I,{\psi }_I$ similar to those in the expression of outflow. This defines a characteristic time for input:

$$W_I\left(t\right)=A_I{\left(\cos \left(2\pi t+{\phi }_I\right)+{\psi }_I\right)}^{b_I}$$

(57)

To find the total input we add the anthropogenic emissions, $I_{\mathrm{A}}\left(t\right)$, which are known as described is Section 3.2:

$$I\left(t\right)=I_{\mathrm{N}}\left(t\right)+I_{\mathrm{A}}\left(t\right)$$

(58)

In this, we have neglected inflows from volcanism and other outgassing sources (e.g. related to El Niño–Southern Oscillation), which may induce some inaccuracies in our modelling.

To apply the differential equation with the observations, which are in discrete time with monthly step, we discretize the time, $t=j\Delta$ and to avoid an implicit numerical scheme, we reduce the time step to half monthly ($\Delta \approx 1/24$ years, but it varies slightly among months because of the different durations of the months) and estimate the values of $S(t)$ at the mid-month as the average of the values in consecutive months. The differential equation is then written in discrete time as:

$$\frac{S\left(\left(j+1\right)\Delta \right)-S\left(j\Delta \right)}{\Delta }+Q\left(j\Delta \right)=I_{\mathrm{N}}\left(j\Delta \right)+I_{\mathrm{A}}\left(j\Delta \right)$$

(59)

from which we find the storage at the step $\left(j+1\right)\Delta$ from the values of the three variables at step $j\Delta$, i.e.,

$$S\left(\left(j+1\right)\Delta \right)=S\left(j\Delta \right)+N\left(j\Delta \right)\Delta$$

(60)

where $N\left(j\Delta \right)$ is the net inflow at step $j\Delta$:

$$N\left(j\Delta \right) :=I_{\mathrm{N}}\left(j\Delta \right)+I_{\mathrm{A}}\left(j\Delta \right)-Q\left(j\Delta \right)$$

(61)

The model has eight parameters, two dimensional exponents $b,b_I$ of the storage-outflow and storage-inflow power laws, two dimensional $A$ (in time units), and four dimensionless for the cosine functions that describe seasonality. We apply the model to two locations (Barrow and Mauna Loa) letting the parameter values be different in the two locations, except for the exponents $b,b_I$, which we assume to be the same at both locations. To estimate these parameters, we make a complete simulation at both locations and perform an optimization (using a typical solver in a common spreadsheet software).

The quantities that provide means for comparing actual data and simulations are $S\left(j\Delta \right)$ and $N\left(j\Delta \right)$. Observed time series of the former are readily available, while the latter are estimated from the former from Equation (60). The simulated values, which we denote $\widehat{S}\left(j\Delta \right)$ and $\widehat{N}\left(j\Delta \right)$, are determined from Equations (55) – (61). After completing a simulation with a specified parameter set, we determine the explained variance for each quantity, $S\left(j\Delta \right)$ and $N\left(j\Delta \right)$, as:

$${\mathrm{E}\mathrm{V}}_S :=1-\frac{\mathrm{v}\mathrm{a}\mathrm{r}\left[\widehat{S}\left(j\Delta \right)-S\left(j\Delta \right)\right]}{\mathrm{v}\mathrm{a}\mathrm{r}\left[S\left(j\Delta \right)\right]}, {\mathrm{E}\mathrm{V}}_N :=1-\frac{\mathrm{v}\mathrm{a}\mathrm{r}\left[\widehat{N}\left(j\Delta \right)-N\left(j\Delta \right)\right]}{\mathrm{v}\mathrm{a}\mathrm{r}\left[N\left(j\Delta \right)\right]}$$

(62)

respectively. These quantities are equivalent to the coefficient of determination ($R^2$) in linear regression and are also known as Nash-Sutcliffe efficiencies. We form an objective function as the sum of these two quantities in both locations and we maximize this sum by changing the parameters.

In a second phase (Phase 2), we focus on determining the exponents $b,b_I$ by numerous simulation runs, by fixing $b$ and optimizing all other parameters. Figure 11 shows the optimized $b_I$ as a function of $b$. Consistently with what was already observed in the discussion of Figure 10, the optimal $b_I$ is always smaller than $b$ by 0.05-0.06. This small difference is essential to keep, as no good fitting is possible if the two are assumed equal.

Figure 12 shows the explained variances as functions of the specified exponent $b$. Those of $S\left(j\Delta \right)$ are extraordinarily high (higher than 0.997), while those of $N\left(j\Delta \right)$ are around 0.85. We observe that the explained variances for Barrow are consistently increasing with the increase of $b$, while those at Mauna Low decrease for $b>1$. From a Pareto optimality point of view, these results suggest that there is no meaning in adopting any value $b<1$. Therefore. we finally choose $b=1$, i.e. a linear reservoir, for the additional reason that it is simpler, exact in its analytical solution, and more intuitive.

3.5 Results of Final Modelling

The next phase (Phase 3) includes the fine tuning of the results of Phase 2, after choosing $b=1$, and their detailed presentation. The final optimized parameters are shown in Table 1.

The evolution of storage, $S(t)$, observed and simulated, is shown in Figure 13, along with a visualization of the seasonal variation of the characteristic times $W,W_I$. The agreement of observed and simulated is impressively good, as visually seen and also indicated in the values of explained variances, marked in the figures.

The evolution of the simulated inflow and outflow ($I(t)$, $Q(t)$) is shown in Figure 14. Here we do not have means for comparison as there are no observations of these quantities. Yet we know that in the last ten years, the average outflow should be close to 104.9 ppm/year (Figure 9). In the optimization we introduced a constraint that the average outflow should not depart more than 5% from the value of 104.9 ppm/year, and this constraint is satisfied. It is informative to compare in this figure the natural emissions to the anthropogenic ones, which are a very small portion of the total.

The evolution of the observed and simulated net inflow ($N\left(t\right)=I\left(t\right)-Q(t)$) is shown in Figure 15. There is a good agreement between the two, reflected in explained variances, of ~86%, as shown in the figures. Yet there is some discrepancy in reproducing the increasing variation of the net inflow with time in Barrow, which is due to the biosphere expansion.

It is not easy to improve the fitting in terms of the latter discrepancy and better represent the biosphere expansion in the last years, unless we sacrifice the parsimony in modelling. An example is given in Figure 16, where we replaced the cosine function for the inflow characteristic time $W_I$ (Equation (57)) with 24 half-monthly values, while not modifying the corresponding expression for the outflow. As seen in the figure, the explained variance was improved from 86% to 92%. Yet there is no perfection in reproducing the biosphere expansion, and, anyhow, the pursuit of parsimony should never be neglected. For these reasons, we may regard our final solution that presented in Figure 13 through Figure 15, rather than that of Figure 16.

3.6. Results for Imaginary Cases

Our next phase (Phase 4) is devoted to examining some imaginary cases to offer additional insights. In this phase we examine the following four cases:

Human emissions are disregarded and only natural processes are considered.

The natural processes are neglected and only the anthropogenic emissions are considered.

In addition to anthropogenic emissions, natural outputs (but no inputs) are also considered.

All processes are considered, but the biosphere expansion is neglected.

The results for case 1 are shown in Figure 18 in terms of comparisons of observed and simulated time series of storage, $S(t)$. The agreement of observed and simulated series as impressively good as that in the complete modelling shown in Figure 13. Hence, the inclusion or omission of the anthropogenic contribution does not offer anything important in modelling, except in altering the model parameters.

Figure 17. Comparison of observed and simulated storage, $S\left(t\right)\equiv \left[\mathrm{C}{\mathrm{O}}_2\right]$, for (upper) Mauna Loa and (lower) Barrow, as in Figure 13, but omitting anthropogenic emissions (imaginary case 1). The insets show the seasonal variation of the characteristic times $W,W_I$.

Figure 17. Comparison of observed and simulated storage, $S\left(t\right)\equiv \left[\mathrm{C}{\mathrm{O}}_2\right]$, for (upper) Mauna Loa and (lower) Barrow, as in Figure 13, but omitting anthropogenic emissions (imaginary case 1). The insets show the seasonal variation of the characteristic times $W,W_I$.

The results for all other three cases of Phase 4 are shown in Figure 18 and are not satisfactory. The worst of all is case 2, in which merely the anthropogenic emissions are considered. The results have no relationship with reality. Case 3 is better, but again neither the overyear trend nor the seasonality are captured. Case 4 is even better as it captures seasonality, but the overyear trend is again not well represented. To implement case 4 case (all processes but without biosphere expansion), $S_0$ was substituted for $S$ in Equations (55) and (56). For this case, Figure 19 shows the observed and simulated net inflow ($N\left(t\right)=I\left(t\right)-Q(t)$), where the inability to capture the observed behaviour, namely the increasing variation of the net inflow with time in Barrow, which is due to the biosphere expansion, is manifest.

Figure 18. Comparison of observed and simulated storage, $S\left(t\right)\equiv \left[\mathrm{C}{\mathrm{O}}_2\right]$, for (upper) Mauna Loa and (lower) Barrow, as in Figure 13 but for the indicated imaginary cases. The insets show the seasonal variation of the characteristic times $W,W_I$ for the case of all emissions without biosphere expansion (imaginary case 4). The explained variance noted is for the same case, while it is smaller in the other indicated cases, decreased to less than –100% for the worst case of human emissions only.

Figure 18. Comparison of observed and simulated storage, $S\left(t\right)\equiv \left[\mathrm{C}{\mathrm{O}}_2\right]$, for (upper) Mauna Loa and (lower) Barrow, as in Figure 13 but for the indicated imaginary cases. The insets show the seasonal variation of the characteristic times $W,W_I$ for the case of all emissions without biosphere expansion (imaginary case 4). The explained variance noted is for the same case, while it is smaller in the other indicated cases, decreased to less than –100% for the worst case of human emissions only.

3.7 Residence Times

The fitted model parameters directly give the reservoir characteristics and their seasonal variation. That variation, which is substantial, is easy to explain as, contrary to the common belief highlighting anthropogenic emissions, the carbon cycle is dominated by the natural emissions and absorptions. It is useful to find an annual average, in the following manner:

$$W_{\mathrm{m}}≔\frac{{\int }_0^{1}S\left(t\right)\mathrm{d}t}{{\int }_0^{1}Q\left(S\left(t\right)\right)\mathrm{d}t}=\frac{{\int }_0^{1}S\left(t\right)\mathrm{d}t}{{\int }_0^1\frac{S_0}{A}{\left(\frac{S(t)}{S_0\left(\cos \left(2\pi t\right)+\psi \right)}\right)}^b\mathrm{d}t}$$

(63)

where for convenience (and without introducing errors) we have changed the time origin so that the phase φ be zero. Assuming much smaller variation of $S(t)$ than $Q(t)$ within a year, we simplify the expression to:

$$W_{\mathrm{m}}\approx \frac{S_0}{\frac{S_0}{A}{\int }_0^1{\left(\frac{1}{\left(\cos \left(2\pi t\right)+\psi \right)}\right)}^b\mathrm{d}t}=\frac{A}{{\int }_0^1{\left(\frac{1}{\left(\cos \left(2\pi t\right)+\psi \right)}\right)}^b\mathrm{d}t}$$

(64)

The integral can be evaluated analytically for any $b$, but here we give its value for $b=1$ which we was resulted from our analysis. The integral in this case is $1/\sqrt{{\psi }^2-1}$

$$W_{\mathrm{m}}\approx A\sqrt{{\psi }^2-1}=\sqrt{\left(A\left(\psi -1\right)\right)\left(A\left(\psi +1\right)\right)}=\sqrt{W_{\mathrm{m}\mathrm{i}\mathrm{n}}{W}_{\mathrm{m}\mathrm{a}\mathrm{x}}}$$

(65)

given that the seasonal minimum and maximum values of $W$ are, respectively, $W_{\mathrm{m}\mathrm{i}\mathrm{n}}=A\left(\psi -1\right),W_{\mathrm{m}ax}=A\left(\psi +1\right)$.

In other words, the annual mean of residence time is the geometric mean of the minimum and maximum values of $W(t)$. The characteristic seasonal and annual mean residence times are shown in Table 2. They vary seasonally from ~1.5 to ~10 years at Barrow, with a narrower range (~2 to ~6 years) at Mauna Loa. On an annual basis, the residence time is ~3.5 to ~4 years.

4. Discussion

In light of the above analyses and results, we can discuss some relevant questions of general interest. The first question is what part of human emissions through the period 1850 to date (the period for which emission data are available) has remained in the current atmosphere.

To answer this question, we observe that, from the mass $\mathrm{d}{m}_I\left(t\right)$ that entered the atmosphere at time $\left[t,t+dt\right]$, there remains a portion equal to $P\left\{\underset{\_}{W}>t_c-t\right\}$, where $t_c$ is the current time, which is equal to $=1-F(t_c-t)$. In other words, the mass remaining is

(66)

By integrating from $t_0=1850$ to $t_c=2023$ we can find the total mass remaining, $M_R$. If $M_A$ is the total mass of anthropogenic emissions through this period, then the proportion remaining is:

(67)

Application with emission data and with $W_0=4$ years results in $M_{\mathrm{R}}=163$ Gt CO₂ or 20.9 ppm, while $M_{\mathrm{A}}=2612$ Gt CO₂ or 334.9 ppm, so that $M_{\mathrm{R}}/M_{\mathrm{A}}=8\%$, comparable to (somewhat smaller than) the estimate ~10% by Stallinga [36] and also slightly smaller than the cumulative emissions of the last 4 years (as is reasonable). This contradicts the IPCC assertion [27 (p. 676, also repeated many times in AR6)] that:

Over the past six decades, the average fraction of anthropogenic CO2 emissions that has accumulated in the atmosphere (referred to as the airborne fraction) has remained nearly constant at approximately 44%.

As a second question, we examine whether or not there is some truth in the above quoted IPCC’s statement that “15 to 40% of an emitted CO₂ pulse will remain in the atmosphere longer than 1,000 years, 10 to 25% will remain about ten thousand years, and the rest will be removed over several hundred thousand years”. We examine it also in connection with the statement that the “turnover time is only about 4 years”, which we deem correct as it agrees with the results of this study. We make the following calculations:

• The probability that a molecule remains after 1000 years is $p=1-F_{\underset{\_}{W}}\left(1000\right)=1-F_{\underset{\_}{w}}\left(1000/4\right)=e^{-250}={10}^{-108.6}$, where we have used Equation (29) to evaluate the $F_{\underset{\_}{w}}\left(w\right)$.
• The probability that out of $N$ molecules none remains after 1000 years is ${(1-p)}^N$ and the probability that at least one molecule remains is ${p_1=1-(1-p)}^N$. Given that as $p\to 0$, $\left({1-\left(1-p\right)}^N\right)/pN\to 1$, for small $p$ (as in our case), we have $p_1=pN$.
• According to IPCC [27 (Figure 5.12)] the atmospheric CO₂ amounts to 850 Pg C = $8.5\times {10}^{17}$g C. Thus, the mass of CO₂ is $8.5\times {10}^{17}\times \left(44/12\right)=3.1\times {10}^{18}$ g (where 44 and 12 are the molecular masses of CO₂ and C). The number of moles is $3.1\times {10}^{18} \mathrm{g} / \left(44\mathrm{g}/\mathrm{m}\mathrm{o}\mathrm{l}\right)=7.1\times {10}^{16} \mathrm{m}\mathrm{o}\mathrm{l}$.
• The Avogadro constant is $6.022\times {10}^{23} \mathrm{m}\mathrm{o}{\mathrm{l}}^{-1}$ and thus the number of CO₂ molecules in the atmosphere is $N=7.1\times {10}^{16} \mathrm{m}\mathrm{o}\mathrm{l}\times 6.022\times {10}^{23} \mathrm{m}\mathrm{o}{\mathrm{l}}^{-1}=4.3\times {10}^{40}={10}^{40.6}$.
• Hence, the probability that, after 1000 years, at least one out of the $N={10}^{40.6}$ molecules remains in the atmosphere is $p_1=pN={10}^{-108.6}\times {10}^{40.6}={10}^{-68}$.
• A probability ${10}^{-68}$ is almost an impossibility. Hence, we can be certain that none of the molecules existing in the atmosphere now, whether due to “emitted CO₂ pulse” or existing before it, will remain after 1000 years—let alone after “ten thousand years” or after “several hundred thousand years”.
• To make this probability a reasonable rarity of 1% (${10}^{-2}$) we need to make $p=p_1/N={10}^{-2}/{10}^{40.6}={10}^{-42.6}$. This would occur at time $t$ such that $1-F_{\underset{\_}{W}}\left(t\right)=1-F_{\underset{\_}{w}}\left(t/4.0\right)=e^{-t/4.0}={10}^{-42.6}$, which yields $t=392$ years.

In other words, the IPCC’s statement “15 to 40% of an emitted CO₂ pulse will remain in the atmosphere longer than 1,000 years” needs to be corrected to “not even one molecule from an emitted CO₂ pulse will remain in the atmosphere longer than 400 years, even if that emitted pulse amounts to the entire current atmospheric CO₂ content”.

The above results are based on the linear reservoir dynamics ($b=1$, according to our optimized model run). If $b$ was smaller, e.g. $b=0.77$ as in the discussion in Section 3.3, the processes would be faster and the probabilities even smaller (as discussed below Equation (46), the residence time $\underset{\_}{w}$, as a stochastic variable, becomes bounded from above and the IRF becomes zero at a finite time).

5. Conclusions

The study offers a comprehensive framework of reservoir routing which is of some usefulness for several problems in hydrology, hydraulics and water management. Additionally it offers some insights on the application of mass balance (continuity equation) with linear or nonlinear dynamics. It defines and clarifies the relevant quantities, which are often confused in literature. It provides an exact solution for the instantaneous response function and the response time, whether in the linear or nonlinear case, and a working approximation of the outflow and the residence time (including its probability distribution and statistical characteristics). These are obtained by theoretical analyses and are readily applicable in practical problems.

On the theoretical aspect, our analyses provide a case where the instantaneous response function results directly from the system dynamics, rather than from stochastic, data-based means, thus complementing the recent causality framework by Koutsoyiannis et al. [28,29,30]. It further provides an extension of this framework for nonlinear dynamics, which deserves further pursuit. Additionally, it confirms the importance of the feature of this framework to include a nonnegative value at zero time lag (a value which actually in the reservoir case is the global maximum of the function), contrary to the Granger causality scheme [57] which excludes the zero time lag (see additional discussion on that issue in [28]).

While, our framework is fairly general and comprehensive, it cannot represent every problem related to storage systems. In particular, the paper’s scope leaves out prediction problems, which may require advanced stochastic methodologies. These may be dealt with in future research.

The application of the reservoir routing framework to the atmospheric CO₂ gives useful insights, in terms of characteristic residence times, which have been an issue of controversy. The theoretical framework results are in excellent agreement with real world data of carbon dioxide concentration. The atmosphere appears to behave as a linear reservoir in terms of the atmospheric CO₂, whose exchange is clearly dominated by the biosphere processes, with the human emissions playing a minor role. The quantification of the atmospheric CO₂ exchange with our framework yields reliable and intuitive results, complying with observations, in contrast to the disputable results of complex climate models. The mean residence time of atmospheric CO₂ is no more than four years and the response time is smaller than that, thus contradicting the mainstream estimates which suggest times of hundreds or thousands of years.

Undoubtedly, numerous natural processes are involved in the carbon cycle, which operate on widely ranging time scales. Indeed, we have rapid processes (photosynthesis, respiration), which occur over days to years and slower processes (e.g. ocean-atmosphere exchange), which operate over timescales of decades to centuries. There are also very slow processes (e.g. carbonate formation) that operate over millennial timescales. But this is irrelevant, as rapid processes remove the CO₂ molecules at their pertinent scales, without waiting for the slow or the very slow processes to act.

Clearly, the atmospheric CO₂ observational data are not consistent with the climate narrative. They rather contradict it. In this, the present study complements earlier studies that (a) causality direction between temperature and atmospheric CO₂ is opposite to commonly assumed [29,30], (b) climate models misrepresent the causality direction that is identified by the data [30], (c) there are no discernible signs of anthropogenic CO₂ emissions on the greenhouse effect, which is dominated by water vapour and clouds [58], or on the isotopic synthesis of atmospheric CO₂, which is determined by the biosphere processes [i].