## 1. Introduction

One of the biggest issues facing humanity today is the observed ongoing global climate change. Consequently, the prediction of future climate as well as changes in climate due to changes in natural processes and human-caused factors (e.g. greenhouse gas emissions) are issues that have deservedly received significant attention. The essential, powerful and effective methodology for solving this class of problems is computational simulation of the Earth’s climate system (ECS) and its components with the use of mathematical climate models that can range from relatively simple to fairly complex. Over the past several decades, the use of climate models as an aid in the understanding past, present and future climates has been substantial. The ECS is a natural, extremely complex, large-scale physical system that includes the atmosphere (the Earth’s gaseous envelope), hydrosphere (oceans, rivers and lakes), land surface, cryosphere (ice and snow) and finally the biosphere together with lots of natural and anthropogenic cycles (e.g. water cycle, carbon and nitrogen cycles). It is important that the constituent elements of the ECS are characterized by their own specific physical, dynamical and chemical properties [1–6]. Dynamics of the ECS has turbulent nature and displays wave-like fluctuations within a broad time-space spectrum and, therefore, are characterized by high non-linearity [3, 7, 8]. From mathematical viewpoint, the ECS is an extremely sophisticated, interactive, multi-scale, non-linear dynamical system. In this context, climate simulation represents one of the most complex and important applications of dynamical systems theory, its concepts and methods. The instantaneous state of certain characteristics of the ECS (such as temperature, humidity, atmospheric pressure, wind and precipitation) is referred simply to as the weather. Climate, meanwhile, represents the ‘average weather’ and is characterized by a statistical ensemble of states through which the ECS travels for decades (usually over ~30 years, according to the World Meteorological Organization’s definition).

State-of-the-art mathematical climate models used in variety of applications represent systems of multi-dimensional, non-linear differential equations in partial derivatives that are the mathematical statements of basic physical laws, primarily the conservation laws for momentum, mass and energy. Such models also include a variety of empirical and semi-empirical relationships and equations that are based on observations and experience rather than theories. Mathematical climate models are mostly deterministic with a large-phase space dimension, containing a vast number of various parameters. Equations that describe the evolution of the ECS and its processes are quite complicated. Therefore, in the majority of situations, we, unfortunately, cannot solve them analytically with an arbitrary set of initial conditions, even for very simple cases. We can only find an approximate solution using numerical methods such as, for example, Galerkin projection or finite-difference technique. Consequently, climate models have finite space and time resolutions. Due to the limited resolutions of climate models, many physical processes those are very important for climate dynamics cannot be adequately resolved by the model space-time grid and, therefore, should be parameterized, i.e. described parametrically. As a result, the number of model parameters increases significantly. A large number of them have various meanings and are not well-known *a priori* and, hence, are not well defined. Parameter errors and their variabilities in time and space generate parametric uncertainty in mathematical climate models and, undoubtedly, affect the output results. Assessment of the potential impact of variations in climate model parameters on the model behaviour represents essential element of model building and quality assurance. Sensitivity analysis in dynamical systems is a powerful tool that allows us to estimate the influence of model parameters and their variations on the results of computer simulations. It is important to ensure that some model parameters describe external forcing that can strongly influence the climate model behaviour. All human-caused impacts on the ECS (greenhouse gas emissions due to combustion of fossil fuels and industrial processes, aerosol emissions due to biomass burning, changes in albedo due to deforestation, soil tillage and land degradation) can be considered as small external perturbations that are described in climate models via parameter variations. Hence, both equilibrium and transient climate system sensitivity to external forcing can be examined within the framework of sensitivity analysis in dynamical system.

## 2. Elements of dynamical systems theory

Dynamical systems theory [3, 9–11] serves as a very powerful and reliable framework for modelling, studying and predicting the temporal-spatial behaviour of the ECS and its constituent elements. Generally, a certain abstract dynamical system represents a pair (*X*, *S*_{t}) , where *X* is the phase space of a system, and *S*_{t} : *X* → *X* is a family of evolution functions that is parameterized by a real variable *t* ∈ *T* . Commonly this variable performs the role of time. It is assumed that the phase space is a complete metric or Banach space that can be either finite- or infinite-dimensional. The set *γ*_{x} = {*x*(*t*) : *t* ∈ *T*} is called trajectory (or orbit), where *x*(*t*) is continuous function with values in *X* such that *S*_{τ}*x*(*t*) = *x*(*t* + *τ*) for all *τ* ∈ *T*_{+} and *t* ∈ *T*. In climate studies, the semi-dynamical systems are of prime interest. Semi-dynamical systems are dynamical systems whose evolution is considered only for *t* ≥ 0 . For semi-dynamical systems, a family of mappings *S*_{t}, *t* ≥ 0, forms a semi-group that satisfies the following conditions [12]:

Continuous-time dynamical system is commonly generated by the set of autonomous ordinary differential equations (ODEs) assuming that *t* ∈ ℝ_{+}

Here *x*_{0} ∈ *X* is a system state specified at *t* = 0 and *f* is continuous vector-valued function. The solution of a system (1) *x* = *x*(*x*, *t*) is determined for all *t* ≥ 0 and can be represented as *x*(*t*) = *x*(*x*_{0}, *t*) = *S*_{t}*x*_{0}. However, in climate simulations, we ordinarily deal with discrete dynamical systems. To convert the set of infinite-dimensional differential equations (1) into finite-dimensional discrete form, a few methods can be applied (e.g. finite-difference approximation or spectral projection technique). Consequently, instead of continuous-time dynamical system (1) we can obtain a discrete in time and space dynamical systems that can be solved numerically for given initial conditions:

Given the system state *x*_{0} at the initial time *t* = 0, we define the trajectory of *x*_{0} under *f* to be the sequence of points *x*_{k} = *f*^{ k}(*x*_{0}), where *f*^{ k} denotes the *k*-fold composition of *f* with itself. Note that *f*^{ 0}(*x*)≡*x* . Eq. (2) uniquely specifies the trajectory of discrete dynamical system if the map *f* : *X* → *X* and the initial state *x*_{0} are given. The state of dynamical system *x*_{k} at time *t*_{k} is defined by the system state *x*_{k − 1} at the previous time *t*_{k − 1} as *f*(*x*_{k − 1}).

It is important to highlight that climate dynamical systems possess a number of generic properties [3, 9, 12–14]:

First, these systems are non-linear and dissipative. This means that the divergence of corresponding vector field ∇ ⋅ *f*(*x*(*t*)) is strictly negative and the system’s phase volume contracts.

Secondly, from a certain moment of time *t**, the norm of the solution for any initial conditions stays bounded: ‖*x*(*t*)‖ < *V*_{0} at *t* > *t**, where *V*_{0} is the so-called absorbing set in the system phase space. All trajectories of climate dynamical system will ultimately enter to the ball of radius *V*_{0}. This property guaranties the existence inside *V*_{0} of finite-dimensional invariant compact attracting set, which is called the attractor. Attractor is, therefore, a set towards which a system tends to evolve for a wide variety of initial conditions of the system. If starting states of a system are chosen on the attractor, then the corresponding orbits will remain on the attractor. All other trajectories will be attracted to this set fairly fast.

Thirdly, trajectories of climate dynamical systems are generally unstable, exhibiting Lyapunov instability. This means that the *n*_{λ}-dimensional part of phase volume of a system increases along certain directions correspond to *n*_{λ} positive Lyapunov exponents (note that *n*_{λ} < *n*).

Finally, a certain unstable trajectory enclosed in a bounded phase volume (attractor) generates a deterministic dynamical chaos, which means that over time, under certain conditions, the behaviour of simulated system begins to resemble a random process, despite the fact that the system is defined by deterministic laws and described by deterministic equations. This phenomenon of deterministic chaos was first uncovered by E. Lorenz as he observed the sensitive dependence of atmospheric convection model output on initial conditions [15]. Consequently, all trajectories that start close enough will diverge from one another, but will never depart from the attractor. The rate of separation of infinitesimally close trajectories is characterized by the positive Lyapunov exponents. The number of directions along which the orbit is unstable is defined by the number of positive Lyapunov exponents *n*_{λ}.

In general, the dynamics of climate system can be conditionally divided into two phases that are correspondingly the motion towards the attractor and the motion along the attractor. When system orbit travels toward the attractor, the system phase volume contracts to the finite-dimensional volume of the attractor. In many theoretical and practical applications, the ECS evolution is considered on its attractor assuming that the system possesses the properties of ergodicity. Then, statistical characteristics of a system (e.g. first *n*_{λ} directions, the dimension of the attractor cannot be less than the dimension of phase space, and the attractor is nested inside a bounded absorbing set. Consequently, the attractor represents a fractal set of dimension *n*_{A} nested in the absorbing ball, where *n*_{λ} ≤ *n*_{A} < *n*. The fractal dimension of attractors of dissipative dynamical systems can be determined by the Kaplan-Yorke conjecture [16]. The so-called Kaplan-Yorke dimension is defined as follows:

where *j* is the maximum integer such that the sum of the *j* largest exponents is still non-negative, i.e. *e*-folding time. Note that for chaotic dynamics the characteristic time is finite, while for regular motion it is infinite. As an example, let us consider the barotropic vorticity equation on the rotating Earth:

where *ψ* is dimensionless geostrophic stream function, *l* is the Coriolis parameter, *h* is an orography, *ϕ* is an external forcing, *α* and *ν* are friction coefficients, *J*( , ) is a Jacobian and ∇ is the Laplacian. Barotropic vorticity equation has a finite-dimensional global attractor. The formula for the estimated fractal dimension of this attractor is [12]:

where *Gr* = ‖*ϕ*‖_{2}/(|*λ*_{1}|*ν*^{2}) is the generalized Grashof number, *λ*_{1} is the largest (*λ*_{1} < 0) eigenvalue of the Laplacian of a sphere.

It is necessary to underline that one of the most essential properties of climate dynamical systems is the stability. Generally, the stability of dynamical systems refers to their response to external forcing and other inputs. Dynamical system is considered to be stable if it remains in an ‘invariable’ state as long as there is no external forcing, and if the system returns to an ‘invariable’ state when the external forcing is eliminated. The climate system stability may be determined in terms of bounded inputs. In other words, the climate system is stable if every bounded input generates a bounded output. Climate system modelling is concerned not only with the stability of a system but also the degree of its stability. Climate models have a large number of various parameters that may affect the results of computational simulations. Assessment of the climate system response to variations in the model parameters is extremely important and is among the main problems with regards to system stability.

Other important property of climate dynamical system is its predictability. Generally speaking, predictability is the degree to which an accurate prediction of a system’s state can be made either qualitatively or quantitatively. It is important to note that pretty much the same models are applied for both weather prediction and climate simulation. Taking into account the chaotic nature of the atmosphere-ocean system, two types of problems associated with predictability can be defined. Predictability of the first kind is associated with numerical weather prediction, which is a Cauchy (initial value) problem. Numerical weather prediction aims to predict, as precisely as possible, the future state of the atmosphere. This problem requires an accurate knowledge of the current state of the atmosphere-ocean system. However, the evolution of this system is highly sensitive to small errors in the initial conditions [15]. Due to the intrinsic limits with regards to the predictability of atmospheric processes, inaccuracies and missing information in the inputs, as well as inadequacies of forecasting models, detailed and useful weather forecasts are limited to about 2 weeks [18].

By contrast, predictability of the second kind focuses on the prediction of the statistical properties of a system with respect to external forcing. Climate simulation belongs to this class of problems and actually represents a boundary value problem that focuses on much longer time scales than numerical weather prediction (typically several months or even years). Boundary conditions (e.g. the energy that reaches the Earth from the Sun, the energy that goes from the Earth back out to space and so on) constrain climate over a long period of time. If these boundary conditions are imposed correctly then we can simulate the Earth’s climate in the long run, without paying attention to what the initial conditions are. To control initial instability caused by initial conditions, the spin-up period is usually eliminated from the climate simulation, prior to analyzing the results. Thus, predictability of climate systems involves the study of stability of climate model attractors with respect to external forcing. Predictability of dissipative dynamical systems applied to climate simulations deteriorates with time. To quantify predictability, we can use the rate of divergence of system orbits in phase space (i.e. Kolmogorov-Sinai entropy, Lyapunov exponents).

## 3. A formal dynamic model of climate system

Let us consider the ECS in a bounded spatial-temporal domain Ω_{t} = Ω × [0, *τ*]. Let us denote by *φ* ∈ *Q*(Ω_{t}) the state vector that characterizes the ECS in the domain Ω_{t}. Note that *Q*(Ω_{t}) is the infinite real space of sufficiently smooth state functions that satisfy certain boundary conditions at the boundary ∂Ω of the spatial domain Ω, which is usually the Earth’s sphere. Let be the vector of spatial variables and *t* ∈ [0, *τ*] be the time. Mathematically, the temporal evolution of the ECS in the domain Ω_{t} is expressed via the set of partial differential equations, which reflect the specific dynamical, physical and other properties of the ECS:

where is a non-linear, multi-dimensional differential operator that describes the dynamics, dissipation and external forcing of the system, *λ* ∈ *G*(Ω_{t}) is the parameter vector, *G*(Ω_{t}) is the domain of admissible values of the parameters and *φ*_{0} is a given vector valued function (the initial state estimate). The model state vector *φ* includes temperature, pressure, density, humidity, wind velocity and other physical variables. Note that the system (4) characterizes a *continuous medium for which the state vector φ is infinite-dimensional*:

*φ*∈ Φ, where Φ is the infinite-dimensional Hilbert space. The vector of parameters, in its turn, contains any input of the system (4) such as classical parameters, initial and boundary conditions and so on. Essentially, the solution of such complex infinite-dimensional system cannot be found analytically. In order to obtain numerical solution, the original system of infinite-dimensional equations should be transformed, using an appropriate method, into a system with a finite number of degrees of freedom. For example, Eq. (4) can be projected onto the sub-space spanned by the orthogonal base

*φ*can be introduced in the form of normally convergent series:

Substituting (5) into (4) and then applying the Galerkin method, we can obtain, instead of infinite-dimensional distributed parameter system (4), the finite-dimensional lumped system that is formally described by the following set of ordinary differential equations:

Here is the state vector with dimension *n* representing a set of spectral coefficients, *F* is the Galerkin projection of the operator on the base Ψ and *m*-dimensional parameter vector.

It is important to note that space-time spectrum of processes occurring in the climate system is extremely broad. Consequently, the state-of-the-art mathematical climate models due to their spatial-temporal limited resolutions are unable to simulate correctly all of these processes. Physical processes that are too small-scale to be explicitly represented in the model due to its discrete spatial-temporal structure are parameterized, i.e. replaced by simplified parametric schemes generating additional model parameters. Let *τ*_{c} be a characteristic timescale of a certain physical process. Physical processes with timescales smaller than *τ*_{c} (the so-called sub-grid-scale processes) should be parameterized in climate models. Examples of these processes include radiative transfer (short-wave solar radiation and outgoing long-wave radiation emitted by the Earth), cloud formation processes, microphysical processes within clouds that lead to the formation of precipitation, the deep convection in the tropics, land-surface processes, photochemical processes, carbon cycle, etc.

Using finite difference method to approximate time derivatives in Eq. (6), we can obtain the generalized numerical climate model that can be used for computer simulations:

where *m*_{0,k} is non-linear operator that indirectly contains model parameters and propagates the state vector from time *t*_{0} (the initial conditions) to time *t*_{k}, and *K* is the number of time steps. Generally, all climate models may be arranged in several classes, based on various principles, e.g. the complexity of the models or the description and representation of physical processes [19]. However, there is no best or general-purpose climate model. Each particular model is characterized by inherent properties, and has specific advantages and disadvantages. Selecting the ‘best’ model depends on many various factors, including the objectives of simulation and what performance measures are used.

## 4. Forward and adjoint sensitivity analysis

To estimate the impact of model parameter variations on the model performance and state variables, one can use a sensitivity function (or coefficient) that is the partial derivative of a given element of state vector *x*_{i} with respect to a certain model parameter *α*_{j} [20, 21]:

where *δα*_{j} is the infinitesimal perturbation of parameter *α*_{j} around some fixed point *x*(*α*^{0} + *δα*) around *x*(*α*^{0}) by Taylor expansion, one can obtain the following linear equation:

where is a sensitivity matrix. Let us rewrite the vector equation (6) in component form:

Differentiating Eq. (8) with respect to *α*_{j}, we obtain the set of non-homogeneous ODEs, the so-called sensitivity equations:

Sensitivity equations describe the evolution of sensitivity functions along a given trajectory, and, therefore, allow tracing the sensitivity dynamics in time. A system of sensitivity equations (9) can be expressed in a matrix form (see below). Thus, to calculate sensitivity functions with respect to parameter *α*_{j} one should be able to solve the following set of differential equations with given initial conditions:

where *S*_{j} = (∂*x*/∂*α*_{j}) = (*S*_{1j}, ⋯, *S*_{nj})^{T} is the sensitivity vector with respect to parameter *α*_{j}, *D*_{j} = (∂*F*_{1}/∂*α*_{j}, ∂*F*_{2}/∂*α*_{j}, …, ∂*F*_{n}/∂*α*_{j}, )^{T}. Once we have solved Eq. (10), it is possible to analyse the sensitivity of system (8) with respect to the parameter *α*_{j}. Since the model parameter vector *α* has *m* components, to evaluate the model response to variations in the parameter vector *δα*, the set of Eq. (10) must be solved *m* times. Therefore, this approach is acceptable for low-order models. The use of sensitivity functions requites the differentiation of model equations with respect to parameters. However, this is not always possible. Fairly often in sensitivity analysis of complex dynamical systems, the model response to variations in its parameters represents a generic response function that characterizes the dynamical system [22, 23]:

where Φ is a non-linear function of state variables *x* and model parameters *α*. The gradient of functional *R* with respect to vector of parameters *α* around the unperturbed parameter vector *α*^{0} and corresponding unperturbed state vector *x*^{0}:

quantifies the influence of parameters on the model performance. In particular, the effect of the *j*th parameter can be estimated as follows:

where *δα*_{j} is the variation in parameter

The accuracy of sensitivity estimates strongly depends on the choice of perturbation *δα*_{j}. By introducing the Gâteaux differential, the sensitivity analysis problem can be considered in the differential formulation eliminating the need to set the value of *δα*_{j} [22, 23]. The Gâteaux differential for the response function (11) has the following form:

Here, *δx* is the state vector variation due to the variation in the parameter vector in the direction *δα*. Linearizing the non-linear model (Eq. (8)) around an unperturbed trajectory *x*^{0}(*t*), we obtain the following system of variational equations, the so-called tangent linear model, for calculating *δx*:

Then, using Eq. (12), we can calculate the variation in the response functional *δR*. Taking into account that where is a dot-product, the model sensitivity to variations in the parameters is estimated via the gradient of the response functional ∇_{α}*R*. This approach is convenient, however, computationally very expensive since climate models involve a large number of parameters. Adjoint approach allows the calculation of sensitivities within a single numerical experiment. Using adjoint equations, one can calculate the sensitivity gradient ∇_{α}*R* as follows [22]:

where the vector function *x*^{∗} is the solution of adjoint model, which is numerically integrated in the inverse time direction:

One of the key problems in climate simulation is the identification of parameters in mathematical models that describe the climate system. This problem is quite difficult due to both the huge number of state variables and parameters, and the argument that the governing finite-difference equations are non-linear grid functions of these states and parameters. Using the adjoint approach, we can solve the identification parameter problem if the observations are available. This problem is mathematically formulated as an optimal control problem in which model parameters play the role of control variables, and model equations are considered as constraints. Let *y*^{obs} be the set of observations and *H* be observation operator mapping from solution space of model to observation space. Therefore, *y*^{obs} = *H*(*x*) + *ε*^{obs}, where *ε*^{obs} is the vector of observation errors. It is usually assumed that these errors are serially uncorrelated and normally distributed with known covariance matrix *W*. The parameter identification problem seeks to minimize, with respect to *α*, a certain objective function
expressing the ‘distance’ between observations and corresponding model state using the model equations as constraints:

where *α** is a specified parameter vector. The objective function is written as

For illustrative purposes, let us consider the following example. Let *x*_{i}(*t*_{k}, *α*) be, respectively, the observation and model prediction of the *i*th component of state vector at time *t*_{k}, *H* the identity operator. Then the objective function can be written as

where are a *weighted coefficient reflecting the accuracy of* observations (in our case, *α*_{j} is defined as

where

## 5. Application of conventional methods of sensitivity analysis

Both forward and adjoint methods allow us to analyse transient and equilibrium sensitivities of dynamical systems. The exploration of sensitivity of complex dynamical systems requires considerable computational resources. For simple enough low-order models, the computational cost is minor and, for that reason, models of this class are widely used as simple test instruments to emulate more complex systems. We will illustrate the above-described conventional methods of sensitivity analysis based on a coupled non-linear dynamical system, which is composed of fast (the ‘atmosphere’) and slow (the ‘ocean’) versions of the well-known Lorenz [15] model (L63). This model allows us to mimic the atmosphere-ocean system and therefore serves as a key element of a theoretical and computational framework for the study of various aspects of sensitivity analysis. Recall that under certain conditions the Lorenz model exhibits a chaotic behaviour. As mentioned above, the system is obtained by coupling of two versions of the original Lorenz model with distinct timescales that differ by factor *ε* [24, 25]:

where the fast and slow sub-systems are represented, respectively, by lower case and capital letters. The following notations are used in Eqs. (16)–(21): *σ*, *r* and *b* are the parameters of the original L63 model, *a* is a parameter representing the amplitude scale factor, *k* is ‘uncentring’ parameter, *c* is a coupling strength parameter for *x* − *X* and *y* − *Y* variables, and *c*_{z} is a coupling strength parameter for *z* − *Z* variables. One can assume that *a* = 1, *k* = 0 and *c* = *c*_{z} without loss of generality. Therefore, the vector of state variables of coupled model (18) and (19) is *x* = (*x*, *y*, *z*, *X*, *Y*, *Z*)^{T} and the vector of model parameters is *α* = (*σ*, *r*, *b*, *c*, *ε*)^{T}. In the operator form, the set of Eqs. (16)–(21) can be rewritten as follows:

where the non-linear uncoupled operator *L* and linear coupled operator *Q* are represented by the following matrices:

The unperturbed parameter values are selected to be as follows:

Chosen values of *σ*, *r* and *b* correspond to chaotic behaviour of the L63 model. For *σ* = 10 and *b* = 8/3, the critical value of parameter *r* is 24.74, which means that any value of *r* larger than 24.74 induces deterministic chaos [15]. The parameter *ε* = 0.1 indicates that the slow system is 10 times slower than the fast system. Basic dynamical, correlation and spectral properties of this system, and the influence of the coupling strength on power spectrum densities, spectrograms and autocorrelation functions were explored in detail in [26]. Here, we will briefly mention the main features of the system (19).

The coupling strength parameter *c* plays a very important role in qualitative changes in the system dynamics since this parameter controls the interactions between fast and slow sub-systems. Qualitative changes in the dynamical properties of a system can be revealed by determining and analyzing the system’s spectrum of Lyapunov exponents that characterize the average rate of exponential divergence (or convergence) of nearby trajectories in the phase space. In the analysis of coupled dynamical systems, we are dealing with conditional Lyapunov exponents that are normally used to characterize the synchronization with coupled systems. The dynamical system (21) has six distinct Lyapunov exponents. If the parameter *c* tends to zero, then this system has two positive, two zero and two negative Lyapunov exponents. Let us examine the impact of coupling strength parameter on the two largest Lyapunov exponents (**Figure 1**). These exponents, being initially positive, monotonically decrease when the coupling strength parameter increases approaching the *x*-axis at about *c* ≈ 0.8 and become negative at *c* ≈ 0.95. Therefore, phase-synchronous regime for fast and slow sub-systems is observed when *c* > 0.95 [27]. However, when *c* > 1.0 , all six Lyapunov exponents become negative causing a limit cycle dynamics.

Apart from Lyapunov exponents, autocorrelation functions enable one to distinguish between regular and chaotic processes and to detect transition from order to chaos. In particular, for chaotic motions, autocorrelation functions decrease in time, in many cases exponentially, while for regular motions, autocorrelation functions are unchanged or oscillating. In general, however, the behaviour of autocorrelation functions for chaotic motions is frequently very complicated and depends on many factors (e.g. [28]). Knowing autocorrelation functions, one can determine a typical timescale (typical time memory) of the process [29]. Moreover, if autocorrelation functions are positive, the dynamical system may have the persistence property (an intention of a system to remain in the similar state from one time moment to the following). For a given discrete dynamic variable

Signal analysis commonly uses the normalized autocorrelation function (ACF), defined as . Results of numerical experiments show that for relatively small parameter *c* (*c* < 0.4), the ACFs and their envelopes for all variables decrease fairly rapidly to zero, consistently with the chaotic behaviour of the coupled system. However, as expected, the rates of decay of the ACFs of the slow variables are less than that of the fast variables. For coupling strength parameter on the interval 0.4 < *c* < 0.6, the ACFs of the fast variables become smooth and converge to zero. As the parameter *c* increases, the ACFs become periodic and their envelopes decay slowly with time, indicating transition to regularity. For *c* > 0.8, calculated ACFs show periodic signal components. In order to explore the sensitivity of system (19) with respect to coupling strength parameter, let us introduce the following sensitivity functions:

The corresponding sensitivity equations can be written as

Sensitivity functions can be introduced for any particular model parameter. Since the parameter vector *α* consists of five components, five sets of sensitivity equations can be derived from the model (Eq. (19)). The dynamics of sensitivity functions can be traced by solving the sensitivity equations along with the non-linear model. Sensitivity functions, calculated on the time interval [0, 20], are shown in **Figure 2**. Envelopes of these functions grow over time while sensitivity functions themselves oscillate. Sensitivity function is a measure of the change in state variable due to the variation in the estimated parameter. Unfortunately, obtained sensitivity functions are inherently uninformative and misleading. We cannot make a clear conclusion from them about system sensitivity to variations in the parameter *c*. Similar results were obtained when we considered the influence of variations in the parameter *r* on the system dynamics. This parameter plays an important role in the formation of system’s dynamical structure and transition to chaos. Let us define the following sensitivity functions:

The associated system of sensitivity equations can be written as

Envelopes of calculated sensitivity functions grow over time and sensitivity functions demonstrate the oscillating behaviour (**Figure 3**). Obtained functions are also uninformative and inconclusive. Thus, using conventional methods of sensitivity analysis can be questionable in terms of interpretation of the obtained results for chaotic dynamics. As discussed in [21], general solutions of sensitivity equations for oscillatory non-linear dynamical systems grow unbounded as time tends to infinity; therefore, sensitivity functions calculated by conventional approaches have a high degree of uncertainty, quickly becoming uninformative and inconclusive as time increases. In this regard, in climate simulation, the average values of sensitivity functions over a certain period of time can be considered as one of the most important measures of sensitivity, where *R* is a generic response functional (Eq. (11)). However, the gradient cannot be correctly estimated within the framework of conventional methods of sensitivity analysis since for chaotic systems it is observed [30–32] that

does not possess uniform convergence and two limits (*τ* → ∞ and *δα* → 0) would not commute.

## 6. Sensitivity analysis based on shadowing property

The new approach of sensitivity analysis known as shadowing method was introduced in [31] to analyse the sensitivity of highly non-linear and/or chaotic dynamical systems with respect to variations in their parameters. This method is theoretically based on the pseudo-orbit tracing (or shadowing) property for discrete and continuous dynamical systems [33, 34]. Using shadowing method, one can properly estimate the time-average sensitivities making the conclusion on the system sensitivity to variations in its parameters. The pseudo-orbit tracing property means that around an approximate (or pseudo) trajectory of dynamical system under consideration, there exists the exact trajectory, such that it locates evenly close to a pseudo-trajectory. Pseudo-trajectories arise by virtue of various errors of computer simulation (e.g. round-off errors and numerical technique errors). Thus, by making numerical simulations, we actually cannot obtain an exact trajectory of a system, but only an approximate trajectory known as a pseudo-trajectory. The exploration of shadowing property in dynamical systems was originated by Anosov [35] and Boven [36]. To date, the shadowing theory is well-established for the so-called hyperbolic dynamics that distinguished by the existence of both expanding and contracting lines for derivatives. Let us make some basic notes on the shadowing theory.

Let *f* : *X* → *X* be a homeomorphism of a compact metric space (*X*, dist). We will consider the discrete dynamical system *f*:

Let us denote by , the trajectory of a point *x* ∈ *X* with respect to *f*, i.e. Set any metric (distance function) *dist*(⋅,⋅) for *X*. It is said that a set of points is called a *d*-pseudo-trajectory (*d* > 0) of *f* if *dist*(*x*_{k + 1}, *f*(*x*_{k})) < *d* for *k* ∈ *ℤ*. Note that a distance function *dist*(⋅,⋅) defines an interval between each pair of geometric objects inside the brackets. It is said that the dynamical system *f* possesses the shadowing property if for every *ε* > 0 , there is *d* > 0 such that for every *d*-pseudo-trajectory *ζ* , there exists *y* ∈ *X* satisfying *dist*(*f*^{k}(*y*), *x*_{k}) < *ε* for all *k* ∈ *ℤ*. According to the so-called discrete shadowing lemma [33], for any *ε* > 0 , there exists *d* > 0, such that any *d*-pseudo-trajectory can be *ε*-shadowed. If any dynamical system has the shadowing property, its trajectories calculated numerically reflect the reality.

The shadowing lemma for continuous dynamical systems (flows) is more sophisticated than for discrete systems since in this instance re-parameterization of shadowing trajectories is required because for flows, close points of the true trajectory and the pseudo orbit do not correspond to the same time moments [33, 34]. Let *X* be the phase space of continuous dynamical system *d*-pseudo-orbit of *f*^{t} if the inequality *dist*(*f*^{t}(*t*, *g*(*τ*)), *g*(*t* + *τ*)) < *d* holds for any *t* ∈ [−1, 1]. A re-parameterization is actually a monotonically increasing homeomorphism *h* of the line ℝ such that *h*(0) = 0. The set of re-parameterizations *h* denoted by *Rep*(*ε*), where *ε* > 0, is defined as [33]:

To illustrate the applicability of this method, let us consider a continuous one-parameter generic dynamical system

The shadowing sensitivity analysis method is based on the ‘continuous’ shadowing lemma and the following assumptions: (a) the dynamical system is ergodic and (b) model state variables are considered over long time interval *t* ∈ [0, *T*], where *T* → ∞, and an averaged performance measure *R*(*α*). For example, the arbitrary trajectory *x*(*t*) can be chosen as a solution of the model equations, such that it is located nearby a certain reference trajectory *x*_{r}(*t*). According to the shadowing lemma, the closest orbit *x*(*t*) to *x*_{r}(*t*) satisfies the following constrained minimization problem [37]:

where *η* is the parameter that provides the same order of magnitude of the two terms in the integrand and *τ*(*t*) is a time transformation. The second term in the integrand describes re-parameterization. The problem (22) is called the non-linear least square shadowing (LSS) problem, and its solution denoted by *x*_{r}(*t*). The performance measure 〈*R*(*α*)〉 averaged over the time interval *t* ∈ [0, *T*] can be then approximated as:

since *α*. The desired sensitivity estimate

The solutions of this problem *S*(*t*) and *μ*(*t*) relate to the solutions of the non-linear LSS problem (22) as follows:

The time-dependent parameter *μ* is called a time-dilation variable, and it corresponds to the time transformation from the shadowing lemma. Using *S* and *μ*, we can estimate the desired sensitivity (the derivative of the linearized performance measure (22) with respect to the parameter *α*):

Several numerical algorithms can be used to solve the linearized LSS problem (24). One such method is based on variational calculus, which is used to derive optimality conditions representing a system of linear differential equations that are then discretized and solved numerically to calculate variables *S* and *μ* [37].

Here, we consider two sets of numerical experiments: weak coupling (*c* = 0.015) and strong coupling (*c* = 0.8) between fast and slow systems. As an example, the original and pseudo-trajectories for the fast *z* and slow *Z* variables are shown in **Figures 4** (weak coupling) and **5** (strong coupling). Pseudo-trajectories were calculated using the LSS method by adding a small perturbation to the parameter *r*. The deviation that is a difference between the ‘true’ and the approximate variables of the slow and fast sub-systems for weak coupling are illustrated in **Figure 6**. Obtained pseudo-trajectories lie very close to the associated true trajectories, proving the existence of shadowing property. Therefore, one can analyse the sensitivity of coupled dynamical system (20) with respect to parameters by averaging the computed sensitivity functions along the trajectory.

c
| ∂Z/∂r | ∂Y/∂r | ∂X/∂r | ∂z/∂r | ∂y/∂r | ∂x/∂r |
---|---|---|---|---|---|---|

1.0 | 1.10 | 0.05 | 0.01 | 1.08 | 0.04 | 0.03 |

0.8 | 0.69 | 0.08 | 0.03 | 1.02 | 0.07 | 0.07 |

0.4 | 0.95 | 0.03 | −0.01 | 1.03 | 0.09 | 0.09 |

0.15 | 0.91 | −0.08 | −0.09 | 1.01 | −0.01 | −0.01 |

10^{−4} | 1.04 | −0.02 | −0.03 | 1.02 | −0.01 | −0.01 |

The strong coupling does not introduce significant qualitative and quantitative changes in the behaviour of pseudo-trajectories with respect to the true orbits. Sensitivity estimates with respect to the parameter *r* calculated over the time interval [0, 20] for different values of coupling strength parameter are shown in **Table 1**. The most sensitive variables are *z* and *Z*. The sensitivity of variables *x*, *y*, *X* and *Y* with respect to *r* is significantly less than that of variables *z* and *Z*.

The use of shadowing method of sensitivity analysis allows the calculation of average sensitivity functions that can be easily interpreted. Let us recall that sensitivity functions obtained via conventional methods are uninformative and inconclusive because their envelopes grow exponentially over time and functions themselves oscillate.

## 7. Fluctuation dissipation relation in climate sensitivity

To estimate the ensemble-averaged response of the climate system to small external forcing, Leith [38] has proposed using the fluctuation-dissipation theorem (FDT) from statistical physics that quantifies the relation between the fluctuations in a system at thermal equilibrium and the response of the system to applied external perturbations. The FDT was initially obtained by Nyquist [39] in terms of the relation between fluctuations in voltage appearing across a resistor and its impedance. Afterwards, the FDT was formulated and proved for systems modelled by a Langevin equation [40] and for general non-linear dynamical systems in thermodynamic equilibrium [41, 42]. According to the FDT, under certain assumptions, the response of stochastic dynamical system to infinitesimal external perturbations is described by the covariance matrix of the unperturbed system. The two key assumptions commonly considered are: (a) the system state is close to equilibrium, and (b) the probability density function of the unperturbed system is Gaussian. Meanwhile, for a climate system, the standard hypotheses of equilibrium statistical mechanics do not hold since the climate system is highly non-linear and dissipative, and is affected by strong external periodic and stochastic forcing.

In spite of strong non-linearity in the climate system, the linear approximation and time-invariant hypothesis are still extensively applied in climate research [43]. A basic assumption in a linear approximation is that the different external perturbations are acting independently and additively on the system’s response. Thus, the response of some climate variable *x* to external forcing can be represented as where *S* is a sensitivity function (coefficient). For example, to estimate the change in equilibrium surface temperature ∆Θ_{s} due to the increase in radiative forcing, the following simplified relationship between carbon dioxide and radiative forcing can be used [44, 45]:

where *c*(*t*) is the concentration of carbon dioxide in parts per million by volume at time *t* and *c*(0) is the reference concentration. Therefore, *S* = 0.8 *K* ⋅ *W*^{− 1}m^{2}. For doubling of CO_{2} concentration, this gives the warming of ~ 3 K.

The response of non-linear systems to external perturbations is fundamentally different from the reaction of linear systems [45]. This difference is mainly due to the wider involvement of the system’s inherent characteristics and irregularity, and various ways of taking them into account. Consequently, system’s fluctuations represent the integration of external forcing and internal feedbacks. The FDT allows us to clearly identify external forcing and separate them from the system’s natural oscillations. Let us consider a finite-dimensional dynamical system described by state vector *x*:

Since we are interested in statistical characteristics of this system, let us introduce its average state:

Together with the system (27), let us consider the perturbed system by adding some external forcing to the right-hand side of (27):

If Eqs. (27) and (29) have stationary probability density functions *ρ* and *ρ* ′ respectively, then

It is important to consider that the average state of the perturbed system (29) *δx* =

where is a certain non-linear function. For small perturbations, we can expect that the relation between and *δx* is nearly linear. If is differentiable at a certain point , then can be represented by a Taylor series. Then if we omit all the non-linear terms and *leave the first order, linear terms*, we obtain:

where is a linear response operator whose properties are not known *a priori*. Thus, the climate sensitivity problem is to find an operator . However, this problem is not trivial due to the complicated fractal structure of the attractors (as the set of states) of chaotic systems. Attractors are structurally stable with respect to small external perturbations only for hyperbolic dynamical systems [35, 36], and the response of these systems is linear with respect to small enough external forcing [47]. Since the climate system is not hyperbolic and its attractors depend on external perturbations, some special approaches should be considered. The use of ∈-regularization of fractal attractor of chaotic dynamical systems allows us to ensure that the attractor ceases to be fractal [48]. The main idea of this method is to add random noise to the right-hand side of the model Eq. (27) [46, 49, 50]:

where * is a small positive constant and **μ*(*t*) is Gaussian stochastic process. This idea has a reasonable physical base since physical mechanisms that are responsible for energy transformation and dissipation in the climate system are never accurately known. In addition, computer simulation of the system (27) will generate pseudorandom noise due to round-off errors. All of these effects may be considered as white Gaussian noise. Resulting, we can get the Fokker-Plank equation corresponding to (26) for studying the evolution of the probability measure (probability density function) [46, 50]:

The function satisfies the following conditions: *x* ∈ *X*, where *X* is a compact manifold without a boundary, then (i) the stationary solution of the Eq. (31) exists, (ii) this solution is unique and (iii) this solution is asymptotically stable [46]. In a general case of the phase space

The problem now is to find the relation between and

where . If the system dynamic is a stationary random process and is normally distributed, then the response operator can be represented as

and we can obtain the following relation between and

where is a *τ*-lagged covariance matrix of *x*. From this equation, it follows that the response operator for ergodic dynamical system can be calculated from a single trajectory. However, the state vector is dimensionally large; therefore, some procedure should be used to reduce its dimensionality in order to calculate the covariance matrix effectively. For the atmosphere, Eq. (36) was first obtained in [38] under very strict assumptions [42]: (a) the system has at least one quadratic invariant (energy); (b) it has an incompressible phase volume; and (c) it is forced by a weak source. Under these conditions, the system probability density function is Gaussian. Climate systems do not satisfy these assumptions since they have a contracting phase space and do not have any exact quadratic invariant.

The derivation of Eq. (36) represented in [46] is based on more realistic and less restrictive assumptions that make the true validity of Eq. (36) for any non-normal distribution . It can be shown that the relation (36) holds under: (a) weak stochastic regularization; (b) weak forcing anomaly and (c) the system has a Fokker-Planck equation with a unique stationary solution.

As a simple example, consider one-dimensional stochastic dynamical system, i.e. system having only one variable Θ generated by the Langevin equation

Here *α* = 1/*τ*_{0}, where *τ*_{0} is the relaxation time for Θ, and

Angular brackets denote ensemble averages and *δ* is the Dirac delta function. Eq. (37) represents a simple version of radiative balance climate model, which is forced by an external heating whose time dependence is white noise [51, 52]. The variable Θ is the departure from steady state Θ_{0}, i.e. Θ is a climate anomaly. The relaxation time *τ*_{0} in Eq. (36) is *τ*_{0} ≈ 58 days [52]. The autocorrelation function for system (37) is

Solving Eq. (37) yields

Since the ensemble average of vanishes, one can average Eq. (37) to obtain

Thus, the average anomaly decays exponentially to zero in the relaxation time *τ*_{0}. Let the climate system be in equilibrium. Then at *t* = *t* ′, the system is perturbed by a sudden infinitesimal delta function impulse *ϑδ*(*t* − *t* ′) to the temperature. Consequently, Θ(*t*) at *t* = *t* ′ is changed on the value of *ϑ*. Then the temperature anomaly begins to decrease back to zero in accordance with Eq. (39). The mean system response to any infinitesimal change in the forcing can be calculated by the Eq. (37). In particular, for being a staircase function, i.e. a constant ∆_{s} that is activated at *t* = 0, the system response is given by

The validity of relation (36) was verified for different models of the atmosphere [46]. This relation holds with a high accuracy for both barotropic and two-layer baroclinic global atmospheric models and with a satisfactory accuracy for global general circulation models of the atmosphere.