Optimal strategy for removal of greenhouse gas in the atmosphere to avert global climate crisis

1.   Introduction

In recent years, we have witnessed many destructive climatic events disrupting people's lives in countries around the globe, which can be attributed to global warming. According to climate scientists, human activities through the emission of greenhouse gases has caused global warming, with the global surface temperature reaching $1.1$ $^{\circ}C$ above that of the pre-industrialized era 1850–1900 in 2011, and possibly exceeding $1.5$ $^{\circ}C$ before the year 2030 ^[1]. In general, the composition of greenhouse gases includes carbon dioxide $(CO_2)$, nitrous oxide $(N_2O)$ and methane $(CH_4)$ contributing 97%, with the remaining 3% due to various fluorinated gases ^[2]. While a delicate balance of greenhouse gases in the atmosphere is necessary to maintain the global surface temperature in the habitable range, any excess concentration contributes to global warming. Greenhouse gases form a blanket around the planet, trapping heat energy and thereby warming the planet, causing melting of the polar icecap and glaciers and raising the ocean water temperature, leading to adverse climatic changes. The primary source of $CO_2$ emissions from human activities is from burning fossil fuels for electricity, heat and transportation ^[3]. To a large extent, emissions of methane and nitrous oxide are related to large scale farming and food production. Therefore, feeding the growing world population and maintaining their health and well being, and at the same time controlling greenhouse gas emissions, is a great challenge. There is an urgent need to strategically reduce greenhouse gas emissions and develop sustainable green energy and food production techniques.

Scientific literature is rich with extensive studies on global warming, from modeling of climate change, to estimation of greenhouse gases released from various sources, to development of advanced technologies for absorbing greenhouse gases. This paper focuses on dynamic modeling of greenhouse gas in the atmosphere and its removal. In general, climate modeling is a very complex subject as it involves many interacting geophysical processes, such as radiative heat transfer, diffusion, convection and transport, hydrological cycle, gravity, turbulence and ocean circulation. Accordingly, there are many different types of mathematical models that can be found in the literature ^[4], such as simple heat transfer based models, probabilistic and time series models, atmospheric circulation models and coupled atmospheric and oceanic circulation models ^[5,6]. Balance of heat energy on the earth's surface has been used in ^[7] to develop a climate model considering the heat energy received from the sun and reduced by reflection, diffusion and transport. Sensitivity to climate due to the atmospheric carbon dioxide level has been investigated based on parametric variations of a probability density function ^[8] and the effects of radiative feedback associated with water vapor, cloud, snow and other physical mechanisms ^[9]. Attempts have also been made to quantify greenhouse gas emissions from ground ^[10] and dairy farms ^[11]. Reference ^[12] develops a stochastic model of the earth's climate system considering ocean-atmosphere interaction subject to solar radiation, lunar gravity and mass and energy balance, as well as their uncertainties.

Changes in global vegetation patterns can also be used as a predictor of climate change ^[13], as vegetation growth is affected by climate change ^[14,15,16] and at the same time vegetation helps maintain a critical balance of temperature, precipitation and carbon dioxide in the atmosphere. The survey paper ^[13] reviews dynamic interactions of environmental factors and vegetation growth as a way of predicting climate change; in particular, various reaction-diffusion-advection models ^[17,18] have been used to describe interactions between temperature, water and vegetation density. Reference ^[19] considers temporal nonlocal interactions of vegetation and water in arid or semi-arid regions to understand formation of vegetation patterns. Regression of local climate variations in Ireland onto annual global mean surface temperature has been used in ^[20] to track emergence of global climate change due to human activities. Human activities also affect global water quality as evaluated in ^[21] based on nitrogen and phosphorous concentration in runoff water. For a comprehensive overview of climate change and its effects and mitigation strategies, we refer the reader to the report Climate Change 2023 ^[1] of the Intergovernmental Panel on Climate Change (IPCC) and the websites of the United States Environmental Protection Agency ^[22] and the European Environment Agency ^[23].

This paper is an attempt to quantify greenhouse gas concentration in the atmosphere and its optimal removal strategy based on an infinite dimensional framework. We use the diffusion-advection model to express the greenhouse gas concentration in the atmosphere including greenhouse sources and absorbers in a very general setting. We prove the existence of optimal policies subject to control constraints and present necessary conditions for optimal removal of greenhouse gases over a given time period along with a convergence theorem for its numerical computation. Simulation results are presented that illustrate removal of greenhouse gases in a prescribed manner achieving a desired pre-determined level at the end of the control period. Currently, various passive strategies ^[22] have been advocated to minimize greenhouse gas emissions, such as increasing energy efficiency for industrial and residential applications, switching to renewable energy, increasing fuel efficiency for transportation, and replacing fossil fuel based automobiles by electric vehicles. Active carbon capture, utilization and storage (CCUS) technologies have also gained momentum in recent years with the installation of thirty facilities that are operational worldwide and an additional 164 facilities at various stages of development ^[24,25]. This paper provides a conceptual framework for the determination of the optimal carbon capture rate that may be considered for systematic removal of greenhouse gases over a desired control time period.

The rest of the paper is organized as follows: Section 2 presents a diffusion-advection model of greenhouse gas concentration in the atmosphere considering sources and absorbers in a very general setting, along with the formulation of optimal removal strategies. The existence of a solution of the dynamic model is presented in Section 3 followed by the existence of an optimal control in Section 4. Necessary conditions for optimal control (removal policy) and convergence theorem for its computation are presented in Section 5. Section 6 extends the results that include removal of greenhouse gases through absorption by plants and vegetation in the region. Simulation results are presented in Section 7.

2.   Dynamic model of greenhouse gas concentration

In this paper, we present a dynamic model that describes the space-time concentration of greenhouse gases in the atmosphere. Let us consider a spatial region $\Omega,$ an open bounded subset of $R^3$ with smooth boundary $\partial \Omega$ and let $I \equiv [0, T]$ be a closed bounded time interval. Let $\rho = \rho(t, x), (t, x) \in (I \times \Omega)$ denote the concentration (or spatial density) of greenhouse gas at time $t$ and position $x$. The spread of the pollutant can occur due to diffusion and convection. Let $\mu$ denote the diffusion coefficient and $v = (v_1, v_2, v_3)'$ the wind velocity in the region possibly both dependent on time and space. Then the concentration $\rho$ is given by the solution of the following partial differential equation,

where $f = f(t, x)$ denotes the source of greenhouse gas, $u = u(t, x)$ denotes the control efforts (removal rate) actuated by an appropriate operator $B$ (clarified in the next section), $\rho_0$ denotes the initial density (state) and $\nabla_{n}\rho$ denotes outward normal (spatial) derivative of $\rho$ at any point on the boundary $\partial \Omega$. For any vector valued function $y(x) \in R^3, \; x \in \Omega$, we have used the notation $\bigtriangledown \cdot y$ for the $div(y)$ and for any scalar valued function $z(x) \in R, \; x \in \Omega$, we have used $\triangledown z$ to denote the gradient vector of $z.$

The main objective is to reduce the greenhouse gas concentration over a plan period $I \equiv [0, T]$ by cutting down the emission and physical removal through absorption and other available means. Suppose we want the greenhouse gas concentration reduced by a certain percentage of the current level subject to limited resources while avoiding significant economic impact. A reasonable objective functional is given by

where $\{w_1, w_2\}$ are nonnegative bounded measurable functions on $I\times \Omega$ representing the weights (or cost) assigned to each of the measures of performance and cost of control. Similarly, $w_3$ is a nonnegative bounded measurable function on $\Omega$ penalizing any mismatch between the desired (goal) state and the actual final state reached. The function $\beta_1$ can be chosen as $\beta_1(t) = e^{-\lambda t}, \; \lambda > 0$ and the constant $\beta_2$ as $0 < \beta_2 < 1,$ depending on the percentage reduction desired at the end of the plan period. The problem is to find a control policy that minimizes the cost functional (2.2).

3.   Compact formulation of the control problem

Control theory for distributed parameter systems or equivalently partial differential equations is well developed. Here, we can use the well-known theoretical results to solve the problem stated above. Let $H = L_2(\Omega)$ denote the real Hilbert space with the usual norm $\parallel \varphi\parallel_{H} = \bigl(\int_{\Omega} |\varphi(x)|^2 dx \bigr)^{1/2}$ and $H^1$ the Sobolev space given by

with the norm $\parallel \varphi\parallel_{H^1} = \bigl (\parallel \varphi\parallel_{H}^2 + \parallel \bigtriangledown \varphi \parallel_{H}^2)\bigr)^{1/2}$. Note that $\bigtriangledown \varphi = \{\partial_i\varphi, i = 1, 2, 3\}$. It is clear that $H^1$ is a Hilbert space and that $H^1 \subset H$. It is well known that the topological (continuous) dual of $H_0^1$ is given by $H^{-1}$ and that $H \subset H^{-1}$.

For example, for any $\eta \in H$ and every $\{i = 1, 2, 3\}$, the element $\partial_i\eta \in H^{-1}$. This follows readily from the following inequality,

Also from this inequality, we observe that the embedding $H\hookrightarrow H^{-1}$ is continuous. We denote the space $H_0^1$ by $V$, and its dual $H^{-1}$ by $V^*$. Identifying $H$ with its own dual, this gives us the triple $\{V, H, V^*\}$, known as the Gelfand triple, satisfying $V\hookrightarrow H \hookrightarrow V^*$ with the embeddings being continuous and dense. Now, we consider the time dependent operator

and define the associated operator $A$ subject to the homogeneous Neumann boundary condition as indicated in Eq (2.1). This is given by

For each $t\in I$, $A(t)$ is a bounded linear operator from $V$ to $V^*$, denoted $A \in {\cal L}(V, V^*)$, while it is an unbounded operator on $H$ and, as seen later in Remark 3.2, as a function of $t$, it is a family of bounded operators, denoted $A \in L_{\infty}(I, {\cal L}(V, V^*))$. We prove the following lemma:

Lemma 3.1 Suppose $0 < d \leq \mu(t, x) \leq D, \; (t, x) \in I\times \Omega$ and the wind velocity $v$ and its divergence $\bigtriangledown\cdot v$ are uniformly bounded on $I\times \Omega$. Then, for all $t \in I$, the operator $A(t)$ is coercive in the sense that there exist nonnegative numbers $\alpha, \; \gamma > 0$ such that for every $\varphi \in V$,

and there exists a positive number $c$ such that

Proof The proof is based on integration by parts and the Hölder inequality. Take any $\varphi \in V$ and consider the bilinear form

Integrating by parts the first term of the above expression and using the lower bound of the diffusion coefficient, we obtain

Considering the second term in Eq (3.1), let us first note that

Now taking the scalar product of this equation with $\varphi \in H$, we obtain

Since, by assumption, the wind velocity and its divergence are uniformly bounded on $I\times \Omega$, there exist nonnegative numbers $a, b > 0$ such that

Thus, it follows from the expression (3.3) and Schwartz inequality that

Using the Cauchy inequality, it follows from the above expression that

for any $\varepsilon > 0.$ Using the expressions (3.1), (3.2) and (3.5), we arrive at the following inequality

Since $\varepsilon > 0$ is otherwise arbitrary, we can choose $\varepsilon = \varepsilon_0 < 1$ sufficiently small so that $d > ((a\varepsilon_0)/2).$ With this choice, we define

and, hence, it follows from (3.6) that

This proves the inequality (ⅰ) and hence the coercivity of the family of operators $\{A(t), t \in I\}$. To prove the second inequality (ⅱ), let us take any $\varphi, \; \psi \in V$ and note that

Using the Schwartz inequality and the upper bounds of the diffusion coefficient, the wind velocity and its divergence, one can verify that

For any set of positive numbers $\{c_1, c_2, c_3\},$ it follows from the elementary inequality $(c_1+c_2 +c_3)^2 \leq\; 4(c_1^2 + c_2^2 + c_3^2)$ that

We note that, for our purpose, it is not necessary to find a tight upper bound. Hence, by simply adding the missing terms to make up the $V$-norm in each of the components of the above inequality and taking $c^2\equiv 4(D^2+b^2+a^2)$, we obtain

Hence, for all $t\in I$, we have $| < A(t)\varphi, \psi > |\; \leq \; c \parallel \varphi\parallel_{V} \parallel \psi \parallel_{V}$, proving the bound (ⅱ). This completes the proof.

Remark 3.2 As indicated in the introduction of Section 2, the family of operators $A(t), t \in I$ is unbounded on the Hilbert space $H$. However, since it satisfies the inequality (ⅱ) for arbitrary $\varphi, \psi \in V$, it is clear that for each $t \in I,$ $A(t) \in {\cal L}(V, V^*)$ and, hence, it is a family of bounded linear operators from $V$ to $V^*$. Further, it follows from our assumption that $|v|_{R^3}$, $|div (v)|$ and the diffusion coefficient $\mu$, are all uniformly bounded measurable functions defined on $I\times \Omega$. Thus, we conclude that $A \in L_{\infty}(I, {\cal L}(V, V^*))$.

Using Lemma 3.1, we can rewrite system (2.1) described by the partial differential equation as an abstract ordinary differential equation on the Hilbert space $H$ involving the Gelfand triple $\{V, H, V^*\}$ as follows

Any real (or complex) valued function $h$ defined on $I$ $\times$ $\Omega$ satisfying $\int_{\Omega} |h(t, x)|^2 dx < \infty$, can be described as an $H$ valued function for any $t\in I$. If the integral is essentially bounded on the interval $I$, we express this by stating that $h \in L_{\infty}(I, H).$ Similarly, if the function is square integrable in both the variables, we may express it by stating that $h \in L_2(I, H)$ and write $\int_{I} \parallel h(t)\parallel_{H}^2 dt < \infty.$

For admissible controls, let $E$ denote another Hilbert space and $U\subset E$ a closed bounded convex set, and consider the Hilbert space $L_{2}(I, E).$ Let

denote the set of admissible control policies. This is a closed bounded convex subset of $L_{2}(I, E)$ and, hence, is weakly compact.

We are now prepared to present a theorem stating the existence and uniqueness of solution of the differential Eq (3.11).

Theorem 3.3 Consider the system (3.11) with the operator valued function $A = \{ A(t), t\in I\}$ satisfying the assumptions (ⅰ) and (ⅱ) of Lemma 3.1, and $f \in L_2(I, V^*)$, $B \in L_{\infty}(I, {\cal L}(E, V^*))$ and $u \in L_2(I, E).$ Then, for every initial state $\rho_0 \in H$, the system (3.11) has a unique solution $\rho \in L_{\infty}(I, H) \cap L_2(I, V)$ with $\dot \rho \in L_2(I, V^*).$

Proof. The proof of existence and uniqueness of solution follows from Ahmed and Teo (^[26], Theorem 5.1.1, p278, see also Lions ^[27]). Here, we present a brief account on the regularity of solution, such as $\rho \in L_{\infty}(I, H) \cap L_2(I, V)$ and $\dot \rho \in L_2(I, V^*).$ In general, one can admit $g\equiv (f+ Bu) \in L_2(I, V^*).$ In order to verify the regularity, we scalar multiply the Eq (3.11) by $\rho$ and note that, in the sense of distribution, $d/dt (\rho, \rho)_{H} = 2 < \dot \rho, \rho > _{V^*, V}$ is well defined. Using this fact, we have

Integrating this equation over the interval $[0, t]$, we obtain

Using the coercivity property (ⅰ) of the operator $A$ in the above expression, one can easily verify that

It follows from the Cauchy inequality that

Choosing $\varepsilon = \alpha$ and using the above expression in inequality (3.14), we arrive at the following inequality

By virtue of the Gronwall inequality, it follows from the above expression that

It is clear from the above inequalities that $\rho \in L_{\infty}(I, H) \cap L_2(I, V).$ Considering the last statement of the theorem, we note that, by virtue of the property (ⅱ) (see Lemma 3.1), $A(t)\rho(t) \in V^*$ for almost all $t \in I.$ Thus, it follows from Eq (3.11) that

Hence, we have the following inequality

Since $\rho \in L_2(I, V)$ and $g \in L_2(I, V^*),$ it follows from the above inequality that $\dot \rho \in L_2(I, V^*).$ This proves the existence, uniqueness and the regularity properties of solutions as stated.

In the following corollary, we show that the solution $\rho$ has a stronger and more interesting regularity property, such as $\rho \in C(I, H).$

Corollary 3.4 Under the assumptions of Theorem 3.3, Eq (3.11) has a unique solution $\rho \in C(I, H).$

Proof Let $X$ denote the linear vector space given by

and endow this with the norm topology

With respect to the above norm topology, $X$ is a Hilbert space. This is a special case of (^[26], Theorem 1.2.15, p27), and it follows from this theorem that $X \subset C(I, H).$ Clearly, according to Theorem 3.3, the system (3.11) has a unique solution $\rho \in L_{\infty}(I, H) \cap L_2(I, V)$ and $\dot \rho \in L_2(I, V^*).$ Thus, $\rho \in X$ and, hence, we conclude that $\rho \in C(I, H).$

Corollary 3.5 Under the assumptions of Theorem 3.3, the operator valued function, $A(t), t \in I,$ generates an evolution operator $\{G(t, s), 0 \leq s \leq t \leq T\}$ in the Hilbert space $H$ satisfying

and for each $\rho_0 \in H,$ $f \in L_2(I, V^*)$ and $u \in L_2(I, E),$ the solution $\rho$ of Eq (3.11) has the representation

Proof Existence of an evolution operator $G$ satisfying the properties (a)–(c) follows from the fact that, for any starting time $\tau \in I$ and for any initial state $\xi \in H$, Eq (3.11), with $f\equiv 0$ and $u = 0$, has a unique solution $\rho(t), t\geq \tau$ in $H$ giving $\rho(t) = G(t, \tau)\xi.$ This holds for every $\xi \in H$ and, hence, $G(t, \tau) \in {\cal L}(H)$ for all $\tau \leq t$ and $t \in I.$ Thus $\dot \rho + A(t)\rho = 0$, for $\rho(\tau) = \xi$, $0\leq \tau < t \leq T$. The property (d) follows from the fact that for any $f \in L_2(I, V^*)$ and $Bu \in L_2(I, H) \subset L_2(I, V^*)$, Eq (3.11) has a unique solution $\rho \in C(I, H)$ with $\dot \rho \in L_2(I, V^*).$ This implies that, for any $g \in L_2(I, V^*),$ the convolution

is well defined and $h$ is an element of $C(I, H).$ For this to hold, it is necessary that $G(t, s) \in {\cal L}(V^*, H)$ for all $0\leq s \leq t \leq T$ and

for all $t \in I.$ This completes the proof.

4.   Existence of optimal control

We consider the cost functional given by the expression (2.2) reproduced below as follows,

where the functions $\{w_1, w_2, w_3\}$ and $\{\beta_1, \beta_2\}$ are as described following Eq (2.2). The objective is to reduce the greenhouse gas concentration by a certain percentage of the initial level $\rho_0$ over the plan period $I.$ Denoting $\rho_d(t, x) = \beta_1(t) \rho_0(x)$ and suppressing the space variable, we may express this as $\rho_d(t)\equiv \rho_d(t, \cdot) \in H$ for all $t \in I.$ Thus, it follows from our assumption on $\beta_1$ and $\rho_0$ that $\rho_d \in L_2(I, H).$ Similarly, the terminal cost $\rho_{\tau} = \beta_2 \rho_0(\cdot) \in H.$ Using this notation, we can rewrite the cost functional (4.1) as follows:

Since the set of weights $\{w_1, w_2, w_3\}$ are nonnegative, it is clear that the operators $\{Q_1, Q_2, Q_3\}$ are all nonnegative self adjoint operators in Hilbert spaces $H$, $E$ and $H$, respectively.

Theorem 4.1 Consider the system (3.11) with the cost functional (4.2), and suppose the assumptions of Theorem 3.3 hold with the admissible controls ${\cal U}_{ad}.$ Suppose the operators $\{Q_1, Q_2, Q_3\}$ are nonnegative self adjoint. Then there exists an optimal control.

Proof We present only a brief outline of the proof. For more details on the question of existence of optimal controls, we refer the reader to (^[26], Theorems 5.1.2–5.1.4, p285). Since the admissible set of controls ${\cal U}_{ad}$ is a bounded subset of $L_2(I, E)$, $B \in L_{\infty}(I, {\cal L}(E, V^*))$ and $f \in L_2(I, V^*),$ it is clear that the set $G_{ad} \equiv \{ g: g = f+Bu, u \in {\cal U}_{ad} \}$ is contained in a bounded subset $L_2(I, V^*).$ Hence, it follows from the inequalities (3.15), (3.16) and (3.18) that the set of solutions ${\cal S}_{ad} \equiv \{ \rho(u), u \in {\cal U}_{ad} \}$ is contained in a bounded subset of the vector space $X$ as described in Corollary 3.4. Since $L_2(I, V^*)$ is a Hilbert space (so a reflexive Banach space), the set $G_{ad}$ is relatively weakly sequentially compact. The set of admissible controls ${\cal U}_{ad}$ is compact in the weak topology. Thus, $G_{ad}$ is weakly closed and, hence, a weakly sequentially compact subset of $L_2(I, V^*).$ The system (3.11) is linear and so the control to solution map $u \longrightarrow \rho(u)$ is affine continuous and, hence, weakly continuous. Thus, the first component of the cost functional, being quadratic in $\rho$, is weakly lower semicontinuous. By virtue of Corollary 3.4, we have $\rho \in C(I, H)$ and, hence, $\rho(T) \in H.$ Thus, the third component of the cost functional is a well defined quadratic functional and, hence, it is also weakly lower continuous. The second component of the cost functional is quadratic in control, and so weakly lower semicontinuous on ${\cal U}_{ad}.$ Hence, the cost functional $u \longrightarrow J(u)$ is weakly lower semicontinuous on ${\cal U}_{ad}$. Thus, it follows from weak compactness of the set ${\cal U}_{ad}$ that $J$ attains its minimum on it proving existence of an optimal control. This completes the outline of our proof.

In theorem 4.1, we have used the fact that controls take values in a bounded set $U\subset E.$ This is not essential; we may consider $U$ to be an unbounded set or the entire space $E$ provided an additional condition is satisfied. This is stated in the following Corollary:

Corollary 4.2 Suppose the assumptions of Theorem 4.1 hold with the exception that ${\cal U}_{ad} = L_2(I, E)$ and the operator $Q_2$ is positive self adjoint in the sense that there exists a positive number $\lambda > 0$ such that $\int_I < Q_2(t)u(t), u(t) > _{E}dt \geq \lambda.$ Then there exists a unique optimal control.

Proof It is clear that $J(u)\geq 0$ for $u \in {\cal U}_{ad},$ and it follows from the assumption on the operator valued function $Q_2$ that $\lim_{ \parallel u \parallel_{L_2(I, E)}\rightarrow \infty} J(u) = \infty.$ Hence, there exists a nonnegative number $M$ such that $\inf_{u \in {\cal U}_{ad}}J(u) = M$ is well defined with $M < \infty.$ Let $\{u^n\} \in {\cal U}_{ad}$ be a minimizing sequence so that $\lim_{n\rightarrow \infty}J(u^n) = M.$ Clearly, this implies that $\{u^n\}$ is a bounded sequence in the Hilbert space ${\cal U}_{ad} = L_2(I, E).$ Thus, there exists a subsequence of the sequence $\{u^n\}$, relabeled as the original sequence, and an element $u^o\in {\cal U}_{ad}$ such that $u^n \buildrel w \over\longrightarrow u^o$ (converges in the weak topology). Since $J$ is weakly lower semicontinuous, we have

Thus, we conclude that the following inequalities hold,

Since $u^o \in {\cal U}_{ad}$, it is clear that $M \leq J(u^o).$ Thus, we have $J(u^o) = M$ proving the existence of an optimal control.

Remark 4.3 It is not essential to have a unique optimal control. It is enough if the set of optimal controls is closed. In that case, $w_2$ is not required to be strictly positive having a positive lower bound.

5.   Necessary conditions of optimality and algorithm

It follows from Theorem 4.1 that an optimal control policy exists. Given that an optimal control exists, it is now reasonable to look for and find ways to construct it. In this section, we develop necessary conditions of optimality using which one can construct the optimal control policies.

Theorem 5.1 Consider the system (3.11) with the set of admissible control ${\cal U}_{ad}$ and the cost functional (4.2). Suppose the assumptions of Theorems 3.3 and 4.1 hold. Then, in order for a control $u^o \in {\cal U}_{ad}$ and the corresponding solution $\rho^o \in {\cal S}_{ad}$ to be optimal, it is necessary that there exists a $\psi \in X$ such that the following inequality and the adjoint and state equations are satisfied:

and

where $A^*(t)$ and $B^*(t)$ denote the adjoints to the operators $A(t)$ and $B(t)$, respectively.

Proof Let $u^o \in {\cal U}_{ad}$ denote the optimal control and $u\in {\cal U}_{ad}$ any other control. For any $\varepsilon \in [0, 1],$ let us construct a control of the form $u^{\varepsilon} \equiv u^o+ \varepsilon (u-u^o).$ Since ${\cal U}_{ad}$ is a closed convex set, $u^{\varepsilon} \in {\cal U}_{ad}$. Thus, by virtue of optimality of $u^o$, we have

Let $\rho^{\varepsilon}$ and $\rho^o$ denote the solutions of Eq (3.11) corresponding to controls $u^{\varepsilon}$ and $u^o$, respectively. Letting $\vartheta$ denote the limit

one can easily verify that $\vartheta$ satisfies the following variational equation

Since $B(u-u^o) \in L_2(I, H) \subset L_2(I, V^*),$ it follows from Theorem 3.3 and Corollary 3.4 that this equation has a unique solution $\vartheta \in X$ (see Eq (3.19)). Hence, the control to solution map, $B(u-u^o)\longrightarrow \vartheta,$ is a continuous linear operator from $L_2(I, V^*)$ to the Hilbert space $X$ and, hence, bounded. Next, computing the difference quotient $\lim_{\varepsilon \longrightarrow 0} (1/ \varepsilon) (J(u^{\varepsilon})- J(u^o))$, we obtain the Gâteaux differential of $J$ evaluated at $u^o$ in the direction $(u-u^o)$ satisfying the inequality,

For convenience of presentation, let us introduce the functional

Since $\vartheta \in X,$ it follows from Corollary 3.4 that $\vartheta \in C(I, H) \subset L_2(I, H).$ Thus, recalling that the element $Q_1(\rho^o-\rho_d) \in L_2(I, H)$ in the first component and the element $Q_3(\rho^o(T)-\rho_{\tau})\in H$ in the second component, it is clear that $\vartheta \longrightarrow L(\vartheta)$ is a continuous linear functional on $L_2(I, H).$ Thus, the composition map $\tilde L$ given by

is a continuous linear functional on $L_2(I, V^*)$. Note that the spaces $\{L_2(I, V), \; L_2(I, V^*)\}$ are reflexive and, hence, by duality, there exists a $\psi \in (L_2(I, V^*))^* = L_2(I, V)$ such that

Thus, it follows from the expressions (5.6)–(5.9) that

This proves the necessary condition (5.1). To prove the necessary condition (5.2), we follow the following steps. Considering the scalar product $< \vartheta(t)$, $\psi(t) >$ and differentiating this (in the sense of distribution), we obtain

Using the variational equation in the above expression, we arrive at the following identity

Setting

in the expression (5.12) and integrating, we obtain

It is clear that the expression on the left hand side of the above equation coincides with the functional $L(\vartheta)$ given by Eq (5.7), and the expression on the right coincides with $\tilde L(B(u-u^o))$ given by Eq (5.9) as necessary. Hence, it follows from Eq (5.13) that $\psi$ must satisfy the adjoint Eq (5.2), which is a necessary condition. Equation (5.3) is the system equation corresponding to the optimal control $u^o \in {\cal U}_{ad}$, so there is nothing to prove. This completes the proof of all the necessary conditions of optimality as stated in the theorem.

In order to determine the optimal control policy, we need an algorithm. Here, we present the algorithm including a proof of its convergence.

Theorem 5.2 Suppose the assumptions of Theorem 5.1 hold. Then there exists a sequence of controls $\{u^n\} \subset {\cal U}_{ad}$ (which can be constructed step by step) along which the sequence $\{ J(u^n)\}$ converges monotonically to a minimum.

Proof We follow the following steps:

Step 1: Choose any control $u^1 \in {\cal U}_{ad}$ and solve Eq (5.3) by replacing $u^o$ by $u^1$ giving the solution $\rho^1.$

Step 2: Use the solution $\rho^1$ (from the previous step $\{u^1, \rho^1\}$) in the adjoint Eq (5.2) in place of $\rho^o$ and solve giving $\psi^1.$

Step 3: Use the pair $\{u^1, \psi^1\}$ in place of $\{u^o, \psi\}$ in the expression (5.1) giving

If the above inequality holds, $u^1$ qualifies to be optimal. This is a rare possibility with very little hope to hit the optimal in one shot. So we move to the next step.

Step 4: We choose $\varepsilon > 0$ sufficiently small and construct a new control from the available data $\{u^1, \psi^1\}$ as follows

so that $u^2 \in {\cal U}_{ad}.$ Choosing $u^2$ for $u$ in the expression (5.15), we obtain

Using the Lagrange formula, one can approximate the cost functional corresponding to $u^2$ as

For $\varepsilon > 0$ sufficiently small, it is clear from the above expression that $J(u^1) > J(u^2).$

Step 5: Using $u^2$ and returning to step 1, and repeating the process, we can construct a sequence of controls $\{u^n\}$ and a corresponding sequence of cost functionals $\{J(u^n)\}$ satisfying the following series of inequalities $J(u^1) > J(u^2) > J(u^3) > \cdots J(u^k) > J(u^{k+1} > \cdots.$ Since $J(u) \geq 0, \; \forall\; u \in {\cal U}_{ad}$ and $J(u^n)$ is a monotone decreasing sequence, it converges to a minimum, say, $\lim_{n\rightarrow \infty} J(u^n) = m_0\geq 0.$ This completes the proof.

6.   A simple extension

The model presented in Section 2 does not include the natural absorption of components of greenhouse gas like $CO_2$ by plants and vegetation in the region. We include this factor in the following model,

where $\delta$ is a nonnegative bounded measurable function defined on $I \times \Omega$ denoting the coefficient of natural absorption by the surrounding vegetation in the geographical region. Defining

and replacing $\gamma$ by $\tilde \gamma \equiv \gamma + \hat \delta,$ the inequality (ⅰ) in Lemma 3.1 can be set as,

The inequality (ⅱ) remains unchanged. The system (3.11) is modified as follows,

With these modifications, all the conclusions of Theorem 3.3, Corollary 3.4, Theorem 4.1 and Corollary 4.2 remain valid. As an extension of Theorem 5.1, we have the following necessary conditions of optimality.

Theorem 6.1 Consider the system (6.2) with the set of admissible control ${\cal U}_{ad}$ and the cost functional (4.2). Suppose the assumptions of Theorems 3.3 and 4.1 hold. Then, in order for a control $u^o \in {\cal U}_{ad}$ and the corresponding solution $\rho^o \in X$ to be optimal, it is necessary that there exists a $\psi \in X$ such that the following inequality, the adjoint and the state equation are satisfied:

and

Proof The proof is identical to that of Theorem 5.1

Remark 6.2 A more realistic model for absorption is given by a continuous nonlinear nondecreasing function $h(\rho)$, which equals zero over $(-\infty, 0]$ and positive over $[0, \infty)$ and flattens out beyond a finite interval $[0, \beta] \subset [0, \infty).$ In other words, the absorption rate saturates beyond a level of greenhouse gas concentration. The saturation level depends on the vegetation type and its density in the geographical region. Thus, by improving and increasing plantation, the absorption rate can be improved and the climate crisis can be moderated and possibly avoided.

7.   Simulation results

This section presents results of numerical simulation on optimal removal of greenhouse gases (GHG) based on the necessary conditions of optimality presented in Theorem 5.1. We consider the following system model in $\Omega = \{-L \le x \le L, \, -L \le y \le L, \, 0 \le z \le L_z\}$:

where the diffusion coefficient $\mu$ is assumed to be constant, $\Delta$ is the Laplacian operator, and $\nabla_n \rho$ is the outward normal derivative of $\rho$ on the boundary $\partial\Omega$. The wind velocity is given by $v = [v_x, v_y, 0]$, signifying constant wind velocity only in the $x$ and $y$ directions. The greenhouse gas sources and absorbers may be arbitrarily located in $\Omega$ defined by the functions

where $\xi_k$ and $\zeta_k$ are indicator functions, respectively, defined on the spatial domain. Furthermore, it is assumed that the $k$-th source emits greenhouse gas at the rate $s_k(t)$ and the $k$-th absorber absorbs at the rate $u_k(t)$. Equation (7.2) also implies that all sources and absorbers of greenhouse gases are located near the ground level.

It is assumed that there is an initial concentration of greenhouse gas near the ground level at the center of $\Omega$ and is rapidly decreasing with altitude:

For the cost function, we take

where $\beta_1(t) = e^{-\lambda t}$ defines the desired temporal profile of reduction of greenhouse gas from its initial value, and $\beta_2$ is the desired reduction of greenhouse gas at the final time which is taken as $\beta_2 = \beta_1(T)$.

For simulation, we consider the normalized equation with $L = 1$, $L_z = 0.1$, $\mu = 1$, $v_x = 2$, $v_y = 1$ (and $v_z = 0$), all in normalized units. For the initial gas distribution given in Eq (7.3), consider $\alpha_x = \alpha_y = 1000$ and $\alpha_z = 200$. The time interval is taken as 2,000 time steps in normalized time unit. Various weights in the cost function are taken as $w_1 = 10^8, \, w_2 = 1$, $w_3 = 10^7$ and $\beta_2 = 0.5$, which signifies a desired 50% reduction in greenhouse gas concentration at the end of the control period from which one can easily compute $\lambda$ and $\beta_1(t)$ for the temporal profile of reduction of greenhouse gas.

Without any loss of generality, for a proof of concept we consider only one source of greenhouse gas located at the center of $\Omega$ defined by $g(x, y, z) = \xi_{x, y}e^{-50z}, \, (x, y, z) \in [-0.125L, 0.125L]\times[-0.125L, 0.125L]\times[0, L_z]$ with the source $s_1(t) = 1$ for all $t\in[0, T]$. For the absorber of the greenhouse gas, we consider three cases:

● Case Ⅰ: No control, i.e., no absorption of greenhouse gas. This provides the worst case scenario of greenhouse gas buildup in the atmosphere if no absorption takes place.

● Case Ⅱ: Greenhouse gas is absorbed uniformly over the entire region $\Omega$, which may be implemented by extensive vegetation growth.

● Case Ⅲ: Greenhouse gas is absorbed over a small region at the center of $\Omega$ defined by $h(x, y, z) = \zeta_{x, y}e^{-50z}, \, (x, y, z) \in [-0.25L, 0.25L]\times[-0.25L, 0.25L]\times[0, L_z]$. Carbon sequestration technologies ^[24,25] could be used for absorption of $CO_2$, which constitutes the largest component of greenhouse gases in the atmosphere.

Figure 1 summarizes the effects of optimal absorption of greenhouse gas and the corresponding absorption policy. The left side figure shows that total greenhouse gas (i.e., $\int_{\Omega} \rho(x, y, z) \, d\Omega)$ in the region increases with time (Case Ⅰ), whereas it reduces as desired with optimal absorption for Case Ⅱ and Case Ⅲ. The right side figure shows the corresponding optimal absorption policies. Note that the required absorption rate for Case Ⅲ is higher than that for Case Ⅱ which is expected since in Case Ⅲ, the greenhouse gas is absorbed at a higher rate albeit over a small region located at the center of $\Omega$. It is worthwhile to note that the goal of optimal absorption is to reduce the total greenhouse gas concentration by about 50% at the end of the control period which was achieved.

Figure 2 shows temporal variations of greenhouse gas density at the ground level, i.e., $z = 0$, for several time instants for Case Ⅱ. This figure clearly shows both spatial diffusion of the greenhouse gas as well as absorption with time. For higher altitudes, results on diffusion and absorption are very similar to those in Figure 2 except for the magnitude, which are omitted here for the sake of brevity. Reduction of greenhouse gas at higher altitudes due to optimal absorption can be better observed in Figure 3. Note that the initial greenhouse gas density (i.e., at $t = 0$) is higher at the ground level with decreased density at higher altitude. With optimal absorption, the gas density is reduced at all altitudes with time. There is only minor reduction of GHG after time $t = 1000\Delta t$ as expected in light of Figure 1. For Case Ⅲ, the results are similar to that of Case Ⅱ except for small differences and are omitted.

Figure 4 is a two-dimensional rendering of greenhouse gas density at the ground level at the final time $t = T$ that illustrates the effects of wind velocity. The left side figure shows the gas density variation with $x$ at $y = 0$ at the end of the control period. Similarly, the right side figure is for the gas density variation with $y$ at $x = 0$ at the end of the control period. Compared to the initial peak density of $\rho = 1$ at the center of $\Omega$ at ground level, Case Ⅰ shows uncontrolled diffusion of GHG whereas with optimal absorption (Cases Ⅱ and Ⅲ) gas density is further reduced. This figure also shows spatial convection of gas due to wind, with a higher shift in the $x$-direction since the wind velocity is assumed to be higher compared to that in the $y$-direction. Note also that for Case Ⅲ, the GHG density at the center of $\Omega$ is lower than that for Case Ⅱ, which is expected since Case Ⅱ has a higher absorption rate.

Simulation results presented in Figures 1–4 were obtained following the algorithm discussed in Theorem 5.2. The initial guess for the control policy $u(t)$ in Eq (7.1) was taken as zero. Following the gradient algorithm, the control policy was updated at every iteration resulting in corresponding monotonic reduction in cost as shown in Figure 5. Overall, the iterative process converged in 5 iterations as shown leading to the optimal solution presented above in Figures 1–4.

8.   Conclusions

This paper presents a conceptual framework for optimal removal of greenhouse gases from the atmosphere over a prescribed control time period. The optimal removal policy is derived based on a diffusion-advection model of greenhouse gas in the atmosphere. The system model as presented allows for a very general setup of source and absorber locations. For the simulation, we assume that all sources and absorbers of greenhouse gas are located near the ground level which is more natural. The numerical results illustrate the key concepts of temporal evolution of diffusion, transport, and optimal absorption. Absorption of $CO_2$, which constitutes 79% of greenhouse gases in the atmosphere, could be achieved both locally using carbon sequestration technologies as well as through extensive reforestation over the region. Progress has been made in recent years in carbon capture and storage with the installation of thirty operating facilities around the world and several more under development. Our results show a proof-of-concept of greenhouse gas removal that may be used for the determination of the optimal removal rate over a plan period. For $CO_2$ reduction by reforestation, one may also consider the rate of reforestation as the control variable which will require optimization of a coupled system of ecological model for vegetation growth and a diffusion-advection model for greenhouse gas, which is beyond the scope of this research. Overall, global warming is a large scale scientific and engineering problem, and combating it will require an equally large engineering endeavor and international collaboration.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

Publication of this article was funded in part by the Temple University Libraries Open Access Publishing Fund.

Conflict of interest

The authors declare that there are no conflicts of interest.