Reservoir Routing and Its Application to Atmospheric Carbon Dioxide Balance

Preprint

Article

Reservoir Routing and Its Application to Atmospheric Carbon Dioxide Balance

Altmetrics

Downloads

142

Views

260

Comments

A peer-reviewed article of this preprint also exists.

Demetris Koutsoyiannis^*

This version is not peer-reviewed

Submitted:

06 May 2024

Posted:

08 May 2024

You are already at the latest version

Alerts

Abstract

Reservoir routing has been a routine procedure in hydrology, hydraulics and water management. It is typically based on the mass balance (continuity equation) and a conceptual equation relating storage and outflow. If the latter is linear, 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 streamline the reservoir routing framework and extend it to find approximate solutions for nonlinear cases. The proposed framework can also be useful for climatic tasks, such as describing the mass balance of atmospheric carbon dioxide and 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. The mean residence time of atmospheric carbon dioxide turns out to be no more than four years and the response time is smaller than that, thus opposing the much longer mainstream estimates.

Keywords:

Subject: Environmental and Earth Sciences - Water Science and Technology

What is more I loved, and still do love, mathematics for itself as not allowing room for hypocrisy or vagueness, my two pet aversions.

Stendhal [1 (p. 111)].

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{d S (t)}{d t} + 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 (S) = Q_{0} {(\frac{S}{S_{0}})}^{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}

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{d S (t)}{d t} + \frac{1}{W_{0}} S (t) = I (t), 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 (t) = S_{0} ⅇ^{- t / W_{0}} + \frac{1}{W_{0}} \int_{0}^{t} ⅇ^{- (t - s) / W_{0}} I (s) ⅆ s, Q (t) = Q_{0} ⅇ^{- t / W_{0}} + \frac{1}{W_{0}} \int_{0}^{t} ⅇ^{- (t - s) / W_{0}} I (s) ⅆ s

(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 (S)

is sublinear, meaning that

\lim_{S \to \infty} Q (S) / 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 \neq 1

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

\frac{d s (τ)}{d τ} + {(s (τ))}^{b} = \frac{i (τ)}{q_{0}}

(5)

with initial conditions

s (0) = q (0) = 1

(6)

where

τ = \frac{t}{W_{0}}, s (τ) = \frac{S (τ W_{0})}{S_{0}}, q (τ) = \frac{Q (τ W_{0})}{Q_{0}}, i (τ) = \frac{I (τ W_{0})}{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 = τ = 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

{(s (τ))}^{b}

s = 1

(similar to Basha [24]), as

{(s (τ))}^{b} = b s (τ) + 1 - b + O ({(s - 1)}^{2})

(8)

and neglecting the high order terms we get

\frac{d s (τ)}{d τ} + b s (τ) = \frac{i (τ)}{q_{0}} + b - 1

(9)

whose solution is

s (τ) = ⅇ^{- τ} + \int_{0}^{τ} ⅇ^{- (τ - s)} (\frac{i (τ)}{q_{0}} + b - 1) ⅆ s, q (τ) = b ⅇ^{- τ} - b + 1 + b \int_{0}^{τ} ⅇ^{- (τ - s)} (\frac{i (τ)}{q_{0}} + b - 1) ⅆ s

(10)

If the dimensionless inflow is constant (i.e.

i (τ) = 1

, or dimensional inflow

I (t) = I_{0} = Q_{0} / q_{0}

), this simplifies to

s (τ) = 1 - \frac{1}{b} (1 - \frac{1}{q_{0}}) (1 - ⅇ^{- b τ}), q (τ) = 1 - (1 - \frac{1}{q_{0}}) (1 - ⅇ^{- b τ}) = \frac{1}{q_{0}} + (1 - \frac{1}{q_{0}}) ⅇ^{- b τ}

(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

\int_{0}^{t} Q_{0} q (\frac{u}{W_{0}}) d u - I_{0} t = \frac{(q_{0} - 1) Q_{0} W_{0}}{b q_{0}} (1 - ⅇ^{- \frac{b t}{W_{0}}})

(12)

while the change in storage is

S_{0} - S_{0} s (\frac{t}{W_{0}}) = S_{0} \frac{1}{b} (1 - \frac{1}{q_{0}}) (1 - ⅇ^{- \frac{b t}{W_{0}}})

(13)

and upon substituting

W_{0} Q_{0} 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 (τ) = q (τ) = 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 (τ) = 1 - \frac{1}{b} (1 - ⅇ^{- b τ}), q (τ) = ⅇ^{- b τ}

(14)

but it is not exact; the exact one will be discussed below. For

b > 1

and

τ \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) \leq 1 - 1 / b

. For

b < 1

, the resevoir empties at

τ = - \ln (1 - b) / b

and beyond that time

q (τ) = 0

For inflow linearly varying with time, i.e.

i (τ) = 1 + a τ

, the solution becomes

s (τ) = 1 + \frac{a τ}{b q_{0}} - \frac{1}{b} (1 - \frac{1}{q_{0}} + \frac{a}{b q_{0}}) (1 - ⅇ^{- b τ}), q (τ) = \frac{1}{q_{0}} - \frac{a}{b q_{0}} + \frac{a τ}{q_{0}} + (1 - \frac{1}{q_{0}} + \frac{a}{b q_{0}}) ⅇ^{- b τ}

(15)

where the rightmost equation can also be written as

q (τ) = \frac{i (τ - 1 / b)}{q_{0}} + (1 - \frac{i (- 1 / b)}{q_{0}}) ⅇ^{- b τ}

(16)

a \geq 0

, the solution is valid for any τ; if

a < 0

, it is valid for

τ \leq - 1 / a

Some special cases of

b

and

i (τ)

also admit exact analytical solutions, like in the linear reservoir problem. A first case is when the inflow is zeroed at time

τ = 0

, and continues to be zero (

i (τ) = 0

) at any

τ > 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 (τ) = ((b - 1) τ + 1)}^{\frac{1}{1 - b}}, {q (τ) = ((b - 1) τ + 1)}^{\frac{b}{1 - b}}

(17)

and it apparently holds for

(1 - b) τ \leq 1

, that is, for any

τ

b \geq 1

, but for

τ \leq 1 / (1 - b)

b < 1

; beyond that time

s (τ) = q (τ) = 0

. As

b \to 1

, we have the limiting case:

s (τ) = q (τ) = e^{- τ}

(18)

We also examine the following two special cases, which are useful as benchmarks. The first benchmark case is superlinear,

q (τ) = {(s (τ))}^{2}, i (τ) = 1

, for which:

q (τ) = \frac{1}{q_{0}} {(1 - \frac{2 (1 - \sqrt{q_{0}})}{(1 + \sqrt{q_{0}}) ⅇ^{\frac{2 τ}{\sqrt{q_{0}}}} + (1 - \sqrt{q_{0}})})}^{2}, s (τ) = \sqrt{q (τ)}

(19)

The second is sublinear,

q (τ) = {(s (τ))}^{1 / 2}, i (τ) = 1

, with

q (t) = \frac{1}{q_{0}} (W ((q_{0} - 1) ⅇ^{- \frac{q_{0} τ}{2} + q_{0} - 1}) + 1), s (τ) = {(q (τ))}^{2}

(20)

where

W (z)

denotes the Lambert W function, defined as the inverse of the function

z = f (w) ≔ w e^{w}

. The case

q (τ) = {(s (τ))}^{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 [0, L]

, 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

V (x) = Q_{0} / A (x)

. 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

d x

in time

d t = d x / V (x) = A (x) 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 (x) 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

\underline{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 + d t

the mass is

S (t) + I (t) d t - Q (t) d t

. Of the molecules contained in

S (t)

at time

t

, the proportion of those removed at time

t + d t

is found, by neglecting second order differentials, to be

Q (t) 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

{\underline{W} > W

}. The event that it leaves the reservoir in the next elementary interval

d W

{\underline{W} \leq W + d W}

and its probability, conditional on

{\underline{W} > W}

, is equal to the ratio of the mass leaving during this interval,

Q (W) d W

, to the total mass:

P \{\underline{W} \leq W + d W | \underline{W} > W\} = \frac{Q (W) d W}{S (W)}

(22)

Applying the definition of conditional probability, we write

\frac{P \{W < \underline{W} < W + d W\}}{P \{\underline{W} > W\}} = \frac{Q (W) d W}{S (W)}

(23)

The numerator is the probability density multiplied by

d W

, i.e. the derivative of the distribution function

F_{\underline{W}} (W)

. Hence

\frac{F_{\underline{W}}^{'} (W)}{1 - F_{\underline{W}} (W)} = \frac{Q (W)}{S (W)}

(24)

The solution of this differential equation is

F_{\underline{W}} (W) = 1 - \exp (- \int_{0}^{W} \frac{Q (t)}{S (t)} ⅆ t)

(25)

and in dimensionless form, with

\underline{w} : = \frac{\underline{W}}{W_{0}}

(26)

the distribution function of the dimensionless residence time is

F_{\underline{w}} (w) = 1 - \exp (- \int_{0}^{w} \frac{q (τ)}{s (τ)} ⅆ τ)

(27)

F_{\underline{w}} (w) = 1 - \exp (- \int_{0}^{w} {(q (τ))}^{1 - \frac{1}{b}} ⅆ τ)

(28)

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

Linear reservoir (in which $q (τ) / s (τ) = 1$ ), any inflow:

$F_{\underline{w}} (w) = 1 - e^{- w}$

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

$F_{\underline{w}} (w) = 1 - \frac{2 ⅇ^{\frac{w}{\sqrt{q_{0}}}}}{(2 + (ⅇ^{\frac{2 w}{\sqrt{q_{0}}}} - 1) (1 + \sqrt{q_{0}}))}$

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

$F_{\underline{w}} (w) = 1 - \frac{{W ((q_{0} - 1) ⅇ^{- \frac{q_{0} w}{2} + q_{0} - 1})}^{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_{\underline{w}} (w) / (1 - F_{\underline{w}} (w))

on a logarithmic axis, instead of the distribution function

F_{\underline{w}} (w)

, 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 (τ) / s (t)

appearing in Equation (27) we propose the approximation:

r (τ) ≔ \frac{q (τ)}{s (τ)} = {(q (τ))}^{1 - \frac{1}{b}} \approx d - \frac{A c}{1 + c τ} = : r_{A} (ρ)

(32)

where

A, c, d

are parameters. Then the integral in Equation (27) becomes

\int_{0}^{w} \frac{q (τ)}{s (τ)} ⅆ τ \approx \int_{0}^{w} (d - \frac{A c}{1 + c τ}) ⅆ τ = d w - \ln {(1 + c w)}^{A}

(33)

so that Equation (27) becomes

F_{\underline{w}} (w) = 1 - e^{- d w} {(1 + c w)}^{A}

(34)

The distribution is interesting, and its domain is

w \geq 0

c \geq 0

and

0 \leq w < - 1 / c

c \geq 0

. In both cases it has all its moments finite if

d \neq 0

. For the special case

d = 0

, the distribution becomes

F_{\underline{w}} (w) = 1 - {(1 + c w)}^{A}

(35)

To find the parameters, we match the true values of

{(q (τ))}^{1 - 1 / b}

τ = 0

and

\infty

and its derivative at

τ = 0

. The latter is found by approximating

q (τ)

from Equation (11). The obtained equations are

1 = d - A c, {(\frac{1}{q_{0}})}^{1 - \frac{1}{b}} = d, (b - 1) (\frac{1}{q_{0}} - 1) = A c^{2}

(36)

and their solution is

d = q_{0}^{\frac{1}{b} - 1}, c = \frac{b - 1}{d - 1} (\frac{1}{q_{0}} - 1), 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_{\underline{w}}^{'} (w)

we find that the mean is

μ_{w} = \{\begin{array}{l} \frac{1}{d}, & A = 0 \\ \frac{ⅇ^{d / c} Γ (A + 1, d / c)}{c {(d / c)}^{A + 1}}, & c > 0, d > 0 \\ {d ⅇ}^{- d} & c = 0, d > 0 \\ \frac{ⅇ^{d / c} (Γ (A + 1, d / c) - Γ (A + 1))}{c {(d / c)}^{A + 1}}, & c < 0, d > 0, A > 0 \\ \frac{- 1}{(A + 1) c}, & d = 0 a n d \{\begin{array}{l} c > 0, A < - 1 \\ o r c < 0, A > 0 \end{array} \end{array}

(38)

while in subcases not included above (e.g.

c > 0, A > - 1

), the mean diverges to infinity. From the equation

F_{\underline{w}} (w_{1 / 2}) = 1 / 2

the median is found to be

w_{1 / 2} = \{\begin{array}{l} \frac{\ln 2}{d}, & c = 0 o r A = 0 \\ \frac{2^{- \frac{1}{A}} - 1}{c}, & d = 0 \\ - \frac{A}{d} W (- \frac{2^{- \frac{1}{A}} d ⅇ^{- \frac{d}{A c}}}{A c}) - \frac{1}{c}, & o t h e r w i s e \end{array}

(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

μ_{w} = 1, w_{1 / 2} = \ln 2

. In the exponential distribution we have

F_{\underline{w}} (μ_{w}) = 1 - 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}, i . e . μ_{W} = W_{0} μ_{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:

μ_{w} = b^{0.45 \ln (q_{0} - a)}, w_{1 / 2} = \ln 2 b^{γ \ln q_{0}}, γ = \{\begin{matrix} 0.2, & b \leq 1 \\ 0.25, & b \geq 1 \end{matrix}

(40)

where the approximation of

μ_{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

μ_{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 (η)

for dimensionless time lag

η = h / W_{0}

is identical to the output function

q (τ)

for

η = τ

, which results from the impulse of an otherwise empty reservoir; and (c) the system is causal, which means that

g (η) = 0

for

η < 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

μ_{η}

and the median

η_{1 / 2}

of the function g(η), namely:

μ_{η} ≔ \int_{0}^{\infty} η g (η) d η, \int_{0}^{η_{1 / 2}} g (η) d η = \frac{1}{2}

(41)

Notice that the median

η_{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 (η) = ((b - 1) η + 1)}^{\frac{b}{1 - b}}

(42)

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

μ_{η} = \frac{1}{2 - b}, η_{1 / 2} = \frac{2^{b - 1} - 1}{b - 1}

(43)

Notice that, for

b \geq 2

, the mean

μ_{η}

diverges to infinity. It is stressed that the quantity

η_{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

ζ_{1 / 2}

and implicitly defined as

q (ζ_{1 / 2}) = 1 / 2

). For completeness, the resulting expression for the latter is

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

(44)

While the median

η_{1 / 2}

(similar to the mean

μ_{η}

) is an increasing function of

b

ζ_{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:

μ_{η} = 1, η_{1 / 2} = ζ_{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 (η)

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

τ = w = 0

. Using Equations (27) and (17) we find that its exact distribution function is

F_{\underline{w}} (w) = 1 - {(1 + (b - 1) w)}^{\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

\underline{w}

is also infinite.

Furthermore, it can be observed that

F_{\underline{w}} (w) = 1 - s (w) = \int_{0}^{w} q (τ) d t

(47)

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

f_{\underline{w}} (w) = q (w)

(48)

This is not a coincidence, but it can be easily shown that it holds for any function

q (τ) = f (s (τ))

. 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 (τ) = f (s (τ))

). 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:

Ω : = \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

Ω

are [

M^{b - 1} 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 (h) = a_{0} + \sum_{i = 1}^{3} a_{i} 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}^{(T_{t} - T_{0}) / (10 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 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 (t)

Q (S) = \frac{S_{0}}{W (t)} {(\frac{S}{S_{0}})}^{b}

(53)

In a first phase (Phase 1), we made several preliminary model runs with different parameterizations, including mathematical expressions of

Q (S)

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 (t)

W (t) = A {(\cos (2 π t + φ) + ψ)}^{b}

(54)

where the parameter

A

has dimensions of time (as does

W (t)

) and the parameters

φ

(phase) and

ψ

are dimensionless. Hence:

Q (t) = \frac{S_{0}}{A} {(\frac{S (t)}{S_{0} (\cos (2 π t + φ) + ψ)})}^{b}

(55)

For two times

t_{1}

and

t_{2}

differing by

2 k π t

, where

k

is an integer (meaning: referring to the same date in different years), we have

Q (t_{1}) / S {(t_{1})}^{b} = Q (t_{2}) / S {(t_{2})}^{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_{N} (t)

, which is not measured:

I_{N} (t) = \frac{S_{0}}{A_{I}} {(\frac{S (t)}{S_{0} (\cos (2 π t + φ_{I}) + ψ_{I})})}^{b_{I}}

(56)

with parameters

b_{I}, A_{I}, φ_{I}, ψ_{I}

similar to those in the expression of outflow. This defines a characteristic time for input:

W_{I} (t) = A_{I} {(\cos (2 π t + φ_{I}) + ψ_{I})}^{b_{I}}

(57)

To find the total input we add the anthropogenic emissions,

I_{A} (t)

, which are known as described is Section 3.2:

I (t) = I_{N} (t) + I_{A} (t)

(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 Δ

and to avoid an implicit numerical scheme, we reduce the time step to half monthly (

Δ \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 ((j + 1) Δ) - S (j Δ)}{Δ} + Q (j Δ) = I_{N} (j Δ) + I_{A} (j Δ)

(59)

from which we find the storage at the step

(j + 1) Δ

from the values of the three variables at step

j Δ

, i.e.,

S ((j + 1) Δ) = S (j Δ) + N (j Δ) Δ

(60)

where

N (j Δ)

is the net inflow at step

j Δ

N (j Δ) : = I_{N} (j Δ) + I_{A} (j Δ) - Q (j Δ)

(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 (j Δ)

and

N (j Δ)

. 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

\hat{S} (j Δ)

and

\hat{N} (j Δ)

, are determined from Equations (55) – (61). After completing a simulation with a specified parameter set, we determine the explained variance for each quantity,

S (j Δ)

and

N (j Δ)

, as:

{E V}_{S} : = 1 - \frac{v a r [\hat{S} (j Δ) - S (j Δ)]}{v a r [S (j Δ)]}, {E V}_{N} : = 1 - \frac{v a r [\hat{N} (j Δ) - N (j Δ)]}{v a r [N (j Δ)]}

(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 (j Δ)

are extraordinarily high (higher than 0.997), while those of

N (j Δ)

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 (t) = I (t) - 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 (t) \equiv [C O_{2}]

, 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 (t) \equiv [C O_{2}]

, 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 (t) = I (t) - 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 (t) \equiv [C O_{2}]

, 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 (t) \equiv [C O_{2}]

, 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}

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_{m} ≔ \frac{\int_{0}^{1} S (t) d t}{\int_{0}^{1} Q (S (t)) d t} = \frac{\int_{0}^{1} S (t) d t}{\int_{0}^{1} \frac{S_{0}}{A} {(\frac{S (t)}{S_{0} (\cos (2 π t) + ψ)})}^{b} 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_{m} \approx \frac{S_{0}}{\frac{S_{0}}{A} \int_{0}^{1} {(\frac{1}{(\cos (2 π t) + ψ)})}^{b} d t} = \frac{A}{\int_{0}^{1} {(\frac{1}{(\cos (2 π t) + ψ)})}^{b} 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{ψ^{2} - 1}

W_{m} \approx A \sqrt{ψ^{2} - 1} = \sqrt{(A (ψ - 1)) (A (ψ + 1))} = \sqrt{W_{m i n} W_{m a x}}

(65)

given that the seasonal minimum and maximum values of

W

are, respectively,

W_{m i n} = A (ψ - 1), W_{m a x} = A (ψ + 1)

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

d m_{I} (t)

that entered the atmosphere at time

[t, t + d t]

, there remains a portion equal to

P \{\underline{W} > t_{c} - t\}

, 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

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_{R} = 163

Gt CO₂ or 20.9 ppm, while

M_{A} = 2612

Gt CO₂ or 334.9 ppm, so that

M_{R} / M_{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_{\underline{W}} (1000) = 1 - F_{\underline{w}} (1000 / 4) = e^{- 250} = 10^{- 108.6}$ , where we have used Equation (29) to evaluate the $F_{\underline{w}} (w)$ .
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$ , $({1 - (1 - p)}^{N}) / p N \to 1$ , for small $p$ (as in our case), we have $p_{1} = p N$ .
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 (44 / 12) = 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} g / (44 g / m o l) = 7.1 \times 10^{16} m o l$ .
The Avogadro constant is $6.022 \times 10^{23} m o l^{- 1}$ and thus the number of CO₂ molecules in the atmosphere is $N = 7.1 \times 10^{16} m o l \times 6.022 \times 10^{23} m o 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} = p N = 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_{\underline{W}} (t) = 1 - F_{\underline{w}} (t / 4.0) = 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

\underline{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].

Supplementary Materials

Not applicable.

Author Contributions

Not applicable.

Funding

This research received no external funding, but was conducted out of scientific curiosity.

Data Availability Statement

This research uses no new data. The data sets used have been retrieved from the sources described in detail in the text.

Acknowledgments

I thank Ioannis Benekos and G.-Fivos Sargentis for their discussions during the preparation of this study.

Conflicts of Interest

The author declares no conflicts of interest.

Appendix A.1: Alternative approximations of a sublinear or superlinear reservoir

With reference to the dimensionless form of the nonlinear reservoir dynamics, here we offer four alternative approximations, which are more involved than the first order approximation discussed and used in the body of this paper. The first is a second order approximation, which works for

b > 1 .

Its form preserves the exact values of

q (s) = s^{b}

for

s = 0

and

s = 1

, and the exact value of its derivative for

s = 1

. The expression is

q (s) = (1 + (b - 1) (s - 1)) s

(A1)

and its solution for constant

i (τ) = 1

q (t) = \frac{1}{q_{0}} - \frac{b^{2} q_{0} - 4 b q_{0} + 4 b + 4 q_{0} - 4}{4 (b - 1) q_{0}} {s e c}^{2} (\frac{t \sqrt{{- q_{0} (b - 2)}^{2} - 4 b + 4}}{2 \sqrt{q_{0}}} - {t a n}^{- 1} (\frac{b \sqrt{q_{0}}}{\sqrt{{- q}_{0} {(b - 2)}^{2} - 4 b + 4}}))

(A2)

The second approximation is appropriate for

b < 1

and also preserves the same exact values. It has the fractional form

q (s) = \frac{s}{b + (1 - b) s}

(A3)

and its solution for constant

i (τ) = 1

q (t) = \frac{b + q_{0} - 1}{(b - 1) q_{0} (W (\frac{- q_{0} b + b + q_{0} - 1}{b q_{0}} \exp (- \frac{t {(b + q_{0} - 1)}^{2} + (b - 1) (q_{0} - 1) q_{0}}{b q_{0}^{2}})) + 1)} - \frac{1}{b - 1}

(A4)

An obvious advantage of both these approximations is that they result in

q (s) = 0

when

s = 0

, while the linear approximation in Equation (8) does not have this property. The third approximation does not have this property either (i.e.,

q (0) \neq 0

while it preserves the other two values,

q (1) = 1, q^{'} (1) = b

), but it has the advantage of being easier to handle. Its expression is

q (s) = 1 - b + b ⅇ^{- 1 + s}

(A5)

and its solution for constant

i (τ) = 1

q (t) = 1 + \frac{(1 - \frac{1}{q_{0}}) (1 - ⅇ^{(b - 1 + \frac{1}{q_{0}}) t})}{ⅇ^{(b - 1 + \frac{1}{q_{0}}) t} - \frac{1}{b} (1 - \frac{1}{q_{0}})}

(A6)

Next we study an approximation in which the reservoir is decomposed into two linear reservoirs with characteristic times

W_{1}, W_{2}

, storages

S_{1} (t), S_{2} (t),

and outflows

Q_{1} (t), Q_{2} (t)

, so that the totals equal the quantities of the real reservoir:

S_{1} (t) + S_{2} (t) = S (t), Q_{1} (t) + Q_{2} (t) = Q (t), a_{1} + a_{2} = 1, W_{1} \frac{d Q_{1} (t)}{d t} + Q_{1} (t) = a_{1} I (t), W_{2} \frac{d Q_{2} (t)}{d t} + Q_{2} (t) = a_{2} I (t)

(A7)

To specify the constituent reservoirs, we study them for the case that the inflow is zero

W_{1} \frac{d Q_{1} (t)}{d t} + Q_{1} (t) = 0, W_{2} \frac{d Q_{2} (t)}{d t} + Q_{2} (t) = 0

(A8)

with solutions

Q_{1} (t) = a_{1} Q_{0} e^{- \frac{t}{W_{1}}}, Q_{2} (t) = a_{2} Q_{0} e^{- \frac{t}{W_{2}}}

(A9)

where we have assumed that the initial conditions in each of the two reservoirs are

Q_{1} (0) = a_{1} Q_{0}, Q_{2} (0) = a_{2} Q_{0}, Q_{0} = Q (0)

On the other hand, the exact solution is this case (homogenous differential equation exists (see Equation (17)) and is

{Q (t) = Q_{0} ((b - 1) \frac{t}{W} + 1)}^{\frac{b}{1 - b}}

(A10)

Therefore, the problem is to specify the weights

a_{1}, a_{2}, a_{1} + a_{2} = 1

and the characteristic times

W_{1}, W_{2}

that best approximate

Q (t)

. To this aim we determine the first three derivatives of the exact solution and the approximation at

t = 0

and equate them. After the algebraic manipulations we find

a_{1,2} = \frac{1}{2} \pm \frac{\sqrt{(b - 1) (2 b - 1)}}{2 (2 b - 1)}, W_{1,2} = \frac{2 a_{1,2}}{b}

(A11)

The approximation works for

b > 1

, so that the quantity in the square root be positive. In the case of

b < 1

there is no meaning in using exponential functions for the approximation, as the domain for which

Q (t) > 0

is finite.

In a similar manner, we can make approximations with more than two linear reservoirs (see also Appendix A.2), but we should have in mind that these are approximations and not exact solutions. In contrast, a differential equation of a reservoir, which is a first order equation, generally does not have an exact solution that is a sum of exponential functions.

Appendix A.2: Notes on the Sum of Exponential Functions as a Response Function

To illustrate the case that the response function is a sum of exponential functions, we consider the IRF of Equation (50) as given by Joos et al. [37], with their original coefficients, which are reproduced in Table A1 . Here we note that the coefficients are given in [37 (Table 5)] without units, with the four

a_{i}

adding up to 1 (except that the table caption contains some instructions with coefficients to multiply and get certain units). This vagueness does not affect our calculations below, as we are only interested in temporal characteristics of their proposed IRF.

Table A1. Parameters of the IRF of Joos et al. [37] as given for CO₂ in their table 5 and as used in both the IPCC AR5 and IR6 [46].

Term	$i = 0$	$i = 1$	$i = 2$	$i = 3$
$Coefficient a_{i}$	0.2173	0.224	0.2824	0.2763
$Coefficient W_{i}$ (years)	∞	394.4	36.54	4.304

With the parameter values in Table A1, the IRF is

Q (h) = 0.2173 + 0.224 ⅇ^{- h / 394.4} + 0.2824 ⅇ^{- h / 36.54} + 0.2763 ⅇ^{- h / 4.304}

(A12)

With some algebra, we may find that, if we erase the constant term

a_{0}

, Equation (A12) is a solution of the linear differential equation

\frac{d^{3} Q (t)}{d t^{3}} - 0.2622 \frac{d^{2} Q (t)}{d t^{2}} + 0.007017 \frac{d Q (t)}{d t} - 0.00001612 (t) = 0

(A13)

This is a third-order linear equation and does not typically represent the dynamics of a reservoir. Instead, this solution can be the output of a cascade of three linear reservoirs, where the first has characteristic time 4.304 years and input 0, the second has characteristic time 36.54 years and its input is the output of the first reservoir, and the third has characteristic time 394.4 years and its input is the output of the second reservoir. While the notion of a cascade has been useful in hydrological applications, i.e. in representing a basin as a cascade of reservoirs, it can hardly be imagined as a proper model for the atmosphere in the carbon balance case of the atmosphere as a whole.

It is stressed that, in the above case, the reservoirs of the cascade are arranged in series and not in parallel (which would make more sense). If different sinks operate in parallel at a single linear reservoir, with each one having a characteristic time

W_{i}

, which means that each of the outputs is proportional to the storage:

Q_{i} = \frac{S}{W_{i}}

(A14)

then the overall reservoir dynamics will be

\frac{d S (t)}{d t} + \sum_{i} \frac{S (t)}{W_{i}} = I (t)

(A15)

\frac{d S (t)}{d t} + \frac{1}{W_{0}} S (t) = I (t), \frac{1}{W_{0}} = \sum_{i} \frac{1}{W_{i}}

(A16)

Obviously, this is the dynamics of a simple linear reservoir with characteristic time equal to the harmonic mean of

W_{i}

. For the case examined with the coefficients in Table A1, the harmonic mean is 11.4 years.

On the other hand, in the cascade of reservoir in series, characterized by the differential equation (A13) and its solution in Equation (A12), if we omit the constant term

a_{0}

, following Equation (41), the mean is given by

μ_{h} ≔ \int_{0}^{\infty} h (Q (h) - a_{0}) d h / \int_{0}^{\infty} (Q (h) - a_{0}) d h = \frac{35 225}{99.9} = 353 y e a r s

(A17)

Now, if we take into account the constant term

a_{0} =

0.2173 in Equation (A12) (i.e. do not subtract it from

Q (h)

), then apparently the mean response time becomes infinite. However, this is not the only problem caused by the constant term. An IRF is the response from an instantaneous impulse, which has infinite mass flow rate, but finite (unit) total mass. On the other hand, the total outflow integral due to this constant term only is

\int_{0}^{\infty} 0.2173 d h = \infty

. Therefore, there is a violation of the mass balance, with an imbalance equal to infinity. A third problem is that if we solve the reservoir differential equation (1) in its homogenous form, there will appear a linear term in

S (t)

, namely

- 0.2173 t

, which at some time (depending on the initial condition of storage) will make the storage negative.

Since Joos et al. [37] state that they fitted their model for lags < 1000 years, it is useful to repeat the above calculation including their constant

a_{0}

, but replacing infinity with the upper limit of their fitting, i.e., 1000 years. This gives

μ_{h} ≔ \int_{0}^{1000} h Q (h) d h / \int_{0}^{1000} Q (h) d h = \frac{134 116}{310.2} = 432 y e a r s

(A18)

This can be regarded as the lowest possible value implied by Joos et al. model, assuming that their response function vanishes off beyond their fitting limit of 1000 years (which however is not stated in their paper).

References

Stendhal [penname of Henri Marie Beyle]. The Life of Henry Brulard [Vie de Henry Brulard, translated to English by J. Sturrock], ISBN 978-1-681371-22-1, The New York Review of Books, New York, NY, USA, 1995. Available online, https://archive.org/details/stendhal-life-of-henry-brulard/ (accessed on 30 April 2024).
Müller-Neuhof, B.; Betts, A.V.G.; Wilcox, G. Jawa, Northeastern Jordan: the first 14C dates for the early occupation phase. Zeitschrift für Orient-Archäologie 2015, 8, 124–131. [Google Scholar]
Fahlbusch, H. Early dams. Proceedings of the ICE - Engineering History and Heritage 2009, 162, 13–18. [Google Scholar] [CrossRef]
Makarewicz, C.A.; Finlayson, B. Constructing community in the Neolithic of southern Jordan: Quotidian practice in communal architecture. PLOS ONE 2018, 13(6), e0193712. [Google Scholar] [CrossRef]
Bible, Genesis , 41, https://www.bible.com/el/bible/2503/GEN.41.GRCBRENT?parallel=392. (accessed on 30 April 2024).
Bell, B. Climate and the History of Egypt: The Middle Kingdom. American Journal of Archaeology 1975, 79, 223–269. [Google Scholar] [CrossRef]
Said, R. The River Nile: Geology, Hydrology and Utilization; Pergamon Press: Oxford, UK, 1993. [Google Scholar]
Hurst, H.E. Long term storage capacities of reservoirs. Trans. Am. Soc. Civil Eng. 1951, 116, 776–808. [Google Scholar] [CrossRef]
Mandelbrot, B.B.; Wallis, J.R. Noah, Joseph, and operational hydrology. Water Resources Research 1968, 4(5), 909–918. [Google Scholar] [CrossRef]
Koutsoyiannis, D. Stochastics of Hydroclimatic Extremes - A Cool Look at Risk, Edition 3, ISBN: 978-618-85370-0-2, 391 pp., Kallipos Open Academic Editions, Athens, 2023. [CrossRef]
Klemeš, V. Watershed as semiinfinite storage reservoir. Journal of the Irrigation and Drainage Division 1973, 99(4), 477–491. [Google Scholar] [CrossRef]
Liu, Z.; Todini, E. Assessing the TOPKAPI non-linear reservoir cascade approximation by means of a characteristic lines solution. Hydrological Processes 2005, 19, 1983–2006. [Google Scholar] [CrossRef]
Fenicia, F.; Savenije, H.H.G.; Matgen, P.; Pfister, L. Is the groundwater reservoir linear? Learning from data in hydrological modelling, Hydrol. Earth Syst. Sci. 2006, 10, 139–150. [Google Scholar] [CrossRef]
Koutsoyiannis, D. Revisiting the global hydrological cycle: is it intensifying? Hydrol. Earth Syst. Sci. 2020, 24, 3899–3932. [Google Scholar] [CrossRef]
Nash, J.E.; Farrel, J.P. A graphical solution of the flood-routing equation for linear storage-discharge relation. Eos, Transactions American Geophysical Union 1955, 36, 319–320. [Google Scholar]
Langbein, W.B. Queuing theory and water storage. Journal of the Hydraulics Division 1958, 84, 1–24. [Google Scholar] [CrossRef]
Nash, J.E. Systematic determination of unit hydrograph parameters. Journal of Geophysical Research 1959, 64, 111–115. [Google Scholar] [CrossRef]
Dooge, J. Linear Theory of Hydrologic Systems. (No. 1468) Agricultural Research Service, US Department of Agriculture, Washington DC, 1973, https://books.google.com/books?id=iVgTfUhBi2gC (accessed on 20 April 2024).
Klemeš, V. Probability distribution of outflow from a linear reservoir. Journal of Hydrology 1974, 21, 305–314. [Google Scholar] [CrossRef]
Klemeš, V.; Borůvka, L. Output from a cascade of discrete linear reservoirs with stochastic input. Journal of Hydrology 1975, 27, 1–13. [Google Scholar] [CrossRef]
Basha, H.A. Nonlinear reservoir routing: particular analytical solution. Journal of Hydraulic Engineering 1994, 120(5), 624–632. [Google Scholar] [CrossRef]
Paik, K. Analytical derivation of reservoir routing and hydrological risk evaluation of detention basins. Journal of Hydrology 2008, 352, 191–201. [Google Scholar] [CrossRef]
Nematollahi, B.; Niazkar, M.; Talebbeydokhti, N. Analytical and numerical solutions to level pool routing equations for simplified shapes of inflow hydrographs. Iranian Journal of Science and Technology, Transactions of Civil Engineering 2021, 1–15. [Google Scholar]
Basha, H.A. Routing equations for detention reservoirs. Journal of Hydraulic Engineering 1995, 121, 885–888. [Google Scholar] [CrossRef]
Berry, E.X. Human CO2 emissions have little effect on atmospheric CO2. International Journal of Atmospheric and Oceanic Sciences 2019, 3, 13–26. [Google Scholar] [CrossRef]
IPCC. Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change [Stocker, T.F., Qin, D., Plattner, G.-K., Tignor, M., Allen, S.K., Boschung, J., Nauels, A., Xia, Y., Bex V., Midgley, P.M., Eds.]. Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, 1535 pp. 2013.
IPCC. Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change [Masson-Delmotte, V., Zhai, P., Pirani, A., Connors, S.L., Péan, C., Berger, S., Caud, N., Chen, Y., Goldfarb, L., Gomis, M.I., Huang, M., Leitzell, K., Lonnoy, E., Matthews, J.B.R., Maycock, T.K., Waterfield, T., Yelekçi, O., Yu, R., Zhou, B. Eds.]. Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, 2391 pp., 2021. [CrossRef]
Koutsoyiannis, D.; Onof, C.; Christofides, A.; Kundzewicz, Z.W. Revisiting causality using stochastics: 1. Theory. Proc. R. Soc. A 2022, 478, 20210836. [Google Scholar] [CrossRef]
Koutsoyiannis, D.; Onof, C.; Christofides, A.; Kundzewicz, Z.W. Revisiting causality using stochastics: 2. Applications, Proc. R. Soc. A, 2022, 478, 20210835. [Google Scholar] [CrossRef]
Koutsoyiannis, D.; Onof, C.; Kundzewicz, Z.W.; Christofides, A. On Hens, Eggs, Temperatures and CO₂: Causal Links in Earth’s Atmosphere. Sci 2023, 5, 35. [Google Scholar] [CrossRef]
Salby, M.L. Physics of the Atmosphere and Climate; Cambridge University Press: New York, NY, USA, 2012. [Google Scholar]
Humlum, O.; Stordahl, K.; Solheim, J.E. The phase relation between atmospheric carbon dioxide and global temperature. Global and Planetary Change 2013, 100, 51–69. [Google Scholar] [CrossRef]
Harde, H. Scrutinizing the carbon cycle and CO2 residence time in the atmosphere. Global and Planetary Change 2017, 152, 19–26. [Google Scholar] [CrossRef]
Harde, H. What humans contribute to atmospheric CO2: Comparison of carbon cycle models with observations. Earth Sciences 2019, 8, 139–159. [Google Scholar] [CrossRef]
Poyet, P. The Rational Climate e-Book. 2nd Edition, ISBN: 978-99957-1-929-62022, Available online, https://www.researchgate.net/publication/347150306 (accessed on 30 April 2024).
Stallinga, P. Residence time vs. adjustment time of carbon dioxide in the atmosphere. Entropy 2023, 25, 384. [Google Scholar] [CrossRef] [PubMed]
Joos, F.; Roth, R.; Fuglestvedt, J.S.; Peters, G.P.; Enting, I.G.; von Bloh, W.; Brovkin, V.; Burke, E.J.; Eby, M.; Edwards, N.R.; Friedrich, T.; Frölicher, T.L.; Halloran, P.R.; Holden, P.B.; Jones, C.; Kleinen, T.; Mackenzie, F.T.; Matsumoto, K.; Meinshausen, M.; Plattner, G.-K.; Reisinger, A.; Segschneider, J.; Shaffer, G.; Steinacher, M.; Strassmann, K.; Tanaka, K.; Timmermann, A.; Weaver, A.J. Carbon dioxide and climate impulse response functions for the computation of greenhouse gas metrics: a multi-model analysis. Atmos. Chem. Phys. 2013, 13, 2793–2825. [Google Scholar] [CrossRef]
Joos, F.; Bruno, M.; Fink, R.; Siegenthaler, U.; Stocker, T.F.; Le Quere, C.; Sarmiento, J.L. An efficient and accurate representation of complex oceanic and biospheric models of anthropogenic carbon uptake. Tellus B 1996, 48, 397–417. [Google Scholar] [CrossRef]
Archer, D.; Brovkin, V. The millennial atmospheric lifetime of anthropogenic CO2. Climatic Change 2008, 90, 283–297. [Google Scholar] [CrossRef]
Archer, D.; Eby, M.; Brovkin, V.; Ridgwell, A.; Cao, L.; Mikolajewicz, U.; Caldeira, K.; Matsumoto, K.; Munhoven, G.; Montenegro, A.; Tokos, K. ,. Atmospheric lifetime of fossil fuel carbon dioxide. Annual Review of Earth and Planetary Sciences 2009, 37, 117–134. [Google Scholar] [CrossRef]
MIT Climate Portal Writing Team featuring guest expert Ed Boyle, How Do We Know How Long Carbon Dioxide Remains in the Atmosphere? MIT Climate Portal, 2023. Available online, https://climate.mit.edu/ask-mit/how-do-we-know-how-long-carbon-dioxide-remains-atmosphere (accessed on 30 April 2024).
Buis, A. The Atmosphere: Getting a Handle on Carbon Dioxide, NASA's Jet Propulsion Laboratory, 2019. Available online, https://science.nasa.gov/earth/climate-change/greenhouse-gases/the-atmosphere-getting-a-handle-on-carbon-dioxide/ (accessed on 30 April 2024).
Berry, E. Why Climate Change Is a Fraud, Available online, https://edberry.com/why-climate-change-is-a-fraud/ (accessed on 30 April 2024).
Myhre, G.; Shindell, D.; Bréon, F.-M.; Collins, W.; Fuglestvedt, J.; Huang, J.; Koch, D.; Lamarque, J.-F.; Lee, D.; Mendoza, B.; Nakajima, T.; Robock, A.; Stephens, G.; Takemura, T.; Zhang H. Anthropogenic and natural radiative forcing supplementary material. In: Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change [Stocker, T.F., Qin, D., Plattner, G.-K., Tignor, M., Allen, S.K., Boschung, J., Nauels, A., Xia, Y., Bex V., Midgley, P.M., Eds.]. Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, 1535 pp. 2013. Available from https://www.ipcc.ch/report/ar5/wg1/chapter-8sm-anthropogenic-and-natural-radiative-forcing-supplementary-material/ (accessed on 30 April 2024).
Strassmann, K.M.; Joos, F. The Bern Simple Climate Model (BernSCM) v1.0: an extensible and fully documented open-source re-implementation of the Bern reduced-form model for global carbon cycle–climate simulations. Geosci. Model Dev. 2018, 11, 1887–1908. [Google Scholar] [CrossRef]
Luo, X.; Xia, T.; Huang, J.; Xiong, D.; Ridoutt, B. Radiative forcing climate footprints in the agricultural sector: Comparison of models from the IPCC 5th and 6th Assessment Reports. Farming System 2023, 1, 100057. [Google Scholar] [CrossRef]
Keeling, C.D.; Piper, S.C.; Whorf, T.P.; Keeling, R.F. Evolution of natural and anthropogenic fluxes of atmospheric CO2 from 1957 to 2003. Tellus B: Chemical and Physical Meteorology 2011, 63, 1–22. [Google Scholar] [CrossRef]
Scripps CO2 Program, Sampling Station Records, Available online, https://scrippsco2.ucsd.edu/data/atmospheric_co2/sampling_stations.html (accessed on 30 April 2024).
Keeling, C.D.; Piper, S.C.; Bacastow, R.B.; Wahlen, M;. Whorf, T.P.; Heimann, M.; Meijer, H.A. Exchanges of atmospheric CO2 and 13CO2 with the terrestrial biosphere and oceans from 1978 to 2000. I. Global aspects, SIO Reference Series, No. 01-06, Scripps Institution of Oceanography, San Diego, 88 pp., 2001.
Keeling, C.D.; Piper, S.C.; Bacastow, R.B.; Wahlen, M;. Whorf, T.P.; Heimann, M.; Meijer, H.A. Atmospheric CO2 and 13CO2 exchange with the terrestrial biosphere and oceans from 1978 to 2000: observations and carbon cycle implications. In A History of Atmospheric CO2 and its effects on Plants, eds, Ehleringer, J.R.; Cerling, T.E.; Dearing, M.D., 83-113, Springer Verlag, New York, 2005.
Ritchie, H.; Roser, M. CO₂ emissions. OurWorldInData.org, 2020. Available online, https://ourworldindata.org/co2-emissions (accessed on 30 April 2024).
Global Carbon Budget (2023) – with major processing by Our World in Data. “Annual CO₂ emissions – GCB” [dataset]. Global Carbon Project, “Global Carbon Budget” [original data]. Available online, https://ourworldindata.org/co2-and-greenhouse-gas-emissions (accessed on 30 April 2024).
International Energy Agency. CO2 Emissions in 2023, 2024. Available online, https://iea.blob.core.windows.net/assets/33e2badc-b839-4c18-84ce-f6387b3c008f/CO2Emissionsin2023.pdf (accessed on 30 April 2024).
Koutsoyiannis, D.; Kundzewicz, Z.W. Atmospheric temperature and CO2: hen-or-egg causality? Sci 2020, 2, 83. [Google Scholar] [CrossRef]
Patel, K.F.; Bond-Lamberty, B.; Jian, J. l Morris, K.A.; McKever, S.A.; Norris, C.G.; Zheng, J.; Bailey, V.L.. Carbon flux estimates are sensitive to data source: a comparison of field and lab temperature sensitivity data. Environmental Research Letters 2022, 17, 113003. [Google Scholar] [CrossRef]
Robinson, C. Microbial respiration, the engine of ocean deoxygenation. Frontiers in Marine Science 2019, 5, 533. [Google Scholar] [CrossRef]
Granger, C.W. Investigating causal relations by econometric models and cross-spectral methods. Econometrica 1969, 37, 424–438. [Google Scholar] [CrossRef]
Koutsoyiannis, D.; Vournas, C. Revisiting the greenhouse effect—a hydrological perspective. Hydrological Sciences Journal 2024, 69, 151–164. [Google Scholar] [CrossRef]

Figure 1. Comparison of the exact (lines) and approximate (dots) solutions for the indicated cases and for constant inflow. The exact solutions are derived from (left) Equation (19) (superlinear reservoir,

b = 2

) and (right) Equation (20) (sublinear reservoir,

b = 1 / 2

). The approximate solutions are derived from Equation (11).

b = 2

) and (right) Equation (20) (sublinear reservoir,

b = 1 / 2

). The approximate solutions are derived from Equation (11).

Figure 2. Comparison of the exact (lines) and approximate (points) solutions for

b = 2

and (left)

q_{0} = 0.8

and (right)

q_{0} = 1.25

, and for the indicated cases. In the case of the constant inflow (

a = 0

), the exact solution is derived from Equation (19) and the approximate solution from Equation (11). The cases of increasing and decreasing flow correspond to

a = 0.1

and

a = - 0.1

, respectively, and the exact curves were determined by numerical integration, while the approximate solutions were determined by Equation (15).

Figure 2. Comparison of the exact (lines) and approximate (points) solutions for

b = 2

and (left)

q_{0} = 0.8

and (right)

q_{0} = 1.25

, and for the indicated cases. In the case of the constant inflow (

a = 0

), the exact solution is derived from Equation (19) and the approximate solution from Equation (11). The cases of increasing and decreasing flow correspond to

a = 0.1

and

a = - 0.1

, respectively, and the exact curves were determined by numerical integration, while the approximate solutions were determined by Equation (15).

Figure 3. Probability plots of the distribution function of the dimensionless residence time for the indicated cases and constant inflow. The exact distributions (lines) are determined by numerical integration of Equations (29)-(31) and are compared to approximate analytical solutions (dots) determined by Equation (34).

Figure 4. Comparison of the mean dimensionless residence time from exact benchmark solutions (Equations (30) and (31) with numerical integration to find

μ_{w}

) and approximations 1 and 2 (Equations (38) and (40), respectively).

Figure 4. Comparison of the mean dimensionless residence time from exact benchmark solutions (Equations (30) and (31) with numerical integration to find

μ_{w}

) and approximations 1 and 2 (Equations (38) and (40), respectively).

Figure 5. Comparison of the mean dimensionless residence time for the indicated cases. Continuous lines for constant inflow (

a = 0

) were obtained from Equations (30) and (31) (with numerical integration to find

μ_{w}

), and those for increasing (

a = 0.1

) and decreasing (

a = - 0.1

) inflow were found by numerical integration of the exact differential equations (considering a finite upper limit of integration, determined such as to avoid numerical instabilities). Beside each exact curve, that of approximation 2 (Equation (40)) is also plotted.

Figure 5. Comparison of the mean dimensionless residence time for the indicated cases. Continuous lines for constant inflow (

a = 0

) were obtained from Equations (30) and (31) (with numerical integration to find

μ_{w}

), and those for increasing (

a = 0.1

) and decreasing (

a = - 0.1

Figure 6. Comparison of characteristic dimensionless time lags for the indicated cases. Mean and median residence times are determined from approximation 2 (Equation (40)). Response times are determined from the exact solution (Equation (43)). The mean response time in the bottom left graph (

b = 2

) is infinite.

b = 2

) is infinite.

Figure 7. Dimensionless outflow for zero inflow in the indicated cases, representing an impulse response function, with comparison of exact solutions (Equation (17)) and approximate ones (Equation (14)). In the right panel (sublinear reservoir,

b = 1 / 2

), according to the exact solution the reservoir empties at time

τ = 1 / (1 - b) = 2

, while according to the linear reservoir approximation the emptying time is

τ = - \ln (1 - b) / b = 1.39

; beyond emptying time the outflow is zero.

b = 1 / 2

), according to the exact solution the reservoir empties at time

τ = 1 / (1 - b) = 2

, while according to the linear reservoir approximation the emptying time is

τ = - \ln (1 - b) / b = 1.39

; beyond emptying time the outflow is zero.

Figure 8. Probability plots of the distribution function of dimensionless residence time for the indicated cases and zero inflow, with comparison of the exact solution (Equation (46)) with the linear reservoir approximation (Equation (14), which yields

F (w) = \min (1, (1 - ⅇ^{- b w}) / b)

F (w) = \min (1, (1 - ⅇ^{- b w}) / b)

Figure 9. Annual carbon balance in the Earth’s atmosphere, in ppm CO₂/year, based on the IPCC estimates [27 (Figure 5.12)]. The balance of 2.4 ppm CO₂/year is the annual CO₂ accumulation in the atmosphere. The modern natural additions (64.5 + 36.6 – (52.2 + 25.6)) = 23.3 ppm is 4.5 times larger than the human emissions (4.4 + 0.8 = 5.2 ppm). (Adapted from [30].)

Figure 10. Atmospheric carbon dioxide inflows and outflows as a function of carbon dioxide concentration. Large circles correspond to IPCC estimates for 1750 and the recent decade (assigned to 2017) and the line joining them is a power law with slope

b = 0.7

. The other circles have been determined using the

Q_{10}

method as in Koutsoyiannis et al. [30] and the triangles by the same method but using 30% higher values of

Q_{10}

. A power law for output (not drawn in the graph) has a slope of 0.77.

b = 0.7

. The other circles have been determined using the

Q_{10}

method as in Koutsoyiannis et al. [30] and the triangles by the same method but using 30% higher values of

Q_{10}

. A power law for output (not drawn in the graph) has a slope of 0.77.

Figure 11. Optimized values of

b_{I}

(the exponent of the power law of inflow) as a function of the values of

b

(the exponent of the power law of outflow). The optimization was made by fixing

b

to a specified value and letting all other parameters free.

Figure 11. Optimized values of

b_{I}

(the exponent of the power law of inflow) as a function of the values of

b

(the exponent of the power law of outflow). The optimization was made by fixing

b

to a specified value and letting all other parameters free.

Figure 12. Achieved values of the indicated explained variances, after maximizing their sum by fixing

b

to a specified value and letting all other parameters free.

Figure 12. Achieved values of the indicated explained variances, after maximizing their sum by fixing

b

to a specified value and letting all other parameters free.

Figure 13. Comparison of observed and simulated storage,

S (t) \equiv [C O_{2}]

, for (upper) Mauna Loa and (lower) Barrow. The insets show the seasonal variation of the characteristic times

W, W_{I}

Figure 13. Comparison of observed and simulated storage,

S (t) \equiv [C O_{2}]

, for (upper) Mauna Loa and (lower) Barrow. The insets show the seasonal variation of the characteristic times

W, W_{I}

Figure 14. Simulated inflows and outflows for (upper) Mauna Loa and (lower) Barrow.

Figure 15. Comparison of observed and simulated net inflows for (upper) Mauna Loa and (lower) Barrow.

Figure 16. Comparison of observed and simulated (upper) storages and(lower) net inflows for Barrow, in the case that the cosine function of inflow is replaced by an arbitrary function, defined through 24 coordinates.

Figure 19. Comparison of observed and simulated net inflows for (upper) Mauna Loa and (lower) Barrow as in Figure 15, but without biosphere expansion (imaginary case 4).

Table 1. Fitted parameters of Equations (55) and (56).

Site	$b$	$φ$	$A$ (years)	$ψ$	$b_{I}$	$φ_{I}$	$A_{I}$ (years)	$ψ_{I}$
Mauna Loa	1	5.448	1.964	2.117	0.945	5.253	1.454	2.858
Barrow	1	5.758	4.182	1.369	0.945	5.151	3.081	1.633

Table 2. Mean residence times, seasonal (

W_{m i n}, W_{m a x}

) and annual

W_{m}

(in years).

Table 2. Mean residence times, seasonal (

W_{m i n}, W_{m a x}

) and annual

W_{m}

(in years).

Site	$Minimum,$ $W_{m i n} = A (ψ - 1)$	$Maximum W_{m a x} = A (ψ + 1)$	$Arithmetic average, A ψ$	$Theoretical mean, W_{m} = \sqrt{W_{m i n} W_{m a x}}$	$Empirical mean W_{m}$ , beginning year	$Empirical mean W_{m}$ , ending year
Mauna Loa	2.21	6.13	4.17	3.68	3.69	3.69
Barrow	1.55	9.90	5.72	3.91	3.95	3.95

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.

MDPI Initiatives

Important Links

Choose an area of interest and we will send you notifications of new preprints at your preferred frequency.

Disclaimer

Reservoir Routing and Its Application to Atmospheric Carbon Dioxide Balance

Abstract

1. Introduction

2. Theoretical Analysis

2.1 System Components and Determination of Their Temporal Evolution

2.2 Residence Time

2.3 Response Time

2.4 Parameters and their Estimation

3. Application to the Carbon Cycle

3.1 A Summary of the Established Approach

3.1.1 Improper Terminology, Obscureness, Ambiguity and Vagueness

3.1.2 Separate Treatment of CO₂ Depending on its Origin

3.1.2 Inappropriate Modelling and Inconsistent Results

3.2 Data

3.3 Premises of the Application

3.4 Model and its Fitting Methodology

3.5 Results of Final Modelling

3.6. Results for Imaginary Cases

3.7 Residence Times

4. Discussion

5. Conclusions

Supplementary Materials

Author Contributions

Funding

Data Availability Statement

Acknowledgments

Conflicts of Interest

Appendix A.1: Alternative approximations of a sublinear or superlinear reservoir

Appendix A.2: Notes on the Sum of Exponential Functions as a Response Function

References

MDPI Initiatives

Important Links

Subscribe