Mathematical Modelling of Climate Change and Variability in the Context of Outdoor Ergonomics

: The current climate change, unlike previous ones, is caused by human activity and is characterized by an unprecedented rate of increase in the near-surface temperature and an increase in the frequency and intensity of hazardous weather and climate events. To survive, society must be prepared to implement adaptation strategies and measures to mitigate the negative effects of climate change. This requires, ﬁrst of all, knowledge of how the climate will change in the future. To date, mathematical modelling remains the only method and effective tool that is used to predict the climate system’s evolution under the inﬂuence of natural and anthropogenic perturbations. It is important that mathematics and its methods and approaches have played a vital role in climate research for several decades. In this study, we examined some mathematical methods and approaches, primarily, mathematical modelling and sensitivity analysis, for studying the Earth’s climate system, taking into account the dependence of human health on environmental conditions. The essential features of stochastic climate models and their application for the exploration of climate variability are examined in detail. As an illustrative example, we looked at the application of a low-order energy balance model to study climate variability. The effects of variations in feedbacks and the climate system’s inertia on the power spectrum of global mean surface temperature ﬂuctuations that characterized the distribution of temperature variance over frequencies were estimated using a sensitivity analysis approach. Our conﬁdence in the obtained results was based on the satisfactory agreement between the theoretical power spectrum that was derived from the energy balance model and the power spectrum that was obtained from observations and coupled climate models, including historical runs of the CMIP5 models.


Introduction
Mathematics represents a very effective and powerful instrument for comprehending the world and solving complex problems in various sciences, engineering and technologies (e.g., [1][2][3][4]).In this aspect, one cannot underestimate the essential role of mathematics in solving planetary problems, among which, the problem of the interaction between society and nature stands out [5,6].Humans live and execute their diverse and multifaceted activities in interaction with the environment.On the one hand, humans affect the environment, changing its properties; on the other hand, environmental conditions affect humans, in particular, their health and even their ability to survive.At the turn of the millennium, humankind is clearly faced with a new pressing challenge posed by climate change.Current climate change, unlike previous ones, is human-induced [7,8] and characterized by an unprecedented rate of increase in the global mean surface temperature (GMST).A report [9] stated that since 1880, the GMST has increased by an average of 0.07 • C per decade.Meanwhile, the growth rate of the GMST in the first two decades of the XXI century reached a value of 0.17 • C per decade, i.e., more than double the average rate.Among the most significant signs of human activities is a substantial increase in the concentration of carbon dioxide (CO 2 ) in the atmosphere.This greenhouse gas (GHG), which is a by-product of fossil fuel combustion, is the main contributor to global warming.
Global warming has already affected climatic conditions in almost all regions of Earth [7,8], but this effect varies significantly from region to region.Obviously, global and regional climate changes have some impact on people and their activities [10,11].Assessing the impact of changing climatic conditions on humans involves assessing climate risks [12].For this purpose, special mathematical models are used.Their input parameters are climatic variables, such as temperature, humidity and precipitation (e.g., [13]).In turn, climate models of varying degrees of complexity are used to predict climate change and variability on global and regional scales [14].It should be emphasized that the existing estimates of global and regional climate changes are characterized by a high degree of uncertainty [15,16].The main causes of these uncertainties are considered, for example, in [17][18][19].The need of stakeholders for a deeper and clear understanding of what is happening with the climate of our planet makes it necessary to obtain new knowledge about the behaviour of the Earth's climate system (ECS) under the influence of external perturbations, both natural and human-made.This is very important for preparing society for the development and implementation of the so-called adaptation strategies and measures that are required to reduce the effects of climate change.By adaptation, we mean the process of society's adjusting to current and projected changes in the climate system and their effects [7].
Mathematics and its methods and approaches already play a vital role in climate research (e.g., [20,21]).One can say without exaggeration that only mathematics allows us to quantify and predict the effects of external natural and anthropogenic perturbations on the climate system and its dynamics.In fact, quantitative research is based on mathematical models of systems and objects under examination.In this regard, we note that reproducing climate and projecting its change under the influence of external forcing, in contrast to conventional problems of physics, have the essential peculiarity that is associated with the impossibility of carrying out full-scale direct physical experiments.Due to several specific properties of the climate system [22], the implementation of laboratory experiments looks very problematic as well.It should be added that the time series of climatic variables, which are obtained from instrumental observations, are short and contain data for only a few decades, which makes it difficult to obtain statistically significant estimates of the state of the climate system.This evidence suggests that mathematical modelling remains the only method for projecting the evolution of the climate system under the influence of natural and anthropogenic perturbations.
In this study, in the context of outdoor environmental ergonomics, we examined some mathematical methods and approaches, primarily mathematical modelling, that are applied to climate research.Stochastic dynamical systems that are applied to study climate variability and sensitivity analysis are explored in detail.

Notes on Outdoor Environmental Ergonomics through the Prism of Climate Change
According to the Merriam-Webster Dictionary [23], "ergonomics is an applied science concerned with designing and arranging things people use so that the people and things interact most efficiently and safety".It is also called "human engineering" or "human factors engineering".Ergonomics' areas of specialization include environmental ergonomics, both indoor and outdoor, which explores how humans interact with their physical environment, as characterized by climatic conditions and several other environmental parameters.Since humans are compelled to perform many types of work outdoors, their vital activity, working capacity and health are largely determined by the weather and climate conditions.Weather and climate can affect workers in different ways, causing a variety of health problems and changing their working efficiency.In this sense, ergonomics is of great importance in reducing health risks that arise from adverse weather and climate conditions.In order to use ergonomic techniques for assessing climate-related risk factors, one should know, first, the current climate conditions in a particular geographical area of interest and, second, how the climate is likely to change in the near term.It is clear that current human-induced global warming affects climate around the world.Humans working outdoors can be expected to be among the first to be affected by climate change [24][25][26].
Climate affects outdoor workers in several ways: via air temperature, wind, humidity, precipitation and some other natural hazardous events, the frequency and intensity of which change with global warming.Climate change apparently could exacerbate health risks via the increase in frequency and intensity of weather and climate-related events, such as cold and heat waves, heavy rains and runoff, severe storms, droughts, floods, biochemical hazards and air pollution.Climate change also contributes to depletion in the atmospheric ozone content, increasing the ultraviolet (UV) radiation at the Earth's surface and thereby causing outdoor workers more intense exposure to UV radiation.In addition, climate conditions significantly determine the requirements for protective equipment for outdoor workers.This makes it necessary to predictably take into account the peculiarities of the climate in various geographical regions regarding protective equipment usage, taking into consideration the long-term use of the type of protective equipment that is adopted.The intensity of many potentially dangerous physical factors (e.g., acoustic, electromagnetic, vibration) depends significantly on the conditions of their propagation, which are determined by the climate variables in a given geographical region of professional activity.Difficulties in taking into account the current and prognostic characteristics of the climate when conducting ergonomic studies are determined both by the insufficient development of mathematical models of the combined effect of physical factors on a human (primarily due to the lack of experimental information) and the insufficient development of mathematical models that simultaneously take into account the combined effects of physical and psychophysiological factors.The development of such mathematical models is currently being intensively pursued by mathematicians together with engineers, as well as physiologists.In the first case, the priorities of developing models are associated with the implementation of technologies of digital twins of the human body, which make it possible to predict its reactions to factors of the conditions of activity for representatives of various socio-professional groups of workers.In the second case, the priorities of the development of models are associated with the study of the mechanisms of the influence of physical factors on a person.In particular, the Nobel Prize in Physiology or Medicine 2021 was awarded for the study of receptors for temperature and touch, the results of which are undoubtedly important for the successful solution of the problems that are mentioned in this article.Climate is defined as the average weather over 30 or more years (World Meteorological Organization recommended using 30 years).Climate variables, such as temperature, humidity and wind speed, are used as inputs in methods and models of environmental ergonomics.In turn, climate variables can be derived from climate projections, which are simulations of the planetary climate in future years and decades based on suggested scenarios for the concentrations of atmospheric greenhouse gases, aerosols and radiatively active gases.Climate projection results are based on a hierarchy of climate models of varying complexity, ranging from coupled atmosphere-ocean general circulation models and Earth system models to low-order conceptual climate models.To obtain climate projections, all these models are exposed to radiative forcing, the level of which is determined by emission scenarios of greenhouse gases and other radiatively active components.Thus, climate models remain the only viable tool for quantifying the climate change that is caused by natural and human-induced perturbations.Meanwhile, climate change projections for the XXI century that are obtained from different models are inherently uncertain.The main sources of uncertainty, which were analysed and summarized in several studies (e.g., [15][16][17]19,27,28]), are as follows: uncertainty in greenhouse gas emission scenarios, natural climate variability and climate models' uncertainties.Since we focus on climate system simulations, the climate models' uncertainties are of primary concern.Climate models, in general, qualitatively agree in their response to external forcing but quantitatively differ in the projected climate change for a given increase in concentrations of atmospheric greenhouse gases, aerosols and radiatively active gases.In addition, each climate model has an individual "internal variability", which is significantly different from the real world climatic variability [29].This is exactly why climate models are constantly being revised and improved based on advances in knowledge.

Climate as a Complex Dynamical System
A system is understood as a set of interrelated and interdependent elements that work together to achieve certain goals [30].A system is called complex if it has properties such as emergent behaviour, self-organization, a hierarchical structure, nonlinearity and spontaneous order.The reasons for the emergence of complexity in systems are quite numerous.For example, nonlinear space-time interactions between system components lead to the emergence of new dynamical properties that cannot be observed by exploring individual elements of the system.In this context, the ECS is a complex large-scale physical system that consists of five basic and interacting subsystems [31]: the atmosphere, hydrosphere, cryosphere, lithosphere and biosphere.Each of these subsystems is characterized by a finite set of distributed variables, the values of which at a given time determine the state of the subsystem.The atmosphere is the most rapidly fluctuating and unstable ECS element.The climate system of our planet is complex for the following reasons: - As mentioned earlier, it is a five-element large-scale physical system with several global hydrological and biochemical cycles.Its elements, being heterogeneous thermodynamical systems, have significant differences in their structure, dynamics, physics and chemistry.Dynamical and physical processes that occur in the ECS subsystems differ in their scales, both spatial and temporal.Elements of the ECS link together through numerous physical coupling mechanisms, both weak and strong, including feedback mechanisms.In turn, each ECS subsystem can be viewed as complex, consisting of subsystems, which themselves are composed of low-order subsystems.
The atmosphere, for instance, can be divided into several vertical layers depending on its thermal stratification: the troposphere, stratosphere, mesosphere and thermosphere.Each ECS component has a specific response time, which must be considered when building ECS models.For example, the atmosphere can be considered the only component of the ECS model for dynamical processes with timescales from days to weeks since the tropospheric response time is about one month, while oceans, land surface and ice cover can be used to specify boundary conditions and/or external forcing.To explore the dynamical processes with timescales ranging from months to years, the ECS model must include the atmosphere and ocean, along with sea ice.Thus, ECS models are built from a hierarchy of models that ultimately form a complex integrated model.- The ECS is an oscillating system that is characterized by fluctuations that are under the influence of internal factors (natural oscillations), as well as due to external perturbations (forced oscillations).Naturally occurring fluctuations result from internal instabilities (e.g., hydrodynamic instabilities, such as barotropic and baroclinic instabilities) and heat transport within the climate system that is caused by the interaction of its components [7].Intentional and unintentional human impacts belong to the class of external forcings.-Since the ECS exchanges energy with the environment, in this sense, it is a thermodynamically open and non-isolated system.However, the ECS is a closed system with regard to the exchange of matter with the environment.The main energy source that drives the ECS is solar energy.Climate is affected by variations in external driving forces that imply natural causes, such as fluctuations in solar and volcanic activities, as well as changes in the gaseous and chemical composition of the atmosphere due to anthropogenic factors.The impact of the ECS on outer space is nonessential.At present, climate change is most influenced by variations in the composition of atmospheric particles and gases.Carbon dioxide (CO 2 ), the concentration of which in the atmosphere has been continuously increasing since the Industrial Revolution, has the greatest impact on current climate change.- The ECS and its subsystems possess emergent properties, examples of which are atmospheric emerging phenomena, such as clouds, large-scale baroclinic and barotropic eddies (cyclones, hurricanes) and small-scale eddies (tornadoes).An example of an emerging climate event is the El Niño-Southern Oscillation, which is a sporadic quasiperiodic variation in ocean surface temperatures across the tropical Pacific Ocean that affects the global atmospheric circulation and ocean circulation patterns.Natural emergent phenomena occur suddenly under some favorable conditions.
The ECS has several other important and specific properties that should be considered when exploring climatic processes.However, even the properties listed above allow us to conclude that only mathematical modelling can serve as the main tool for studying the ECS and, in particular, climate projections.The experience of recent decades shows that the main results of climate theory were obtained using global climate models.Ideally, the models of the ECS should allow for reproducing the modern climate, exploring its sensitivity to external perturbations and predicting the state of the ECS for many years and even decades to come.It should be emphasized that methods of the dynamical systems theory constitute the basis of methods of the mathematical theory of climate [32].To apply this theory to the study of a real climate system, it must be assigned with some abstract mathematical object that represents the idealization of the ECS and can be called its "perfect" model.It is assumed that such a perfect model exists and the observed dynamics of the ECS represents a realization of the trajectory that is generated by this model.It is usually believed that a perfect ECS model is a deterministic semi-dynamical system that is dissipative, ergodic and possesses a global attractor, while any trajectory of it is unstable [30].In a formalized form, a model of the ECS is a set of 3D nonlinear differential equations in partial derivatives that generates a deterministic finite-dimensional semi-dynamical system [32][33][34]: where x is the vector of state variables characterizing the system state at a given time t, x 0 is a given initial state of a system, n is the dynamical system dimension, α ∈ R m is the m-dimensional parameter vector and f is an external forcing.Since analytical solutions for the climate system models are commonly not available, we used numerical methods to solve Equation (1).To this end, the original system of partial differentials in Equation (1) was replaced by discrete space-time approximations using any suitable technique (e.g., Galerkin approach, finite-difference method).Thus, in climate modelling, we consider discrete dynamical systems.The numerical solution of the discrete model equations is only possible with the use of high-performance computers.Due to the discrete structure of ECS models, a large number of physical processes and cycles cannot be explicitly represented on the model grid.Theoretically, discrete ECS models are unable to realistically reproduce the physical process on spatial scales of the order of two-fold resolution of the model grid [35].Such processes are usually parameterized, i.e., are expressed parametrically in a simplified manner.The generalized discrete climate model that is used in computer simulations can be formally written as follows: where M 0,k is a nonlinear model operator that propagates state variables from the initial time t 0 to time t k (this operator indirectly contains model parameters, including external forcing) and K is the number of time steps.
For the sake of completeness, it makes sense to list some of the generic properties of dynamical systems that are used in climate research [32,33]: -These systems belong to the class of dissipative dynamical systems that possess (strange) attractors.-Trajectories of climate dynamical systems are generally unstable in the Lyapunov sense.-Some unstable trajectory, enclosed in an attractor, generates deterministic dynamical chaos.
In climate research, the ECS evolution is generally considered on its attractor, assuming that the ECS is an ergodic dynamical system.This allows for calculating the statistical characteristics of a system (e.g., first, x = x , and second, var(x) = x 2 − x 2 moments) by averaging along a certain trajectory [32].It should be noted that attractors of dissipative dynamical systems, which are commonly called strange attractors, are characterized by a very complex fractal structure.To determine their fractal dimension, the Kaplan-Yorke conjecture is used [36].
Deterministic climate models are very useful for exploring and projecting long-term climate change trends.To study the natural climate variability against the background of human-induced global warming, it is advisable to use stochastic climate models, which were first suggested by Hasselmann [37]: where x ∈ R n is the multidimensional random process satisfying the initial conditions x| t=0 = x 0 with probability 1, W ∈ R s is a vector of independent Wiener processes and g(x, t) is a matrix that describes the dependence of the sub-grid noise on the state variables.
In climate research, noise can play the role of fluctuations of solar insolation at the top of the atmosphere, sea surface temperature oscillations, weather processes, etc. Stochastic models are a useful tool for studying the climate system's response to random external perturbations.In this case, the Fokker-Planck equation that is associated with Equation ( 3) can be used to trace the evolution of the probability density function p(x, t).

Climate System Sensitivity to External Forcing
External forcing in climate models is described via some model parameters.Consequently, the climate system's response to external perturbations can be estimated on the basis of a sensitivity analysis of dynamical systems [38], allowing us to ultimately assess the climate change that is caused by natural factors and human activities.The choice of sensitivity analysis method largely depends on the objectives of the study and the complexity of the model that is used in computational experiments.Conventional sensitivity analysis methods, both forward and adjoint, are very useful for simple (low-order, conceptual) climate models.Using these methods, the effect of model parameter uncertainties and/or variations in model dynamics is estimated via sensitivity functions that are defined as the derivatives of model state variables with respect to the parameters.For a state variable x i (i = 1, . . ., n), a sensitivity function with respect to the parameter α j (j = 1, . . ., m) is defined as S ij = ∂x i /∂α j α 0 j , where α 0 j is the parameter value around which the sensitivity function is estimated.If δα j is the uncertainty (or small change) in the parameter α j , then the corresponding change δx i in the state variable x i is estimated as follows: δx i δα j ≈ δα j S ij .
In some cases, when we use simple climate models, the sensitivity functions can be calculated from a system of inhomogeneous linear equations called sensitivity equations, which are obtained by differentiating the model equations with respect to the parameters: where S α = ∂x i /∂α j ∈ R n×m is the sensitivity matrix, J x (F ) ∈ R n×n and J α (F ) ∈ R n×m are the Jacobians given by J x (F ) = ∂F i /∂x j and J α (F ) = ∂F i /∂α j .
To calculate the sensitivity functions, it is required to solve the model equations, together with the sensitivity equations, with given initial conditions.Since the model parameter vector is m-dimensional, to estimate the model's sensitivity to its parameters, we need to solve m times the model equations together with sensitivity equations.This technique is suitable for simple climate models with few parameters.In addition, it is not always possible to differentiate the model equations with respect to all parameters.In climate research, the model response to changes in parameters can be represented by some generic response functional that characterizes the climate model [39,40]: where ϕ is a nonlinear function of model state variables x and parameters α.
The gradient of R with respect to parameter vector α, calculated in the vicinity of the point x 0 , α 0 , is used to quantify the impact of parameter variations on the model dynamics: where dR/dα j x 0 ,α 0 = ∑ n i=1 (∂R/∂x i ) ∂x i /∂α j + ∂R/∂α j .
Using complex climate models, the sensitivity to the parameter α j can be estimated numerically as follows: In this case, however, the accuracy of the sensitivity estimate is very dependent on the magnitude of the perturbation δα j .To overcome this problem, we can introduce the Gâteaux differential [39], which allows us to consider the sensitivity analysis in a differential formulation.The Gâteaux differential for the response function R is of the form: where δx is the variation in the state vector due to the parameter vector variation in the direction δα.
If we linearize the model equations around the unperturbed trajectory, then we derive a set of variational equations for calculating δx: Then, the variation in the response function δR is calculated using Equation ( 9).This approach is ineffective from the computational viewpoint if the model has a large number of parameters.Using the adjoint approach, we can estimate the sensitivities in just one numerical experiment since the gradient of the response function can be computed from the following equation [40]: where x * is the solution of the adjoint equations integrated numerically in the inverse time direction:

Sensitivity Analysis of a Chaotic Dynamical System
Sensitivity analysis based on the forward and adjoint conventional approaches discussed above does not yield positive results when we explore deterministic dynamical systems that exhibit chaotic behaviour [41][42][43].This is because general solutions of sensitivity equations for dynamical systems of this kind grow unbounded as time tends to infinity.As a result, the sensitivity functions obtained via conventional approaches are highly uncertain.A deterministic dynamical system that exhibits chaotic behaviour is very sensitive to initial conditions [44]; consequently, the rate of separation in the phase space of the two trajectories that are initially infinitely close to each other exponentially grow with time as δx(t) ≈ δx(0)e λt , where λ is the leading Lyapunov exponent and δx(0) is the initial separation, resulting in large errors in the calculation of sensitivity functions.The largest Lyapunov exponent for the chaotic dynamics is positive, i.e., λ > 0; therefore, the two orbits diverge exponentially in time.Sensitivity functions that are calculated using this manner are extremely inaccurate and inherently uninformative such that it does not allow us to correctly estimate the system's sensitivity to variations in its parameters.This is primarily because the time-averaged sensitivities ∇ α R(↑ α) cannot be correctly estimated using conventional sensitivity analysis methods since for chaotic dynamics does not converge uniformly and the two limits would not commute [39][40][41].In this case, the method that is based on the theory of shadowing pseudo-orbits (approximate trajectories) in dynamical systems comes to the rescue [45].Using this approach, one can estimate the time-averaged sensitivities ∇ α R(α) and make an unambiguous conclusion regarding the system's sensitivity to its parameters.Note that the shadowing property (or pseudo-orbit property) implies that near an approximate system's trajectory, there exists the exact trajectory such that it lies uniformly close to a pseudo-trajectory.In this case, the sensitivity analysis problem is reduced to finding a pseudo-trajectory, assuming that the exact trajectory is a solution to the unperturbed problem.A detailed description of the "shadowing" method is rather cumbersome and can be found in [42,43]; therefore, it will not be discussed here.To illustrate the method and its accuracy and efficiency, we applied a conceptual model that was obtained by coupling two variants of the chaotic Lorenz model (L63) [44], the timescales of which differ by a factor ε [46]: where lowercase and capital letters denote fast and slow subsystems, respectively; σ, r and b are the parameters of the L63 model, which are proportional to the Prandtl number, Rayleigh number and dimension of the considered atmospheric layer, respectively (the L63 system exhibits chaotic behaviour for the following parameter values: σ = 10, r = 28 and b = 8/3); c is the coupling strength parameter for variables x and y; c z is the coupling strength parameter for variable z; k is the "uncentering" parameter; and a is the parameter that represents the amplitude scale factor.Without loss of generality, one can assume that and c = c z .
As an example, let us consider the sensitivity of the state variables of the model in Equation ( 14) to the parameter r, which plays a very important role in shaping the model dynamics.The system of sensitivity equations is derived by differentiating the system in Equation ( 14) with respect to parameter r: .
where the sensitivity functions are defined as follows: To find how the sensitivity functions evolve over time, we needed to solve the sensitivity system in Equation (15), along with the model in Equation ( 14).The calculated sensitivity functions S 3,r and S 6,r are shown in Figure 1 (note that we have chosen the functions S 3,r and S 6,r arbitrarily).These sensitivity functions exhibited fluctuating behaviour and their envelopes grew exponentially over time.Consequently, it is very difficult to draw a clear conclusion about the system's sensitivity to variations in the parameter r.Similar behaviour was inherent in the other sensitivity functions.To apply the shadowing sensitivity analysis method, first, we needed to calculate both the original and pseudo trajectories (Figure 2).Then, the sensitivity of the model state variables to the parameter r could be estimated (Table 1).To apply the shadowing sensitivity analysis method, first, we needed to calculate both the original and pseudo trajectories (Figure 2).Then, the sensitivity of the model state variables to the parameter r could be estimated (Table 1).
Table 1 shows that the variables z and Z were most sensitive to variations in the parameter r, while the sensitivities of the fast x and y and slow X and Y variables to the parameter r were much less than that of the variables z and Z.To apply the shadowing sensitivity analysis method, first, we needed to calculate both the original and pseudo trajectories (Figure 2).Then, the sensitivity of the model state variables to the parameter r could be estimated (Table 1).   1 shows that the variables z and Z were most sensitive to variations in the parameter r, while the sensitivities of the fast x and y and slow X and Y variables to the parameter r were much less than that of the variables z and Z.

Climate System Response to Small External Forcing Based on the Fluctuation Dissipation Relation
The fluctuation dissipation theorem (FDT), which is a powerful instrument in statistical physics, was first proposed in [47] and can be applied to estimate the climate system response to a small external forcing.Generally speaking, this theorem establishes the quantitative relationship between the fluctuations in a stochastic system at thermal equilibrium and the system's response to an applied external perturbation.On the basis of the  The fluctuation dissipation theorem (FDT), which is a powerful instrument in statistical physics, was first proposed in [47] and can be applied to estimate the climate system response to a small external forcing.Generally speaking, this theorem establishes the quantitative relationship between the fluctuations in a stochastic system at thermal equilibrium and the system's response to an applied external perturbation.On the basis of the FDT for a system with certain properties, it is possible to find a linear operator that characterizes the system's response to a small external forcing by using the covariance matrix of the unperturbed system.However, the climate system is a nonlinear dissipative dynamical system for which the standard assumptions of equilibrium statistical mechanics are not fulfilled.
Despite the strong nonlinearity of the climate system, linear approximation and the time-invariant hypothesis are widely used in climate research [7].The fundamental assumption of the linear approximation is that various external perturbations act independently and additively on the system's response.In this case, the response of the climate system, which is characterized by some climatic variable x, to the external forcing ∆ f can be represented as ∆x = S∆ f , where S is the sensitivity coefficient.For instance, to estimate the change in the equilibrium surface temperature ∆Θ s due to the radiative forcing generated by an increase in atmospheric CO 2 concentration, a simple relationship between the CO 2 concentration and the corresponding radiative forcing can be used [48]: ∆ f CO 2 (t) = 5.35 × ln(c(t)/c(0)), where c(t) is the CO 2 concentration (in parts per million volume) at time t and c(0) is the reference CO 2 concentration.Using the formula ∆Θ s = S∆ f CO 2 , where the sensitivity coefficient has a value of 0.8 K W −1 m 2 [48], for a doubling of pre-industrial CO 2 level, we obtain ∆Θ s ≈ 3 • C. Radiative forcing of 3.7 W m −2 for a doubling of CO 2 was used in this calculation.
The reaction of nonlinear systems to external forcing is fundamentally different from the response of linear systems.This difference is mostly due to the wider involvement of the system's inherent characteristics and irregularities and the various ways of taking them into consideration.Hence, the fluctuations of nonlinear systems are the integration of external perturbations and internal feedback mechanisms.By using the FDT, it becomes possible to clearly identify external perturbations and then separate them from the natural system's fluctuations.Let the evolution of some finite-dimensional deterministic dynamical system, the state of which is characterized by the state vector x, be represented by the following equation: The time-averaged state of this system (the system's "climate") is determined by the expression x = lim τ→∞ (1/τ) τ 0 x(t)dt.Let the system in Equation ( 17) be perturbed by some external forcing: Obviously, the time-averaged state of this system x differs from x and, thus, the difference δx = x − x depends on δ f : δx = U (δ f ), where U is some nonlinear function.For a sufficiently small (infinitesimal) δ f , one can expect that δx and δ f are linearly dependent.If a function U is differentiable with respect to δ f at some point, then U can be expanded into a Fourier series.Then, by omitting high-order terms and leaving only first-order linear terms, one can obtain an approximate linear relationship between δx and δ f : δx ≈ L(δ f ), where L = ∂U /∂δ f | δ f 0 is a linear response operator.Finding the response operator for dynamical systems that are used in climate research is a nontrivial problem since these systems are not hyperbolic, for which a linear relation between δx and δ f was proved [49].The use of the -regularization technique [50] makes it possible to solve the problem of finding the response operator for a climate dynamical system.For a "quasi-Gaussian" formulation of the FDT, the response operator L is of the form [51,52]: where C(τ) is the covariance matrix with time lag τ: C(τ) = x(t + τ)x T (t) .Thus, the relationship between the climate response δx and the external forcing δ f is as follows: As shown in several papers (e.g., [52]), Equation ( 20) is satisfied with sufficiently high accuracy for global atmospheric general circulation models.For illustrative purposes, consider the one-dimensional stochastic dynamical system that has only one state variable x and is generated by the Langevin equation: .
where the dot over the letter indicates a derivative with respect to time; α = 1/τ e , where τ e is the relaxation time for x; and f (t) is an external forcing with zero mean f (t) = 0 and a correlation function that is given by f (t) f (t ) = ϑδ(t − t ), where δ is the Dirac delta function.
Note that the Langevin equation (Equation ( 21)) is a very simple version of a stochastic energy balance climate model in which the external forcing is Gaussian white noise and the variable x is the global mean surface temperature anomaly.The autocorrelation function for the state variable x is x(t)x(t ) / x 2 = e −α|t−t | .The solution to Equation ( 21) is given by The ensemble average of f (t) is zero; therefore: .
Thus, for example, if the system in Equation ( 21), being initially in an equilibrium state, is affected by an impulse perturbation that is described by the delta function ϑδ(t − t ) at time t = t , which leads to an instant increase in temperature anomaly at t = t by ϑ, then the temperature anomaly begins to decay, exponentially decreasing back to zero.The rate of decay is determined by the relaxation time.
The mean system response to any infinitesimal external forcing F (t) can be calculated using Equation (20).In particular, for a staircase function F (t), i.e., a constant ∆ s that is activated at t = 0, the system response is given by

General Notes on Climate Modelling
The Earth's climate system, as we discussed above, is a multicomponent, multiscale and fully coupled system.Mathematical modelling is the most effective and, perhaps, the only powerful tool for studying the climate and projecting its future state.Mathematical modelling involves the development of extremely complex and sophisticated mathematical models in which simulated processes are represented as a dynamical system and described by a set of ordinary or partial differential equations.In general, these equations are based on the fundamental laws of physics, such as the conservation of mass, momentum and energy.
The basis for the development of climate models is the mathematical theory of climate, which can be conditionally divided into the following directions: - The mathematical formulation of the problem, i.e., the translation of the real-world problem into the form of mathematical equations to be solved.When developing climate models, it is usually assumed that the Navier-Stokes equations for a compressible fluid are valid for describing the atmospheric and ocean dynamics, and the equations of classical equilibrium thermodynamics are locally valid as well.De facto, in climate models, the Reynolds-averaged Navier-Stokes equations are used in which the effects of sub-grid processes with scales smaller than the averaging scale are expressed via the characteristics of large-scale processes.The main sub-grid processes include the following: -Atmospheric radiative transfer processes (transfer of shortwave (solar) and longwave (terrestrial) radiation through the atmosphere Among the variety of climate models that are used in climate research, four main types of models can be distinguished [53,54]: simple climate models (two-dimensional, one-dimensional or even zero-dimensional); intermediate complexity models; general circulation models of the atmosphere with a simplified description of the upper mixed ocean layer and sea ice; and fully coupled atmosphere-ocean general circulation models (AOGCMs), which are at the top of the climate model hierarchy.At first glance, the more complex the model that is used in climate research, the more accurate and reliable the results that are obtained.To a large extent, this is correct.Nevertheless, in many problems of climate theory, models of intermediate complexity and simple climate models are used with great success.Choosing a suitable model is a creative process.This choice is determined by the objectives of the study, taking into account numerous factors, including the available computational resources.Idealized, well-understood simple climate models that sit at the bottom of the climate model hierarchy are very useful tools for improving our understanding of the climate system, analyzing possible climate change due to current and future greenhouse gas emissions, estimating the effect of volcano activity on the cooling of our planet and similar problems regarding climate theory.Climate system models of intermediate complexity are mainly used to study the climate system over long timescales at reduced computational costs.This class of models is also useful in the exploration of individual climate-forming processes, feedbacks and the interactions between them.The area of the application of atmospheric general circulation models is close to the area of research of intermediate complexity climate models.Note that the gap between the models of these two classes is not so large since some models of intermediate complexity are derived directly from atmospheric general circulation models.At present, the use of simple models, models of intermediate complexity and atmospheric general circulation models in climate research is of an auxiliary nature.The most advanced tools that are currently available for simulating the current climate and the climate response to external natural and human-made perturbations are AOGCMs, which describe all components of the climate system and their interactions in sufficient detail.In recent decades, the development of AOGCMs has been marked by significant progress due to both advances in the study of the climate system itself and an increase in computing resources, which provide ever greater detail and completeness of the descriptions of climatically significant processes.Presentday developments in computer technology provide the ability to numerically integrate AOGCMs for many hundreds of years.

Governing Equations
AOGCMs are developed and constantly updated by different modelling groups around the world, primarily to understand how the Earth's climate has changed in the past and how it might change in the future.AOGCMs provide the ability to model physical, chemical and biological processes in the atmosphere, ocean, land and cryosphere in great detail, requiring the most powerful supercomputers to produce climate projections.Research groups that are involved in climate simulations coordinate their activity via the Coupled Model Intercomparison Projects (CMIP), and the obtained results are published in the Intergovernmental Panel on Climate Change (IPCC) assessment reports.In the preparation of the last, Sixth Assessment Report (AR6), published this year, the results that were obtained using 100 different AOGCMs that were developed by 46 different modelling groups were used.All AOGCMs are constructed from mathematical equations that describe the evolution of climate system components.These equations, as we mentioned earlier, express the basic laws of physics, such as conservation of momentum, mass, energy, water vapor and other gases and aerosols of other substances.Although the system of differential equations that is used to build climate models varies from model to model, any AOGCM includes horizontal motion equations, an equation for vertical velocity (or hydrostatic equation), a thermodynamic equation, a continuity equation and a moisture equation.To close this system of differential equations, some additional diagnostic equations are usually introduced, representing the relationships between climatic variables.One of the AOGCMs that has been involved in the CMIP over the years is the INMCM, which was developed in the Institute for Numerical Mathematics of the Russian Academy of Sciences to simulate climate processes on a wide range of temporal scales, from synoptic timescales to centennial climate change [55][56][57][58].The system of equations that is used in the INMCM to describe the atmospheric general circulation is written in spherical coordinates (λ, ϕ) as the horizontal coordinates (here λ and ϕ are the longitude and latitude, respectively), and the dimensionless pressure σ = p/p s (here p is the pressure and p s is the surface pressure) is used as the vertical coordinate.The system of INMCM equations comprises the following: (a) Two momentum equations: (b) The thermodynamic equation: (c) The equation for specific humidity that describes the hydrological cycle: (d) The continuity equation: (e) Hydrostatic equation: (f) The equation of state: Here u and v are the zonal and meridional wind components; a is the Earth's mean radius; . σ = dσ/dt is the analogue of the vertical velocity in the σ coordinate system; T is the temperature; Φ is the geopotential; Ω = ζ + f is the absolute vorticity, where ζ is the relative vorticity; f is the Coriolis parameter; K is the kinetic energy; R is the gas constant for dry air; q is the specific humidity; c p is the specific heat of dry air at constant pressure; T v is the virtual temperature; and ρ is the air density.The terms F u and F v describe the vertical friction and the horizontal diffusion processes; is the diabatic heating rate; F T and F q describe the vertical and horizontal diffusion of heat and water vapor, respectively; and C and E describe the source and sink processes for water vapor, respectively.The absolute vorticity and kinetic energy are given by The boundary conditions are periodic around a latitude circle.At the poles, the solution is assumed to be bounded.The vertical boundary conditions that are imposed on the bottom (σ = 1) and the top (σ = 0) of the atmosphere are .σ = 0, meaning the vanishing of vertical velocity.The equations that are used to simulate the ocean dynamics are written in a spherical coordinate system in the Boussinesq hydrostatic approximation with a rigid lid as follows: (a) Two momentum equations: (b) The thermodynamic equation: (c) The equation for the mass continuity of salinity: (d) The continuity equation: (e) Hydrostatic equation: (f) The equation of state: where Here σ = z/H, where z is the depth that is measured from the undisturbed ocean surface; H = H(λ, ϕ) is the topography of the ocean bottom; u, v and .σ are the components of the current velocity vector along the longitude, latitude and vertical coordinates, respectively; T, S, p and ρ are, respectively, the temperature, salinity, pressure and departure of density from the value of ρ 0 = 1200 g m −3 (in this case, ρ(t, p, S) is a known nonlinear function that defines the state of seawater, including its compression with increasing depth); the right-hand side terms F u , F v , F T and F S have the same meaning as the corresponding terms in atmospheric circulation equations and describe the turbulent dissipation and turbulent exchange of heat and salinity.
The system in Equations ( 32)-( 38) is written for a cylindrical non-simply-connected domain that is bounded by the unperturbed ocean surface (σ = 0) and its bottom (σ = 1).The boundary conditions for an analog of the vertical velocity at these horizons are of the form . σ = 0.The model equations are solved numerically on the high-performance supercomputer Lomonosov-2, which is deployed at the Moscow State University.
The INMCM is a participant of the CMIP, which provides a standard protocol for the assessment of AOGCMs.Validation of the INMCM against observational data and the reanalysis products shows that the model reproduces the current climate with reasonable accuracy, comparable to that of other similar coupled climate models [7,53,54].That is why this model has been successfully applied for many years for projecting climate using various scenarios of human-made impact on the climate system of our planet.In particular, the analysis of the results of numerous numerical experiments with the INMCM allows us to conclude that the model response operator in Equation ( 20) reproduces the response of the climate system to heating sources of various spatial structures with high accuracy.The accuracy of reproducing the spatial structure of the response and its amplitude is, on average, 72 and 85%, respectively, for the near-surface temperature.For some other components of the response (for example, for the stream function and temperature in the middle troposphere), the values of the correlations turn out to be much higher.
It is very important to emphasize that although climate models that are developed in different research centers predict human-made warming in the future quite consistently, the indicators that characterize the climate system response to a given radiative perturbation differ significantly from model to model (e.g., [59]).Accordingly, global and regional climate change projections have a rather high degree of uncertainty for given scenarios of greenhouse gas emissions.For instance, the range of the projected global mean surface temperature (GMST) for the late 21st century relative to the reference period of 1986-2005 is likely to be 2.6 to 4.8 • C under the RCP8.5 greenhouse gas concentration scenario.Note that RCP8.5 corresponds to the concentration pathway with the highest greenhouse gas emissions.The range of the projected GMST under other emission scenarios (RCP2.6,RCP4.5, RCP6.0) is also quite large (0.3 to 1.7 • C, 1.1 to 2.6 • C, and 1.4 to 3.1 • C, respectively) [7].This uncertainty arises mainly from inter-model differences in the description of radiative feedback mechanisms.Thus, the ongoing improvements of the AOGCMs and their ability to represent the current climate are an extremely pressing challenge for modern climate science.

Using Low-Order Simple Models to Study Climate Variability
Despite the fact that climate projections are based on ensembles of highly complex and computationally expensive AOGCMs, low-order simple models are still very popular in exploring various problems of climate theory.One of these very important problems is the study of natural climate variability against the background of human-induced global warming.Climate variability characterizes the variations of climate variables relative to the mean state.Climate variability arises from natural quasi-periodic processes that are inherent in the climate system, some of which are not fully understood (e.g., El Niño-Southern Oscillation).External forcing due to volcanism or fluctuations in solar activity can also cause climate variability.It is important that climate variability occurs over different temporal and spatial scales.From the viewpoint of socio-economic development, it is vital to quantify the variability of climate variables over a wide range of timescales, from years to decades.Climate variability that is estimated on the basis of AOGCMs shows a wide statistical range across models.For instance, the decadal variance of the GMST in CMIP5 models differs from each other by a factor of more than four [29].Some studies have shown that changes in the variance of the GMST that is simulated by CMIP5 models differ from each other significantly not only on decadal timescales but also on inter-annual and multi-decadal timescales (e.g., [60]).However, the reasons for such a large inter-model spread of the GMST variability are not entirely clear.Simple low-order climate models in a stochastic formulation are useful tools for exploring how various physical mechanisms (e.g., feedbacks) affect the formation of climate variability over a wide range of timescales.
Energy balance models (EBMs) are among the simple but quite effective models of the climate system.These models were introduced almost simultaneously about 50 years ago by M. Budyko [61] and W.D. Sellers [62].The purpose of EBMs is to obtain a better understanding of the Earth's climate system, why the current state of the climate system is what it is, how sensitive the climate is to external forcing, how feedbacks and climate system inertia affect the climate system, etc. Simple climate models, such as EBMs, are usually calibrated to global scale observations in order to be practically useful tools.De-spite the intensive development of computer capabilities, simple climate models have not lost their relevance.On the contrary, for decades, these models serve as a link between theoretical concepts and the results obtained from complex climate models.Over the years, EBMs have been successfully applied to explore various aspects of the general climate theory (e.g., [63][64][65][66][67][68][69][70][71]), including feedbacks (e.g., [72,73]), climate-vegetation interaction (e.g., [74,75]) and climate-biosphere interaction (e.g., [76][77][78][79]).Here, we used a two-box EBM representing the entire climate system with only two variables, the anomaly of the GMST T and the anomaly of global mean deep ocean temperature T D [80,81].Since the two-box EBM is linear, it allows for the construction of an analytical solution that describes the climate system's behaviour near an equilibrium state.The model employs a deep ocean box.Consequently, within the framework of this model, the climate system response to radiative forcing is characterized by the "fast" τ f ≈ 3.9 yr and "slow" τ s ≈ 242 yr relaxation times, allowing one to study the climate system evolution on timescales that exceed the interannual cycle [82].As shown in several studies (e.g., [19,37,65,[83][84][85][86][87][88][89][90]), simple linear systems, such as a two-box EBM with additive stochastic forcing that is parameterized by Gaussian white noise, were shown to be useful for exploring climate variability on timescales from years to decades.These models are able to explain much of the observed climate variability on timescales that vary from annual through to several decades.Since we were interested in climate variability on interannual-interdecadal timescales, we applied a two-box EBM.
The model equations were written as follows: where C is the effective heat capacity of the upper box, consisting of the atmosphere and the oceanic surface mixed layer and characterized the climate system inertia; λ, which is called the climate feedback parameter, is the proportionality coefficient between the radiative response and the GMST; γ is the deep ocean heat uptake parameter; C D is the effective heat capacity of the deep ocean; and F s is the stochastic forcing.The patterns of real radiative stochastic forcing are rather complex and de facto undefined.Hence, F s in climate simulations is usually considered additive Gaussian white noise, which is delta-correlated in time random process with a zero mean F s = 0 and correlation function R s (τ) = F s (t)F s (t + τ) = 2D s δ(τ), where D s is the diffusion coefficient and δ(τ) is the delta function.The diffusion coefficient D s is determined by the variance of the random white noise process σ 2 s and its correlation time τ s as D s = σ 2 s τ s [91].As a measure of climatic variability, we used the variance of the GMST anomaly σ 2 T , which can be obtained in two different ways: via the power spectrum density (PSD) of the GMST fluctuations, and via the Fokker-Planck (forward Kolmogorov) equation that is associated with Equations ( 39) and (40).Let us consider the first approach to find the variance σ 2 T .The PSD of stochastic forcing F s is, according to the Wiener-Khinchin theorem [92], the Fourier transform of its autocorrelation function: If stochastic forcing F s (t) passes through the EBM linear operator to get T(t), then the PSD of T(t) is as follows: where H(ω) is the Fourier transform of the EBM transfer function and ω is the real angular frequency.Using the notation q 2 s = 2D s = 2σ 2 s τ s and substituting Equation (41) into Equation (42), we obtain The two-box EBM transfer function is the Fourier transform of the differential Equations (38) and (39) with zero initial conditions under the assumption that F s (t) = δ(t), where δ(t) is the Dirac delta function: Then, the PSD of GMST fluctuations can be obtained using Equation ( 43): The integral of S T (ω) over positive frequencies is the variance of the GMST anomaly σ 2 T : The climate response to (human-made) radiative forcing has a wide range of timescales.For example, the response timescale of the atmosphere is just a few weeks.The response timescale of a system that involves the atmosphere and the ocean mixing layer is a few years, and the system that also involves the deep ocean is characterized by timescales from decades to centuries and even longer.Within the framework of a two-box EBM, the response timescale is determined by the parameter C, which characterizes the climate system inertia since its value depends on the depth of the upper ocean layer that is included in the model.In turn, the response strength depends on the feedback parameter λ.Since these mechanisms play a critical role in the formation of the climate system's response to radiative forcing, we focused on analyzing the climate sensitivity to climate system inertia and feedbacks.To estimate the effect of variations in the feedbacks and climate system inertia on the climate variability, we applied sensitivity analysis using absolute and relative sensitivity functions.It should be emphasized that since the signal power spectrum displays the signal's variance as a function of frequency, allowing one to estimate the distribution of the GMST variance over varying timescales, to assess the effect of climate feedbacks and inertia on its variability, we explored the sensitivity of the GMST power spectrum to variations in the model parameters λ and C. Sequentially differentiating Equation (45) with respect to the feedback parameter λ and climate inertia parameter C, we derived the absolute sensitivity functions (ASF) ψ λ and ψ C that characterize the effect of feedbacks and thermal inertia on the PSD of GMST fluctuations: The corresponding relative sensitivity functions (RSFs) ψ R λ and ψ R C are: The reference values of EBM parameters that were used in the calculation were taken from [93]: Values for these parameters were estimated based on a CMIP5 AOGCM analysis and are rounded multi-model means.The ranges of feedback λ and climate system inertia C, which characterize the uncertainty in the corresponding parameter, were as follows: 0.61 ≤ λ ≤ 1.70 W m −2 K −1 and 4.7 ≤ C ≤ 8.6 W yr m −2 K −1 [94].
For the sake of convenience, we introduce the dimensionless feedback coefficient f, which is usually used in the system's theory to characterize the feedbacks that are inherent in the system.This coefficient is given by f = 1 − λ/λ 0 , where λ 0 ≈ 3.3 W m −2 K −1 is the Planck feedback parameter, which characterizes the increase (decrease) in outgoing terrestrial radiation per unit of warming (cooling) of our planet.The range of the feedback coefficient f that corresponds to the range of feedback parameter λ is 0.49 ≤ f ≤ 0.82 with an ensemble average value of 0.66.Note that higher values of the feedback parameter λ correspond to smaller values of the feedback coefficient f .The parameter q s , which characterizes the intensity of stochastic radiative forcing, was estimated using the asymptotic relation q 2 s = σ 2 s τ yr , where σ 2 s is the radiative forcing variance averaged over the period τ yr [19].Since in the model, a one-year averaging period is used, τ yr = 1 yr, and time is also measured in years, then q 2 s = σ 2 s .The standard deviation range of the global mean radiative forcing is 0.16 − 0.40 W m 2 , with a multi-model average value of 0.26 W m 2 [94].
The ASFs to parameters λ and C that were calculated for different values of the feedback coefficient f are displayed in Figure 3.The coefficient f is varied over a wide range of values since feedback uncertainty is one of the main sources of uncertainty in the climate system's response to radiative forcing.Analysis of Figure 3 shows that the absolute value of ψ λ increased as the period of the GMST fluctuations increased.The rate of change in ψ λ was strongly dependent on the value of the feedback coefficient.The larger the coefficient f, the greater the rate of change in ψ λ .An important finding was that radiative feedbacks were strong across the entire frequency range, and continued to become stronger at very long timescales.The ASF  exhibited more complex behaviour, which was shaped like an inverted bell.The point at which the function  had an extremum corresponded to the frequency  * = ( + ) 2 ⁄ .For the reference parameter values, this frequency corresponded to an oscillation period of 25 years.It is interesting that the frequency  * defined the bending point at which the spectrum mode changed, i.e.,  * =  (Figure 4).At frequencies larger than  (i.e., for fluctuations with  >  ), the power spectrum corresponded to red noise, while at low frequencies ( <  ), the power spectrum reached the plateau region in which the PSD was constant.Let us briefly discuss how the theoretical power spectrum obtained with the two-box EBM fits the observational data.The PSD of the surface temperature fluctuations had a power-law dependence on frequency  () ∝  , where  was the scaling exponent (e.g., [95][96][97][98]).In the frequency range, corresponding to timescales from months to decades,  = 0.80 ± 0.15 [98].The scaling exponent value that was obtained from the CMIP5 model results [19] for frequencies that corresponded to timescales from years to decades was  = 0.82 ± 0.10 with 95% confidence.Dividing this power spectrum into a low-frequency range with decadal and interdecadal scales, and a The ASF ψ C exhibited more complex behaviour, which was shaped like an inverted bell.The point at which the function ψ C had an extremum corresponded to the frequency ν * = (λ + γ)/2πC.For the reference parameter values, this frequency corresponded to an oscillation period of 25 years.It is interesting that the frequency ν * defined the bending point at which the spectrum mode changed, i.e., ν * = ν c (Figure 4).At frequencies larger than ν c (i.e., for fluctuations with ν > ν c ), the power spectrum corresponded to red noise, while at low frequencies (ν < ν c ), the power spectrum reached the plateau region in which the PSD was constant.Let us briefly discuss how the theoretical power spectrum obtained with the two-box EBM fits the observational data.The PSD of the surface temperature fluctuations had a power-law dependence on frequency S TT (ν) ∝ ν −γ , where γ was the scaling exponent (e.g., [95][96][97][98]).In the frequency range, corresponding to timescales from months to decades, γ = 0.80 ± 0.15 [98].The scaling exponent value that was obtained from the CMIP5 model results [19] for frequencies that corresponded to timescales from years to decades was γ = 0.82 ± 0.10 with 95% confidence.Dividing this power spectrum into a low-frequency range with decadal and interdecadal scales, and a high-frequency range with an interannual scale, it was found that for these ranges, the values of the scaling exponent were 0.40 and 1.53, respectively.The scaling exponent for the theoretical power spectrum of the GMST fluctuations, as estimated via the two-box EBM with decadal and longer timescales, was 0.2-0.4,depending on the value of the feedback factor; meanwhile, the value of γ for the interannual timescales was closer to 2. The RSFs ψ R λ and ψ R C that allowed us to rank the parameters λ and C in accordance with their importance in the formation of the power spectrum are presented in Figure 5.This figure shows that the RSFs were intrinsically nonlinear, varying significantly with frequency, but their behaviour was monotonic.The intersection point of ψ R λ and ψ R C corresponded to the frequency of the bending point ν c , which was determined earlier.In the high-frequency range of the power spectrum, the influence of the climate system's inertia on the power spectrum was more significant than the influence of the feedbacks.In contrast, in the low-frequency range, the influence of the feedbacks on the power spectrum exceeded the influence of the climate system's inertia.The effect of the parameter uncertainty on the PSD was determined as follows: ( ) ≈ ±| |||, where the variation  characterizes the absolute uncertainty in The effect of the parameter uncertainty on the PSD was determined as follows: δ(S TT ) δα ≈ ±|ψ α ||δα|, where the variation δα characterizes the absolute uncertainty in the parameter α = (λ, C), and the variation δ(S TT ) δα characterized the uncertainty in the power spectrum caused by δα.The fractional (percentage) uncertainty in the PSD that was caused by uncertainties in the parameters λ and C was estimated as follows: [δ(S TT )/S TT ] α ≈ ±ψ R α δα/α 0 × 100%.The absolute and relative uncertainties in the power spectrum of the GMST fluctuations were estimated with the assumption that the parameter uncertainties corresponded to one standard deviation of the mean, i.e., we assumed that δλ = 0.31 W m 2 K −1 and δC = 1.1 W yr m −2 K −1 (Table 2).Table 2 shows that the behaviour of δ(S TT ) δα was similar to the behaviour of the corresponding ASF due to the linear relationship between them.The calculated fractional uncertainties presented in Table 2 clearly confirmed the results of assessing the relative role of the parameters λ and C in the formation of the power spectrum of the GMST fluctuations in different frequency ranges.The effect of the parameter uncertainty on the PSD was determined as follows: ( ) ≈ ±| |||, where the variation  characterizes the absolute uncertainty in the parameter  = (, ), and the variation ( ) characterized the uncertainty in the power spectrum caused by .The fractional (percentage) uncertainty in the PSD that was caused by uncertainties in the parameters  and  was estimated as follows: ( )  ⁄ ≈ ± |  ⁄ | × 100% .The absolute and relative uncertainties in the power spectrum of the GMST fluctuations were estimated with the assumption that the parameter uncertainties corresponded to one standard deviation of the mean, i.e., we assumed that  = 0.31 W m 2 K 1 and  = 1.1 W yr m 2 K 1 (Table 2).Table 2 shows that the behaviour of ( ) was similar to the behaviour of the corresponding ASF due to the linear relationship between them.The calculated fractional uncertainties presented in Table 2 clearly confirmed the results of assessing the relative role of the parameters λ and C in the formation of the power spectrum of the GMST fluctuations in different frequency ranges.The role of the anthropogenic effect was not considered in this stochastic EBM.Meanwhile, in [19,82,94], the deterministic version of a two-box EBM was applied to analyze the effect of human-induced radiative forcing on climate variability.Similarly, anthropogenic forcing can be taken into account in the stochastic model.However, we should keep in mind that the spectrum's shape is invariant under rescaling.
In the future, we intend to consider a more complex model, namely, a 1D EBM in stochastic formulation in which random forcing is a function of latitude.

Concluding Remarks
Mathematical theories, models and tools help us to solve some of the most challenging problems in the physical, biological and life sciences, including climate sciences.Mathematics and its methods and approaches have played a vital role in climate research for several decades.To date, mathematical modelling remains the only method and effective instrument that is used to predict the climate system's evolution under the influence of natural and anthropogenic perturbations.Current climate change is one of the greatest threats to humans in various aspects.Understanding the climate formation processes and projecting the future state of the climate system are the main issues in climate science, which cannot be solved without the use of mathematical methods for studying complex systems.In this study, we considered some mathematical methods and approaches, primarily mathematical modelling and sensitivity analysis, in the context of their application to climate system exploration.
The current climate is characterized not only by a long-term trend but also by variability against the background of human-induced global warming.The exploration of climate variability requires the development of stochastic climate models, which is not a trivial problem.Here we discussed the essential features of stochastic climate models and their application to climate variability exploration.As an illustrative example, we considered the application of a low-order energy balance model to study climate variability.

Figure 1 .
Figure 1.Sensitivity functions of the variables z and Z to the parameter r, which were calculated for strong coupling between fast and slow systems ( = 0.8).

Figure 1 .
Figure 1.Sensitivity functions of the variables z and Z to the parameter r, which were calculated for strong coupling between fast and slow systems (c = 0.8).

Figure 1 .
Figure 1.Sensitivity functions of the variables z and Z to the parameter r, which were calculated for strong coupling between fast and slow systems ( = 0.8).

Figure 2 .
Figure 2. Original (red) and pseudo (blue) trajectories for the fast z and slow Z variables, calculated for strong coupling between fast and slow systems ( = 0.8).

Figure 2 .
Figure 2. Original (red) and pseudo (blue) trajectories for the fast z and slow Z variables, calculated for strong coupling between fast and slow systems (c = 0.8).

Figure 3 .
Figure 3. Absolute sensitivity functions  (a) and  (b), which were calculated for different values of the feedback coefficient f.

Figure 3 .
Figure 3. Absolute sensitivity functions ψ λ (a) and ψ C (b), which were calculated for different values of the feedback coefficient f.

Figure 4 .
Figure 4. Power spectra of the global mean surface temperature fluctuations that were obtained from the two-box random energy balance model for different values of the feedback coefficient f.The dashed line shows the characteristic 1  ⁄ slope.

Figure 5 .
Figure 5. Relative sensitivity functions  and  , which were calculated for different values of the feedback coefficient .

Figure 4 .
Figure 4. Power spectra of the global mean surface temperature fluctuations that were obtained from the two-box random energy balance model for different values of the feedback coefficient f.The dashed line shows the characteristic 1/ν 2 slope.

Figure 4 .
Figure 4. Power spectra of the global mean surface temperature fluctuations that were obtained from the two-box random energy balance model for different values of the feedback coefficient f.The dashed line shows the characteristic 1  ⁄ slope.

Figure 5 .
Figure 5. Relative sensitivity functions  and  , which were calculated for different values of the feedback coefficient .

Figure 5 .
Figure 5. Relative sensitivity functions ψ R λ and ψ R C , which were calculated for different values of the feedback coefficient f .

Table 1 .
Sensitivity functions for fast and slow variables for strong ( = 0.8) and weak ( = 0.15) coupling between systems.

Table 1 .
Sensitivity functions for fast and slow variables for strong (c = 0.8) and weak (c = 0.15) coupling between systems.