*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–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 *TD*, respectively, the equations that govern these perturbations can be written as follows:

$$
\mathcal{L}\_{lI}\frac{dT}{dt} = -\lambda T - \gamma (T - T\_D) + \Delta R\_{\text{CO}\_2} + (1 - \mathfrak{a}\_0)\Delta R\_A \tag{2}
$$

$$
\mathcal{L}\_D \frac{dT\_D}{dt} = \gamma (T - T\_D) \tag{3}
$$

Here, *CU* and *CD* are the effective heat capacities of the upper and lower models, respectively (note that *CU CD*); *λ* is a climate radiative feedback parameter; *γ* is a coupling strength parameter that describes the rate of heat loss by the upper layer; Δ*R*CO2 is the radiative forcing caused by global increase in the atmospheric CO2 concentration; Δ*RA* is the negative radiative forcing generated by the artificial aerosols at the top of the atmosphere; and *α*<sup>0</sup> is Earth's planetary albedo. We will assume that the temperature anomaly *T* is identified with the global mean surface temperature change *Tsfc* [53,54].

Despite its simplicity, this model imitates climate changes under external radiative forcing with reasonable accuracy [56–59]. We have chosen values of 7.34 W yr m−<sup>2</sup> K−1, 105.5 W yr m−<sup>2</sup> K−1, 1.13 W m−<sup>2</sup> K<sup>−</sup>1, and 0.7 W m−<sup>2</sup> K−<sup>1</sup> for parameters *CU*, *CD*, *λ*, and *γ*, 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 + bT\_D + \frac{\Delta R\_{\text{CO}\_2}}{\text{C}\_{U}} + \frac{(1 - \alpha\_0)\Delta R\_{\text{A}}}{\text{C}\_{U}}\tag{4}$$

$$\frac{dT\_D}{dt} = pT - pT\_D \tag{5}$$

where

$$a = \frac{\lambda + \gamma}{\mathbb{C}\_{II}}, \quad b = \frac{\gamma}{\mathbb{C}\_{II}}, \quad p = \frac{\gamma}{\mathbb{C}\_{D}}.\tag{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–71]:

$$
\Delta \mathcal{R}\_A = -\mathfrak{a}\_A \mathcal{Q}\_0 \tag{7}
$$

Here, *α<sup>A</sup>* is the instant albedo of the global aerosol layer; *Q*<sup>0</sup> is the global average incoming solar radiation on the top of the atmosphere defined as *Q*<sup>0</sup> = *I*0/4, where *I*<sup>0</sup> = 1368 W m<sup>2</sup> is a solar constant [72,73]. Thus, to estimate the radiative effect of stratospheric aerosol, we need to calculate the albedo *α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 *EA*. To determine *EA* from the known *αA*, the mass balance equation is used:

$$\frac{dM\_A}{dt} = E\_A - \frac{M\_A}{\tau\_A} \tag{8}$$

where *τ<sup>A</sup>* is the residence time of stratospheric aerosol particles; *MA* is the global mass of the stratospheric aerosols, which is linearly related to the albedo *α<sup>A</sup>* [69]:

$$
\mathfrak{a}\_A = M\_A(\beta\_A k\_A / Q\_0 \mathbb{S}\_\mathfrak{e}) \tag{9}
$$

where the coefficient *β<sup>A</sup>* = 24 W m−<sup>2</sup> [70,71]; *kA* = 7.6 m2g−<sup>1</sup> is the mass extinction coefficient [69]; *Se* is Earth's area determined as *Se* = 4*πR*<sup>2</sup> *<sup>e</sup>*, where *Re* = 6371 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 SO2, hydrogen sulfide H2S, 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 *ES* (in Tg S yr<sup>−</sup>1) and *MS* (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 *ES* ≈ *EA*/4 and *MS* ≈ *MA*/4. As the relationship between *MA* and *α<sup>A</sup>* is linear, the following predictive equation for *α<sup>A</sup>* can be derived from Equation (8):

$$\frac{d\mathfrak{a}\_A}{dt} = \chi^{-1} E\_S - \frac{\mathfrak{a}\_A}{\mathfrak{r}\_A} \tag{10}$$

where *<sup>χ</sup>* = *<sup>Q</sup>*0*Se*/(4*βAkA*) ≈ 2.39 × <sup>10</sup><sup>2</sup> Tg S.

Thus, solving the OCP, we can find the optimal control law *α*∗ *<sup>A</sup>*(*t*) and then calculate the optimal aerosol emission rate *E*∗ *<sup>S</sup>*(*t*) 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 CO2 can be parameterized as a function of CO2 only [75,76]: <sup>Δ</sup>*R*CO2 <sup>=</sup> *<sup>κ</sup>* <sup>×</sup> ln [*C*CO2 (*t*)/*C*(0) CO2 ], where *κ* (W m−2) is the empirical coefficient; *C*CO2 (*t*) is the CO2 concentration at time *<sup>t</sup>*; and *<sup>C</sup>*(0) CO2 is the reference CO2 concentration level. A typical

value for the parameter *κ* is near 5.35 W m2 [75,76]. In our model, we have taken the total global mean anthropogenic and natural radiative forcing Δ*RN* as prescribed by the different scenarios and approximated by a linear function of time:

$$
\Delta \mathcal{R}\_N = \eta t \tag{11}
$$

where *η* is the annual rate of forcing (see Table 1).

**Table 1.** Annual radiative forcing rate *η*.


#### *2.4. Optimal Control Problem Formulation*

We let [t0, tf] be a finite and fixed time interval. The OCP is defined as follows:

We find the control function *αA*(*t*) generating the corresponding temperature anomalies *T*(*t*) and *TD*(*t*) that minimizes the objective function:

$$\mathcal{J} = \frac{1}{2} \int\_{t\_0}^{t\_f} a\_A^2(t)dt\tag{12}$$

*subject to the dynamics* (4) *and* (5) *and given initial T*(*t*0) = 0 *and TD*(*t*0) = 0*, as well as final (terminal) T* - *tf* = *T<sup>f</sup>* conditions.

In this formulation, the terminal condition *T<sup>f</sup>* is interpreted as a target change in the global mean surface temperature at *t* = *tf* , and the performance index (Equation (12)) characterizes the aerosol consumption for SRM operations (recall that *α<sup>A</sup>* and *MA* 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 *tf* 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 *α<sup>A</sup>* ∈ [0, *U*], where *U* is the maximum value of technically feasible and affordable albedo *αA*.
