Estimating the Impact of Artificially Injected Stratospheric Aerosols on the Global Mean Surface Temperature in the 21th Century

Abstract

In this paper, we apply the optimal control theory to obtain the analytic solutions of the two-component globally averaged energy balance model in order to estimate the influence of solar radiation management (SRM) operations on the global mean surface temperature in the 21st century. It is assumed that SRM is executed via injection of sulfur aerosols into the stratosphere to limit the global temperature increase in the year 2100 by 1.5 °C and keeping global temperature over the specified period (2020–2100) within 2 °C as required by the Paris climate agreement. The radiative forcing produced by the rise in the atmospheric concentrations of greenhouse gases is defined by the Representative Concentration Pathways and the 1pctCO₂ (1% per year CO₂ increase) scenario. The goal of SRM is formulated in terms of extremal problem, which entails finding a control function (the albedo of aerosol layer) that minimizes the amount of aerosols injected into the upper atmosphere to satisfy the Paris climate target. For each climate change scenario, the optimal albedo of the aerosol layer and the corresponding global mean surface temperature changes were obtained. In addition, the aerosol emission rates required to create an aerosol cloud with optimal optical properties were calculated.

1. Introduction

Climate change is among the most significant threats to human civilization in the 21st century and beyond [1]. The Paris Climate Accord proposed to hold average temperature increase “to well below 2 °C above pre-industrial levels” and to pursue efforts to keep warming “below 1.5 °C above pre-industrial levels” [2]. To reach these goals, eight countries have already presented long-term low-emission strategies, which aims to reduce greenhouse gas emissions; several countries are currently in the process of preparing such strategies [3]. Meanwhile, the World Meteorological Organization’s (WMO) “Statement of the State of the Global Climate in 2017” released in January 2018 said, “The global mean temperature in 2017 was approximately 1.1 °C above the pre-industrial era” [4]. There is high confidence that planetary warming will continue throughout the 21st century even if we immediately stopped emitting greenhouse gases into the atmosphere (e.g., References [5,6,7,8,9]). Some resent studies (e.g., References [10,11,12,13]) suggest that geoengineering technologies can serve as a supplementary measure to stabilize climate as “in the absence of external cooling influence” [14], it is hard to achieve the Paris Agreement climate goals.

Solar radiation management (SRM) by injection of sulfur aerosols into the stratosphere [15,16] is one of the most feasible and promising solutions for inducing negative radiative forcing (RF) from aerosols in order to at least partially compensate the positive RF from atmospheric greenhouse gases. The current state of understanding of climate engineering technologies, including SRM, has been discussed in References [17,18,19,20,21,22,23,24,25,26]. Over the years, climate models have played a key role in exploring geoengineering techniques and predicting and quantifying their potential effects on Earth’s climate (e.g., References [27,28,29,30,31,32,33,34,35,36]). Due to the uncertainties inherent in climate models that could not be sufficiently reduced over the last decade [37], the resulting range of possible outcomes of hypothetical geoengineering efforts remains quite vague. To handle the climate response uncertainties, some studies (e.g., References [38,39,40,41,42,43,44,45]) have suggested modeling the Earth’s climate as a control system with feedbacks, which allows planning scenarios for geoengineering using the so-called “design model”. This formulation makes it possible to design the control law and calculate the amount of SRM forcing as a function of time needed to offset the rise in global mean surface temperature due to human-caused positive RF. Meanwhile, exploring Earth’s global climate as controlled dynamical system, we can approach geoengineering from the perspective of optimal control theory [46,47,48,49]. Within the optimal control framework, the goal of geoengineering can be formulated in terms of extremal problem, which involves finding control functions and the corresponding climate system trajectory that minimize or maximize a certain objective functional (also referred to as performance measure or index) subject to various constraints (e.g., References [50,51]. If $x$ is the state vector of climate system and $u$ is the vector of control variables, then the abstract extremal problem can be formulated as follows:

$$\mathcal{J}\left(x,u\right)\to extr, \mathcal{F}\left(x, u\right)=0, \left(x, u\right)\in \mathcal{M}\subset \mathcal{X}\times \mathcal{U}$$

(1)

The statement of this problem includes a set $\mathcal{X}\times \mathcal{U}$ on which the (real) functional $\mathcal{J}\left(x,u\right)$ is defined and constraints imposed on state and control variables given by the model of control object $\mathcal{F}\left(x, u\right)=0$ (dynamic constraints) and by the subset $\mathcal{M}$ in $\mathcal{X}\times \mathcal{U}$. The solution to the extremal problem (Equation (1)) is the optimal process $\left(x^{\ast },u^{\ast }\right)$. Thus, by solving the optimal control problem (OCP), we can obtain the mathematically rigorous control law and the corresponding system’s trajectory that are relevant for the specified performance measure $\mathcal{J}\left(x,u\right)$.

This paper deals with a simple mathematical model for controlling the global mean surface temperature $T_{sfc}$ in the 21st century by the injection of sulfur aerosols into the stratosphere to limit the global temperature increase in the year 2100 by 1.5 °C above pre-industrial level and keeping global temperature over the period of 2020–2100 within 2 °C as required by the Paris climate agreement. The objective is to minimize resources (the total mass of aerosols) required to achieve the desired final state of the climate system. In the model, the positive RF produced by the rise in the atmospheric concentrations of greenhouse gases is specified in accordance with the Representative Concentration Pathways [52] and the 1pctCO₂ (1% per year CO₂ increase) scenario.

The mathematical statement of OCP is the collection of the following key elements: objective function defined to judge the effectiveness of control process, mathematical model of the controlled object, equality and inequality constraints to be satisfied by state and control variables, and boundary and initial conditions (if any) for state variables. To imitate the behavior of the climate system, we applied a two-component energy balance model [53,54,55] in which the global mean surface temperature anomaly (perturbation) represents the variable that interests us the most, and the albedo of the global aerosol layer is designated as the control variable. We derived analytical expressions for both the optimal albedo of the global aerosol layer and the corresponding change in the global mean surface temperature.

The results of illustrative calculations are presented for the period 2020–2100. For each climate change scenario, the optimal albedo of the aerosol layer—and therefore the aerosol emission rates—as well as the associated global mean surface temperature changes were found.

We need to emphasize that the main reason for using such a model is that similar two-layer models have been considered and analyzed in a number of papers considering the response to forced climate change. For example, Geoffroy et al. [56,57] obtained and discussed the general analytical solutions of the two-layer model for different hypothetical climate forcing scenarios and suggested the approach of calibrating the model parameters to imitate the time response of coupled general circulation models (CGCMs) from CMIP5 to radiative forcing. Gregory et al. [58] analyzed the two-layer model and discussed the transient climate response, the global mean surface air temperature change under two scenarios: one with a step forcing (the abrupt 4xCO₂ experiment) and one with the 1pctCO₂ scenario. Despite the fact that the two-layer model is one of the simplest tools to mimic climate dynamics under external radiative forcing, it was able to simulate the evolution of average global surface temperature over time in response to both abrupt and time-dependent forcing with reasonable accuracy (e.g., References [56,59]).

In the two-layer model, climate control is carried out via changing Earth’s planetary albedo by injection of sulfur aerosols into the stratosphere (“albedo modification”). Sulfur aerosols increase the amount of sunlight that is scattered back to space, thereby reducing the amount of sunlight absorbed by Earth. Inherently, the planetary albedo is an average of the local albedo, averaged over the entire globe. The local albedo, in turn, is a highly variable dimensionless parameter that depends on a number of local factors, such as the composition of the atmosphere and in particular the presence of aerosols, the cloud amount and properties [60], the sea ice cover [61], the land use [62,63,64], the snow cover [65], etc. A typical value of Earth’s planetary albedo is about 0.3 [66].

As change in the albedo of our planet is a powerful driver of climate (indeed, a 1% change in the Earth’s planetary albedo generates the radiative effect of 3.42 Wm⁻², which is commensurate with radiative forcing due to a doubling of CO₂ concentrations in the atmosphere), scientists have proposed “albedo modification” as a powerful tool to deal with global warming (e.g., References [16,20,23,29,31,33,67]).

2. Materials and Methods

2.1. The Model of Control Object

The control object is Earth’s climate system. To simulate the climate system dynamics under the influence of external radiative forcing, we have applied the mathematical model consisting of two subsystems: One is the upper layer subsystem, which combines the atmosphere, the land surface, and the upper ocean; the other is the lower layer subsystem, which represents the deep ocean [53,54,55]. The state of each subsystem is characterized by the corresponding temperature perturbation (anomaly) with respect to initial climate “equilibrium” state. Denoting temperature anomalies for upper and lower subsystems by $T$ and $T_D$, respectively, the equations that govern these perturbations can be written as follows:

$$C_U\frac{dT}{dt}=-\lambda T-\gamma \left(T-T_D\right)+\mathsf{\Delta }{R}_{{\mathrm{CO}}_2}+\left(1-{\alpha }_0\right)\mathsf{\Delta }{R}_A$$

(2)

$$C_D\frac{d{T}_D}{dt}=\gamma \left(T-T_D\right)$$

(3)

Here, $C_U$ and $C_D$ are the effective heat capacities of the upper and lower models, respectively (note that $C_U\ll C_D$); $\lambda$ is a climate radiative feedback parameter; $\gamma$ is a coupling strength parameter that describes the rate of heat loss by the upper layer; $\mathsf{\Delta }{R}_{{\mathrm{CO}}_2}$ is the radiative forcing caused by global increase in the atmospheric CO₂ concentration; $\mathsf{\Delta }{R}_A$ is the negative radiative forcing generated by the artificial aerosols at the top of the atmosphere; and ${\alpha }_0$ is Earth’s planetary albedo. We will assume that the temperature anomaly $T$ is identified with the global mean surface temperature change $T_{sfc}$ [53,54].

Despite its simplicity, this model imitates climate changes under external radiative forcing with reasonable accuracy [56,57,58,59]. We have chosen values of $7.34 \mathrm{W} \mathrm{y}\mathrm{r} {\mathrm{m}}^{-2} {\mathrm{K}}^{-1}$, $105.5 \mathrm{W} \mathrm{y}\mathrm{r} {\mathrm{m}}^{-2} {\mathrm{K}}^{-1}$, $1.13 \mathrm{W} {\mathrm{m}}^{-2} {\mathrm{K}}^{-1}$, and $0.7 \mathrm{W} {\mathrm{m}}^{-2} {\mathrm{K}}^{-1}$ for parameters $C_U, C_D, \lambda$, and $\gamma$, respectively. These values are taken in accordance with values consistent with the CMIP5 multimodel mean under climate change derived in Reference [56].

For convenience sake, we have rewritten the model Equations (2) and (3) as follows:

$$\frac{dT}{dt}=-aT+b{T}_D+\frac{\mathsf{\Delta }{R}_{{\mathrm{CO}}_2}}{C_U}+\frac{\left(1-{\alpha }_0\right)\mathsf{\Delta }{R}_{\mathrm{A}}}{C_U}$$

(4)

$$\frac{d{T}_D}{dt}=pT-p{T}_D$$

(5)

where

$$a=\frac{\lambda +\gamma }{C_U}, b=\frac{\gamma }{C_U}, p=\frac{\gamma }{C_D}.$$

(6)

Ordinary Differential Equations (4) and (5) represent a mathematical model of the object to be controlled.

2.2. Parameterization of the Aerosols’ Radiative Effect

The climate control is assumed to be executed through the injection of nonabsorptive sulfate aerosols into the stratosphere. Injected aerosol particles scatter shortwave solar radiation back to the outer space and consequently change the radiative balance of our planet, increasing Earth’s planetary albedo and therefore causing the negative RF at the top of the atmosphere [68,69,70,71]:

$$\mathsf{\Delta }{R}_A=-{\alpha }_A{Q}_0$$

(7)

Here, ${\alpha }_A$ is the instant albedo of the global aerosol layer; $Q_0$ is the global average incoming solar radiation on the top of the atmosphere defined as $Q_0=I_0/4$, where $I_0=1368 {\mathrm{W}\mathrm{m}}^2$ is a solar constant [72,73]. Thus, to estimate the radiative effect of stratospheric aerosol, we need to calculate the albedo ${\alpha }_A$, which is considered as the control variable. However, in reality, we have the ability to manipulate the emission rate of aerosols injected into the stratosphere $E_A$. To determine $E_A$ from the known ${\alpha }_A$, the mass balance equation is used:

$$\frac{d{M}_A}{dt}=E_A-\frac{M_A}{{\tau }_A}$$

(8)

where ${\tau }_A$ is the residence time of stratospheric aerosol particles; $M_A$ is the global mass of the stratospheric aerosols, which is linearly related to the albedo ${\alpha }_A$ [69]:

$${\alpha }_A=M_A\left({\beta }_A{k}_A/Q_0{S}_e\right)$$

(9)

where the coefficient ${\beta }_A=24$ W m⁻² [70,71]; $k_A=7.6{\mathrm{m}}^2{\mathrm{g}}^{-1}$ is the mass extinction coefficient [69]; $S_e$ is Earth’s area determined as $S_e=4\pi R_e^2$, where $R_e=6371 \mathrm{km}$ is Earth’s radius.

In geoengineering, sulfate aerosol particles are not directly injected into the stratosphere but can be formed from gaseous precursors, such as sulfur dioxide SO₂, hydrogen sulfide H₂S, carbonyl sulfide OCS, or dimethyl sulfide (DMS), which then convert into aerosols. We will express the emission rate of aerosol precursors as well as the mass of sulfate aerosols in units of sulfur, denoting them by $E_S$ (in Tg S yr⁻¹) and $M_S$ (in Tg S), respectively. Assuming that 1 Tg of sulfur injected into the stratosphere forms approximately 4 Tg of aerosol particles [74], we obtain that $E_S\approx E_A/4$ and $M_S\approx M_A/4$. As the relationship between $M_A$ and ${\alpha }_A$ is linear, the following predictive equation for ${\alpha }_A$ can be derived from Equation (8):

$$\frac{d{\alpha }_A}{dt}={\chi }^{-1}{E}_S-\frac{{\alpha }_A}{{\tau }_A}$$

(10)

where $\chi =Q_0{S}_e/\left(4{\beta }_A{k}_A\right)\approx 2.39\times {10}^2 \mathrm{Tg} \mathrm{S}$.

Thus, solving the OCP, we can find the optimal control law ${\alpha }_A^{\ast }\left(t\right)$ and then calculate the optimal aerosol emission rate $E_S^{\ast }\left(t\right)$ using Equation (10).

2.3. Parameterization of the Anthropogenic Radiative Forcing

In energy balance models, simple empirical expressions are generally used to calculate radiative forcing due to the increase in atmospheric greenhouse gases. For example, the radiative forcing caused by a perturbation of the atmospheric burden of CO₂ can be parameterized as a function of CO₂ only [75,76]: $\mathsf{\Delta }{R}_{{\mathrm{CO}}_2}=\kappa \times \ln [C_{{\mathrm{CO}}_2}\left(t\right)/C_{{\mathrm{CO}}_2}^{\left(0\right)}]$, where $\kappa$ (W m⁻²) is the empirical coefficient; $C_{{\mathrm{CO}}_2}\left(t\right)$ is the CO₂ concentration at time t; and $C_{{\mathrm{CO}}_2}^{\left(0\right)}$ is the reference CO₂ concentration level. A typical value for the parameter $\kappa$ is near 5.35 W m² [75,76]. In our model, we have taken the total global mean anthropogenic and natural radiative forcing $\mathsf{\Delta }{R}_N$ as prescribed by the different scenarios and approximated by a linear function of time:

$$\Delta R_N=\eta t$$

(11)

where $\eta$ is the annual rate of forcing (see

Table 1. Annual radiative forcing rate η.

Scenario	RCP8.5	1pctCO2	RCP6	RCP4.5
η (W m−2 yr−1)	$7.14\times$ 10−2	$5.29\times$ 10−2	$3.814\times$ 10−2	$2.17\times$ 10−2

).

2.4. Optimal Control Problem Formulation

We let [${\mathrm{t}}_0,{\mathrm{t}}_{\mathrm{f}}$] be a finite and fixed time interval. The OCP is defined as follows:

We find the control function ${\alpha }_A\left(t\right)$ generating the corresponding temperature anomalies $T\left(t\right)$ and $T_D\left(t\right)$ that minimizes the objective function:

$$\mathcal{J}=\frac{1}{2}\underset{t_0}{\overset{t_f}{{\displaystyle \int }}}{\alpha }_A^2\left(t\right)dt$$

(12)

subject to the dynamics (4) and (5) and given initial $T\left(t_0\right)=0$ and $T_D\left(t_0\right)=0$, as well as final (terminal) $T\left(t_f\right)=T^f$ conditions.

In this formulation, the terminal condition $T^f$ is interpreted as a target change in the global mean surface temperature at $t=t_f$, and the performance index (Equation (12)) characterizes the aerosol consumption for SRM operations (recall that ${\alpha }_A$ and $M_A$ are linearly dependent functions). Thus, we wish to minimize the mass of aerosols required to reach the target surface temperature change at the final time. The global mean deep ocean temperature anomaly at the final time $t_f$ is not defined because changes in the global mean surface temperature are of primary concern, while changes in the deep ocean temperature are only of secondary concern. The total amount of aerosols annually emitted to the stratosphere can be limited by the available technical equipment. In this case, the minimization problem (Equation (12)) should be considered within the framework of control-constrained OCP. The set of admissible controls is given formally by ${\alpha }_A\in \left[0, U\right]$, where $U$ is the maximum value of technically feasible and affordable albedo ${\alpha }_A$.

2.5. Method for Solving the Optimal Control Problem

Before proceeding further, we rewrite the model Equations (4) and (5) by replacing $\Delta R_{C{O}_2}$ with $\eta t$ (11) and $\Delta R_A$ with $-{\alpha }_A{Q}_0$ (7):

$$\frac{dT}{dt}=-aT+b{T}_D+ct-q{\alpha }_A$$

(13)

$$\frac{d{T}_D}{dt}=pT-p{T}_D$$

(14)

where $c=\mathsf{\eta }/C_U$ and $q=\left(1-{\alpha }_0\right)Q_0/C_U$.

We solve the formulated OCP using the Pontryagin’s maximum principle (PMP) [47]. The Hamiltonian function for the problem (12) is defined as follows:

$$H=-\frac{1}{2}{\alpha }_A^2+{\psi }_1\left(-aT+b{T}_D+ct-q{\alpha }_A\right)+{\psi }_2\left(pT-p{T}_D\right)$$

(15)

where ${\psi }_1$ and ${\psi }_2$ are time-varying Lagrange multipliers, also known as costate or adjoint variables, which satisfy the adjoint system:

$$\frac{d{\psi }_1}{dt}=-\frac{ \partial H}{\partial T}=a{\psi }_1-p{\psi }_2$$

(16)

$$\frac{d{\psi }_2}{dt}=-\frac{\partial H}{\partial T_D}=-b{\psi }_1+p{\psi }_2$$

(17)

The PMP states that the optimal control ${\alpha }_A^{\ast }\left(t\right)\in \left[0, U\right]$ is one that would maximize the Hamiltonian (Equation (13)) at each fixed time $t\in \left[t_0, t_f\right]$:

$${\alpha }_A^{\ast }=\underset{{\alpha }_A\in \left[0,U\right]}{\arg \max }H\left({\alpha }_A\right)$$

(18)

Therefore, to find the optimal control ${\alpha }_A^{\ast }$, we must maximize $H$ with respect to ${\alpha }_A$, where the control belongs to the admissible control region ${\alpha }_A\in \left[0,U\right]$. The Hamiltonian maximization condition is as follows:

$$\frac{\partial H}{\partial {\alpha }_A}=-{\alpha }_A-q{\psi }_1=0$$

(19)

Thus, to find the optimal control and the corresponding climate system’s trajectory, we need to solve the set of four ordinary Differential Equations (13), (14), (16), and (17) in four unknowns $T$, $T_D$, ${\psi }_1$, and ${\psi }_2$ with given initial and terminal conditions. As the variable $T_D$ is not defined at $t_f$, the following transversality condition for costate variable ${\psi }_2$ applies: ${\psi }_2\left(t_f\right)=0$ [48,49]. The analytic expressions derived for the control variable ${\alpha }_A$ and temperature anomalies $T$ and $T_D$ can be written as follows:

$${\alpha }_A\left(t\right)=-C_{1}q\left[v_{11}{e}^{{\lambda }_{1}t}+e^{\left({\lambda }_1-{\lambda }_2\right)t_f}{v}_{21}{e}^{{\lambda }_{2}t}\right]$$

(20)

$$T\left(t\right)=C_1{\alpha }_1\left(e^{{\lambda }_{1}t}-e^{{\lambda }_{2}t}\right)+C_3{e}^{-{\lambda }_{1}t}+C_4{e}^{-{\lambda }_{2}t}++w_{2}t+w_1$$

(21)

$$T_D\left(t\right)=C_3\frac{a-{\lambda }_1}{b}{e}^{-{\lambda }_{1}t}+C_4\frac{a-{\lambda }_2}{b}{e}^{-{\lambda }_{2}t}+$$

(22)

$$+C_1\left[\frac{{\alpha }_1\left(a+{\alpha }_1\right)-q^2{v}_{11}}{b}{e}^{{\lambda }_{1}t}-\frac{{\alpha }_2\left(a+{\alpha }_2\right)-q^2{v}_{21}{e}^{\left({\lambda }_1-{\lambda }_2\right)t_f}}{b}{e}^{{\lambda }_{2}t}\right]+$$

$$+\frac{a{w}_2-c}{b}t+\frac{a{w}_1+w_2}{b}$$

where $C_1$, $C_3$, and $C_4$ are arbitrary integration constants (note that the integration constant $C_2=-C_1{e}^{\left({\lambda }_1-{\lambda }_2\right)t_f}$); ${\lambda }_1$ and ${\lambda }_2$ are the eigenvalues of the coefficient matrix of the adjoint system, Equations (14) and (15); $v_{11}$ and $v_{21}$ are the components of the corresponding eigenvectors.

$${\alpha }_1=\frac{q^2{v}_{11}\left({\lambda }_1+p\right)}{{\lambda }_1^2+{\lambda }_1\left(a+p\right)+\left(ap-pb\right)}$$

$${\alpha }_2=\frac{q^2{v}_{21}\left({\lambda }_2+p\right)e^{\left({\lambda }_1-{\lambda }_2\right)t_f}}{{\lambda }_1^2+{\lambda }_1\left(a+p\right)+\left(ap-pb\right)}$$

$$w_1=\frac{c\left[\left(ap-pb\right)-p\left(a+p\right)\right]}{{\left(ap-pb\right)}^2}$$

$$w_2=\frac{pc}{ap-pb}$$

The constants of integration are determined by applying the boundary conditions.

If we consider climate engineering as a state-constrained OCP with constraints on the state variables, then additional necessary conditions for optimality, known as the complementary slackness conditions, should be specified [77]. In this study, we express the OCP with the following state constraint:

$$T\left(t\right)\le C_T \forall t\in \left[t_0, t_f\right]$$

(23)

where $C_T$ is the threshold parameter whose value should be set. The meaning of the condition (Equation (23)) known as a path constraint is that throughout the geoengineering project, the global mean surface temperature change should not exceed a certain value $C_T$, which is determined a priori. We should highlight that state constraints add a great deal of complexity to the OCP [77,78].

3. Results and Discussion

In the calculations, we took calendar years 2020 and 2100 as the initial $t_0$ and the final (terminal) $t_f$ time, respectively, which meant that we were examining the climate control problem on the finite time interval 2020–2100. To formulate the boundary conditions and impose a constraint on change in the global mean surface temperature $T_{sfc}$, we assumed the following:

-

The temperature anomalies $T$ and $T_D$ were calculated relative to 2020, i.e., the boundary conditions for $T$ and $T_D$ at $t=t_0$ were $T_{2020}=0$ and $T_{D, 2020}=0$, respectively, where the numerical subscript referred to the year 2020.

-

By 2020, $T_{sfc}$ would exceed the pre-industrial level by 1.1 °C, i.e., $\Delta T_{2020}=1.1 ℃$.

-

By 2100, $T_{sfc}$ would exceed the pre-industrial level by 1.5 °C, i.e., $\Delta T_{2100}=1.5 ℃$.

-

For the 2020 to 2100 period, the rise in $T_{sfc}$ should not exceed 2 °C above the pre-industrial level.

Then, the permissible increase in the temperature anomaly $T_{2100}$ by year 2100 relative to 2020 would be $T_{2100}=\Delta T_{2100}-\Delta T_{2020}=0.4 ℃$. This value represents the boundary condition for T at $t=t_f$. The threshold parameter, which defines a path constraint (Equation (23)), is $C_T=2 ℃-\Delta T_{2020}=0.9 ℃$.

Changes in both global mean surface temperature and deep ocean temperature calculated for different climate change scenarios in the absence of climate engineering interventions are illustrated in Figure 1. The corresponding temperature changes in the year 2100 are shown in

Table 2. Calculated temperature changes T and T_D from 2020 to 2100 (changes relative to the pre-industrial level are shown in brackets).

Scenario	RCP8.5	1pctCO2	RCP6.0	RCP4.5
T (K)	3.16 (4.26)	2.34 (3.44)	1.70 (2.80)	0.96 (2.06)

. According to Reference [79], without additional measures to reduce GHG emissions (RCP8.5 scenario), increases in global mean surface temperatures are expected to be between 3.7 and 4.8 °C by the year 2100 versus pre-industrial levels (this range is based on median climate response). As seen in Table 2, by year 2100, the model outlined here projects globally averaged surface temperature increases of 4.26, 3.44, and 2.80 °C for the RCP8.5, 1pctCO₂, and RCP6.0 scenarios, respectively (relative to pre-industrial period). Thus, geoengineering can be regarded as one of supplementary measures needed to achieve the climate targets of the Paris Agreement.

We considered results of calculations for the RCP8.5 (the worst-case) scenario in more detail. Figure 2 shows (a) the optimal albedo of the global stratospheric aerosol layer, (b) the corresponding surface temperature anomaly, (c) the mass of the global aerosol layer, and (d) the optimal emission rate of aerosol particles calculated for RCP8.5 pathway with and without constraint on the global mean surface temperature increase. In the absence of state constraint, the optimal albedo ${\alpha }_A^{\ast }$ and, accordingly, the optimal emission rate of aerosol particles $E_S^{\ast }$ would increase exponentially. This optimal aerosols emission rate ensures that the target temperature anomaly $T_{2100}=0.4 ℃$ is satisfied. However, within the given time interval 2020–2100, a temperature rise would exceed the set point $C_T$, i.e., $T\left(t\right)>C_T$ (the “overshooting” phenomenon [80]). The maximum increases in global mean surface temperature for different climate change scenarios are presented in

Table 3. Maximum global mean surface temperature anomaly T calculated without state constraint.

Scenario	RCP8.5	1pctCO2	RCP6	RCP4.5
T (K)	2.48	1.84	1.34	0.78

. The use of the constraint (Equation (21)) allows us to avoid overshoot; however, compared to the unconstrained case, keeping the increase in global mean surface temperature below the target constrained level $C_T$ would require additional amount of aerosols (see

Table 4. The mass of aerosols $M_{S, tot}$ (Tg S) injected into the stratosphere from 2020 to 2100.

Scenario	RCP8.5	1pctCO2	RCP6	RCP4.5
$M_{S,tot}$ without state constraint	36.5	25.6	17.0	7.7
$M_{S,tot}$ with state constraint	73.6	44.46	23.3	-

). For example, for the RCP8.5 scenario, the total mass of aerosol particles injected in the stratosphere from the year 2020 to 2100 is about 73.6 Tg S, which is about 2 times larger than $M_{S, tot}^{\ast }$, calculated by solving an unconstrained OCP.

Results obtained for 1pctCO₂, RCP6.0, and RCP4.5 scenarios are represented in Figure 3, Figure 4 and Figure 5, respectively. These figures show that the overshooting phenomenon is also observed for the 1pctCO₂ and RCP6.0 scenarios. The only exception is the RCP4.5 scenario.

If the optimal control problem is considered with control variable constraint ${\mathsf{\alpha }}_{\mathrm{A}}\left(\mathrm{t}\right)\le U$, then the target value for the temperature anomaly in the year 2100 may not necessarily be achieved (this depends on the value of the constraint U and the scenario in question). Assuming, for example, that $U=0.02$, then the corresponding instant mass of aerosols is estimated to be 4.8 Tg S. In such a case, for the RCP8.5 scenario, the calculated temperature anomaly in the year 2100 relative to 2020 would exceed the target value by 0.3 °C, which is equivalent to exceeding the pre-industrial level by 1.8 °C. It needs to be recalled that constraint on the control variable is associated with a possible limitation on resources required to implement the project, namely, the amount of aerosols available to the project executors.

We emphasize that the results of calculations discussed above are for illustration purposes only. The primary outcome presented in this paper is the optimal control-based approach that can be used to design projects targeting purposeful manipulation of climate and weather.

4. Concluding Remarks

The use of fine aerosol particles, artificially injected into the stratosphere, is considered to be one of the most effective and feasible measures to counter global warming in the 21st century and beyond. Computer simulation using mathematical climate models of various degrees of sophistication and complexity is the most popular and reliable technique for exploring and estimating the effectiveness of stratospheric aerosol climate engineering and climate and weather manipulation. Numerical simulation of climate engineering requires the design of fairly realistic scenarios for aerosol injections. This paper introduced the optimal-control-based method for designing climate engineering scenarios. Considering Earth’s climate as controlled dynamical system, we proposed to approach geoengineering from the standpoint of the optimal control theory, thereby formulating the goal of geoengineering projects in terms of extremal problem. The capability to apply this technique was illustrated using the two-layer energy balance model in which the global mean surface temperature anomaly and the deep ocean temperature perturbation were the state variables, and the emission rate of aerosol precursors was the control variable. Solutions to the unconstrained as well as state and control constrained problems were obtained on the basis of classical Pontryagin’s maximum principle. The proposed method will provide additional useful insights for the development of optimal climate manipulation strategies to counter global warming in the 21st century.