InTechOpen uses cookies to offer you the best online experience. By continuing to use our site, you agree to our Privacy Policy.

Earth and Planetary Sciences » Oceanography and Atmospheric Sciences » "Topics in Climate Modeling", book edited by Theodore Hromadka and Prasada Rao, ISBN 978-953-51-2661-4, Print ISBN 978-953-51-2660-7, Published: October 5, 2016 under CC BY 3.0 license. © The Author(s).

# Climate Model Sensitivity with Respect to Parameters and External Forcing

By Sergei A. Soldatenko and Denis Chichkine
DOI: 10.5772/64710

Article top

## Overview

Figure 1. Two largest conditional Lyapunov exponents as functions of coupling strength parameter.

Figure 2. Time dynamics of sensitivity functions with respect to parameter c on the time interval [0, 20] for c0 = 0.9.

Figure 3. Time dynamics of sensitivity functions with respect to parameter r on the time interval [0, 25] for c0 = 0.9.

Figure 4. Original (in red) and pseudo (in blue) orbits for fast z and slow Z variables for c = 0.015.

Figure 5. Original (in red) and pseudo (in blue) orbits for fast z and slow Z variables for c = 0.8.

Figure 6. Differences between variables that correspond to the original trajectories and pseudo orbits for c = 0.01.

# Climate Model Sensitivity with Respect to Parameters and External Forcing

Sergei A. Soldatenko1 and Denis Chichkine2
Show details

## Abstract

Mathematical modelling is one of the most powerful methods for the study and understanding of the Earth’s climate system and its components. Modern climate models used in variety of applications are derived from a set of multi-dimensional non-linear differential equations in partial derivatives, which describe dynamical, physical and chemical processes and cycles in the climate system. Climate models are mostly deterministic with a large-phase space dimension containing a vast number of parameters that have various meanings. Most of them are not well-known a priori and, hence, are not well defined. Parameter errors and their time and space variabilities generate parametric uncertainty. Some model parameters describe external forcing that can strongly influence the climate model behaviour. It is, therefore, very important to estimate the influence of variations in parameters on the model behaviour and results of simulations. Questions like these can be answered by performing sensitivity analysis. This chapter considers various methods of sensitivity analysis that can be used: first, to estimate the influence of model parameter variations on its behaviour; second, to identify parameters of climate models and third, to study the model response to external forcing.

Keywords: climate, dynamical systems, sensitivity analysis

## 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 [16]. 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, 911] 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 (XSt) , where X is the phase space of a system, and St : 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 St, t ≥ 0, forms a semi-group that satisfies the following conditions [12]:

1. S0 ≡ I, where I is the identity operator

2. Stx is continuous in both t and x ∈ X

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

 ẋ=fx,x0=x0,x∈X (1)

Here x0 ∈ X is a system state specified at t = 0 and f is continuous vector-valued function. The solution of a system (1) x = x(xt) is determined for all t ≥ 0 and can be represented as x(t) = x(x0t) = Stx0. 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:

 xk+1=fxk,k∈Z+ (2)

Given the system state x0 at the initial time t = 0, we define the trajectory of x0 under f to be the sequence of points xkX:k+ such that xk = f k(x0), 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 x0 are given. The state of dynamical system xk at time tk is defined by the system state xk − 1 at the previous time tk − 1 as f(xk − 1).

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

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)‖ < V0 at t > t*, where V0 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 V0. This property guaranties the existence inside V0 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 x¯=x and second varx=x2-x¯2 moments) can be calculated by averaging along a certain trajectory. However, attractors of dissipative dynamical systems have highly complex fractal structure. Such attractors are commonly referred to strange attractors. Since the phase volume of dissipative dynamical systems contracts continuously in the limit of large time, the dimension of attractor is found to be smaller than the dimension of system phase space. Since the phase volume of a system is expanded in 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 nA nested in the absorbing ball, where nλ ≤ nA < 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:

 DKY=j+∑i=1jλiλi+1 (3)

where j is the maximum integer such that the sum of the j largest exponents is still non-negative, i.e. i=1jλi>0. The sum of all positive Lyapunov exponents, according to the theorem [17], gives an estimate of Kolmogorov-Sinai entropy, i.e. the value showing the mean divergence of the trajectories on the attractors. The arrangement of the Lyapunov exponents in Eq. (3) is as follows: λ1λ2λnd. The multiplicative inverse (reciprocal) of the largest Lyapunov exponent is referred as characteristic 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:

ψt+Jψ,ψ+l+h=ϕ-αψ+ν2ψ,ψ0=ψ0

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]:

DF12π23Gr2312+ln32π+lnG13

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 φ ∈ Qt) the state vector that characterizes the ECS in the domain Ωt. Note that Qt) 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:

 ∂tφrt=Lφrt,λrt,φr0=φ0r (4)

where is a non-linear, multi-dimensional differential operator that describes the dynamics, dissipation and external forcing of the system, λ ∈ Gt) is the parameter vector, Gt) 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 Ψ=ψii=1n that is defined to represent state vector on the domain Ω. Thus, state vector φ can be introduced in the form of normally convergent series:

 φrt≈∑i=1nxitψir (5)

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:

 ẋ=Fx,α,t∈0τ,x0=x0 (6)

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 is the 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:

xk=m0,kx0,k=1,,K

where m0,k is non-linear operator that indirectly contains model parameters and propagates the state vector from time t0 (the initial conditions) to time tk, 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 xi with respect to a certain model parameter αj [20, 21]:

 Sijtα≡∂xitα∂αjα0=limδαj→0xit,αj0+δαj-xitαj0δαj (7)

where δαj is the infinitesimal perturbation of parameter αj around some fixed point αj0. Approximating the state vector x(α0 + δα) around x(α0) by Taylor expansion, one can obtain the following linear equation:

xα0+δα=xα0+xαα0δα+H.O.T.

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

 ẋi=Fixα,t∈0τ,xi0=xi0,i∈1n (8)

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

 ddt∂xi∂αj=∑k=1n∂Fi∂xk∂xk∂αj+∂Fi∂αj,i∈1n,j∈1m (9)

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:

 ẋi=Fixα,xi0=xi0Ṡj=M∙Sj+Dj,Sj0=Sj0 (10)

where Sj = (∂x/∂αj) = (S1j, ⋯, Snj)T is the sensitivity vector with respect to parameter αj, MMij=Fi/xj is a Jacobian matrix and Dj = (∂F1/∂αj, ∂F2/∂αj, …, ∂Fn/∂α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]:

 Rxα=∫0τΦtxαdt (11)

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 x0:

αRx0α0=dRdα1dRdαmTx0,α0

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

dRdαjαj0Rx0+δx;α10,,αj0+δαj,,αm0-Rx0α0δαj

where δαj is the variation in parameter αj0. Note that

dRdαj=i=1nRxixiαj+Rαj=i=1nSijxiαj+Rαj

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:

 δRx0α0=∫0τ∂Φ∂xx0,α0∙δx+∂Φ∂αx0,α0∙δαdt (12)

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 x0(t), we obtain the following system of variational equations, the so-called tangent linear model, for calculating δx:

δxt=Fxx0,α0δx+Fαx0,α0δα,t0τ,δx0=δx0

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]:

 ∇αRx0α0=∫0τ∂Φ∂αx0,α0-∂F∂αx0,α0T∙x*dt (13)

where the vector function x is the solution of adjoint model, which is numerically integrated in the inverse time direction:

 -∂x*∂t-∂F∂xx0,α0T∙x*=-∂Φ∂xx0,α0,t∈0τ,x*0=0 (14)

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 yobs be the set of observations and H be observation operator mapping from solution space of model to observation space. Therefore, yobs = 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:

α*=argminJxα

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

 Jα=12Hx-yobsW2 (15)

For illustrative purposes, let us consider the following example. Let yiobstk and xi(tkα) be, respectively, the observation and model prediction of the ith component of state vector at time tk, σi2 the variance of yiobstk, and H the identity operator. Then the objective function can be written as

 Jxα=12∑i∑kwixitkα-yiobstk2 (16)

where are a weighted coefficient reflecting the accuracy of observations (in our case, wi=1/σi2). Many optimization algorithms rely on descent methods that require the computation of the gradient of the objective function. The gradient of Eq. (16) with respect to parameter αj is defined as

 ∂J∂αj=-∑i∑kwi∆xitk∂xitkα∂αj (17)

where xitk=xitkα-yiobstk. The right-hand side of Eq. (17) shows that sensitivity functions play a critical role in determining the corrected values of model parameters.

## 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]:

 x.=σy−x−caX+ky.=rx−y−xz+caY+kz.=xy−bz+czZ (18)
 X.=εσY−X−cx+kY.=εrX−Y−aXZ+cy+kZ.=εaXY−bZ−czz (19)

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 cz is a coupling strength parameter for z − Z variables. One can assume that a = 1, k = 0 and c = cz without loss of generality. Therefore, the vector of state variables of coupled model (18) and (19) is x = (xyzXYZ)T and the vector of model parameters is α = (σrbcε)T. In the operator form, the set of Eqs. (16)(21) can be rewritten as follows:

 ẋ=L+Qx (20)

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

L=σσ0r1x00xbεσεσ00εrεεX0εXεb,Q=c0000c000cc000c0000c

The unperturbed parameter values are selected to be as follows:

σ0=10,r0=28,b0=8/3,ε0=0.1,c00.11.2

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.

#### Figure 1.

Two largest conditional Lyapunov exponents as functions of coupling strength parameter.

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 xii=0n, an autocorrelation function is determined as where the brackets denote ensemble averaging. Assuming time series originates from a stationary and ergodic process, ensemble averaging can be replaced by time averaging over a single normal realization .

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:

S1c=x/c,S2c=y/c,S3c=z/c
S4c=X/c,S5c=Y/c,S6c=Z/c

The corresponding sensitivity equations can be written as

Ṡ1c=σS2c-S1c-cS4c-X
Ṡ2c=rS1c-S2c-xS3c-zS1c+cS5c+Y
Ṡ3c=xS2c+yS1c-bS3c+cS6c+Z
Ṡ4c=εσS5c-S4c-cS1c-x
Ṡ5c=εrS4c-S5c-XS6c-ZS4c+cS2c+y
Ṡ6c=εXS5c+YS4c-bS6c-cS3c-z

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:

#### Figure 2.

Time dynamics of sensitivity functions with respect to parameter c on the time interval [0, 20] for c0 = 0.9.

S1r=x/r,S2r=y/r,S3r=z/r
S4r=X/r,S5r=Y/r,S6r=Z/r

The associated system of sensitivity equations can be written as

Ṡ1r=σS2r-S1r-cS4r
Ṡ2r=x+rS1r-S2r-xS3r+zS1r+cS5r
Ṡ3r=xS2r+yS1r-bS3r+cS6r
Ṡ4r=εσS5r-S5r-cS1r
Ṡ5r=εX+rS4r-S5r-XS6r+ZS4r+cS2r
Ṡ6r=εXS5r+YS4r-bS6r-cS3r

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 [3032] that αR(α)αR(α) . This is because the integral

I=limτ0τlimδα0Rα+δα-Rαδαdt

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

#### Figure 3.

Time dynamics of sensitivity functions with respect to parameter r on the time interval [0, 25] for c0 = 0.9.

## 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 :×XX generated by homeomorphism f:

Fkx=fkx,kZandxX

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(xk + 1f(xk)) < 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(fk(y), xk) < ε 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 ft:×XX generated by a set of ordinary differential equations. A function g:X is a d-pseudo-orbit of ft if the inequality dist(ft(tg(τ)), g(t + τ)) < d holds for any τ and 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]:

Repε=hRepht1-ht2t1-t2ε,foranydifferentt1,t2R.

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

 ẋ=fx,α,x∈ℝn,α∈R (21)

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α=limT1T0TRxtα,αdt is of the most interest for us. With these assumptions, we can use the arbitrarily chosen trajectory of the system to trace the state variables along the orbit and calculate 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 xr(t). According to the shadowing lemma, the closest orbit x(t) to xr(t) satisfies the following constrained minimization problem [37]:

 minx,τ1T∫0Txτt-xrt2+η2dτdt-12dtsuchthatẋ=fx,α (22)

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 xsTtα and τsTtα is a solution of the model equations and time transformation that provides the trajectory xsTtα to be close to xr(t). The performance measure 〈R(α)〉 averaged over the time interval t ∈ [0, T] can be then approximated as:

 Rα≈RsTα=1τT-τ0∫τ0τTRxsTtα,αdt (23)

since xsTtα satisfies the model equation at a different α. The desired sensitivity estimate αRsTα can be computed by solving the following linearized LLS problem [37]:

 minS,μ1T∫0TS2+η2μ2dtsuchthatdSdt=∂f∂xS+∂f∂α+μfxrα (24)

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

 St=ddαxsTτsTtα,α,μt=ddαdτsTtαdt (25)

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 α):

 ∇αRsTα≈1T∫0T∂R∂xS+∂R∂α+μR-R¯dt,whereR¯=1T∫0TRdt (26)

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].

#### Figure 4.

Original (in red) and pseudo (in blue) orbits for fast z and slow Z variables for c = 0.015.

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.

#### Figure 5.

Original (in red) and pseudo (in blue) orbits for fast z and slow Z variables for c = 0.8.

#### Figure 6.

Differences between variables that correspond to the original trajectories and pseudo orbits for c = 0.01.

c Z/∂rY/∂rX/∂rz/∂ry/∂rx/∂r
1.01.100.050.011.080.040.03
0.80.690.080.031.020.070.07
0.40.950.03−0.011.030.090.09
0.150.91−0.08−0.091.01−0.01−0.01
10−41.04−0.02−0.031.02−0.01−0.01

#### Table 1.

Sensitivity estimates of fast and slow variables with respect to parameter r.

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]:

FCO2t=5.35×lnct/c0

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=SΔFCO2, where S = 0.8 K ⋅ W− 1m2. For doubling of CO2 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:

 ẋ=fx,x0=x0,x∈ℝn (27)

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

 x=limT→∞1T∫0Txtdt (28)

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

 ẋ′=fx′+δF,x′0=x0′ (29)

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

 x=∫xρxdx,x′=∫x′ρ′x′dx′ (30)

It is important to consider that the average state of the perturbed system (29) x' can be different from that of the unperturbed system (27) and, therefore, the difference δx = xx' depends on the external forcing [46]:

δx=NδF

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:

δx=LδF+H.O.T

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]:

 ẋ=fx+ϵμt,x0=x0 (31)

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]:

 ∂ρϵ∂t=d∇ρϵ-divfxρϵ (32)

The function satisfies the following conditions: ρ0 and . If 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 n, the problem, unfortunately, remains unsolved. In a similar way, we can obtain the Fokker-Plank equation for the perturbed system:

 ∂ρ′∂t=ϵ2∇ρ′-divfx′+δFρϵ′ (33)

The problem now is to find the relation between and δF. If the increment δF is small enough then

 δx=LδF=∫0t∫xt+τBxtTρϵdxdτδF (34)

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

 Lt=∫0tCτC0-1dτ (35)

and we can obtain the following relation between and δF:

 δxt=∫0tCτC0-1dτ⋅δFt (36)

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

 Θ̇+αΘ=Ft (37)

Here α = 1/τ0, where τ0 is the relaxation time for Θ, and Ft is a broadband noise forcing such that

 Ft=0,FtFt′=ϑδt-t′ (38)

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

ΘtΘtΘ2=e-αt-t

Solving Eq. (37) yields

 Θt=Θ0e-αt+∫0te-αt′-tFt′dt′ (39)

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

Θ̇+αΘ=0

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

Θ=slimt0te-ατdτ=slimt1-e-αtα=sα

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.

## References

1 - Robinson WA (2001) Modeling Dynamic Climate System. New York: Springer-Verlag, 213 p.
2 - Marshall J, Plumb RA (2007) Atmosphere, Ocean, and Climate dynamics: An Introductory Text. Burlington: Elsevier Academic Press, 344 p.
3 - Dijkstra HA (2013) Nonlinear Climate Dynamics. New York: Cambridge University Press, 367 p.
4 - McGuffie K, Henderson-Sellers A (2014) The Climate Modeling Premier. 4th edition. New York: Wiley-Blackwell, 456 p.
5 - Barry RG, Hall-McKim EA (2014) Essentials of the Earth’s Climate System. New York: Cambridge University Press, 271 p.
6 - Goosse H (2015) Climate System Dynamics and Modeling. New York: Cambridge University Press, 378 p.
7 - Pedlosky J (1987) Geophysical Fluid Dynamics. New York: Springer-Verlag, 710 p.
8 - Pedlosky J (2003) Waves in the Ocean and Atmosphere. Berlin-Heidelberg: Springer-Verlag, 264 p.
9 - Katok A, Hasselbatt B (1987) Introduction to the Modern Theory of Dynamical Systems. New York: Cambridge University Press, 802 p.
10 - Tchuenche JM (2013) Dynamical Systems: Theory, Applications & Future Directions. New York: Nova Publ., 312 p.
11 - Brin M, Stuck G (2015) Introduction to Dynamical Systems. New York: Cambridge University Press, 256 p.
12 - Dymnikov VP, Filatov AN (1997) Mathematics of Climate Modeling. Boston: Birkhäuser, 264 p.
13 - Gill M, Childress S (1987) Topics in Geophysical Fluid Dynamics: Atmospheric Dynamics, Dynamo Theory, and Climate Dynamics. New York: Springer-Verlag, 512 p.
14 - Buzzi J (2012) Chaos and Ergodic Theory. In: Meyers RA, ed., Mathematics of Complexity and Dynamical Systems. New York: Springer-Verlag, pp. 63–87.
15 - Lorenz E (1963) Deterministic nonperiodic flow. Journal of the Atmospheric Sciences 20: 130–141.
16 - Kaplan JL, Yorke JA (1979) Chaotic behavior of multidimensional difference equations. In: Peitgen HO and Walter HO, eds., Functional Differential Equations and Approximations of Fixed Points. Lecture Notes in Mathematics. Berlin: Springer-Verlag, pp. 228–237.
17 - Pesin BJ (1977) Characteristic Lyapunov exponents and smooth ergodic theory. Russian Mathematical Surveys 32: 55–114.
18 - Palmer TN (2014) Predictability of Weather and Climate: From Theory to Practice. In: Palmer T. and R. Hagedorn, eds., Predictability of Weather and Climate. New York: Cambridge University Press, pp. 1–29.
19 - Stocker T (2011) Introduction to Climate Modelling. Advances in Geophysical and Environmental Mechanics and Mathematics.. Berlin-Heidelberg: Springer-Verlag, 182 p.
20 - Eslami M (1994) Theory of Sensitivity in Dynamical Systems: An Introduction. Berlin-Heidelberg: Springer-Verlag, 601 p.
21 - Rosenwasser E, Yusupov R (2000) Sensitivity of Automatic Control Systems. Boca Raton: CRC Press, 456 p.
22 - Cacuci DG (2003) Sensitivity and Uncertainty Analysis, Volume I: Theory. Boca Raton: CRC Press, 285 p.
23 - Cacuci DG, Ionesku-Bujor M, Navon IM (2005) Sensitivity and Uncertainty Analysis, Volume II: Applications to Large-Scale Systems. Boca Raton: CRC Press, 353 p.
24 - Bofetta G, Crisanti A, Paparella F, Provenzale A, Vulpiani A (1998) Slow and fast dynamics in coupled systems: A time series analysis view. Physica D 116: 301–312.
25 - Peña M, Kalnay E (2004) Separating fast and slow models in coupled chaotic systems. Nonlinear Processes in Geophysics 11: 319–327.
26 - Soldatenko S, Chichkine D (2014) Correlation and spectral properties of a coupled nonlinear dynamical system in the context of numerical weather prediction and climate modeling. Discrete Dynamics in Nature and Society 2014: 498184.
27 - Rosenblum M, Pikovsky A, Kurth J (1996) Phase synchronization of chaotic oscillators. Physical Review Letters 76: 1804–1807.
28 - Liverani C (1995) Decay of correlation. Annals of Mathematics 142: 239–301.
29 - Panchev S, Spassova T (2005) Simple general circulation and climate models with memory. Advances in Atmospheric Sciences 22: 765–769.
30 - Lea D, Allen M, Haine T (2000) Sensitivity analysis of the climate of a chaotic system. Tellus 52A: 523–532.
31 - Wang Q (2013) Forward and adjoint sensitivity computation for chaotic dynamical systems. Journal of Computational Physics 235: 1–13.
32 - Soldatenko S, Steinle P, Tingwell C, Chichkine D (2015) Some aspects of sensitivity analysis in variational data assimilation for coupled dynamical systems. Advances in Meteorology 2015: 1–22.
33 - Pilyugin S (1999) Shadowing in Dynamical Systems, Lecture Notes in Mathematics. Berlin: Springer-Verlag, 276 p.
34 - Palmer HJ (2000) Shadowing in Dynamical Systems. Theory and Applications.. Dordrecht: Kluwer, 299 p.
35 - Anosov DV (1967) Geodesic flows on closed Riemann manifolds with negative curvature. Proceedings of the Steklov Mathematical Institute 90: 1–209.
36 - Bowen R (1975) Equilibrium States and the Ergodic Theory of Anosov Diffeomorphisms. Lecture Notes in Mathematics.. Berlin: Springer-Verlag, 77 p.
37 - Wang Q, Hu R, Blonigan P (2014) Least squares shadowing sensitivity analysis of chaotic limit cycle oscillations. Journal of Computational Physics 267: 210–224.
38 - Leith CE (1975) Climate response and fluctuation dissipation. Journal of the Atmospheric Sciences 32: 2022–2026.
39 - Nyquist H (1928) Thermal adaptation of electric charge in conductors. Physical Review 32: 110–113.
40 - Callen HB, Welton TA (1951) Irreversibility and generalized noise. Physical Review 83: 34–41.
41 - Kubo R (1957) Statistical-mechanical theory of irreversible processes. I. General theory and simple applications to magnetic and conduction problems. Journal of Physical Society of Japan 12: 570–586.
42 - Kraichnan RH (1959) Classical fluctuation-relaxation theorem. Physical Review 113: 1181–1182.
43 - IPCC, 2013: Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. Stocker TF, Qin D, Plattner GK, Tignor M, Allen SK, Boschung J, Nauels A, Xia Y, Bex V and Midgley PM, eds., Cambridge: Cambridge University Press, 1535 p.
44 - Myhre G, Highwood EJ, Shine KP, Strodal F (1998) New estimates of radiative forcing due to well mixed greenhouse gases. Geophysical Research Letters 25: 2715–2718.
45 - Ruzmaikin A (2014) Climate as a game of chance. Physics Uspekhi 184: 297–311.
46 - Gritsun A, Branstator G (2007) Climate response using a three-dimensional operator based on the fluctuation-dissipation theorem. Journal of the Atmospheric Sciences 64: 2558–2575.
47 - Ruelle D (1998) General linear response formula in statistical mechanics and the fluctuation-dissipation theorem far from equilibrium. Physics Letters A 245: 220–224.
48 - Zeeman EC (1998) Stability of dynamical systems. Nonlinearity 1: 115–135.
49 - Risken H. (1984) The Fokker-Plank Equation: Methods of Solution and Applications. New York: Springer-Verlag, 454 p.
50 - Dymnikov VP (2002) Fluctuation-dissipation relations for dynamic stochastic equations with coefficients periodically dependent on time for dissipative systems with random forcing. Atmospheric and Oceanic Physics 38: 658–663.
51 - Budyko MI (1969) The effect of solar radiation variations on the climate of the earth. Tellus 21: 611–619.
52 - North GR, Cahalan RF (1981) Predictability of a solvable climate model. Journal of the Atmospheric Sciences 38: 504–513.