Predictability in Deterministic Dynamical Systems with Application to Weather Forecasting and Climate Modelling

Abstract

Climate system consisting of the atmosphere, ocean, cryosphere, land and biota is considered as a complex adaptive dynamical system along with its essential physical properties. Since climate system is a nonlinear dissipative dynamical system that possesses a global attractor and its dynamics on the attractor are chaotic, the prediction of weather and climate change has a finite time horizon. There are two kinds of predictability of climate system: one is generated by uncertainties in the initial conditions (predictability of the first kind) and another is produced by uncertainties in parameters that describe the external forcing (predictability of the second kind). Using the concept of the ‘perfect’ climate model, two kinds of predictability are considered from the standpoint of the mathematical theory of climate.

1. Introduction

High‐complexity computational models that simulate earth's climate system (ECS) have earned well‐deserved recognition as the indispensable and primary instrument for numerical weather prediction (NWP) as well as for the study of climate change and variability caused by both natural processes and human activities [1–4]. In spite of dramatic progress achieved over the past few decades in weather forecasting and climate simulation thanks to the advances in computing hardware and algorithms and to a substantial increase in the volume of climatological data, contemporary computational climate models can reconstruct the real world only with a certain degree of validity [3]. There are several major sources of discrepancy between climate model simulation results and reality. First of all, climate models remain an ideal mathematical abstraction of a real physical system, namely the ECS. These models ignore some physical, dynamical and chemical processes or, at least, represent them in a simplified fashion. As a result, various physical simplifications in the formulation of climate models substantially influence their adequacy [5]. Second, the NWP and climate simulation are mathematically an initial‐value (Cauchy) and/or a boundary‐value (Dirichlet or von Neumann) problem, which is solved numerically using finite‐difference, spectral or another appropriate method. Consequently, uncertainties emerging in the initial and boundary conditions as well as in the climate model parameters and external forcing, approximation, truncation and round‐off errors lead to distinctions between the model output and the observed real state of the ECS. Third, let us suppose that we have the ‘perfect’ model of the ECS. It means that exact governing equations are known exactly and can be solved. However, even in this, hypothetically ideal, case the ability of climate models to predict the future remains limited. This can be explained by the fact that the atmosphere, which is the most rapidly changing component of the ECS, is strongly nonlinear and exhibits irregular (chaotic) spatial‐temporal oscillations on all scales ranging from millimetre seconds (turbulent fluctuations) to thousands of kilometres and several years (climate variability). This phenomenon known as deterministic chaos was first discovered by Lorenz [6]. The chaotic nature of the atmosphere significantly limits our ability to successfully predict the weather and climate since the predicted trajectory of the ECS is unstable with respect to both the infinitesimal errors in initial conditions and external forcing [7]. Even with a perfect atmospheric model and accurate initial condition, we cannot predict the weather beyond approximately two weeks.

For further discussion, we need to clarify that terms ‘weather’ and ‘climate’ have different meanings. Weather is defined as the daily conditions of the atmosphere in terms of such atmospheric variables as temperature, humidity, wind direction and velocity, surface pressure, cloud cover and precipitation. In turn, the climate represents an ensemble of states traversed by climate system over a sufficiently long temporal interval (about 30 years, according to the World Meteorological Organization). Here, the ensemble includes not only a set of system states but also the probability measure defined on this set. Therefore, climate, roughly speaking, can be considered as the ‘average’ weather, in terms of mean and variance, in a certain geographical location over many years.

Time horizon of a forecast's usefulness and validity can be characterized by the specific measure known as predictability. Predictability is commonly understood as the degree to which it is possible to make an accurate qualitative or quantitative forecast of the future system's state. The study of atmospheric predictability was initiated by Thompson [8] and Lorenz [6, 9] more than 50 years ago and was extensively explored theoretically using various numerical and statistical models since then (e.g. [10–17]). One of the obvious measures of predictability that can be used to verify a weather forecast is the mean‐squared error (the average of the squared differences between forecasts and observations). This measure increases over time and asymptotically approaches some finite value known as the saturation value. Therefore, predictability is lost when the forecast errors become comparable to the saturation value in magnitude. If this happens, the forecast result is not better than any randomly selected trajectory of the system. However, for a number of reasons, mean‐squared error and other weather forecast verification metrics (e.g. mean absolute error and mean error) are rarely used to estimate the climate system predictability in practice (for details, see Ref. [18]).

Predictability characterizes both the physical system itself and the model of this system that is used to make a forecast. However, in atmospheric and climate studies we are interested in the predictability of real dynamical processes rather than the predictability of the model used in simulations.

According to Lorenz [19], in weather and climate modelling we are facing the predictability of two kinds reflecting the internal and external variability of the climate system, respectively. The predictability of the first kind relates to the Cauchy (initial value) problem, namely the prediction of sequential states of the ECS for constant values of external parameters and given variations in the initial conditions. In contrast, the predictability of the second kind refers to a boundary‐value problem, specifically to the prediction of response of the climate system in asymptotical equilibrium to perturbations in external parameters (forcing).

This chapter considers both the predictability of atmospheric and climate processes with respect to the initial data errors (predictability of the first kind) as well as the predictability with respect to external perturbations (predictability of the second kind). The stability of dynamical system is also discussed since stability is a key problem related to predictability in dynamical systems.

2. Climate system as a complex adaptive dynamical system

Let us begin with some preliminary notes and definitions which will be used in this chapter.

The term ‘system’ generally refers to a goal‐oriented set of interconnected and interdependent elements that operate together to achieve some objectives [20]. The system is called complex if it possesses such characteristics as emergent behaviour, nonlinearity and high sensitivity to initial conditions and/or to perturbations, self‐organization, chaotic behaviour, feedback loop, spontaneous order, robustness and hierarchical structure. Complexity in systems arises from nonlinear spatio‐temporal interactions between their components. These nonlinear interactions lead to the appearance of new dynamical properties (for example, synchronous oscillations and other structural changes) that cannot be observed by exploring constituent elements individually.

Complex systems include a special class of systems that have the capacity to adapt to system's environment. These systems are known as complex adaptive systems. In a complex adaptive system, parts are linked together in such a way that the entire system as a whole has the capacity to transform fundamentally the interrelations and interdependences between its components, the collective behaviour of a system and also the behaviour of individual components due to the external forcing. Complex adaptive systems are dynamical systems since they evolve and change over time. These systems have a number of properties that include the following [21, 22]: co‐evolution, connectivity, sub‐optimality, requisite variety and iteration, edge of chaos and, certainly, emergence and self‐organization.

The ECS ( $S$ ) is understood as a complex, large‐scale physical system that consists of the following five basic and interacting constituent subsystems [23]:

1. Atmosphere ( $A$ ), the gaseous and aerosol envelope of the earth that propagates from the land, water bodies and ice‐covered surface outwards to space.

2. Hydrosphere ( $H$ ), the ocean and other water bodies on the surface of our planet, and water that is underground and in the atmosphere.

3. Cryosphere ( $C$ ), the sea ice, freshwater ice, snow cover, glaciers, ice caps and ice sheets and permafrost.

4. Lithosphere ( $L$ ), the solid, external part of the earth.

5. Biosphere ( $B$ ), the part of our planet where life exists, i.e. $S=A\cup H\cup C\cup L\cup B$

The ECS components are characterized by a finite set of distributed variables whose values at a given time determine their state. The most unstable and rapidly oscillating component of the ECS is the atmosphere.

The ECS is a large‐scale and unique physical system that possesses a number of specific properties (e.g. [24–29]) making the exploration of this system a high complexity problem. In contradistinction to many problems in physics, the study of the climate system, its change and variability cannot be implemented by a direct physical experiment due to climate system's essential features as a large‐scale physical system. Laboratory experiments and analytical approaches have a very limited applicability to climate exploration by virtue of extreme complexity of the ECS. As a result, in climate studies the computational simulation represents the primary instrument and as such requires the development of appropriate mathematical models and numerical algorithms.

The utilization of mathematical models in climate research involves the development of a specific mathematical theory that allows one to explore the climate system along with its mathematical models. The contemporary mathematical theory of climate is based on methods of the qualitative theory of differential equations that enables us to explore the behaviour of climate system in its phase space [30]. In other words, the dynamical system theory is currently the theoretical foundation of mathematical climate theory. In this context, the ECS can be viewed as a complex adaptive dynamical system [21, 22].

The ECS belongs to the class of complex adaptive systems due to the following factors:

1. The ECS is a complex large‐scale physical system combining the atmosphere, hydrosphere, cryosphere, land and biota together with global biochemical cycles (such as cycles of CO₂, N₂O and CH₄) and aerosols. Components of the climate system are heterogeneous thermo‐dynamical subsystems characterized by specific variables that determine their states. Elements of the ECS have strong differences in their structure, dynamics, physics and chemistry. They cover processes with different temporal and spatial scales, and link together via numerous physical coupling mechanisms, which can be either weak or strong. Each subsystem of the ECS can in turn be viewed as being composed of subsystems, which are themselves composed of subsystems. For example, the atmosphere can be divided into several layers based on its vertical temperature distribution. These layers are respectively the troposphere, stratosphere, mesosphere and thermosphere. The atmosphere can also be divided into surface layer, boundary layer and free atmosphere based on the influence of surface friction.

2. Each component of the ECS is characterized by a specific response time. This fact is very important to building the ECS’ models. The relation of a certain component to some ECS’ model is determined by the ratio between the temporal scale of processes under consideration and its response time. For example, the atmosphere, which has a response time of about one month in the troposphere, can be considered a sole component of the ECS’ model for processes with temporal scales of days to weeks. In this case, oceans, land surface and ice cover are considered as the boundary conditions and/or external forcing. If we study processes which have temporal scales of months to years, the atmosphere and ocean must be included in the ECS’ model together with sea ice. Thus, computational models of the ECS are built up from hierarchy of models, forming finally a complex integrated model.

3. The ECS has a large number of positive and negative feedback mechanisms which control the behaviour of the ECS. Some examples of these mechanisms are ice‐albedo feedback (positive feedback), water vapour feedback (positive feedback), cloud feedback (both positive and negative feedbacks), carbon cycle feedback (negative feedback), feedback due to Arctic methane release (positive feedback) and many others.

4. Physical and dynamical processes in the ECS cover a broad spectrum of temporal and spatial scales. Time scales are varied from seconds to decades, and spatial spectrum of dynamical processes covers molecular to planetary scales. Dynamical processes in the ECS and its components are nonlinear. Subsystems of the ECS interact with one another nonlinearly producing, under certain conditions, a chaotic behaviour of subsystems and the overall climate system.

5. The ECS and its components inherently have emergent properties. Examples of atmospheric emergent phenomena include but are not limited to clouds, large‐scale eddies (cyclones and anticyclones) and small‐scale vortices such as tornados. Examples of climate emergent phenomena are the El Niño–Southern Oscillation, which is a quasi‐periodical irregular variation in the ocean surface temperature over the Pacific in tropics that strongly influences global climate, ocean circulation patterns and glacial‐interglacial cycles. Natural emergent phenomena appear spontaneously under certain favourable conditions.

6. The ECS is a thermodynamically open and non‐isolated system because it exchanges energy with its surroundings. However, the ECS is a closed system with respect to the exchange of matter with its surroundings. The energy that drives the ECS is mainly solar energy. The ECS is affected by changes in external driving forces, which imply natural causes such as solar activity variations and volcanic activities, as well as man‐made changes in chemical composition of the atmosphere. However, the impact of the ECS on the outer space is insignificant. Currently, changes in climate are mostly affected by variations in the atmospheric composition of particles and gases. In the Arctic, the role of changes in albedo (reflection coefficient) is also tangible. The most influential gas component to affect the climate is CO₂, which comprises about 70% points of the global warming potential.

7. The components of the ECS are also non‐isolated systems. They act as cascading systems and interact with each other in various ways including through the transfer of momentum, sensible and latent heat, gases and particles. All together they compose the climate system, which is a unique large‐scale natural system.

8. Dynamical processes in the ECS fluctuate due to both internal factors (natural oscillations) and external forcing (forced oscillations). Natural fluctuations are caused by internal instability (for example, hydrodynamic instability such as barotropic and baroclinic) with respect to stochastic perturbations. Human impacts, both intentional and unintentional, belong to the category of external forcing.

Undoubtedly, there are other specific properties of the ECS that should be taken into account while studying climate as a complex adaptive system and building models of the ECS.

To simulate the ECS, we should assign some mathematical object that is an abstract representation of the real climate system taking into account its essential features mentioned above. This object is known as a perfect model of the ECS. It is usually assumed that a perfect model is deterministic semi‐dynamical system that is dissipative, ergodic and possesses a global attractor. It is also assumed that any trajectory generated by the model is unstable [30].

Formally, an abstract climate system model represents a set of multi‐dimensional nonlinear differential equations in partial derivatives, which generates finite dimensional deterministic semi‐dynamical system of the form [24, 30]

$dx/dt=F\left(x,p,f\right), x\in {\mathbf{R}}^n,x{|}_{t=\text{0}}=x_0,t\ge \text{0},$ E1

where x is the state vector, the components of which characterize the state of a system at a given time t, $x_{\text{0}}$ is a given initial state of a system, n is the dimension of dynamical system, $p\in {\mathbf{R}}^p$ is the vector of model parameters and $f$ is the external forcing. The solution to climate model equations (1) cannot be found analytically and one needs to employ available numerical methods. For that reason, in order to obtain numerical solution, the original set of partial differential equations is replaced with discrete spatio‐temporal approximations using any appropriate technique (e.g. finite‐difference method, Galerkin approach, etc.). Thus, in weather and climate simulation we mainly deal with discrete dynamical systems.

Suppose that the set of n real variables $x_{\text{1}},x_2,\dots ,x_n$ defines the current state of discrete‐time dynamical system representing the ECS. A certain particular state $x=\left(x_{\text{1}},x_2,\dots ,x_n\right)$ corresponds to a point in an n‐dimensional space $Q\subseteq {\mathbf{R}}^n$ , known as the system phase space. Let $t_m\in {\mathbf{Z}}_{+}$ $\left(m=\text{0,1,2,}\dots \right)$ be the discrete time, and let $f=\left(f_{\text{1}},f_2,\dots ,f_n\right)$ be a smooth vector‐valued function defined in the domain $Q\subseteq {\mathbf{R}}^n$ . This function describes the evolution of the system state from one moment to another. Then, a deterministic discrete‐time semi‐dynamical system that approximates the continuous time dynamical system (1) can be specified by the following equation:

$x\left(t_{m+\text{1}}\right)=f\left(x\left(t_m\right)\right),x\left(t_{\text{0}}\right)=x_0,m=\text{0, 1, 2,}\dots ,.$ E2

It is obvious that a family of operators forms a semi‐group:

$f_{s+p}=f_s\circ f_p,f_{\text{0}}=I,\forall s,p\in {\mathbf{Z}}_{+},$ E3

where I is the identity operator. Therefore, the system state $x\left(t_m\right)$ at time $t_m$ can be explicitly expressed via the initial condition $x_{\text{0}}$ :

$x\left(t_m\right)=f^m\left(x_{\text{0}}\right),$ E4

where $f^m$ denotes an m‐folding application of $f$ to $x_{\text{0}}$ . The sequence ${\left\{x\left(t_m\right)\right\}}_{m=\text{0}}^{\infty }$ is a trajectory of system (2) in its phase space, which is uniquely defined by the initial values of state variables $x_{\text{0}}$ .

For reference, let us reproduce a couple of definitions [30].

Definition 1. The solution $x\left(t\right)$ to system (1) is Lyapunov stable if $\forall \epsilon >\text{0},\exists \delta \left(\epsilon \right)>0$ such that

$\Vert x_{\text{0}}-x_{\text{0}}^{*}\Vert <\delta \left(\epsilon \right)\Rightarrow \Vert x\left(t\right)-x^{*}\left(t\right)\Vert <\epsilon ,\forall t\ge \text{0},$ E5

where $x^{*}\left(t\right)$ is the solution to the system

$d{x}^{*}/dt=F\left(x^{*},p,f\right),x^{*}{|}_{t=\text{0}}=x_{0.}^{*}\text{.}$ E6

Definition 2. The solution $x\left(t\right)$ to system (1) is stable with respect to the continuous perturbation $\delta F$ if $\forall \epsilon >\text{0, }\exists \delta \left(\epsilon \right)>0$ such that

$\Vert \delta F\Vert <\delta \left(\epsilon \right)\Rightarrow \Vert x\left(t\right)-x^{*}\left(t\right)\Vert <\epsilon ,\forall t\ge \text{0},$ E7

where $x^{}\left(t\right)$ is the solution to the following perturbed equation:

$d{x}^{*}/dt=F\left(x^{*},p,f\right)+\delta F,x^{*}{|}_{t=\text{0}}=x_0^{*}.$ E8

These definitions are important when considering both kinds of predictability.

The key point for further consideration is the assumption that climate system model described by the set of nonlinear partial differential equations (1) is ‘perfect'. We suppose that system (1) is nonlinear dissipative semi‐dynamical system ( $t\ge \text{0}$ ) that has an absorbing set in the phase space and its solution exists and is unique for any $t\ge \text{0}$ . Next, we assume that the system (1) possesses a global attractor of finite dimension that is a certain set in the system's phase space towards which a system tends to evolve for a wide variety of initial conditions of the system. Global attractor is characterized by the attraction property and invariance [30]. So, the dynamics of system (1) can be formally divided into to two phases: (1) movement towards the attractor and (2) motion on the attractor. When studying the climate system stability and predictability we assume that the system trajectory is on the attractor and its dynamics are chaotic. We also assume that system (1) possesses the property of ergodicity. Thus, statistical characteristics of the climate system (e.g. the first $\overline{x}=⟨x⟩$ and second $var\left(x\right)=⟨x^{\text{2}}⟩-{\overline{x}}^{\text{2}}$ moments) can be calculated by averaging along a certain trajectory.

Structurally, any climate system model represents a set of interacting and interdependent models of lower level (i.e. atmospheric model, model of the ocean, etc.). The number of these lower level models is determined by the objectives of a problem under consideration. For example, to study the large‐scale climate variability the model can include the following major components: tropical, mid‐latitude and polar troposphere, stratosphere, ocean, land ice, ocean and sea ice, surface and boundary layers, hydrological cycle, clouds (e.g. convective and stratiform), precipitation, aerosols, CO₂ and CH₄ cycles, solar radiation, terrestrial emission, etc. Other subsystems of the ECS (e.g. vegetation, land surface and biota) can be considered as the boundary conditions and external forcing. In numerical weather prediction problem, some atmospheric model (either global, regional or local) is the main component of the forecasting system, while ocean, sea ice, land surface are used only to impose boundary conditions. Note that models of general circulation of the atmosphere and the ocean represent main computational instruments for simulating the ECS.

3. Climate model governing equations

The main energy source of the ECS is the Sun. Spatial inhomogeneity and temporal changes of the heat energy that the earth's surface receives from the Sun are the main cause of motions in the atmosphere and ocean. Equations that govern the atmospheric and oceanic circulation represent the mathematical expressions of fundamental laws of physics: conservation of momentum, conservation of mass, conservation of water and conservation of energy (the first law of thermodynamics). Some diagnostic relationships between variables are also used (i.e. the equation of state). Almost every model uses a slightly different set of equations tailored to a specific problem. However, all climate models include the following basic equations: two equations for horizontal motions (or equation for the vorticity and divergence), equation for the vertical velocity (or hydrostatic equation), continuity equation, as well as thermodynamic and moisture equations. Equations of motion are derived from the law of conservation of momentum applicable to a rotating system. These equations describe all types and scales of atmospheric motions that are important for the formation of weather and climate (i.e. large‐scale Rossby waves, planetary waves and gravity waves). Conservation of mass is mathematically expressed in the form of continuity equation, equation for conservation of moisture and equations for conservation of other substances taken into account in a particular climate model.

The set of equations that describes the general circulation of the atmosphere can be written in the spherical co‐ordinate system ( $\lambda ,\varphi$ ) defined by longitude $\lambda$ and latitude $\varphi$ , with normalized pressure as a vertical coordinate $\sigma =p/p_s$ , where p is pressure and $p_s$ is the surface pressure [1, 31]. The set of the model equations includes two momentum equations:

$\frac{\partial u}{\partial t}=\eta v-\frac{\text{1}}{a\text{cos}\varphi }\frac{\partial }{\partial \lambda }\left(\Phi +K\right)-\frac{\mathrm{R}{T}_v}{a\text{cos}\varphi }\frac{\partial \text{ln}{p}_s}{\partial \lambda }-\dot{\sigma }\frac{\partial u}{\partial \sigma }=F_{uV}+F_{uH},$ E9

$\frac{\partial v}{\partial t}=-\eta u-\frac{\text{1}}{a}\frac{\partial }{\partial \varphi }\left(\Phi +K\right)-\frac{\mathrm{R}{T}_v}{a}\frac{\partial \text{ln}{p}_s}{\partial \varphi }-\dot{\sigma }\frac{\partial v}{\partial \sigma }=F_{vV}+F_{vH},$ E10

where $u$ and $v$ are zonal and meridional velocities, $a$ is the earth's average radius, $\sigma =d\sigma /dt$ is the vertical velocity in the $\sigma$ co‐ordinate system, $\Phi$ is geopotential, $T$ is temperature, R is the gas constant for dry air, $K=\left(u^{\text{2}}+v^2\right)/\text{2}$ is the kinetic energy, $\eta =\varsigma +f$ is the absolute vorticity, $f$ is the Coriolis parameter and $\varsigma$ is the relative vorticity that is given by

$\varsigma =\frac{\text{1}}{a\text{cos}\varphi }\left[\frac{\partial v}{\partial \lambda }-\frac{\partial }{\partial \varphi }\left(u\text{cos}\varphi \right)\right]\text{.}$ E11

The virtual temperature $T_v$ is represented as

$T_v=T\left[\text{1}+\left(\frac{{\mathrm{R}}_v}{\mathrm{R}}-\text{1}\right)q\right],$ E12

where $T$ is the temperature, $q$ is the specific humidity and ${\mathrm{R}}_v$ is the gas constant for water vapour. The terms $F_{uV}$ and $F_{vV}$ describe the vertical friction and terms $F_{uH}$ and $F_{vH}$ the horizontal diffusion. Generally, however, the momentum equations are transformed into the equations for the absolute vorticity $\eta$ and the divergence $D$ using new independent variable $\mu =\text{sin}\varphi$ :

$\frac{\partial \eta }{\partial t}=\frac{\text{1}}{a\left(\text{1}-{\mu }^{\text{2}}\right)}\frac{\partial }{\partial \lambda }\left(N_v+\text{cos}\varphi F_{vV}\right)-\frac{\text{1}}{a}\frac{\partial }{\partial \mu }\left(N_u+\text{cos}\varphi F_{uV}\right)+F_{\eta H},$ E13

$\begin{array}{l}\frac{\partial D}{\partial t}=\frac{\text{1}}{a\left(\text{1}-{\mu }^{\text{2}}\right)}\frac{\partial }{\partial \lambda }\left(N_u+\text{cos}\varphi F_{uV}\right)+\frac{\text{1}}{a}\frac{\partial }{\nabla \mu }\left(N_v+\text{cos}\varphi F_{vV}\right)+F_{DH}\\ -{\nabla }^{\text{2}}\left(\Phi +K+\mathrm{R}{T}_0\text{ln}{p}_s\right),\end{array}$ E14

where the horizontal divergence is given by

$D=\frac{\text{1}}{a\text{cos}\varphi }\left[\frac{\partial u}{\partial \lambda }+\frac{\partial }{\partial \varphi }\left(v\text{cos}\varphi \right)\right]\text{.}$ E15

The spherical horizontal Laplacian can be written as

${\nabla }^{\text{2}}=\frac{\text{1}}{a^{\text{2}}\left(\text{1}-{\mu }^{\text{2}}\right)}\frac{{\partial }^{\text{2}}}{\partial {\lambda }^{\text{2}}}+\frac{\text{1}}{a^{\text{2}}}\frac{\partial }{\partial \mu }\left[\left(\text{1}-{\mu }^{\text{2}}\right)\frac{\partial }{\partial \mu }\right]\text{.}$ E16

To provide the computational effectiveness of numerical integration scheme, the virtual temperature is partitioned into two parts, one of which $T_{\text{0}}$ is a function of the vertical coordinate only, i.e. $T_v\left(\lambda ,\mu ,\sigma ,t\right)=T_0\left(\sigma \right)+T{\text{'}}_v\left(\lambda ,\mu ,\sigma ,t\right)$ . Then, the nonlinear dynamical terms $N_u$ and $N_v$ can be represented in the following form:

$N_u=\eta V-\mathrm{R}{T}_v^{\text{'}}\frac{\text{1}}{a}\frac{\partial \text{ln}{p}_s}{\partial \lambda }-\dot{\sigma }\frac{\partial U}{\partial \sigma },$ E17

$N_v=-\eta U-\mathrm{R}{T}_v^{\text{'}}\frac{\left(\text{1}-{\mu }^{\text{2}}\right)}{a}\frac{\partial \text{ln}{p}_s}{\partial \mu }-\dot{\sigma }\frac{\partial V}{\partial \sigma },$ E18

where $U=u\text{cos}\varphi$ and $V=v\text{cos}\varphi$ .

The thermodynamic equation, which represents the mathematical expression of the first law of thermodynamic, is written for a perturbation in temperature $T\text{'}$ calculated with respect to the mean $T_{\text{0}}\left(\sigma \right)$ mentioned above:

$\begin{array}{l}\frac{\partial T^{\text{'}}}{\partial t}=-\frac{\text{1}}{a\left(\text{1}-{\mu }^{\text{2}}\right)}\frac{\partial }{\partial \lambda }\left(U{T}^{\text{'}}\right)-\frac{\text{1}}{a}\frac{\partial }{\partial \mu }\left(V{T}^{\text{'}}\right)+T^{\text{'}}D-\dot{\sigma }\frac{\partial T^{\text{'}}}{\partial \sigma }+\frac{\mathrm{R}{T}_v}{c_p^{*}}\frac{\omega }{p}\\ +Q+F_{TV}+F_{TH}-\frac{\text{1}}{c_p^{*}}\left[u\left(F_{uV}+F_{uH}\right)+v\left(F_{vV}+F_{vH}\right)\right],\end{array}$ E19

where $Q$ is the diabatic heating rate, $\omega$ is the pressure vertical velocity and $c_p^{*}$ is given by

$c_p^{*}=c_p\left[\text{1}+\left(\frac{c_v}{c_p}-\text{1}\right)\right]\text{.}$ E20

Here, $c_p$ is the specific heat of dry air at a constant pressure and $c_v$ is the specific heat of water vapour at a constant pressure.

The equation for specific humidity is used to describe the hydrologic cycle in the atmosphere:

$\frac{\partial q}{\partial t}=-\frac{\text{1}}{a\left(\text{1}-{\mu }^{\text{2}}\right)}\frac{\partial }{\partial \lambda }\left(Uq\right)-\frac{\text{1}}{a}\frac{\partial }{\partial \mu }\left(Vq\right)+qD-\dot{\sigma }\frac{\partial q}{\partial \sigma }+S+F_{qV}+F_{qH},$ E21

where the term $S$ describes the source/sink processes for water vapour, and $F_{qV}$ and $F_{qH}$ are the vertical and horizontal water vapour diffusion.

Let us consider now the continuity equation that represents the conservation of mass law:

$\frac{\partial \text{ln}{p}_s}{\partial t}=-\frac{U}{a(\text{1}-{\mu }^{\text{2}})}\frac{\partial \text{ln}{\mathit{p}}_s}{\partial \lambda }-\frac{V}{a}\frac{\partial \text{ln}{\mathit{p}}_s}{\partial \mu }-D-\frac{\partial \dot{\sigma }}{\partial \sigma }\text{.}$ E22

Integrating this equation from the top ( $\sigma =\text{0}$ ) to the bottom ( $\sigma =\text{1}$ ), with the vertical boundary conditions $\dot{\sigma }=\text{0}$ at $\sigma =\text{1}$ and $\sigma =\text{0}$ , one can obtain the equation for surface pressure prediction:

$\frac{\partial \text{ln}{p}_s}{\partial t}={\displaystyle \underset{\text{0}}{\overset{\text{1}}{\int }}\left[D+\frac{U}{a(\text{1}-{\mu }^{\text{2}})}\frac{\partial \text{ln}{\mathit{p}}_s}{\partial \lambda }+\frac{V}{a}\frac{\partial \text{ln}{\mathit{p}}_s}{\partial \mu }\right]d\sigma }\text{.}$ E23

Combining the continuity equation and the equation for the surface pressure, one can derive the diagnostic equation for the vertical velocity $\dot{\sigma }$ :

$\dot{\sigma }=\sigma {\displaystyle \underset{\text{0}}{\overset{\text{1}}{\int }}\left[D+\frac{U}{a(\text{1}-{\mu }^{\text{2}})}\frac{\partial \text{ln}{\mathit{p}}_s}{\partial \lambda }+\frac{V}{a}\frac{\partial \text{ln}{\mathit{p}}_s}{\partial \mu }\right]d\sigma -{\displaystyle \underset{\text{0}}{\overset{\sigma }{\int }}\left[D+\frac{U}{a(\text{1}-{\mu }^{\text{2}})}\frac{\partial \text{ln}{\mathit{p}}_s}{\partial \lambda }+\frac{V}{a}\frac{\partial \text{ln}{\mathit{p}}_s}{\partial \mu }\right]d\sigma }}\text{.}$ E24

Two diagnostic equations, the hydrostatic equation and the equation of state, are also components of a set of equations that are used to simulate the atmospheric general circulation. The hydrostatic equation is

$\partial \Phi /\partial \text{ln}\mathit{\sigma }=-\mathrm{R}{T}_v\text{.}$ E25

In the integral form, this equation can be written as

$\Phi ={\Phi }_s-{\displaystyle \underset{\text{1}}{\overset{\sigma }{\int }}\mathrm{R}{T}_{v}d\text{ln}\sigma ,}$ E26

where ${\Phi }_s$ is the geopotential at the earth's surface. The equation of state is written as

$p=\rho \mathrm{R}{T}_v,$ E27

where $\rho$ is the air density.

Boundary conditions in the longitudinal direction are periodic, and the solution to the atmospheric model equations is bounded at the north and south poles. Vertical boundary conditions represent the vanishing of vertical velocity both at the bottom and at the top of the atmosphere: $\dot{\sigma }=\text{0}$ at $\sigma =\text{1}$ and $\sigma =\text{0}$ .

Equations used in the ocean model are written in the Boussinesq hydrostatic approximation with a rigid lid in the spherical coordinate system, with depth $z$ as a vertical coordinate defined as negative downwards from $z=\text{0}$ , which denotes the ocean surface [1, 31]. The set of model equations include the following:

1. The horizontal equations of motion:

   $\frac{\partial u}{\partial t}+L\left(u\right)-\left(f+\frac{u}{a}\text{tan}\varphi \right)v+\frac{\text{1}}{a{\rho }_o\text{cos}\varphi }\frac{\partial p}{\partial \lambda }=k_V\frac{{\partial }^{\text{2}}u}{\partial z^{\text{2}}}+F_u,$ E28
   $\frac{\partial v}{\partial t}+L\left(v\right)+\left(f+\frac{u}{a}\text{tan}\varphi \right)u+\frac{\text{1}}{a{\rho }_o}\frac{\partial p}{\partial \varphi }=k_V\frac{{\partial }^{\text{2}}v}{\partial z^{\text{2}}}+F_v,$ E29

   where $k_V$ is the vertical eddy viscosity coefficient, ${\rho }_{\text{0}}$ is the density of sea water and the advection operator, $L\left(\alpha \right)$ , is given by

   $L\left(\alpha \right)=\frac{\text{1}}{a\text{cos}\varphi }\left(\frac{\partial u\alpha }{\partial \lambda }+\frac{\partial v\alpha \text{cos}\varphi }{\partial \varphi }\right)+\frac{\partial w\alpha }{\partial z}\text{.}$ E30

   The horizontal viscosity terms, $F_u$ and $F_v$ , are defined as

   $F_u=k_H\left[{\nabla }^{\text{2}}u+\frac{(\text{1}-{\text{tan}}^{\text{2}}\varphi )u}{a^{\text{2}}}-\frac{\text{2sin}\varphi }{a^{\text{2}}{\text{cos}}^{\text{2}}\varphi }\frac{\partial v}{\partial \lambda }\right],$ E31
   $F_v=k_H\left[{\nabla }^{\text{2}}v+\frac{(\text{1}-{\text{tan}}^{\text{2}}\varphi )v}{a^{\text{2}}}+\frac{\text{2sin}\varphi }{a^{\text{2}}{\text{cos}}^{\text{2}}\varphi }\frac{\partial u}{\partial \lambda }\right],$ E32

   where $k_H$ is the horizontal eddy viscosity coefficient. The given form of the diffusion terms, $F_u$ and $F_v$ , is required for conserving angular momentum property.

2. The hydrostatic equation:

   $\partial p/\partial z=-g\rho .$ E33

3. The thermodynamic equation:

   $\frac{\partial T}{\partial t}+L\left(T\right)={\kappa }_V\frac{{\partial }^{\text{2}}T}{\partial z^{\text{2}}}+{\kappa }_H{\nabla }^{\text{2}}T,$ E34

   where ${\kappa }_V$ and ${\kappa }_H$ are, respectively, the vertical and horizontal eddy diffusivity coefficients.

4. The equation for the mass continuity of salinity:

   $\frac{\partial S}{\partial t}+L\left(S\right)={\kappa }_V\frac{{\partial }^{\text{2}}S}{\partial z^{\text{2}}}+{\kappa }_H{\nabla }^{\text{2}}S\text{.}$ E35

5. The equation of continuity:

   $\frac{\partial w}{\partial z}=-\frac{\text{1}}{a\text{cos}\varphi }\frac{\partial u}{\partial \lambda }-\frac{\text{1}}{a\text{cos}\varphi }\frac{\partial v\text{cos}\varphi }{\partial \varphi }.$ E36

6. The equation of state:

   $\rho =\rho \left(T,S\text{.}p\right)\text{.}$ E37

Due to their extreme complexity, weather and climate models can be implemented on computers only using numerical techniques. Since models are based on partial differential equations, it is necessary, first, to ensure that the problem under consideration is well posed, i.e. it has a unique solution that depends on the boundary and initial conditions. Thus, both the initial and boundary conditions must be properly specified. Next, weather and climate mathematical models should be transformed into numerical models that can be implemented on computers. The most widely used technique for solving differential equations of weather and climate models is the finite‐difference method according to which the derivatives in the partial differential equations are approximated on a certain temporal‐spatial grid. Thus, instead of continuous functions, which describe the state of climate system and its components, we deal with discrete functions defined only for specific times separated by the time step $\Delta t$ and specific space locations separated by spatial (horizontal $\Delta s$ and vertical $\Delta h$ ) steps. As a result, instead of partial differential equation we obtain finite‐difference equations (numerical model). It is very important that numerical schemes used for the discretization of model differential equations must satisfy several fundamental requirements: finite‐difference equations must be consistent with model differential equations, the solution of finite‐difference equations must converge to the solution of differential equations and numerical schemes must be computationally stable. In practice, finite difference is not the only method used to solve weather and climate problems. The most popular among other methods are the family of Galerkin techniques, spectral, finite‐volume and finite element approaches.

In contemporary climate models, due to their discrete spatial and temporal structure, a large number of physical processes and cycles cannot be clearly represented and formulated by model equations. Climate models are theoretically incapable of simulating processes on spatial scales of the order of magnitude that is twice the model grid length [32]. Such thermo‐dynamical, physical and chemical processes and cycles are parameterized, i.e. expressed parametrically using simplified description. Study of the climate system by computer simulation requires extensive computational resources. As a result, the predictability problem is usually studied either on the basis of low‐order models, which possess the main properties of the climate system (nonlinearity, chaos, dissipative, etc.), or on the basis of complex climate models using the ensemble approach or the Monte Carlo method.

4. Predictability of climate system

4.1. Predictability of the first kind

The first kind predictability of climate processes (predictability of climate processes with respect to the initial conditions) will be considered under the assumption that the climate system (1) evolves on its attractor. Since system (1) is a nonlinear dissipative dynamical system, its attractor, known as a strange attractor, has an extremely complex fractal structure and can be characterized by such parameters as dimension, characteristic Lyapunov exponents, invariant measure and asymptotically steady solution and others. If some trajectory of system (1) is enclosed in a bounded phase volume (attractor), then the system's dynamics demonstrate deterministic chaos: the behaviour of simulated system resembles a stochastic process despite the fact that the system is described by deterministic laws and its evolution is governed by deterministic differential equations. So, all orbits of a system that start close enough will diverge from one another, however, will never depart from the attractor. The rate of separation of infinitesimally close orbits is characterized by positive Lyapunov exponents. The number of directions along which the orbit is unstable is equal to the number of positive Lyapunov exponents $n_{\lambda }$ (note that $n_{\lambda }<n$ , where n is a system's dimension). Thus, trajectories of climate dynamical systems are Lyapunov unstable.

To consider the initial growth rates of errors in the initial conditions let us linearize Eq. (1) around some trajectory to obtain the equation in variations:

$d{x}^{\text{'}}/dt=M_{t}x{\text{'}}_0,$ E38

where $M_t=\partial F/\partial x$ is the tangent propagator along the trajectory between the initial state $x{\text{'}}_0$ and the forecast state $x\text{'}$ at a certain time t (actually $M_t$ is a Jacobian matrix). Obviously, one can obtain

${\Vert x^{\text{'}}(t)\Vert }^{\text{2}}=\left(M_{t}x{\text{'}}_0,M_{t}x{\text{'}}_0\right)=\left(M_t^{*}{M}_{t}x{\text{'}}_0,x{\text{'}}_0\right),$ E39

where (·,·) is the inner product in ${\mathbf{R}}^n$ and $M^{*}$ is the transpose of $M$ . Since the operator $M_t^{*}{M}_t$ is self‐adjoint, then for any t one can consider the following eigenvalue problem:

$M_t^{*}{M}_t{\psi }_i={\sigma }_i{\psi }_i,$ E40

where ${\sigma }_i$ is the ith eigenvalue of the matrix $M_t^{*}{M}_t$ and ${\psi }_i$ is the corresponding eigenvector. Representing $x{\text{'}}_0$ in the form of series as $x{\text{'}}_0={\displaystyle \sum _i{\alpha }_i{\psi }_i}$ , one can get ${\Vert x^{\text{'}}(t)\Vert }^{\text{2}}={\displaystyle \sum _i{\sigma }_i{\alpha }_i{}^{\text{2}}}$ . So, the forecast error on temporal interval $\left[\text{0,}t\right]$ depends on errors in the initial distribution of eigenvectors ${\psi }_i$ and singular values of the tangent linear propagator $M_t$ . Since system (1) is ergodic, we can also calculate the Lyapunov exponents ${\lambda }_i$ in accordance with the multiplicative theorem [33]:

${\lambda }_i=\underset{t\to \infty }{{\displaystyle \lim }}\frac{\text{1}}{t}\text{ln}{\sigma }_i\left(M_t^{*}{M}_t\right),i=\text{1,}\dots ,n.$ E41

The Lyapunov exponents define the exponential growth (decay) of linear independent components of x’ at $x^{\text{'}}\to \text{0}$ . The knowledge of the Lyapunov exponent spectrum of a dynamical system allows one to estimate the attractor fractal dimension, the rate of Kolmogorov‐Sinai entropy production and the characteristic e‐folding time. Knowledge of these parameters is very important for the stability and predictability analysis of dynamical systems. The fractal dimension of attractors of dissipative dynamical systems can be determined by applying the Kaplan‐Yorke conjecture [34]:

$D_{KY}=J+{\displaystyle \sum _{i=\text{1}}^J{\lambda }_i/\left|{\lambda }_{i+\text{1}}\right|},$ E42

where J is the maximum integer such that the sum of the J largest exponents is still non‐negative, i.e. ${\displaystyle {\sum }_{i=\text{1}}^J{\lambda }_i}>\text{0}$ . The sum of all positive Lyapunov exponents, according to theorem [35], gives an estimate of the Kolmogorov‐Sinai entropy, i.e. the value showing mean divergence of the trajectories on attractors. The arrangement of the Lyapunov exponents in (42) is as follows: ${\lambda }_{\text{1}}\ge {\lambda }_{\text{2}}\ge \dots \ge {\lambda }_{n_d}$ . The multiplicative inverse (reciprocal) of the largest Lyapunov exponent is referred to as the characteristic $e$ ‐folding time.

Let ${\delta }_{\text{0}}$ be the initial perturbation of $x\text{'}$ used to integrate equation (8). Since the system is Lyapunov unstable, after some sufficiently large temporal interval of integration the distance between two hyper‐points in the phase space reaches the value of ${\delta }_t$ . Let ${\overline{\delta }}_t$ be the accepted error tolerance, then the predictability time of a system can be roughly estimated as

$T_p\approx {\lambda }_{\text{max}}^{-\text{1}}\text{ln}\left({\overline{\delta }}_t/{\delta }_{\text{0}}\right),$ E43

where ${\lambda }_{\text{max}}$ is the leading Lyapunov exponent. The error doubling time can be calculated as $t=\text{ln2}/{\lambda }_{\text{max}}$ . However, Lyapunov exponents are very useful instrument to estimate the predictability of low‐order dynamical systems [36].

Climate data observations are subject to measurement errors. The simplest way to represent the resulting uncertainty is to define the probability density function (PDF) $\rho \left(x,t_0\right)$ or, generally, the set of a finite measure ${\mu }_{\text{0}}$ on which the initial state $x_{\text{0}}$ is concentrated. The time evolution of a system leads to a divergence and mixing of points of this set. Since the initial state $x_{\text{0}}$ is concentrated on a set having the measure ${\mu }_{\text{0}}$ , then after some period of time the measure will become ${\mu }_t$ . Let $\overline{\mu }$ be the invariant ergodic measure. Suppose the convergence theorem $\mu \to \overline{\mu }$ does exist. Hence, at a certain time $t\to t_{\epsilon }$ the measure ${\mu }_t$ falls into the $\epsilon$ ‐neighbourhood of $\overline{\mu }$ . Consequently, the initial data information characterized by ${\mu }_{\text{0}}$ will be completely lost. So, one can say that the time $t_{\epsilon }$ defines the potential predictability of a system under consideration [16]. Thus, a focal point of the predictability problem is to prove the existence of ergodic measure and the existence of convergence theorem. This problem, however, is extremely difficult to solve because the structure of the invariant measure generated on the system attractor is sophisticated and non‐smooth. To avoid this problem, the stochastic regularization can be applied [37]. So, in lieu of system (1), the following stochastic dynamical system will be considered [16]:

$dx/dt=F\left(x\right)+\eta (t),$ E44

where $\eta$ is a Gaussian stochastic process: $⟨{\eta }_i(t){\eta }_j(t^{\text{'}})⟩=\text{2}{d}_{ij}\delta (t-t^{\text{'}}),d_{ij}\ge \text{0}$ . This procedure is correct since our knowledge about the model parameters is always limited, thus real climate models have random errors, which are represented by the term $\eta$ . Under the assumption that $d_{ij}=d$ , one can write the Fokker‐Plank equation with respect to PDF $\rho \left(x,t\right)$ , which describes the evolution of $\rho$ [30]:

$\partial \rho /\partial t+div\left(F\left(x\right)\rho \right)=d\Delta \rho ,\rho \ge \text{0, }{\displaystyle \int \rho dx=\text{1}\text{.}}$ E45

Let $\overline{\rho }$ be a stationary solution to Eq. (45), i.e. $\text{div}\left(F\left(x\right)\rho \right)=d\Delta \rho$ . If $x$ belongs to the compact manifold without boundary, then $\overline{\rho }$ is asymptotically stable [37]. The existence of a stationary solution (i.e. attractor) at infinity has been proved for finite‐dimensional dynamical systems [38].

Suppose that the initial condition $x_{\text{0}}$ is specified then the condition $\rho {|}_{t=\text{0}}=\delta (x-x_{\text{0}})$ is also specified and enable us to solve Eq. (45). The numerical integration of Eq. (45) transforms the PDF $\rho \left(x,t\right)$ , which asymptotically evolves to the stationary solution $\overline{\rho }$ : $\rho \to \overline{\rho }$ at $t\to t_{\epsilon }$ . Thus, at sufficiently large time $t_{\epsilon }$ predictability is finally lost. There is a question: how can we estimate the time $t_{\epsilon }$ ? Let us consider the following one‐variable stochastic dynamical equation [16].

$dx/dt=-\gamma x+\eta ,$ E46

$x{|}_{t=\text{0}}=x_{\text{0}},⟨\eta \text{(}t\text{)}\eta \text{(}{t}^{\text{'}}\text{)}⟩=\text{2}{\eta }^{\text{2}}\delta \text{(}t-t^{\text{'}}\text{)},⟨\eta ⟩=\text{0},$ E47

where $x_{\text{0}}$ is the known initial condition and $\eta$ is the Gaussian $\delta$ ‐correlated process. If we average Eq. (46) we obtain

$d⟨x⟩/dt=-\gamma ⟨x⟩,⟨x⟩{|}_{t=\text{0}}=x_{\text{0}},$ E48

thus $⟨x⟩=x_0{e}^{-\gamma t}$ . For the newly introduced variable $\theta \left(t\right)=⟨x^{\text{2}}⟩$ , we can obtain the following equation:

$d\theta /dt=-\text{2}\gamma \theta +\text{2}⟨\eta \cdot x⟩\text{.}$ E49

Since $x\left(t\right)=x_0{e}^{-\gamma t}+{\displaystyle \underset{\text{0}}{\overset{t}{\int }}{e}^{\gamma \left(t-\tau \right)}\eta \left(\tau \right)d\tau }$ , then

$d\theta /dt=-\text{2}\gamma \theta +\text{4}{\eta }^{\text{2}}\text{.}$ E50

The solution to this equation is

$\theta \left(t\right)=\frac{\text{2}{\eta }^{\text{2}}}{\gamma }\left(\text{1}-e^{-\text{2}\gamma t}\right).$ E51

Equation for the PDF $\rho$ has the following form:

$\partial \rho /\partial t=\partial (\rho x\gamma )/\partial x+{\eta }^{\text{2}}{\partial }^{\text{2}}\rho /\partial x^{\text{2}}.$ E52

The stationary solution to Eq. (52) can be found if we suppose that the left‐hand side is equal to zero. Then, we have $\overline{\rho }=\left(\text{1}/\sqrt{\pi \overline{\theta }}\right)e^{-x^{\text{2}}/\overline{\theta }}$ , where $\overline{\theta }=\text{2}{\eta }^{\text{2}}/\gamma$ . We assume that the solution to Eq. (48) is of the form

$\rho (t)=\frac{\text{1}}{\sqrt{\pi \theta (t)}}{e}^{-{\left(x-⟨x(t)⟩\right)}^{\text{2}}/\theta (t)}.$ E53

By substituting (53) into (52) one can be convinced that if $\theta \left(t\right)$ and $⟨x(t)⟩$ satisfy Eq. (48) and Eq. (50), respectively, then Eq. (53) is the solution to the Fokker‐Planck equation (52). As a result, any initial data that is normally distributed will be attracted to the steady solution of Eq. (52), which is also normally distributed. The dissipation parameter $\gamma$ determines the rate at which PDF $\rho$ approaches $\overline{\rho }$ . The auto‐correlation function for the stationary stochastic process (46) can be written as

$C\left(\tau \right)=\frac{\text{2}{\eta }^{\text{2}}}{\gamma }{e}^{-\gamma \tau }\equiv \overline{\theta }{e}^{-\gamma t}\text{.}$ E54

Thus, the potential predictability of system (46) can be characterized by the auto‐correlation function of the process $x(t)$ and, therefore, the convergence of $\rho (t)$ to $\overline{\rho }$ can be explored based only on function $C\left(\tau \right)$ with time lag $\tau$ . This conclusion is valid for the set of multi‐dimensional differential equations [16]. In this case, however, the covariance matrix is used instead of the auto‐correlation function. It is very important that for climate models the convergence of the covariance matrix $C(t)$ to the covariance matrix of stationary process $\overline{C}$ is defined only by climatological values of climate model variables. As a result, potential predictability is also determined by climatological data.

Generally, the potential predictability can be defined as the convergence time of initial distribution to the equilibrium one. To quantify the rate of convergence of one‐dimensional distributions to the equilibrium ones, the concept of entropy can be used. If the information entropy $S={\displaystyle \int \rho \text{ln}\rho d\alpha }$ is taken as a measure of predictability, then for the Gaussian distribution $\rho =\left(\text{1}/\sqrt{\text{2}\pi {\sigma }^2}\right)e^{-{\left(\alpha -\overline{\alpha }\right)}^{\text{2}}/(\text{2}{\sigma }^{\text{2}})}$ information entropy can be expressed as $S=\text{ln}{\sigma }^{\text{2}}+C$ . It can be shown that the variance and, therefore, the entropy are directly dependent on the Lyapunov exponents [39]. To study the predictability of climate system, the relative entropy $S_r={\displaystyle \int \rho \text{ln}\left(\rho /\overline{\rho }\right)d\alpha }$ , where $\overline{\rho }$ is an equilibrium PDF, is a more suitable measure [40]. Relative entropy is invariant with respect to nonlinear transformations of $\alpha$ and $\rho \to \overline{\rho }$ at $t\to \infty$ .

4.2. Predictability of the second kind

Predictability of the second kind relates to the predictability of changes in climate system caused by infinitesimal perturbations in the parameters that describe the external forcing. Climate prediction does not involve forecasting weather conditions at either a certain geographical region or globally. On the contrary, climate prediction aims to forecast statistics of the climate system averaged over sufficiently long period of time. So, we are interested in how external perturbations affect certain aspects of climate statistic, such as the first $\overline{x}$ (mean) and/or second ${\sigma }_x^{\text{2}}$ (variance) moments. One of the most important problems in the exploration of predictability of the second kind is to distinguish the response signal of the climate system to perturbed external forcing from the noise in the model output results. The signal‐to‐noise ratio can be used to make the conclusion with respect to the usefulness of the obtained climate system response. Thus, the predictability of the second kind is mathematically reduced to finding the response function of the climate system model [39].

Consider the following finite‐dimensional dynamical system that is controlled by some external forcing $f$ (e.g. the concentration of carbon dioxide in the atmosphere):

$dx/dt=F\left(x\right)+f,x{|}_{t=\text{0}}=x_{\text{0}},$ E55

Suppose that system (55) possesses the attractor $A$ and let $\mu$ be its invariant measure. The behaviour of this system will be explored on the attractor $A$ . Since system (55) a priori possesses the property of ergodicity, its statistical characteristics are calculated by averaging along a single, sufficiently long, random trajectory. Thus, the average state $⟨x⟩$ and variance $⟨{\sigma }_x^{\text{2}}⟩$ of system (55) are defined, respectively, as

$⟨x⟩=\underset{T\to \infty }{{\displaystyle \lim }}\frac{\text{1}}{T}{\displaystyle \underset{\text{0}}{\overset{T}{\int }}x\left(t\right)dt}={\displaystyle \underset{A}{\overset{}{\int }}xd\mu ,}⟨{\sigma }_x^{\text{2}}⟩={\displaystyle \underset{A}{\overset{}{\int }}}(x-⟨x⟩{)}^{\text{2}}d\mu .$ E56

Let system (55) be perturbed by an infinitesimal disturbance in the external forcing $\delta f$ such that $\delta f\ll f$ :

$d{x}^{*}/dt=F\left(x^{*}\right)+f+\delta f\text{.}$ E57

For this system $⟨x^{*}⟩={\displaystyle \underset{A}{\overset{}{\int }}{x}^{*}d{\mu }^{*}}$ and $⟨{\sigma }_x^{\text{2}}⟩={\displaystyle \underset{A}{\overset{}{\int }}{\left(x^{*}-⟨x^{*}⟩\right)}^{\text{2}}d{\mu }^{*}}$ . Let us introduce the new variable $x^{\text{'}}(t)=x(t)-x^{*}(t)$ . Assuming that $\Vert x^{\text{'}}\Vert$ is rather small then, combining (55) and (56), one can obtain the following linear equation for variable $x^{\text{'}}$ :

$d{x}^{\text{'}}/dt=J\left(x\right)x^{\text{'}}+\delta f\text{.}$ E58

where $J\left(x\right)=\partial F/\partial x$ is the Jacobian. Let $\delta f$ be a staircase function that is activated at $t=\text{0}$ then the solution to Eq. (58) can be written in terms of the Green's function:

$x^{\text{'}}(t)={\displaystyle \underset{\text{0}}{\overset{t}{\int }}G\left(t,t^{\text{'}}\right)\delta f\left(t^{\text{'}}\right)d{t}^{\text{'}}.}$ E59

The operator $\mathrm{R}={\displaystyle \underset{\text{0}}{\overset{t}{\int }}G\left(t,t\text{'}\right)dt\text{'}}$ is a sought‐for response function (operator). If at $t=\text{0}$ the distribution of initial states is identical for both unperturbed (55) and perturbed (57) systems, then one can calculate the average response operator:

$⟨\mathrm{R}⟩={\displaystyle \underset{\text{0}}{\overset{t}{\int }}⟨G\left(t,t^{\text{'}}\right)⟩dt\text{'}}={\displaystyle \underset{\text{0}}{\overset{t}{\int }}G\left(t-t\text{'}\right)d\left(t-t\text{'}\right)}\text{.}$ E60

By averaging both sides of Eq. (59), one can get the following linear equation to calculate the system's response to the external forcing:

$⟨x^{\text{'}}⟩=⟨\mathrm{R}⟩\delta f.$ E61

Suppose that system (55) is regular, i.e. for this system the quadratic conservation law is valid and system itself satisfies the Liouville equation for incompressibility in the phase space. Assume also that the system is in equilibrium. Taking into consideration the fluctuation dissipation theorem [41], the average impulse response operator of the regular system in equilibrium is expressed via system's statistics:

$⟨G\left(t,t^{\text{'}}\right)⟩=G\left(t-t\text{'}\right)=C(t-t^{\text{'}})C^{-\text{1}}(\text{0}),$ E62

where $C\left(t-t\text{'}\right)=⟨x(t)x^{{\rm T}}(t^{\text{'}})⟩$ is the system's auto‐correlation matrix with time lag $\tau =t-t\text{'}$ . Now we can combine (60) and (62) to get the following well‐known formula [42]:

$⟨x^{\text{'}}⟩={\displaystyle \underset{\text{0}}{\overset{\infty }{\int }}C(t)C^{-\text{1}}(\text{0})dt\cdot \delta f.}$ E63

Thus, the mean response of climate system to external forcing is determined by observations of unperturbed climate oscillation.

5. Concluding remarks

The prediction of climate change caused by natural processes and human‐induced drivers is one of the most critical scientific issues facing the mankind in the 21st century. Computer‐simulated climate models represent a very powerful and, perhaps, the only research instrument for studying climate and its dynamics. One of the key components of climate models, namely the model of the atmospheric general circulation, currently also serves as a primary tool for the numerical weather prediction all around the globe. However, the climate (atmospheric) system's trajectory calculated via numerical integration of multi‐dimensional partial differential equations that describe the climate (atmospheric) system evolution is unstable with respect to both perturbations (errors) in the initial conditions and infinitesimal external forcing expressed by some model parameters and/or boundary conditions. This instability limits the time horizon of the validity of the climate (weather) forecast and leads to predictability problem.

In this chapter, the climate system is considered as a complex adaptive dynamical system that possesses a number of specific properties such as, for example, dissipativity, nonlinearity and chaoticity. From this perspective, the climate predictability problem is best discussed and analysed by formally examine two kinds of predictability. The first kind of predictability refers to the initial value problem (estimating the impact of perturbations in the initial conditions on the forecast skill), while the second kind of predictability relates to the boundary value problem (estimating the impact of external forcing on the system's behaviour).