*2.5. Method for Solving the Optimal Control Problem*

Before proceeding further, we rewrite the model Equations (4) and (5) by replacing Δ*RCO*<sup>2</sup> with *ηt* (11) and Δ*RA* with −*αAQ*<sup>0</sup> (7):

$$\frac{dT}{dt} = -aT + bT\_D + ct - qa\_A \tag{13}$$

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

where *c* = η/*CU* and *q* = (1 − *α*0)*Q*0/*CU*.

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}a\_A^2 + \psi\_1(-aT + bT\_D + ct - qu\_A) + \psi\_2(pT - pT\_D) \tag{15}$$

where *ψ*<sup>1</sup> and *ψ*<sup>2</sup> 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\tag{16}$$

$$\frac{d\psi\_2}{dt} = -\frac{\partial H}{\partial T\_D} = -b\psi\_1 + p\psi\_2\tag{17}$$

The PMP states that the optimal control *α*∗ *<sup>A</sup>*(*t*) ∈ [0, *U*] is one that would maximize the Hamiltonian (Equation (13)) at each fixed time *t* ∈ *t*0, *tf* :

$$\mathfrak{a}\_A^\* = \underset{\mathfrak{a}\_{A} \in [0, lI]}{\text{arg } \max} H(\mathfrak{a}\_A) \tag{18}$$

Therefore, to find the optimal control *α*∗ *<sup>A</sup>*, we must maximize *H* with respect to *αA*, where the control belongs to the admissible control region *α<sup>A</sup>* ∈ [0, *U*]. The Hamiltonian maximization condition is as follows:

$$\frac{\partial H}{\partial \alpha\_A} = -\alpha\_A - q\psi\_1 = 0\tag{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*, *TD*, *ψ*1, and *ψ*<sup>2</sup> with given initial and terminal conditions. As the variable *TD* is not defined at *tf* , the following transversality condition for costate variable *ψ*<sup>2</sup> applies: *ψ*<sup>2</sup> - *tf* = 0 [48,49]. The analytic expressions derived for the control variable *α<sup>A</sup>* and temperature anomalies *T* and *TD* can be written as follows:

$$\boldsymbol{\omega}\_A(t) = -\mathbb{C}\_1 q \left[ \boldsymbol{v}\_{11} \boldsymbol{e}^{\lambda\_1 t} + \boldsymbol{e}^{(\lambda\_1 - \lambda\_2) t\_f} \boldsymbol{v}\_{21} \boldsymbol{e}^{\lambda\_2 t} \right] \tag{20}$$

$$T(t) = \mathbb{C}\_1 u\_1 \left( e^{\lambda\_1 t} - e^{\lambda\_2 t} \right) + \mathbb{C}\_3 e^{-\lambda\_1 t} + \mathbb{C}\_4 e^{-\lambda\_2 t} + \dots + w\_2 t + w\_1 \tag{21}$$

$$T\_D(t) = \mathbb{C}\_3 \frac{a - \lambda\_1}{b} e^{-\lambda\_1 t} + \mathbb{C}\_4 \frac{a - \lambda\_2}{b} e^{-\lambda\_2 t} + \tag{22}$$

$$+ \mathbb{C}\_1 \left[ \frac{a\_1 (a + \alpha\_1) - q^2 \upsilon\_{11}}{b} e^{\lambda\_1 t} - \frac{a\_2 (a + \alpha\_2) - q^2 \upsilon\_{21} e^{(\lambda\_1 - \lambda\_2) t\_f}}{b} e^{\lambda\_2 t} \right] +$$

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

where *C*1, *C*3, and *C*<sup>4</sup> are arbitrary integration constants (note that the integration constant *C*<sup>2</sup> = −*C*1*e* (*λ*1−*λ*2)*tf* ); *<sup>λ</sup>*<sup>1</sup> and *<sup>λ</sup>*<sup>2</sup> are the eigenvalues of the coefficient matrix of the adjoint system, Equations (14) and (15); *v*<sup>11</sup> and *v*<sup>21</sup> are the components of the corresponding eigenvectors.

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

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

$$w\_1 = \frac{c[(ap-pb) - p(a+p)]}{\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(t) \le \mathcal{C}\_T \qquad \forall t \in \left[t\_{0\prime}, t\_f\right] \tag{23}$$

where *CT* 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 *CT*, which is determined a priori. We should highlight that state constraints add a great deal of complexity to the OCP [77,78].
