An Optimal Control Perspective on Weather and Climate Modiﬁcation

: Intentionally altering natural atmospheric processes using various techniques and technologies for changing weather patterns is one of the appropriate human responses to climate change and can be considered a rather drastic adaptation measure. A fundamental understanding of the human ability to modify weather conditions requires collaborative research in various scientiﬁc ﬁelds, including, but not limited to, atmospheric sciences and different branches of mathematics. This article being theoretical and methodological in nature, generalizes and, to some extent, sum-marizes our previous and current research in the ﬁeld of climate and weather modiﬁcation and control. By analyzing the deliberate change in weather and climate from an optimal control and dynamical systems perspective, we get the ability to consider the modiﬁcation of natural atmospheric processes as a dynamic optimization problem with an emphasis on the optimal control problem. Within this conceptual and uniﬁed theoretical framework for developing and synthesizing an optimal control for natural weather phenomena, the atmospheric process in question represents a closed-loop dynamical system described by an appropriate mathematical model or, in other words, by a set of differential equations. In this context, the human control actions can be described by variations of the model parameters selected on the basis of sensitivity analysis as control variables. Application of the proposed approach to the problem of weather and climate modiﬁcation is illustrated using a low-order conceptual model of the Earth’s climate system. For the sake of convenient interpretation, we provide some weather and climate basics, as well as we give a brief glance at control theory and sensitivity analysis of dynamical systems.


Introduction
Weather modification and geoengineering represent the deliberate alteration of atmospheric and climate conditions, locally or globally, by humans using the available assets and resources based on the existing theoretical understanding of weather and climate processes (see [1][2][3][4][5][6] and references herein). Over the years, people have sought to modify the environment, including the atmosphere, in an attempt to adapt to its ever-changing conditions. However, the scientific-based stage of modification of the environment, first of all meteorological processes, only began in the middle of the 20th century when scientists at the General Electric Research Laboratories suggested using dry ice to disperse clouds [7]. Since the early 1960s, dozens of projects have been implemented around the world aimed at modifying various atmospheric phenomena including different types of clouds, precipitation, lighting, hail, tornadoes, thunderstorms, hurricanes, and tropical cyclones. In the current century, the study and assessment of the human ability to modify environmental conditions has taken on a new sound in connection with climate change, which is caused by anthropogenic activities posing a serious threat to all of humanity [8]. To mitigate the impact of climate change on nature and society, scientists have proposed the use of tools and methods of so-called geoengineering (e.g., [9][10][11][12][13]). However, we should make a difference between weather modification and geoengineering. Weather characterizes the current atmospheric conditions in a certain area or point at a particular time, or, in other words, the state of the atmosphere at some period of time is described by such atmospheric variables as temperature, pressure, humidity, wind speed, and direction, etc. Meanwhile, climate is the ensemble of states traversed by the Earth's climate system over a sufficiently long temporal interval (~30 years). In this regard, the atmosphere, one of the five major components of the Earth's climate system, is the most dynamic, unstable, and fastest-responding element of the climate system. The goal of weather modification is to change weather conditions over some limited geographical area or some geographical point. In other words, weather modification technologies are used to affect processes only in the atmosphere, which, as was stated above, is the most variable element of the climate system. In turn, geoengineering, being a planetary-scale process, is aimed at neutralizing anthropogenic radiative forcing and thereby reducing or even preventing human-caused warming of the Earth. Let us note that anthropogenic radiative forcing is produced mainly by greenhouse gases (carbon dioxide CO 2 , methane CH 4 , nitrous oxide N 2 O, and fluorinated gases) entering the atmosphere via burning fossil fuels and industrial and agricultural activities. However, both weather modification and geoengineering have one thing in common: these two procedures are goal-oriented processes implemented by means of external human-produced effects (interventions) to achieve specific objectives under various constraints. Thus, deliberate modification of weather and climate is, in substance, a dynamic optimization problem.
By viewing the atmosphere and/or the climate as a controllable dynamical system, we can approach weather modification and geoengineering from the perspective of optimal control theory. Within this conceptual and unified theoretical framework, the purpose of weather modification or geoengineering is formulated in terms of an extremal problem, which involves finding control functions and the corresponding climate (atmospheric) system trajectory that minimize or maximize a given objective function (also referred to as performance measure or index) subject to various constraints. In this instance, the atmospheric (climate) process in question is considered a closed-loop dynamical system, the evolution of which is described by the appropriate mathematical model, commonly represented by a set of differential equations. In this context, the human control actions can be described by variations in the model parameters selected on the basis of sensitivity analysis as control variables. This multi-disciplinary approach for planning and implementation of weather modification and geoengineering projects is known as geophysical cybernetics, the theoretical foundations of which were laid by the authors of this paper in the 1980s and '90s [14].
This article is theoretical and methodological in nature aiming at generalizing and, to some extent, summarizing our previous and current research related to weather and climate modification and control. In the present paper, we first consider the deliberate change in weather and climate from an optimal control and dynamical systems perspective and, second, illustrate the application of this approach using a low-order conceptual model of the Earth's climate system. For the sake of convenient interpretation, we provide some weather and climate basics, as well as giving a brief glance at control theory and sensitivity analysis of dynamical systems relevant to weather and climate control.

Atmosphere and Climate as Dynamical Systems
The Earth's climate system (ECS) and its components, the atmosphere, hydrosphere, land surface, cryosphere, and biosphere, as control objects have some unique features that are previously discussed in detail in [4]. The ECS is an open large-scale object affected by different external forcing mechanisms including the anthropogenic one, but at the same time, the impact of ECS on the external environment is minor. Each of the ECS subsystems has specific dynamical, physical, and chemical properties. For example, the atmosphere is the most variable and unstable component of the ECS having turbulent nature and showing wave-like oscillations over a wide-range time-space spectrum. ECS simulation, and moreover, its control, is an extremely difficult problem due to a number of objective reasons. First, ECS is an enormously complex natural object with many feedbacks and cycles that are difficult to describe mathematically. Second, the ECS and its elements are not fully identified as control objects; the corresponding mathematical models are not "perfect", and the degree of their reliability and validity is not always sufficiently high. Nevertheless, over the last few decades, models of varying levels of complexity, from conceptual to realistic, have become very powerful tools used for numerical weather prediction and climate simulations [15][16][17].
Typically, deterministic mathematical models are used to numerically forecast weather and study climate trends, while random mathematical models are applied to analyze the variability of the climate system driven by stochastic forcing. For the purposes of weather and climate modification and control, deterministic models are the primary focus. In general terms, atmospheric (climate) models are systems of non-linear differential equations in partial derivatives that are the mathematical statements of basic physical, chemical, and biological laws. Such models also include a variety of empirical and semiempirical relationships and equations that are based on observations and experiences rather than theories, and contain a large number of parameters that are considerably diverse in their physical meaning. The equations describing the atmospheric and climate dynamics are very complicated and, therefore, in most cases can be solved only numerically. To find an approximate solution, various numerical techniques are usually used, for example, the Galerkin projection method or the finite difference approach. As a result, atmospheric and climate models are of finite space and time resolutions. Because of the limited time-spatial resolutions of atmospheric and climate models, some physical, chemical, and biological processes cannot be adequately resolved by the model space-time grid, and therefore can be only parameterized leading to the significant increase in the model parameters. Some model parameters are not well defined, generating parametric uncertainty in atmospheric and climate models affecting the output results. Assessment of parametric uncertainty is a very important component of model building and its quality assurance. Sensitivity analysis in dynamical systems is among the most helpful tool for estimating the influence of model parameters and their variations on the results of numerical simulations. It is necessary to emphasize that both natural external forcing and purposeful and unintentional human-caused forcing on the atmospheric (climate) system are described in atmospheric (climate) models by means of variations in certain parameters. In this context, sensitivity analysis serves as a tool for selecting parameters to be used as control variables. An abstract dynamical system is a pair (X, S t ), where X is the system's phase space and S t : X → X is a family of smooth evolution functions parameterized by a real variable t ∈ T, the time. The set γ s = {x(t) : t ∈ T} is called system's trajectory (orbit), where x(t) is continuous function with values in X such that S τ x(t) = x(t + τ) for all t ∈ T and τ ∈ T + . In atmospheric and climate studies, semi-dynamical systems are of primary interest, for which a family {S t : t ≥ 0} of mappings forms a one-parameter continuous semigroup satisfying the following conditions [18]: S τ x is continuous in both t and x ∈ X.
Let us now consider continuous time finite dimensional deterministic dynamical system with state vector x ∈ R n and vector field f : R n → R n , generated by the following set of autonomous ordinary differential equations (ODEs): We consider an autonomous dynamical system only for convenience's sake, since any non-autonomous system of n unknown variables (x 1 , . . . , x n ) of t can formally be seen as an autonomous system in n + 1 unknown variables (t, x 1 , . . . , x n ). The state variables (components of the state vector) representing a given atmospheric or climate system include such physical quantities as the air and ocean temperatures, barometric pressure, air density, humidity, wind velocity, and some others. It is very important that in atmospheric and climate modelling, we are interested in the so-called "viable" state of a system in question [19]. This means that the state vector belongs to the viable domain in state space known as the viability constraint set, which can be formally defined by the condition x < x • . As mentioned above, in weather prediction and climate modeling, we deal with discrete dynamical systems. Towards that end, using either a projection onto a finite set of basic functions or a discretization of derivatives in time and space, the continuous dynamical system (1) is transformed into discrete dynamical system x i+1 = f (x i ), which approximates the continuous system (1), and therefore can be solved numerically for given initial conditions. Here x i denotes n-dimensional vector of state variables at discrete times l ∈ Z + .
It is important to note that deterministic models of the atmosphere and climate possess some generic properties which must be borne in mind when modeling atmospheric and climate processes, as well as when developing methods for controlling them. Some of these properties are as follows. Firstly, deterministic models used in the numerical weather prediction are very sensitive to initial conditions precluding exact forecasting of weather and leading to chaotic oscillations in solutions of model equations. In other words, over time, a phase space trajectory of the system bears a resemblance to random process, although the system dynamics is determined by deterministic laws and described by deterministic equations. This phenomenon known as deterministic chaos was first uncovered by E. Lorenz [20]. Secondly, dynamical systems used in climate modeling are dissipative, implying that the divergence of vector field ∇· f (x(t)) is negative which leads to contracting the system's phase volume [21]. Thirdly, since climate dynamical systems are dissipative, there exists a bounded absorbing set in the phase space that attracts trajectories of climate system. This property assures the existence of a finite dimensional invariant compact attracting set, called the attractor, toward which a dynamical system tends to evolve over time [21].
Let us remark that small-scale physical processes cannot be always adequately represented on the space-time grid of existing deterministic atmospheric and climate models and are perceived by these models as noise. To take into account the influence of these subgrid physical processes on the atmospheric and climate dynamics, one can use stochastic dynamical systems generated by the following set of stochastic differential equations [22]: where X ∈ R n is the multidimensional stochastic process of interest possessing the initial condition X| t=0 = X 0 with probability one, W ∈ R d is a vector of independent Wiener processes, and g(X(t), t) is a matrix describing the dependence of the (sub-grid) noise on the state vector. Stochastic models can be seen as useful tools for exploring the response of the atmospheric (climate) system to random external forcing. In this case, the evolution of the probability density function p(X, t) is determined using the Fokker-Planck equation associated with Equation (2). However, the probabilistic formulation of the weather and climate optimal control problem is beyond the scope of this article.

Sensitivity Analysis of Atmospheric and Climate Models
Atmospheric and climate models contain hundreds or even thousands of parameters having different physical meanings. Some of them are empirical and semi-empirical, while some others are adjustable and not well defined. Errors in the model parameters generate parametric uncertainty. For convenience of further discussion, let us introduce the k-dimensional parameter vector α ∈ R k and then rewrite the set of autonomous ODEs (1) as follows: dx/dt = f (x, α), where x(0) = x 0 ∈ R n . The space-time spectrum of atmospheric and climatic processes is very wide, therefore it is hardly possible to develop a "univer-sal" model capable of covering a broad range of atmospheric and climatic fluctuations. If the physical (dynamical) process of interest has a characteristic time τ, then all processes whose scales are less than τ cannot be explicitly reproduced by the model and can only be described parametrically. Based on physical considerations and using sensitivity analysis, some of the model parameters can be selected as control variables [23]. Varying control parameters, we can formally control the dynamics of atmospheric and climate processes.
In conventional sensitivity analysis [24], in its attempt to analyze the influence of model parameters (forcing functions) on state variables, sensitivity functions (coefficients) are usually used, defined as partial derivatives of the state vector components with respect to the model parameters. For example, the sensitivity of the state variable x i (i = 1, . . . , n) with respect to the parameter α j (j = 1, . . . , k) is estimated using the following sensitivity function: j is a base value of the parameter α j at which the partial derivative is estimated. Assuming that δα j is an infinitesimal (infinitely small) change in the parameter α j with respect to its unperturbed value α 0 j , then, approximating the state vector around its base value x(α 0 j ) by the Taylor expansion, we get the following expression for the corresponding change in the vector of state variables: However, the absolute changes in the state vector caused by variations in different parameters do not allow the parameters to be ranked in accordance with the degree of their influence on the state vector, since the model parameters have different units. To compare the relative role of model parameters in changing δx, and thereby to rank the parameters, we can use relative sensitivity functions defined as Relative sensitivity functions play a very important role in the justification and selection of those parameters that can be considered as physically feasible control variables [23]. For example, considering atmospheric baroclinic instability as a controllable object, we have identified two fundamental parameters that govern the development of baroclinic instability in the atmosphere: the static stability parameter σ 0 and the vertical wind shear Λ 0 induced by the meridional temperature gradient. Analyzing the relative sensitivity functions, we found the critical value of the wavelength L cr of the growing modes that divides the spectrum of unstable waves into two parts. The growth rates of the amplitudes of unstable waves, the length of which is shorter than L cr , depend mainly on the static stability parameter, while if L < L cr , then the development of baroclinic instability is predominantly affected by the vertical wind shear, that is, by the meridional temperature gradient. For typical values of atmospheric parameters, we found that L er = 3800 km.
Sensitivity functions can be found by solving the set of linear non-homogeneous sensitivity equations obtained by differentiating the model equations with respect to the parameters: where S α (x) = ∂x i /∂α j ∈ R n×k is the sensitivity matrix, J x ( f ) ∈ R n×n and J x ( f ) ∈ R n×k are the Jacobian matrices given by These equations show how sensitivity functions evolve over time along the trajectory of a dynamical system. To estimate the sensitivity of the state vector x with respect to the parameter vector α, we need to solve a set of differential equations that includes the model equations and the sensitivity equations with given initial conditions. However, due to the large dimensionality of the parameter vector of modern atmospheric and climate models, "direct" sensitivity analysis, considered above, is a computationally expensive process. In addition, the derivation of sensitivity equations requires the sequential differentiation of model equations with respect to the desired parameters, which is difficult to perform for complex atmospheric and climate models. To overcome these difficulties, we can introduce some response function R(x,α) and then examine its sensitivity with respect to model parameters in a differential formulation using the adjoint model [25]. Formally, the response function can be written as follows: where Φ is a function that depends on the model state variables and parameters. For example, this function can characterize the total energy norm of a dynamical system Φ = x T Wx, where W is the matrix of weights defining the norm. The gradient of the response function with respect to α at the base point α 0 characterizes the sensitivity of the model state variables to the model parameters. This "sensitivity" gradient can be determined from a single numerical experiment by solving the following equation [25]: where the superscript "T" means "transpose", x 0 is the trajectory of the dynamical system corresponding to the parameter vector α 0 , and the vector-valued function x * is the solution to the adjoint model: The variation in the response function δR caused by the parameter variations is calculated as follows: δR(x 0 , α 0 ) = ∇ α R, δα , where δR(x 0 , α 0 ) = ·, · is a scalar product.
It should be underlined that sensitivity analysis of chaotic dynamical systems used in atmospheric and climate studies usually deserves special examination [26]. However, this topic is beyond the scope of this paper, since it is not related to weather and climate modification and control. A fairly complete and comprehensive survey of sensitivity analysis for nonlinear dynamical systems exhibiting, under certain conditions, chaotic behavior can be found in previously published papers (e.g., [27,28]).

Generic Formulation of the Optimal Control Problem
Intending to study weather and climate manipulations as an optimal control problem, we will consider a continuous time deterministic dynamical system, evolving over a given time interval t 0 , t f , where the terminal time t f > 0 can either be free or fixed. Without loss of generality, we will assume that t f is fixed. The state of a system at any instant of time t ∈ t 0 , t f is characterized by the vector x ∈ R n . We will suppose that the system is controllable. In other words, the system can be moved from some initial state to any other (desired) state in a finite time using certain external manipulations characterized formally by the vector of control variables u ∈ R m . By a system, we mean a set of ordinary differential equations of the form: where f : t 0 , t f × R n × R m → R n is a given vector-valued function. Note that the right-hand side of Equation (10) explicitly depends on the control variables, while uncontrolled parameters are omitted for convenience. It is worth noting that the control variables should be selected based on physical considerations, taking into account the feasibility of technical implementation of deliberate modification of weather and climate [4]. Thus, control variables are formally restricted to some control domain U called as the admissible control set of u: Generally, various physical constraints, expressed mathematically in the form of equalities and/or inequalities, can be imposed on state variables. Formally, this means that the state vector must belong to a certain phase-space domain: where X is a given subset of R n . Equations (10)- (12), considered together, restrict the set of admissible terminal values of state variables, i.e., x(t f ) = x f ∈ X f , where X f is the set of reachable states.
In essence, an optimal control problem is aimed at finding the control law for a given system that ensures the fulfillment of a control objective, which is usually the minimization of some function (functional) J(x, u), the objective functional (or performance index). We should note that the formulation of a performance index is dependent on the problem under study. To be more specific, can we consider the Bolza optimal control problem, assuming that no constraints are imposed on the state and control variables, and the initial and terminal times are fixed: and the boundary conditions An optimal control problem can be solved by Pontryagin's maximum principle, dynamical programming, or classical approaches of the calculus of variations.

Illustrative Example
Holding the increase in global mean surface temperature "to well below 2 • C above pre-industrial levels and pursuing efforts to limit the temperature increase to 1.5 • C above pre-industrial levels" was designated in the 2015 Paris Agreement on Climate as a priority area for combating global warming. This ambitious goal is expected to be achieved through the transition to low-carbon development, which, on the one hand, is a vital necessity, and on the other, a serious challenge. Nevertheless, many countries are already in the transition to a low-carbon economy, developing national strategies to achieve the Paris Agreement goals using, in particular, the Shared Socioeconomic Pathways scenarios of projected social and economic worldwide changes up to 2100 [29]. However, the Earth's climate system possesses significant inertia (e.g., [30]): it takes several decades for the climate system to reach a new equilibrium state in response to low (even zero) emissions of greenhouse gases (GHG). Therefore, after significant reduction of atmospheric GHG concentrations, surface temperature apparently will continue to rise. This phenomenon is known as "lag (time delays) between cause and effect". Thus, geoengineering can be considered as one of the technologically feasible options to stabilize the climate. In this paper, we leave "behind brackets" of such important aspects of geoengineering as physical side effects and environmental risks, ethical, legal and other social issues, recognizing the importance of their consideration and assessment (e.g., [31]).
Geoengineering can be implemented via human intervention in the redistribution of solar radiation flux due, for example, to the injection into the stratosphere of fine aerosol particles, which have the properties to scatter solar radiation in the visible spectral range and weakly absorb in the infrared range. For example, sulfate aerosols have such properties. Controlled emissions of sulfur dioxide or hydrogen sulfide (precursor gases) into the stratosphere ultimately lead to the formation of sulfate aerosol particles. The stratospheric aerosols increase the Earth's planetary albedo α 0 , change the radiation balance, and, as a consequence, decrease the Earth's surface temperature. Sensitivity analysis shows that an increase in α 0 by 1% leads to a decrease in the solar radiation flux at the top of the atmosphere by about 3.4 W/m 2 , which is comparable with the radiation effect of doubling the concentration of atmospheric CO 2 . To assess the effectiveness of geoengineering projects and their consequences for nature and society, numerical modeling is used for given scenarios of anthropogenic GHG emissions and atmospheric concentrations, and heuristically specified geoengineering scenarios. It is obvious that going through all possible options of intentional geoengineering manipulations is an ineffective approach. In contrast to that, we consider climate manipulation within the framework of the optimal control theory. To illustrate this approach, we will apply the two-box low-parametric climate model, taking into account the radiation effects of the sulfate aerosols artificially injected into the stratosphere. The model equations are as follows [32]: where T and T D are temperature anomalies for the upper (the atmosphere and mixing ocean layer) and lower (the deep ocean) boxes, C and C D are the effective heat capacities for the upper and lower boxes, λ is a climate feedback parameter, γ is a coupling strength parameter describing the rate of heat loss by the upper box, ∆R GHG is the radiative forcing produced by GHG, and ∆R A is the radiative forcing generated by stratospheric aerosols.
The two-box model, in spite of its simplicity, is capable of simulating globally averaged climate change caused by human-induced radiative forcing with a reasonable accuracy (e.g., [33,34]). The following values of model parameters have been used in calculations [31]: C = 7.34 (W yr)/ m 2 K , C D = 105.5 (W yr)/ m 2 K , λ = 1.13 W/ m 2 K and γ = 0.7 W/ m 2 K . Radiative forcing ∆R GHG is approximated by a linear function ∆R GHG = ηt, where the parameter η is determined from the Representative Concentration Pathway (RCP) data [35]. For the worst-case emission scenario (RCP8.5 [35]), the annual radiative forcing rate is η = 7.14·10 −2 W/ m 2 yr . Radiative forcing produced by aerosols is calculated by the formula: ∆R A = −α A Q 0 , where α A is the albedo of the aerosol layer (note that α A << 1), Q 0 = 342 W/m 2 is the mean insolation on the top of the Earth atmosphere. This allows albedo α A to be considered as a control variable. However, in reality, we have the ability to control the rate of aerosol emissions E A , which is included in the aerosol mass balance equation: where M A is the total mass of stratospheric aerosols, E A is the aerosol emission rate and τ A is the residence time of stratospheric aerosol particles. The mass of aerosols M A is linearly related to the albedo α A : M A = α A (Q 0 S e /β A k A ), where β A = 24 W/m 2 is the empirical coefficient, k a = 7.6 m 2 /g is the mass extinction coefficient of aerosol particles, S e is the Earth's surface area.
In practice, precursor gases are injected into the stratosphere, therefore the mass of sulfate aerosols and their emission rate are expressed in sulfur units and denoted as E S (TgS/year) and M S (TgS), respectively, taking into account that 1 Tg of sulfur is equivalent to 4 Tg of aerosol particles. Then, Equation (17) can be rewritten as where χ = Q 0 S e /4β A k a ≈ 2.39·10 2 TgS. If the optimal albedo of aerosol layer is determined, then the optimal emission rate E * S (t), which provides the formation of an aerosol layer of the optimal mass M * S (t), is calculated using Equation (18).
The optimal control problem will be considered with regard to the finite time interval t ∈ [t 0 , t e ] on which the dynamics of the control object are described by Equation (16) with given boundary conditions: Thus, the left end of the phase trajectory is fixed, and the right end is fixed only for the variable T, while the variable T D is free. These boundary conditions are chosen since the surface temperature anomaly T is of primary interest. The optimal control problem is formulated as follows: on the finite time interval [t 0 , t e ] find the control variable α * A (t) belonging to an admissible value domain, so that when the dynamic constraints (16) and boundary conditions (19) imposed on the system are satisfied, the given functional characterizing the mass flow rate of aerosols (20) has reached its minimum value.
The terminal condition T f represents a target change in the global mean surface temperature at t = t f . As an example, we assume that T f = 2 • C. In essence, the performance index characterizes the consumption of aerosols for geoengineering manipulations since the albedo of the aerosol layer α A is a linear function of M A . So, we aim at minimizing the mass of aerosols required to achieve the target surface temperature change at final time which is fixed. The amount of aerosols that can be annually delivered to the stratosphere can be limited by the available technical capabilities. Therefore, we will formally assume that the domain of admissible controls is an interval α A ∈ [0, U], where U is the maximum value of technically feasible albedo α A .
To solve the optimal control problem, we use the Pontryagin's maximum principle (PMP), which is the major tool in optimal control. The Hamiltonian function used to solve the problem of optimal control for dynamical system (16) is given by: where a = (λ + γ)/C; b = γ/C; c = η/C; q = (1 − α 0 )Q 0 /C; p = γ/C S ; ψ 1 and ψ 2 are the time-varying Lagrange multipliers that satisfy the adjoint equations: The optimal control α * A (t) ∈ [0, U] at each fixed time t ∈ [t 0 , t f ] must be such that In other words, the optimal control results in a maximum value of H at any t ∈ [t 0 , t f ] The corresponding stationarity conditions for Hamiltonian yields: To find the optimal control α * A and the optimal surface temperature anomaly T * generated by α * A , we must solve the system of four ordinary Equations (16) and (21) in four unknown variables ψ 1 , ψ 2 , T, and T D , with given initial and terminal conditions (18). Since no terminal condition is specified for the variable T D , the following transversality condition ψ 2 t f = 0 is used in calculations. We have derived the following analytical expressions for the optimal albedo of aerosol layer α * A and the corresponding optimal surface temperature anomaly T * : where λ 1 and λ 2 are the eigenvalues of the coefficient matrix of the adjoint system (22), ν 11 and ν 21 are the components of the corresponding eigenvectors, C 1 , C 2 , C 3 , and C 4 , are arbitrary integration constants, while , Considering geoengineering as a state-constrained optimal control problem additional necessary conditions for optimality, the complementary slackness condition, must be specified. The meaning of the path constraint condition (28) is that the global mean surface temperature anomaly should not exceed the value of threshold parameter C T , which is set a priori.
As an example, we consider the results of calculations for the RCP8.5 emission scenario (as we mentioned earlier, this is the most conservative scenario in relation to the growth of atmospheric GHG concentrations). The optimal control problem is examined on the finite time interval 2020-2100. In other words, t 0 = 2020 and t f = 2100. Since the temperature anomalies are calculated relative to 2020, the initial conditions for variables T and T D are as follows: T 2020 = 0 and T D,2020 = 0, where the numerical subscript is referred to the year 2020. We assume that: -By 2020, the surface temperature anomaly would exceed the pre-industrial level by 1.1 • C, i.e., ∆T 2020 = 1.1 • C -By 2100, the surface temperature anomaly would exceed the pre-industrial level by 1.5 • C; -For the 2020 to 2100 period of time, the increase in T should not exceed 2 • C above the pre-industrial level.
Then the allowable temperature growth by year 2100 relative to 2020 would be T 2100 = ∆T 2100 − ∆T 2020 = 0.4 • C. This value is taken as a terminal condition for variable T at t = t f . The threshold parameter, which defines a path constraint, is Calculations show that in the absence of deliberate interventions in the Earth's climate system, the growth of global average surface temperature over the 80-year (2020-2100) interval is about 3.8 • C which is significantly higher than the level established by the Paris Agreement. Figure 1 presents two curves, one for the optimal albedo of the aerosol cloud and the other for the corresponding surface temperature anomaly, as functions of time calculated for the case (a) when no constraints are imposed on the control variable and the surface temperature anomaly, and (b) with constraint imposed on the surface temperature anomaly. As can be seen from Figure 1b, in the absence of limitations on the surface temperature increase, beginning from 2060, some overshoot beyond the threshold C T is observed. However, stratospheric aerosols, artificially injected into the stratosphere, guarantee compliance with the specified terminal condition. Recall that this condition is as follows: T 2100 = 0.4 • C at t = 2100.

Discussion
In response to human-induced climate change, which is one of the main threats natural and anthropogenic systems in the 21st century, the scientific meteorologi community has proposed several technologies known as geoengineering to deceler global warming and stabilize the climate. However, all of these technologies, includi solar radiation management (SRM) considered in this paper, are still theoretical in natu since "Playing with the Earth's climate is a dangerous game with unclear rules" [3 Nevertheless, scientists from different countries continue geoengineering research usi mainly computer simulations expecting that geoengineering techniques could potentia be helpful in the future to enhance ongoing efforts to mitigate global warming by ducing anthropogenic GHG emissions. Meanwhile, in the vast majority of studies heuristic approach is used to developing geoengineering strategies and scenarios (e [11,[37][38][39][40] and references therein). This is because global climate models used to proj climate change are extremely complex. This important circumstance is a serious obsta to the application of optimization methods, and, in particular, optimal control in geoe gineering problems. However, a heuristic technique is not guaranteed to be optimal a The optimal control problem with state constraint requires the consideration of an additional necessary condition termed as a condition of complementary slackness: For the case of C T = 0.9 • C (see above), the optimal albedo and the corresponding surface temperature anomaly are shown in Figure 1 in red color. As can be seen from this figure, by using state constraint, it is possible to prevent an undesirable growth of surface temperature within a given time interval. The amount of aerosols required over an 80-year period to maintain geoengineering is 36.5 TgS for the unconstraint case and 73.6 TgS for the case with state constraint. For RCP6.0 scenario (this scenario is described by the Intergovernmental Panel on Climate Change as an intermediate scenario in which emissions' peak is around 2080, then decline) we obtained that the total amount of aerosols is 17.0 Tg for the unconstrained optimal control problem, and 23.3 Tg for the problem with state constraint.

Discussion
In response to human-induced climate change, which is one of the main threats to natural and anthropogenic systems in the 21st century, the scientific meteorological community has proposed several technologies known as geoengineering to decelerate global warming and stabilize the climate. However, all of these technologies, including solar radiation management (SRM) considered in this paper, are still theoretical in nature since "Playing with the Earth's climate is a dangerous game with unclear rules" [36]. Nevertheless, scientists from different countries continue geoengineering research using mainly computer simulations expecting that geoengineering techniques could potentially be helpful in the future to enhance ongoing efforts to mitigate global warming by reducing anthropogenic GHG emissions. Meanwhile, in the vast majority of studies, a heuristic approach is used to developing geoengineering strategies and scenarios (e.g., [11,[37][38][39][40] and references therein). This is because global climate models used to project climate change are extremely complex. This important circumstance is a serious obstacle to the application of optimization methods, and, in particular, optimal control in geoengineering problems. However, a heuristic technique is not guaranteed to be optimal and rational, but this approach is sufficient for obtaining approximate (satisfactory) solutions in the case when finding an optimal solution is impossible. As a result, theoretical studies of geoengineering operations are usually performed outside the framework of optimization theory. Note that the optimal (best) solution is a solution that is preferable for one reason or another. More specifically, the optimal decision (option, choice, etc.) is the best decision among admissible alternatives if there is a rule of preference for one over the other, known as the optimality criterion. One can make a comparative assessment of possible decisions (alternatives) and choose from the best based on such criterion. In our case, optimality criterion is the total mass of aerosols injected into the stratosphere.
In relation to the problem considered in this paper, the commonly used procedure for determining the albedo of aerosol layer (and, consequently, the total mass of aerosol particles injected into the stratosphere) required to counteract temperature rise is, in a simplified manner, as follows. First, we select the GHG emissions scenario of interest (e.g., RCP8.5) and then calculate the corresponding radiative forcing due to changes in concentrations of the relatively well-mixed GHG, using, for example, the first-order approximation expression [41]: ∆R CO 2 (t) = κ ln(C t /C 0 ), where κ = 5.35 W m 2 is the empirical parameter, C t is the CO 2 -equivalent concentration in parts per million by volume at time t, C 0 is the reference concentration. Second, we assume that stratospheric aerosols, partially reflecting solar radiation back to space, produce radiative forcing ∆R A (t) = −ς∆R CO 2 (t), where the parameter ς defines the portion of anthropogenic radiative forcing that should be neutralized by stratospheric aerosols. Note that if the parameter ς is set to 1, then anthropogenic radiative forcing will be completely neutralized. Next, we calculate the albedo of aerosol layer as follows (see Section 3): α A (t) = −∆R A (t)/Q 0 = ς∆R CO 2 (t)/Q 0 . According to this equation, the albedo α A is linearly related to the radiation forcing ∆R CO 2 . The corresponding aerosol emission rate at time t can be found using the Equation (18). The solution found in such a manner is in a certain sense non-optimal, although the SRM problem is essentially mathematically solved. With the parameter values used in this paper (see Section 3), we found that for the RCP8.5 scenario, assuming that the anthropogenic radiative forcing will be completely compensated by stratospheric aerosols (ς = 1), and also assuming that T 2100 = T 2020 , the albedo α A should increase linearly from zero in 2020 to~0.02 in 2100.
The study of geoengineering operations within the framework of the optimal control theory, which is a branch of mathematical optimization, allows for obtaining, in a certain sense, the optimal solution, considering various constrains, which, of course, must be formalized mathematically in the form of equalities and/or inequalities. We emphasize that the dynamics of the system are defined as optimal only with respect to the selected performance index. Note that the goal-setting problem (formulation of the performance index) is non-trivial and requires special consideration.
In the present study, a very simple climate model, the two-box energy balance model (EBM) is used, allowing for predicting the globally averaged surface temperature from the analysis of the planetary energy balance. The model is linear and does not describe the climate system dynamics. A primary motivation of using such a model is that similar models have been examined in a number of research papers studying the essential features of climate system response to natural and human-induced radiative perturbations that affect the Earth's climate system. Despite its simplicity, the two-box EBM was able to reproduce the evolution of global mean surface temperature over time in response to timedependent radiative forcing with reasonable accuracy [33]. To demonstrate the applicability of new mathematical approaches (in our case, the optimal control theory) to "real world problems", the use of simple modeling tools is a common and useful practice. However, caution should be exercised in evaluating and interpreting the results obtained using this approach. We considered the applicability of optimal control theory to hypothetical deliberate weather and climate modification and, in particular, to the SRM approach of geoengineering. An algorithm for solving the problem of optimal control of the Earth's climate using classical Pontryagin's maximum principle was presented, both with and without constraints imposed on the state variable. In fact, the Earth's climate is a highly complex nonlinear dynamical system with feedbacks and cycles that affect the climate response to forcing generated by GHG [4]. However, nonlinear problems of the Earth's climate system optimal control can be solved mainly numerically using highly complicated coupled general circulation models of the atmosphere and ocean. Undoubtedly, such problems are of extreme complexity. Therefore, we believe that the application of optimal control theory in combination with simple atmospheric and climate models will be very useful in the design of weather and climate control systems, as well as in the development of scenarios for intentional modifications of climate and weather.
It is important to note that two classical mathematical tools for studying optimally control systems are Pontryagin's maximum principle and Bellman's Dynamic programming. However, the application of these methods is very difficult or even impossible if nonlinear models of high complexity are considered. In order to overcome the difficulty, some asymptotic and approximate approaches can be used (e.g., [42][43][44]).
It should be underlined that the results obtained in this paper serve only to illustrate the applicability of optimal control theory to climate modification problems since the use of very simple models of climate systems allows one to solve the optimal control problem analytically. In the next studies, we intend to apply the considered approach to solving the problems of modifying various atmospheric processes and phenomena, as well as to examine a number of geoengineering problems using more complex climate models.

Concluding Remarks
In this paper, we introduce the optimal control-based method for planning and executing weather and climate manipulation projects. The application of this technique is demonstrated using the two-box energy balance model in which the annual emission rate of aerosol precursors is the control variable, while the global mean surface temperature is the main state variable. The optimal control problem in both state unconstrained and constrained formulations is analytically solved using the classical PMP. Our approach provides additional insights for the development of optimal climate manipulation strategies to counter global warming in the 21st century.
The majority of prior geoengineering research has applied some heuristic considerations to set up aerosol emission scenarios. In this paper, we propose considering geoengineering problems within the optimization framework applying methods of the optimal control theory. This approach allows one to obtain a geoengineering scenario by rigorously solving the optimal control problem for a given performance index (objective function). Since the solution was derived in a simplified mathematical formulation, the results obtained should be considered mainly for illustration purposes. In general, the paper serves to demonstrate the capabilities of optimization methods in solving problems of weather and climate modification. We expect this work will attract researchers' attention to explore geoengineering using classical and approximate methods of the optimal control theory.