Other tools from statistical mechanics and nonequilibrium thermodynamics have, therefore, to be considered (Bouchet and Sommeria, 2002; Chavanis and Sommeria, 1996; Farrell and Ioannou, 1996; Lucarini et al., 2014; Lucarini and Sarno, 2011; Majda and Wang, 2006; Robert and Sommeria, 1991; Trefethen et al., 1993), and will be discussed in Sects. IV and V. Combining these tools with those of the successive-bifurcation approach could lead to a more complete physical characterization of gyre modes in realistic models. Preparing the ground for combining dynamical-systems tools and statistical-physics tools in this way is one of the main purposes of our review paper.

### G. Multiple Scales: Stochastic and Memory Effects

In Sect. II.B, we have illustrated in Fig. 12 the multiplicity of time scales present in the climate spectrum. We also pointed out that this multiplicity of scales gives rise to the need for a hierarchy of models that enable the study of separate scales of motion and of the phenomena associated with each, as well as of the interactions between two or more scales, cf. Fig. 11. In this subsection, we discuss several ways in which one can address these issues, using the theory of stochastic processes and taking into account non-markovian effects.

#### 1. Fast Scales and Their Deterministic Parametrization

Let us concentrate, for the sake of definiteness, on variability within a particular range of frequencies $f$ in the climatic power spectrum of Fig. 12a, say seasonal to centennial, i.e., $10^{-2} \text{ yr}^{-1} = f_1 \le f \le f_2 = 10^0 \text{ yr}^{-1}$, and write a model of this variability as

$$ \dot{\mathbf{z}} = \mathbf{H}(\mathbf{z}; \boldsymbol{\mu}). \qquad (23) $$

We have seen that oscillatory modes of both the THC, cf. Sect. III.D, and of the wind-driven circulation, cf. Sec. III.F, lie in this range. How should one then take into account the slower time scales to the left of this range, $f_0 < f < f_1$, and the faster ones to the right, $f_2 < f < f_3$, where $0 \le f_0 < 10^{-2} \text{ yr}^{-1}$ and $1 < f_3 < \infty$?

A time-honored approach in physical modeling is to describe variability to the left as a prescribed (slow) evolution of parameters,

$$ \boldsymbol{\mu} = \boldsymbol{\mu}(\epsilon t), \quad 0 < \epsilon \ll 1, \qquad (24) $$

where $\epsilon$ is small and $\epsilon t$ is, therefore, a slow time. One might, for instance, consider a particular $\mu = \mu(\epsilon t)$ in Eq. (24) above to represent a slow change in the solar constant or in the height of the topography.

To the right, one might approximate the more rapid fluctuations as being very fast with respect to those of main interest or even infinitely fast. There are two distinct approaches based on these ideas: the first one is

purely deterministic, the second one introduces a noise process and thus stochastic considerations.

The standard slow-fast formulation of a system of differential equations — when assuming a large but finite separation of time scales — is given, for $\epsilon \neq 0$, by

$$ x' = F(x, y; \epsilon), \qquad (25a) $$

$$ y' = \epsilon G(x, y; \epsilon). \qquad (25b) $$

Here $\mathbf{z} = (x, y)^T$, with $x$ the fast and $y$ the slow variable, and $\mathbf{H} = (F, G)^T$, while $(\cdot)^T$ designates the transpose and $(\cdot)' = d(\cdot)/dt$ stands for differentiation with respect to the fast time $t$.

As long as $\epsilon \neq 0$, system (25) is equivalent to the so-called slow system

$$ \epsilon \dot{x} = F(x, y; \epsilon), \qquad (26a) $$

$$ \dot{y} = G(x, y; \epsilon), \qquad (26b) $$

in which the dot stands for differentiation $d(\cdot)/d\tau$ with respect to the slow time $\tau = \epsilon t$.

The classical way of dealing with such problems has been matched asymptotic expansions (e.g., Grasman, 1987; Lagerstrom, 1988). This methodology arose originally from dealing with boundary layers in fluid dynamics, with the inner problem refering to the fast variations in the boundary layer, while the outer problem refers to the more slowly varying free flow outside this layer.

More recently, a point of view inspired by dynamical systems theory (Fenichel, 1979) has led to geometric singular perturbation theory (e.g., Jones, 1995). In this approach, one considers the invariant manifolds that arise in the two complementary limits obtained by letting $\epsilon \to 0$ in the fast and the slow system, respectively.

In the fast system (25), the limit is given by

$$ x' = F(x, y; 0), \qquad (27a) $$

$$ y' = 0, \qquad (27b) $$

while in the slow one, it only makes sense if the right-hand side of (26a) is identically zero; if so, the latter limit is given by

$$ 0 = F(x, y; 0), \qquad (28a) $$

$$ \dot{y} = G(x, y; 0). \qquad (28b) $$

The algebraic equation $F(x, y; 0) = 0$ defines the critical manifold $M_c$ on which the solutions of the reduced problem $\dot{y} = G(x, y; 0)$ evolve; here $x = X_F(y)$ are the explicit solutions of the implicit equation (28a).

The splitting of the full, slow-fast system given by either Eqs. (25) or (26) into the two systems (27) and (28) has proven very helpful in the study of relaxation oscillations (e.g., Grasman, 1987). We saw such sawtooth-shaped, slow-fast oscillations arise in either the THC (Sect. III.C.2) or the wind-driven circulation of the oceans (Sect. III.F.2).
Another important application of this methodology is in the reduction of large multiscale problems to much smaller ones. In the systems (25) and (26), we considered both $x$ and $y$ to be scalar variables. We saw in Fig. 11, though, that the characteristic spatial and temporal scales of atmospheric, oceanic and coupled climate phenomena are highly correlated with each other; large-scale motions tend to be slow and the smaller-scale ones faster. Thus, it is much more judicious to consider $\mathbf{z} = (\mathbf{X}^T, \mathbf{Y}^T)^T$, with $\mathbf{X} \in \mathbb{R}^m$, $\mathbf{Y} \in \mathbb{R}^n$, and $m \gg n$, i.e., the number of small and fast degrees of freedom much larger than that of the large and slow ones.

This set-up corresponds conceptually to the parametrization problem, which we defined in Sect. II.B as finding a representation of the unresolved subgridscale processes described by $\mathbf{X} \in \mathbb{R}^m$ in terms of the resolved, larger-scale ones described by $\mathbf{Y} \in \mathbb{R}^n$. A paradigmatic example is that of parametrizing cloud processes, with spatial scales of 1 km and smaller and with temporal scales of one hour and less, given the large-scale fields characterized by lengths of tens and hundreds of kilometers and by durations of substantial fractions of a day and longer. In this case, the critical manifold appears to be $S$-shaped, as for the periodically forced Van der Pol oscillator (e.g., Guckenheimer et al., 2003, Fig. 2.1), with jumps that occur between the branch on which convection, and hence rain, is prevalent, and the one on which the mean vertical stratification is stable, and thus no rain is possible. Next, we discuss specifically the parametrization of convective processes and of clouds in this perspective.

## 2. An Example: Convective Parametrization

Clouds have a dramatic role in climate modeling and in determining the climate's sensitivity to natural and anthropogenic perturbations (e.g., IPCC, 2001, 2014a). A substantial literature exists, therefore, on cloud observations, modeling and simulation (e.g., Emanuel, 1994, and references therein). See our earlier discussion in Sects. II.C.2 and III.G.1. One of the oldest, best known and most widely used cumulus parametrizations is the Arakawa and Schubert (1974) (AS) one. Cumulus convection occurs due to moist convective instability, which converts the potential energy of the large-scale mean state into the kinetic energy of the cumulus clouds. A fundamental parameter in this process is the fractional entrainment rate $\lambda$ of a cumulus updraft. In the AS parametrization of moist atmospheric convection, the key idea is that an ensemble of cumulus clouds is in quasi-equilibrium with the large-scale environment.

The cloud work function $A(\lambda)$ changes in time according to

$$ \dot{A}(\lambda) = J \otimes M_B(\lambda) + F(\lambda); \quad (29) $$

here $M_B(\lambda)$ is the nonnegative mass flux through the cloud base, and $J \otimes M_B$ is a weighted average over cloud types, with $J$ standing for the weights in the averaging integral, while $F$ is the large-scale forcing.

The dot stands for differentiation with respect to the slow time $\tau = \epsilon t$, as in Eq. (26a). The quasi-equilibrium assumption in the AS parametrization corresponds simply to the critical manifold equation Eq. (28a) above, which becomes

$$ 0 = J \otimes M_B(\lambda) + F(\lambda). \quad (30) $$

In this case, the small parameter that corresponds to the $\epsilon$ of the general slow-fast formulation above is the reciprocal of the adjustment time $\tau_{adj}$ of a cloud ensemble to the mean state, $\epsilon \sim 1/\tau_{adj}$.

Pan and Randall (1998) proposed an equation that corresponds to the behavior of a cumulus ensemble off the critical manifold given by (30), which they termed a prognostic closure. In their formulation, one computes a cumulus kinetic energy $K$ from

$$ \dot{K} = B + S - D; \quad (31) $$

here $B$ is the buoyancy production term, $S$ the shear production term, and $D$ the vertically integrated dissipation rate. The main parameters on which the behavior of the slow-fast system given by Eqs. (30, 31) depends are

$$ \alpha = M_B^2/K, \quad \tau_D = K/D. \quad (32) $$

While Pan and Randall (1998) do not determine $\alpha$, $\tau_D$ or $\tau_{adj}$ explicitly, they provide qualitative arguments based on the physics of cumulus convection that make the quasi-equilibrium limit plausible.

## 3. Stochastic Parametrizations

It is of broader interest, though, to consider now the slow-fast formulation of a system of differential equations in the case of infinite separation of time scales, i.e., when the fast motions have infinite frequency or, more precisely, zero decorrelation time. In this case, one has to introduce a white-noise process and the associated stochastic considerations.

The basic idea relies on the Einstein (1905) explanation of Brownian motion, in which a large particle is immersed in a fluid formed of many small ones. Let the large particle move along a straight line with velocity $u = u(t)$, subject to a random force $\eta(t)$ and to linear friction $-\lambda u$, with coefficient $\lambda$. The equation of motion is

$$ du = -\lambda u dt + \eta(t). \quad (33) $$

The random force $\eta(t)$ is assumed to be a “white noise,” i.e., it has mean zero $\mathcal{E}[\eta(t; \omega)] = 0$ and autocorrelation $\mathcal{E}[\eta(t; \omega)\eta(t+s; \omega)] = \sigma^2\delta(s)$, where $\delta(s)$ is a Dirac function, $\sigma^2$ is the variance of the white-noise process, $\omega$
labels the realization of the random process, and $\mathcal{E}$ is the expectation operator, which averages over the realizations $\omega$. Alternative notations for the latter are the overbar, in climate sciences, and the angle brackets, in quantum mechanics, $\mathcal{E}[F] := \bar{F} := \langle F \rangle$.

Equation (33), with $\eta = \sigma dW$, is a linear stochastic differential equation (SDE) of a form that is now referred to as a Langevin equation, where $W(t)$ is a normalized Wiener process, also called Brownian motion. It was introduced into climate dynamics by Hasselmann (1976), who identified slow, “climate” changes with the motion of the large particle and fast, “weather” fluctuations with the motions of the small fluid particles. He also thought of weather as associated with the atmosphere and climate with the ocean, cryosphere and land vegetation.

Specifically, Hasselmann (1976) assumed that, in a system like (23), and without formally introducing the time-scale separation parameter $\epsilon$, one would have $\tau_x \ll \tau_y$, where $\tau_x$ and $\tau_y$ are the characteristic times of the fast x- and slow y-variables, respectively. From this assumption, and relying also upon the results of Taylor (1921), he derived then a linear SDE for the deviations $y'_j$ of the slow variables $\mathbf{Y}$ from a reference state $\mathbf{Y}_0$, and the properties of the corresponding covariance matrix and spectral densities. In particular, the red-noise character of the spectrum $S = S(f)$, with $S \sim f^{-2}$ for many oceanic observed time series, gave considerable credence to the thermal-flywheel role that Hasselmann (1976) attributed to the ocean in the climate system.

In the light of recent mathematical results on very large time scale separation in slow-fast deterministic systems, let us consider a relatively simple — but still sufficiently relevant and instructive — case, in which the reduction to an SDE can be rigorously derived, cf. Pavliotis and Stuart (2008) and Melbourne and Stuart (2011). Their deterministic system of ordinary differential equations (ODEs) is a small modification of Eq. (26), namely

$$ \begin{aligned} \dot{x} &= \epsilon^{-1} f_0(y^{(\epsilon)}) + f_1(x^{(\epsilon)}, y^{(\epsilon)}), \quad x^{(\epsilon)}(0) = x_0, \quad (34a) \\ \dot{y}^{(\epsilon)} &= \epsilon^{-2} g(y^{(\epsilon)}), \quad y^{(\epsilon)}(0) = y_0, \end{aligned} \quad (34b) $$

where $x^{(\epsilon)} \in \mathbb{R}^d$ and $y^{(\epsilon)} \in \mathbb{R}^{d'}$.

The formal difference with respect to the situation studied in the previous subsection is that $F(x, y; \epsilon)$ of (26a) has been expanded in $\epsilon$ as $F(x, y; \epsilon) = f_0(y^{(\epsilon)}) + \epsilon f_1(x^{(\epsilon)}, y^{(\epsilon)})$, while $G(x, y; \epsilon)$ of (26b) has been both simplified, in becoming $x$-independent, and “accelerated,” to read $G(x, y; \epsilon) = \epsilon^{-2}g(y^{(\epsilon)})$. The basic idea is that the chaotic and fast $y^{(\epsilon)}$ induces, as $\epsilon \to 0$, a white-noise driving of the slow $x$. Note that one needs $d' \ge 3$ for the autonomous Eq. (34b) to have chaotic solutions.

Melbourne and Stuart (2011) assume merely that the fast equation (34b) has a compact attractor $\Lambda \in \mathbb{R}^{d'}$, which supports an invariant measure $\mu$, and that $\mathcal{E}_\mu f_0(x) = 0$, along with certain boundedness and regularity conditions on $f(x, y)$. They then show rigorously that $x^{(\epsilon)}(t) \to_p X(t)$ as $\epsilon \to 0$, where the convergence is

with respect to the appropriate probability measure, and $X(t)$ is the solution of the SDE

$$ X(t) = x_0 + \int_0^t \bar{F}(X(s))ds + \sigma W(t). \quad (35) $$

Here $W$ is the Brownian motion with variance $\sigma^2$, such that the white noise in Eq. (33) be given by $\eta = \sigma dW$. Moreover, $\bar{F}(X) = \mathcal{E}_\mu F(x, \cdot)$ with respect to $\mu$.

Fundamental mathematical issues associated with the above diffusive limit of slow-fast systems were explored by Papanicolaou and Kohler (1974) and early results in the physical literature include Beck (1990) and Just et al. (2001). Many aspects of the applications to climate modeling are covered in Palmer and Williams (2009). More specifically, Berner et al. (2017) and Franzke et al. (2015) discuss issues of stochastic parametrization of subgrid-scale processes.

To conclude this subsection, it is of interest to consider, in a broader perspective, the potential for a unified theory of nonautonomous dynamical systems, in which the fast processes may be modeled as either deterministic or stochastic. The theory of purely deterministic, skew-product flows was laid on a solid basis by Sell (1971, and references therein) and, more recently, by Kloeden and Rasmussen (2011). Random dynamical systems are extensively covered by Arnold (1998), with many recent results in an active field.

Berger and Siegmund (2003) point out that “Quite often, results for random dynamical systems and continuous skew product flows are structurally similar,” and thus open the way to a unified theory. They outline both commonalities and distinctions between the two broad classes of nonautonomous dynamical systems, in order to shed further light on existing results, as well as stimulate the development of common concepts and methods.

Caraballo and Han (2017) provide a solid and accessible introduction to random dynamical systems, as well as to deterministically nonautonomous ones, along with several interesting applications. They also consider the two distinct types of formulation of the deterministic ones, the so-called process formulation and the skew-product flow one. System (34) above, for instance, is a particular case of a master-slave system

$$ \dot{x} = f(x, y), \quad \dot{y} = g(y), \quad \text{with } x \in \mathbb{R}^d, y \in \mathbb{R}^{d'} \quad (36) $$

that induces a skew-product flow, where $y(t)$ is the driving force for $x(t)$.

These developments bear following, since the climate sciences offer a rich source of relevant problems and could thus lead to novel and powerful applications of the unified theory. Caraballo and Han (2017) already studied the Lorenz (1984) model, in which seasonal forcing acting on deterministic subseasonal variability can induce interannual variability. Other applications will be discussed in the next sections, especially in Sect. IV.E.
## 4. Modeling Memory Effects

While the study of differential equations goes back to Isaac Newton, the interest for including explicitly delays into evolution equations that govern physical and biological processes is relatively recent. A mid-20th-century reference is Bellman and Cooke (1963), followed by Driver (1977) and Hale (1977). Delay-differential equations (DDEs) were introduced into the climate sciences by Bhattacharya et al. (1982) and have been used widely in studying ENSO (e.g., Ghil et al., 2015; Tziperman et al., 1994a, and references therein).

As we shall see, memory effects can play a key role when there is little separation between scales, in contrast to the assumptions of Hasselmann (1976) and of other authors mentioned in the two immediately preceding subsections. Moreover, when properly incorporated in the mathematical formulation of the climate problem at hand, relying on such effects can lead to highly efficient and accurate model reduction methods.

a. *The Mori-Zwanzig formalism.* In statistical physics, the Mori (1965)-Zwanzig (1961) (MZ) formalism arose from describing the interaction of a Hamiltonian many-particle system with a heat bath. Today, though, it is being used in a large number of applications that include dissipative systems.

The fundamental idea is illustrated by the following very simple, linear system of two ODEs (E and Lu, 2011):

$$ \dot{x} = a_{11}x + a_{12}y, \quad (37a) $$

$$ \dot{y} = a_{21}x + a_{22}y. \quad (37b) $$

The sole assumption is that we are interested in the details of the behavior of $x(t)$ but only in the statistics of $y$; naturally, one thinks of $y$ as fluctuating faster than $x$ but this is not actually required for the formalism outlined below to work. One considers $x$ as a (slowly varying) parameter in solving Eq. (37b) for $y$ by using the variation-of-constants formula

$$ y(t) = e^{a_{22}t}y(0) + \int_{0}^{t} e^{a_{22}(t-s)}a_{21}x(s)ds, $$

and plugs this result back into Eq. (37a), to yield

$$ \dot{x} = a_{11}x + \int_{0}^{t} K(t-s)x(s)ds + f(t). \quad (38) $$

Equation (38) is a generalized Langevin equation (GLE), in which $K(t) = a_{12} \exp(a_{22}t)a_{21}$ is the memory kernel and $f(t) = a_{12} \exp(a_{22}t)y(0)$ is the noise term, since one thinks of $y(0)$ as randomly drawn from the rapidly fluctuating $y(t)$. The essential difference with respect to Eq. (33) is the convolution integral in (38),

which expresses the delayed action of the slow variable $x$ on the fast variable $y$.

The MZ formalism consists — in a general, nonlinear set of Markovian evolution equations with a large or even infinite number of degrees of freedom — in selecting the variables one is interested in via a projection operator, and deriving the generalized form of the GLE above. Examples of Markovian evolution equations are systems of ODEs or PDEs for which an instantaneous initial state carries all the information from the past.

In this general case, the memory term involves repeated convolutions between decaying memory kernels and the resolved modes, and the GLE is therewith a non-Markovian, stochastic integro-differential system that is very difficult to solve without further simplifications. Among the latter, the short-term memory approximation (e.g., Chorin and Stinis, 2007) posits rapidly decaying memory and is equivalent to assuming a relatively large separation of scales, as in Sects. III.G.1 and III.G.3 above.

Fortunately, Kondrashov et al. (2015) have shown that there is a way to approximate the GLE in a very broad setting, efficiently and accurately, by a set of Markovian SDEs without the need for pronounced scale separation and in the presence of so-called intermediate-range memory. This way relies on a methodology that was developed at first quite independently of the MZ formalism, namely empirical model reduction (EMR: Kondrashov et al., 2005b; Kravtsov et al., 2005, 2009).

b. *Empirical model reduction (EMR) methodology.* The purpose of EMR development was deriving relatively simple nonlinear, stochastic-dynamic models from time series of observations or of long simulations with high-end models, such as GCMs. An EMR model can be compactly written as

$$ \dot{x} = -Ax + B(x, x) + L(x, r_t^l, \xi_t, t), \quad 0 \le l \le L-1. \quad (39) $$

Here $x$ typically represents the resolved and most energetic modes. Most often, these are chosen by first selecting a suitable basis of empirical orthogonal functions (EOFs: Preisendorfer, 1988) or other data-adaptive basis (e.g., Kravtsov et al., 2009, and references therein), and retaining a set of principal components that capture a satisfactory fraction of the total variance in the data set.

The terms $-Ax$ and $B(x, x)$ represent, respectively, the linear dissipation and the quadratic self-interactions of these modes, while $L$ is a time-dependent operator that is bilinear in the resolved variables $x$ and the unresolved ones $r_t^l$. These interactions take a prescribed form in the EMR formulation of Eq. (39), and arise by integrating recursively—from the lowest level $L$ to the top level, $l = 0$—the “matrioshka” of linear SDEs

$$ d\mathbf{r}_t^l = \mathbf{M}_l(\mathbf{x}, \mathbf{r}_t^0, \dots, \mathbf{r}_t^l)d + \mathbf{r}_t^{l+1}dt. $$
At each level $l$, the coupling between the variable $\mathbf{r}_l^l$ and the previous-level variables $(\mathbf{x}, \mathbf{r}_l^0, \dots, \mathbf{r}_l^{l-1})$ is modeled by $L-1$ rectangular matrices $\mathbf{M}_l$ of increasing order.

In practice, the matrices $\mathbf{A}, \mathbf{M}_l$, and the quadratic terms $\mathbf{B}$ are estimated by a recursive least-square procedure, which is stopped when the $L$th-level residual noise $\mathbf{r}_{l-1}^L = \xi_l$ has a lag-1 vanishing autocorrelation. The stochastic residuals $\mathbf{r}_l^l$, obtained in this recursive minimization procedure, are ordered in decreasing order of decorrelation time, from $\mathbf{r}_l^0$ to $\xi_l$.

Note that the integral terms arising in the $\mathbf{L}$ operator are convolution integrals between the macro-state variables $\mathbf{x}$ and memory kernels that decay according to the dissipative properties of the matrices $\mathbf{M}_l$. These decay times are not necessarily short and one can thus treat the case of intermediate-range memory, in the MZ terminology. Furthermore, cf. Kondrashov et al. (2005b), note that any external forcing, such as the seasonal cycle, can be typically introduced as a time dependence in the linear part $\mathbf{A}$ of the main level of Eq. (39).

### c. Role of memory effects in EMR.
We propose here a simple analytic example that should help understand the general description of EMR in the previous paragraphs, as well as the connection to the MZ formalism. The model is given by

$$dx = (f(x) + r)dt, \quad (40a)$$

$$dr = (\gamma x - \alpha r)dt + dW_t, \quad (40b)$$

where $f$ is a nonlinear function, $W_t = W_t(\omega)$ is a standard Wiener process, as in Eq. (33), $\alpha > 0$ and $\gamma$ is real. We are interested in $x$, which is slow, and want to parameterize $r$, which is fast.

Proceeding for (40) as we did for the system (37), but now for a fixed realization $\omega \in \Omega$, we get

$$r(t; \omega) = e^{-\alpha t} r_0 + \gamma \int_0^t e^{-\alpha(t-s)} x(s) ds + W_t(\omega),$$

where $r(t_0; \omega) = r_0$. Substituting this back into (40a) yields the following randomly forced integro-differential equation,

$$\dot{x} = f(x) + e^{-\alpha t} r_0 + \gamma \int_0^t e^{-\alpha(t-s)} x(s) ds + W_t(\omega), \quad (41)$$

which is the analog of Eq. (38) in the stochastic-dynamic context of (40).

Kondrashov et al. (2015) proved rigorously that the EMR (39) is equivalent to a suitable generalization of the GLE (41). Pavliotis (2014) defines as quasi-Markovian a stochastic process where the noise in the GLE can be described by adding a finite number of auxiliary variables, as in the case of the EMR approach. Thus, the EMR methodology can be seen as an efficient implementation of the MZ formalism, i.e., an efficient solution of

the associated GLE, even in the absence of large-scale separation. This result explains the remarkable success of EMR in producing reduced models that capture the multimodality as well as the nontrivial power spectrum of phenomena merely known from time series of observations or of high-end model simulations.

In the remainder of this subsection, we give two examples of this success and further references to many more. Alternative approaches to efficient solutions of the GLE can be found, for instance, in Chorin et al. (2002).

### d. EMR applications.
We choose here an EMR model to simulate Northern Hemisphere mid-latitude flow (Kravtsov et al., 2005) and a real-time ENSO prediction model (Kondrashov et al., 2005b). Further examples of successful application of the methodology appear in Kravtsov et al. (2009, 2011) and elsewhere.

Kravtsov et al. (2005) introduced the EMR methodology and illustrated it at first with quadratically nonlinear models of the general form given here in Eq. (39). The applications were to the Lorenz (1963) convection model, the classical double-well potential in one space dimension, and a triple-well potential in two dimensions with an exponential shape for the wells. More challenging was a real-data application to geopotential height data for 44 boreal winters (1 December 1949–31 March 1993). The data set consisted of 44 × 90 = 3960 daily maps of winter data, defined as 90-day sequences starting on 1 December of each year. The best EMR fit for the data required the use of nine principal components and of $L = 3$ levels.

The probability density functions (PDFs) for the observed and the EMR model-generated data sets are plotted in Fig. 32. The EMR clearly captures quite well the three modes obtained with a Gaussian mixture model, cf. Smyth et al. (1999) and Ghil and Robertson (2002). These three modes correspond to three clusters found by Cheng and Wallace (1993) using very different methods on a somewhat different data set. The maps of the corresponding centroids appear as Fig. 1 in Ghil and Robertson (2002) and are discussed therein; they agree quite well with those of Cheng and Wallace (1993), cf. Smyth et al. (1999, Fig. 9).

Kondrashov et al. (2005b) fitted the global SST field between 30°S–60°N over the time interval January 1950–September 2003 by using linear and quadratic EMR models with one and two noise levels, $L = 1, 2$, based on monthly SST anomaly maps and allowing a seasonal dependence of the dissipative terms in Eq. (39). Their results when using $L = 2$ were much better, for either a linear or a quadratic model, which clearly shows the role of memory effects in EMR modeling and the connection with the MZ formalism that was explained above.

The use of the EMR models in prediction was tested by so-called hindcasting or retrospective forecasting, i.e., a protocol — also called “no look-ahead” — in which the
<figure><img src="image_41.png" /><figcaption>FIG. 32: Multimodal probability density function (PDF) for the Northern Hemisphere's geopotential height anomalies of 44 boreal winters; see Smyth et al. (1999) for details of the data set and of the mixture model methodology for computing the PDFs. (a) PDF of the observed height anomalies; and (b) of the anomalies given by the EMR model. Modified from Kravtsov et al. (2005), with the permission of the American Meteorological Society.</figcaption></figure>

data available past a certain time instant are eliminated when constructing the model to be used in the forecast to be evaluated. The results of these tests are plotted in Fig. 33, for $L = 2$ and a linear vs. a quadratic model. The light-black rectangle in the Eastern Tropical Pacific corresponds to the region (5°S–5°N, 150°–90°W) over which SST anomalies are averaged to obtain the Niño-3 index, already mentioned in Sect. III.E.2.

Climate forecast skill is measured mainly via root-mean-square errors and anomaly correlations. The former skill scores are given in Kondrashov et al. (2005b, Fig. 2d) and clearly indicate the superiority of the quadratic model. The latter apeear in Fig. 33 here, too.

Anomaly correlations of roughly 0.5–0.6 (red in the figure) are considered quite useful in climate prediction, and the red area is substantially larger in panel (b) here, covering most of the Tropical Pacific and Indian Oceans. Concerning actual real-time forecasts, Barnston et al. (2012) have found that, over the 2002–2011 interval of their evaluation, the EMR-based forecasts of UCLA’s Theoretical Climate Dynamics (TCD) group were at the top of the eight statistical models being evaluated and exceeded in skill by only a few of the 12 dynamical, high-end models in the group that participated in the IRI ENSO Forecast plume https://iri.columbia.edu/our-expertise/climate/forecasts/enso/current/?enso_tab=enso-sst_table.

e. Explicit derivation of the parametrized equations. The EMR methodology is data-driven, so that it allows to construct bottom-up an effective dynamics able to account for the observed data. Wouters and Lucarini (2012, 2013), showed that the MZ formalism can be used also to derive the parametrizations in a top-down manner in a rather general way. Assume the system of interest is

<figure><img src="image_42.png" /><figcaption>FIG. 33: Validation of forecast skill of the EMR-based ENSO model of Kondrashov et al. (2005b) via anomaly correlation maps of SSTs at 9-month lead time. (a) Using a quadratic EMR model with $L = 1$; and (b) a quadratic model with $L = 2$. Modified from Kondrashov et al. (2005b), with the permission of the American Meteorological Society.</figcaption></figure>

described by the following evolution equations:

$$ \dot{\mathbf{X}} = f_{\mathbf{X}}(\mathbf{X}) + \epsilon \Psi_{\mathbf{X}}(\mathbf{X}, \mathbf{Y}), \quad (42a) $$

$$ \dot{\mathbf{Y}} = f_{\mathbf{Y}}(\mathbf{Y}) + \epsilon \Psi_{\mathbf{Y}}(\mathbf{X}, \mathbf{Y}); \quad (42b) $$

here, again, $\mathbf{X} \in \mathbb{R}^m$, $\mathbf{Y} \in \mathbb{R}^n$, $\mathbf{X}$ is the set of large-scale, energetic variables of interest and, typically, $m \gg n$.

We assume, furthermore, that $\epsilon$ is a small parameter describing the strength of the coupling between the two sets of variables, and that, if $\epsilon = 0$, the dynamics is chaotic for both the $\mathbf{X}$ and $\mathbf{Y}$ variables. By expanding the MZ projection operator, it is possible to derive the following expression for the projected dynamics on the $\mathbf{X}$ variables, which is valid up to order $O(\epsilon^3)$:

$$ \dot{\mathbf{X}} = f_{\mathbf{X}}(\mathbf{X}) + \epsilon M(\mathbf{X}) + \epsilon S(\mathbf{X}) + \epsilon^2 \int K(\mathbf{X}, t-s)ds. \quad (43) $$

This equation provides the explicit expression of the mean-field, deterministic term $M$; the time-correlation properties of the stochastic term $S$ that is, in general, multiplicative; and of the integration kernel $K$, which defines the non-Markovian contribution.

Figure 34 provides a schematic diagram of the three terms of the parametrization: the $M$-term comes from time averaging of the effects of the $\mathbf{Y}$ variables on the $\mathbf{X}$ variables; the $S$-term results from the fluctuations of the forcing of the $\mathbf{Y}$ variables on the $\mathbf{X}$ variables; and the non-Markovian contribution represents the self-interaction of the $\mathbf{X}$ variables on themselves at a later time, mediated by the $\mathbf{Y}$ variables.
These terms are derived using the statistical properties of the uncoupled dynamics of the **Y** variables, at $\epsilon = 0$. In the limit of infinite time scale separation between the **X** and **Y** variables, the non-Markovian term drops out and the stochastic term becomes a white-noise contribution, where one needs to use the Stratonovich definition of the stochastic integral (Pavliotis and Stuart, 2008). If, instead, $\Psi_Y = 0$ in the master-slave system (42), the non-Markovian term is identically zero, as expected.

<figure><img src="image_43.png" /><figcaption>FIG. 34: Schematic representation of the three terms contributing to the parametrization of the fast variables **X** in Eq. (43). (a) Mean-field term **M**; (b) stochastic term **S**; and (c) memory term with kernel **K**. The symbols ($\tau$, $\tau_1$, $\tau_2$) denote the delays involved in these effects. Reproduced with permission from Wouters and Lucarini (2013).</figcaption></figure>

A surprising property of the surrogate dynamics given by Eq. (43) is that the expectation value of any observable $\phi(\mathbf{X})$ is the same as for the full dynamics governed by Eq. (42), up to third order in $\epsilon$ (Wouters and Lucarini, 2013). This property is due to the fact that one can derive Eq. (43) by treating the weak coupling using Ruelle response theory (Ruelle, 1998, 1999, 2009); see Sect. IV.E.3 below. As a result, the average model error due to the use of parametrized dynamics is well under control; see discussion in Hu et al. (2019). Moreover, Vissio and Lucarini (2018a,b) showed by explicit examples that this top-down approach can, in some cases, help derive scale-adaptive parametrizations, while Demaeyer and Vannitsem (2017) demonstrated its effectiveness in an intermediate complexity climate model.

## IV. CLIMATE SENSITIVITY AND RESPONSE

A central goal of the climate sciences is to predict the impact of changes in the system's internal or external parameters — such as the greenhouse gas (GHG) concentration or the solar constant — on its statistical properties. A key concept in doing so is climate sensitivity, which

aims to measure the response of the climate system to external perturbations of Earth's radiative balance. As we shall see below, this measure is being used for projecting, for instance, mean temperature changes over the coming century as a response to increasing concentrations of atmospheric GHGs. While a good start, accurate and flexible predictions of climate changes require, though, more sophisticated concepts and methods.

### A. A Simple Framework for Climate Sensitivity

In order to illustrate the main ideas, let us consider the simple energy balance model (EBM) introduced in Eqs. (13) of Sect. III.A where the net radiation $R = R_i - R_o$ at the top of the atmosphere is related to the corresponding average temperature $T$ near the Earth's surface by $R = R(T)$. This simple, 0-D EBM includes both longwave and shortwave processes, so that $cdT/dt = R(T)$, as in Eq. 13a.

Following Peixoto and Oort (1992) and Zaliapin and Ghil (2010), we assume, furthermore, that there are $N$ climatic variables $\{\alpha_k = \alpha_k(T), k = 1, \dots, N\}$ that are, to a first approximation, directly affected by the temperature change only and which can, in turn, affect the radiative balance. Hence, one can write $R = R(T, \alpha_1(T), \dots, \alpha_N(T))$. Let us assume, furthermore, that, for a certain reference temperature $T = T_0$ one has $R(T_0) = 0$, which corresponds to steady-state conditions.

The simplest framework for climate sensitivity is to think of the difference in global annual mean surface air temperature $\Delta T$ between two statistical steady states, which have distinct $CO_2$ concentration levels. We then assume that changing the $CO_2$ concentration corresponds to applying an extra net radiative forcing $\Delta \tilde{R}$ to the system, and look for the corresponding change $\Delta T$ in the average temperature, so that $R(T_0 + \Delta T) + \Delta \tilde{R} = 0$.

For small $\Delta T$ and smooth $R = R(T)$, the Taylor expansion yields

$$ \begin{aligned} \Delta \tilde{R} &= -\left.\frac{dR}{dT}\right|_{T=T_0} \Delta T + \mathcal{O}\left((\Delta T)^2\right) \\ &= \left.\frac{\partial R}{\partial T}\right|_{T=T_0} \Delta T - \sum_{k=1}^{N} \left.\frac{\partial R}{\partial \alpha_k}\frac{\partial \alpha_k}{\partial T}\right|_{T=T_0} \Delta T \\ &\quad + \mathcal{O}\left((\Delta T)^2\right). \end{aligned} \quad (44) $$

Here, $\mathcal{O}(x)$ is a function such that $\mathcal{O}(x) \le Cx$ as soon as $0 < x < \epsilon$ for some positive constants $C$ and $\epsilon$. While the higher-order terms in $\Delta T$ are usually small, they can become important where the smooth dependence of $R$ on $T$ breaks down. Specifically, rapid climate change may ensue when the system crosses a tipping point, as will be explained further below.
Introducing the notations

$$ \frac{1}{\lambda_0(T_0)} = -\left.\frac{\partial R}{\partial T}\right|_{T=T_0}, \quad (45a) $$

$$ f_k(T_0) = -\lambda_0(T_0) \left. \frac{\partial R}{\partial \alpha_k} \frac{\partial \alpha_k}{\partial T} \right|_{T=T_0}, \quad (45b) $$

for the “reference sensitivity” $\lambda_0$ and the “feedback factors” $f_k$ at the reference state $T = T_0$, we obtain

$$ \Delta \tilde{R} = \frac{1 - \sum_{k=1}^{N} f_k(T_0)}{\lambda_0(T_0)} \Delta T + \mathcal{O}((\Delta T)^2), \quad (46) $$

which readily leads to

$$ \Delta T = \Lambda(T_0) \Delta \tilde{R} + \mathcal{O}((\Delta T)^2). \quad (47) $$

Here

$$ \Lambda(T_0) = -\left.\frac{dR}{dT}\right|_{T=T_0} = \frac{\lambda_0((T_0))}{1 - \sum_{k=1}^{N} f_k(T_0)} \quad (48) $$

is the linear gain factor, which can be defined as long as $\sum_{k=1}^{N} f_k((T_0)) \neq 1$; note that, if the sum of the feedback factors exceeds unity, i.e. $\sum_{k=1}^{N} f_k((T_0)) > 1$, the system is unstable.

Referring again to Eqs. (13), the feedback associated with the dependence of $\alpha$ on the temperature in Eq. (13b) is usually taken to be the ice-albedo feedback, cf. Sect. III.A, while the dependence of the emissivity $m$ on the temperature in Eq. (13c) is associated with the changes in the atmospheric opacity. For the latter one, the standard, reference sensitivity, associated with deviations from Planck’s law for black-body radiation, is $\lambda_0(T_0) = -\partial R/\partial T|_{T=T_0} = 4\sigma m(T_0)T_0^3$.

More specifically, feedbacks that can contribute to changes in reflectivity in Eq. (13b) include, aside from the incremental presence of snow and ice, also the climate-vegetation feedback (e.g., Rombouts and Ghil, 2015; Watson and Lovelock, 1983; Zeng and Neelin, 2000, and references therein). The feedbacks that can affect the sensitivity of emitted radiation in Eq. (13c) include atmospheric alteration in water vapor content and cloud cover, as well as in GHGs and aerosol concentration.

In climate studies, different measures of climate sensitivity are used. The so-called equilibrium climate sensitivity (ECS) denotes the globally and annually averaged surface air temperature increase that would result from sustained doubling of the concentration of carbon dioxide in Earth’s atmosphere vs. that of the reference state, after the climate system had reached a new steady-state equilibrium (Charney et al., 1979). The ECS was used extensively by the IPCC’s first three assessment reports, up to IPCC (2001).

Taking the linear approximation in Eqs. (44)–(47) above, one has:

$$ \begin{aligned} \text{ECS} &= \Lambda(T_0) \Delta \tilde{R}_{2 \times \text{CO}_2} \\ &= \frac{\lambda_0((T_0))}{1 - \sum_{k=1}^{N} f_k(T_0)} \Delta \tilde{R}_{2 \times \text{CO}_2}. \end{aligned} \quad (49) $$

<figure><img src="image_44.png" /><figcaption>FIG. 35: Test of the linear scaling of the long-term climate response with respect to the CO₂ concentration increase in 15 GCMs. On the abscissa: standard ECS, as in Eq. (49); on the ordinate: long-term response of the globally averaged surface air temperature to quadrupling of the CO₂ concentration. Reproduced with permission from Pfister and Stocker (2017).</figcaption></figure>

Note that the radiative forcing is, to a good approximation, proportional to the logarithm of the CO₂ concentration. Hence, in the linear-response regime, and for a given reference state $T_0$, the long-term globally averaged surface air temperature change resulting from a quadrupling of the CO₂ concentration is twice as large as the ECS; see Fig. 35.

The concept of climate sensitivity can be generalized to describe the linear dependence of the long-term average of any climatic observable with respect to the radiative forcing due to changes in CO₂ or in other GHGs, as well as to changes in solar radiation, aerosol concentration or any other sudden changes in the forcing (Ghil, 1976, 2015; von der Heydt et al., 2016; Lucarini et al., 2010b).

## B. Climate Sensitivity: Uncertainties and Ambiguities

The ECS is widely considered to be the most important indicator in understanding climate response to natural and anthropogenic forcings. It is usually estimated from instrumental data coming from the industrial age, from proxy paleoclimatic data, and from climate models of different levels of complexity. In climate models, the ECS results from a nontrivial combination of several model parameters that enter the feedback factors {$f_k$} in Eq. (45b) above, and it requires careful tuning. Despite many years of intense research, major uncertainties still exist in estimating it from past climatic data, as well as substantial discrepancies among different climate models
<figure><img src="image_45.png" /></figure>

FIG. 36: Estimates of the state dependence of the ECS using the Community Earth System Model (CESM). The four curves correspond to CESM simulations with four multiples of the solar constant $Q_0$: in the notation of Eqs. (13) here, $\mu = 0.75, 0.875, 1.0$ and $1.1$ for the blue, green, brown and red curves. (b) ECS as a function of CO$_2$ partial pressure; and (c) ECS as a function of global mean surface air temperature. The temperatures shown in (c) are an average between the base and doubled-CO$_2$ state. The shaded region in (b) and (c) indicates the IPCC estimated range for ECS. Reproduced with permission from Figs. 1(b,c) of Wolf et al. (2018); see also Gómez-Leal et al. (2018).

(IPCC, 2001, 2007, 2014a). In fact, Charney et al. (1979) estimated the ECS uncertainty as 1.5–4.5 K for CO$_2$ doubling and this range of uncertainties has increased rather than decreased over the four intervening decades.

The basic reason for these uncertainties lies in the high sensitivity of the ECS to the strength of the feedbacks $f_k$, as a result of the factor $1/\{1 - \sum_{k=1}^N f_k(T_0)\}$ in Eq. (48). Efforts to reduce the uncertainty in ECS values for the current climate include adopting ultra-high resolution GCMs (Satoh et al., 2018), in which one may better account for feedbacks that act on a larger range of scales, to applying relations that are rigorously valid for simple stochastic models to the climate data (Cox et al., 2018), where one hopes to take advantage of general, and possibly universal, relationships between climatic variables.

In particular, the largest uncertainty in defining the ECS for the current climate state is associated with the difficulty in estimating correctly the strength of the two main feedbacks associated with clouds (Bony et al., 2015;

<figure><img src="image_46.png" /></figure>

FIG. 37: State dependence of the ECS: estimates from proxy data and from climate models. The large number of data sets, models and acronyms is detailed in the von der Heydt et al. (2016) review paper. The $\Delta T$ on the x-axis refers to the difference of the given $T_0$ from the pre-industrial value $T_{00}$, i.e., $T_0 = T_{00} + \Delta T$; the $S_{[CO_2,LI]}$ on the y-axis refers to sensitivity with respect to CO$_2$ concentration corrected for land-ice albedo feedback. Mean values: for data – color-coded circles with shaded probability density functions; for models – squares with error bars. Reproduced with permission from von der Heydt et al. (2016).

Schneider et al., 2017). A warmer climate leads to increased presence of water vapor in the atmosphere, and, in turn, to more clouds. An increased cloud cover leads, on the one side, to an increase in the climate system's albedo (cooling effect) and, on the other side, to a more efficient trapping of longwave radiation near the surface (warming effect). The balance between the two feedbacks changes substantially according to the type of cloud, with the cooling effect dominant for low-lying clouds, while the warming effect is dominant for high-altitude ones. This is a striking example of the multiscale nature of the climate system: an extremely small-scale, short-lived dynamical process — cloud formation — has a substantial effect on the planet's global and long-term energy budget.

Despite the highly simplified description above, it should be clear that the ECS is a state-dependent indicator. This state dependence is further supported by the evidence in Figs. 36 and 37. In particular:

* Both the Planck response and the strength of the feedbacks that determine the gain factor $\Lambda$ in Eq. (48) depend on the reference state $T_0$. As an example, in warmer climates where sea ice is absent, the positive ice-albedo feedback is greatly reduced, thus contributing to a smaller climate sensitivity. On the other hand, in warmer climates the atmosphere is more opaque as a result of the presence of more water vapor, leading to a strong enhancement of the greenhouse effect.
* The radiative forcing is only approximately linear with the logarithm of the CO₂ concentration, so that ΔR̃₂×CO₂ depends on the concentration's reference value. In fact, this dependence is weak across a large range of CO₂ concentrations, but it is greatly strengthened by optical saturation effects in the CO₂ absorption bands.

* Near the moist greenhouse threshold, which corresponds to a tipping point of the Earth system, the ECS is greatly strengthened. Figure 36 shows that for a solar irradiance comparable or stronger than the present one, the peak in the value of the ECS is realised at a surface temperature of about 320 K, which corresponds to a lower CO₂ concentration in the case of weaker irradiance; see the discussion in Gómez-Leal et al. (2018). Note that reaching the moist greenhouse threshold for lower values of the solar irradiance requires exceedingly high CO₂ concentrations.

* Recently, Schneider et al. (2019), using very high-resolution simulations that represent explicitly convective processes, proposed a mechanism of instability of stratocumulus clouds occurring at high CO₂ concentrations that greatly enhances the ECS and eventually leads to an abrupt transition to a much warmer climatic regime.

While useful, the ECS concept faces practical difficulties because its definition assumes that, after the forcing is applied, the climate reaches a new steady state after all transients have died out. Since the climate is multiscale in both time and space, it is extremely non-trivial to define an effective cut-off time scale able to include all transient behavior. Thus, a time scale of 100 years is long compared to atmospheric processes, but short with respect to oceanic ones that involve the deep ocean. While a time scale of 5 000 years is long compared to oceanic processes, but short with respect to cryospheric ones that involve the dynamics of the Antarctic ice sheets. Therefore, one needs to associate each ECS estimate from observational or model data to a reference time scale; see Fig. 38 for an illustrative cartoon, and its discussion by Members et al. (2012) and von der Heydt et al. (2016).

### C. Transient Climate Response (TCR)

Transient climate response (TCR) has recently gained popularity in the study of climate change and climate variability because of its ability to help capture the evolution in time of climate change by addressing the transient, rather than asymptotic, response of the climate system to perturbations in the CO₂ concentration. TCR is defined as the change in the globally averaged surface air temperature recorded at the time at which CO₂

<figure><img src="image_47.png" /><figcaption>FIG. 38: Dependence of the effective ECS on the reference time scale. Consideration of longer time scales entails taking into account a larger set of slow climate processes. Reproduced with permission from http://www.realclimate.org/index.php/archives/2013/01/on-sensitivity-part-i/.</figcaption></figure>

has doubled due to an increase at a 1% annual rate, i.e. roughly after 70 years, having started at a given reference value $T_0$ (Otto et al., 2013). This operational definition agrees rather well with the standard IPCC-like simulation protocols; and TCR is therewith better suited than ECS to test model outputs against observational data sets from the industrial era.

As shown in Fig. 39, the TCR is found to be lower than the ECS for a long time, because of the climate system's thermal inertia, which is dominated by the oceans' heat capacity. A smaller effective heat capacity $c$ in Eq. (13a), and hence a shorter relaxation time, would result in the TCR catching up much faster with the ECS, as is the case in regular diffusion processes; see Fig. 40a. Assuming linearity in the response, a relationship must clearly exist between ECS and TCR. So far, this inference has been based, by-and-large, on heuristic arguments of time scale separation between climate feedbacks rather than being rigorously derived (Otto et al., 2013); it will be derived more systematically in Sect. IV.E herein.

### D. Beyond Climate Sensitivity

The standard viewpoint on climate sensitivity discussed above is associated with the idea that the climate is in equilibrium, in the absence of external perturbations. In the setting of deterministic, autonomous dynamical systems, this view can be described by the change in the position of a fixed point, $X_0 = X_0(\mu)$, as a function of a parameter $\mu$.

We illustrate in Fig. 40 the difference between the ways that a change in a parameter can affect a climate model's behavior in the case of equilibrium solutions, cf. panel (a), vs. more complex dynamical behavior, cf. panels (b) and (c). Assume that the climate state is periodic,
<figure><img src="image_48.png" /><figcaption>FIG. 39: Transient climate response (TCR) estimates from two climate models: (i) a coupled atmosphere-ocean GCM (red curve) and (ii) a simple illustrative model with no energy exchange with the deep ocean (green curve). Time on abscissa from start of CO₂ concentration increase at pre-industrial levels, with change in global mean temperature on the ordinate. The “additional warming commitment” corresponds to temperature stabilization at a given CO₂ level, i.e., at 2 × CO₂ or at 4 × CO₂. Reproduced with permission from IPCC (2001, Fig. 9.1).</figcaption></figure>

i.e., lies on a limit cycle, rather than being a fixed point, as in panel (a). In this case, climate sensitivity can no longer be defined by a single scalar, $\partial \bar{T}/\partial \mu$, but needs four scalars — the sensitivity of the mean temperature along with that of the limit cycle's frequency, amplitude, and phase — or more, e.g., the orbit's ellipticity, too.

But the internal climate variability can be better described in terms of strange attractors than by fixed points or limit cycles. Moreover, the presence of time-dependent forcing, deterministic as well as stochastic, introduces additional complexities into the proper definition of climate sensitivity. It is thus apparent that a rigorous definition of climate sensitivity requires considerably more effort.

Ghil (2015) proposed to measure the change in the overall properties of the attractor before and after the change in forcing by computing the Wasserstein distance $d_W$ between the two invariant measures. The Wasserstein distance or “earth mover’s distance” $d_W$ quantifies the minimum “effort” in morphing one measure into another one of equal mass on a metric space, like an $n$-dimensional Euclidean space (Dobrushin, 1970).

Monge (1781) originally introduced this distance to study a problem of military relevance. Roughly speaking, $d_W$ represents the total work needed to move the “dirt” (i.e., the measure) from a trench you are digging to another one you are filling, over the distance between the two trenches. In general, the shape of the two trenches and the depth along the trench — i.e., the support of the measure and its density — can differ. Robin et al. (2017)

and Vissio and Lucarini (2018b) recently showed the effectiveness of applying this idea to climate problems.

## E. A General Framework for Climate Response

A major use of state-of-the-art climate models is to produce projections of climate change taking into account different possible future scenarios of emission of greenhouse gases and pollutants like aerosols, as well as changes in land use, which has substantial impacts on the terrestrial carbon cycle. Projections are needed not only for quantities like the global average surface air temperature, but on spatially and temporally detailed information is needed for a multitude of practical needs; see, for instance, Fig. 18.

The ECS concept is well suited for describing the properties of equilibrium solutions of heuristically simplified equations of the climate system, like Eq. 13a, and has clear intuitive appeal, as in Fig. 40a. But it also has basic scientific limitations:

* it only addresses long term-climatic changes and no detailed temporal information, an issue only partially addressed by TCR information;

* it only addresses changes in the globally averaged surface air temperature and no spatial information at the regional scale and at different levels of the atmosphere, of the ocean, and of the soil;

* it cannot discriminate between radiative forcings resulting from different physical and chemical processes, e.g. differences resulting from changes in aerosol vs. GHG concentration; the two impact quite differently shortwave and longwave radiation, and different atmospheric levels.

We will thus try to address these shortcomings by taking the complementary points of view of nonequilibrium statistical mechanics and dynamical systems theory. The setting of nonautonomous and of stochastically forced dynamical systems allows one to examine the interaction of internal climate variability with the forcing, whether natural or anthropogenic; it also helps provide a general definition of climate response that takes into account the climate system's nonequilibrium behavior, its time-dependent forcing, and its spatial patterns.

### 1. Pullback Attractors (PBAs)

The climate system experiences forcings that vary on many different time scales (e.g., Saltzman, 2001), and its feedbacks also act on multiple time scales (e.g., Ghil and Childress, 1987). Hence, defining rigorously what climate response to forcing, vs. intrinsic variability, actually is requires some care: Observed variations can be related to
<figure><img src="image_49.png" /><figcaption>FIG. 40: Climate sensitivity (a) for an equilibrium model; (b) for a nonequilibrium, oscillatory model; and (c) for a nonequilibrium model featuring chaotic dynamics and stochastic perturbations. As a forcing (atmospheric CO₂ concentration, say, blue dash-dotted line) changes suddenly, global temperature (red heavy solid) undergoes a transition. (a) Only the mean temperature T̄ changes; (b) the amplitude, frequency, and phase of the oscillation change, too; and (c) all the details of the invariant measure, as well as the correlations at all orders, are affected. Reproduced from Ghil (2017), with permission from the American Institute of Mathematical Sciences.</figcaption></figure>

the presence of natural periodicities — such as the daily and the seasonal cycle, and orbital forcings; to rapid, impulsive forcings, such as volcanic eruptions; or to slow modulations to the parameters of the system, as in the case of anthropogenic climate change.

For starters, consider a dynamical system in continuous time,

$$ \dot{x} = F(x, t) \qquad (50) $$

on a compact manifold $\mathcal{Y} \subset \mathbb{R}^d$; here $x(t) = \phi(t, t_0)x(t_0)$, with initial state $x(t_0) = x_0 \in \mathcal{Y}$. The evolution operator $\phi(t, t_0)$ is assumed to be defined for all $t \ge t_0$, with $\phi(s, s) = \mathbf{1}$, and it thus generates a two-parameter semigroup. In the autonomous case, time-translational invariance reduces the latter to a one-parameter semigroup since, $\forall t \ge s$, $\phi(t, s) = \phi(t - s)$. In the nonautonomous case, in other terms, there is an absolute clock.

We are interested in forced and dissipative systems such that, with probability one, initial states in the remote past are attracted at time t towards $\mathbf{A}(t)$, a time-dependent family of geometrical sets that define the system's pullback attractor (PBA). In the autonomous case, $\mathbf{A}(t) \equiv \mathbf{A}_0$ is the time-independent attractor of the system, and it is known to support, under suitable conditions, a physical measure $\mu(dx)$ (Eckmann and Ruelle, 1985; Ledrappier and Young, 1988).<fnref n="7" />

Such a PBA can also be constructed when random forcing is present (e.g., Arnold, 1988, and references therein),

$$ dx = F(x, t)dt + g(x)d\eta, \qquad (51) $$

where $\eta = \eta(t; \omega)$ is a Wiener process, while $\omega$ labels the particular realization of this random process, and $d\eta(t)$ is commonly referred to as “white noise.” The noise can be multiplicative, and one then uses the Itô calculus for the integration of Eq. (51).

In the random case, the PBA $\mathbf{A}(t; \omega)$ is commonly referred to as a random attractor. A more detailed and mathematically rigorous discussion of these concepts appears in Carvalho et al. (2013); Chekroun et al. (2011) and Ghil et al. (2008, Appendix A). Careful numerical applications of PBAs to explain the wind-driven circulation and the THC are now available (Pierini et al., 2018, 2016; Sevellec and Fedorov, 2015).

In the purely deterministic case, the theory of nonautonomous dynamical systems goes back to the skew-product flows of Sell (1967). A concept that is closely related to PBAs is that of snapshot attractors; it was introduced in a more intuitive and less rigorous manner into the physical literature by Romeiras et al. (1990). Snapshot attractors have also been used for studying time-dependent problems of climatic relevance (Bódai et al., 2013, 2011; Bódai and Tél, 2012; Drótos et al., 2015).

Note that, in the most recent IPCC reports (IPCC, 2001, 2007, 2014a) and according to the standard protocols described in Sect. II.D.2, future climate projections

<footnote n="7">⁷ Among invariant measures, a natural measure is one obtained by flowing a volume forward in time, and a physical measure is one for which the time average equals the ensemble average almost surely with respect to Lebesgue measure; moreover physical implies natural. A particular class of invariant measures of interest are Sinai-Ruelle-Bowen (SRB) measures (e.g., Young, 2002),</footnote>

<footnote>and an important result is that — for a system with no null Lyapunov exponents, except the one corresponding to the flow — an ergodic SRB measure is physical.</footnote>
<figure><img src="image_50.png" /><figcaption>FIG. 41: A cartoon depicting the pullback attractor (PBA) for a nonautonomous system as in Eq. (50). The set **B** of finite Lebesgue measure, initialized at time $t = t_0 - \tau$, evolves towards the set **A**($t_0$). The construction is exact in the limit $\tau \to \infty$; for a random attractor, see Ghil *et al.* (2008, Figs. A.1 and A.2). Adapted from Sevellec and Fedorov (2015), with permission.</figcaption></figure>

are virtually always performed using as initial states the final states of sufficiently long simulations of historical climate conditions. As a result, it is reasonable to assume that the pullback time $\tau$, as defined in Fig. 41, is large enough, and that the covariance properties of the associated $\mathbf{A}(t)$ sets are therefore well approximated.

## 2. Fluctuation Dissipation and Climate Change

The fluctuation-dissipation theorem (FDT) has its roots in the classical theory of many-particle systems in thermodynamic equilibrium. The idea is very simple: the system's return to equilibrium will be the same whether the perturbation that modified its state is due to a small external force or to an internal, random fluctuation. The FDT thus relates natural and forced fluctuations of a system (e.g., Kubo, 1957, 1966); it is a cornerstone of statistical mechanics and has applications in many areas (Marconi *et al.*, 2008, and references therein).

We have emphasized already in Sect. III.B that, even when the climate system is in a steady state, it is not at all in thermodynamic equilibrium (e.g., Ghil, 2019; Lucarini and Ragone, 2011). Still, Leith (1975) showed that FDT applies to a 2-D or QG turbulent flow with two integral invariants, kinetic energy $E$ and enstrophy $Z$, under some additional assumptions of normal distribution of the realizations and stationarity. Soon thereafter, Bell (1980) showed that the FDT still seemed to work for a highly truncated version of such a model, even in the presence of dissipative terms that invalidate the thermodynamic equilibrium assumption.

The FDT has been applied to the output of climate models to predict the climate response to a step-like in-

crease of the solar irradiance (North *et al.*, 1993) as well as to increases in atmospheric CO$_2$ concentration (Cionni *et al.*, 2004; Langen and Alexeev, 2005), while Gritsun and Branstator (2007) and Gritsun *et al.* (2008) used it to predict the response of an atmospheric model to localised heating anomalies. Most recently, Cox *et al.* (2018) tried to reduce the uncertainty in the ECS discussed in Sect. IV.B by a systematic FDT application to an ensemble of model outputs, as well as to the observed instrumental climate variability.

The FDT-based response of the system to perturbations in the above examples reproduces the actual changes at a good qualitative, rather than strictly quantitative level. An important limitation of these insightful studies is the use of a severely simplified version of the FDT that is heuristically constructed by taking a gaussian approximation for the invariant measure of the unperturbed system. This approximation amounts to treating the climate as being in thermodynamic equilibrium; see also Majda *et al.* (2009), who specifically address the applicability of the FDT in a reduced phase space.

In order to address the problem of gaussian approximation, Cooper and Haynes (2011) proposed to construct a kernel-based approximation of the actual invariant measure of the unperturbed system and then use Eq. (56) below to construct the system's Green function. This approach has been applied successfully in a very low-dimensional system, but its robustness may be limited by the kernel's arbitrary cut-off not being able to account for the smaller scales of the invariant measure's fine structure.

## 3. Ruelle Response Theory

FDT generalizations to systems out of equilibrium have been developed since the early 1950s (e.g., Kubo, 1966). But a particularly fruitful change in point of view was provided by D. Ruelle (Ruelle, 1998, 1999, 2009), who considered the problem in the setting of dynamical systems theory, rather than that of statistical mechanics. The former point of view is justified in this context by the so-called chaotic hypothesis (Gallavotti and Cohen, 1995), which states, roughly, that chaotic systems with many degrees of freedom possess a physically relevant invariant measure, as discussed in Sect. IV.E.1.

It is common in practice to assume that a time-dependent measure $\mu_t(dx)$ associated with the evolution of the dynamical system given by Eq. (50) does exist. Still, computing the expectation value of measurable observables with respect to this measure is in general far from trivial and requires setting up a large ensemble of initial states in the Lebesgue measurable set **B** mentioned before. Moreover, PBAs and the physical measures they support only set the stage for predicting the system's sensitivity to small changes in the forcing or the dynamics.
Ruelle's response theory allows one to compute the change in the measure $\mu(dx)$ of an autonomous Axiom A system due to weak perturbations of intensity $\epsilon$ applied to the dynamics, in terms of the unperturbed system's properties. The basic idea behind it is that the invariant measure of such a system, even though supported on a strange attractor, is differentiable with respect to $\epsilon$. Baladi (2008) reviewed extensions of response theory and provided a complementary approach to it, while Lucarini (2016) estimated its radius of convergence in $\epsilon$.

The nonautonomous version of the Ruelle response theory allows one to calculate the time-dependent measure $\mu_t(dx)$ on the PBA by computing the time-dependent corrections to it with respect to a reference state $x(t) = \tilde{x}(t)$. This provides a fairly general formulation for the climate system's response to perturbations.

Let us assume that we can write

$$ \dot{x} = F(x,t) = F(x) + \epsilon X(x,t) \quad (52) $$

where, $\forall t \in \mathbb{R}$ and $\forall x \in \mathcal{Y} \subset \mathbb{R}^d$, $|\epsilon X(x,t)| \ll |F(x)|$. Hence, we can take $F(x)$ as the background dynamics and $\epsilon X(x,t)$ as a perturbation. We only treat here the case of deterministic dynamics; see Lucarini (2012) for stochastic perturbations. As shown by Lucarini et al. (2017), we can restrict our analysis without loss of generality to the separable case of $F(x,t) = F(x) + \epsilon X(x)T(t)$.

To evaluate the expectation value $\langle\Psi\rangle^\epsilon(t)$ of a measurable observable $\Psi(x)$ with respect to the measure $\mu_t(dx)$ of the system governed by Eq. (50), one writes:

$$ \langle\Psi\rangle^\epsilon(t) = \int \Psi(x) \mu_t(dx) = \langle\Psi\rangle_0 + \sum_{j=1}^{\infty} \epsilon^j \langle\Psi\rangle_0^{(j)}(t); \quad (53) $$

here $\langle\Psi\rangle_0 = \int \Psi(x) \bar{\mu}(dx)$ is the expectation value of $\Psi$ with respect to the SRB invariant measure $\bar{\mu}(dx)$ of the autonomous dynamical system $\dot{x} = F(x)$. We will restrict ourselves here to the linear correction term, which can be written as:

$$ \begin{aligned} \langle\Psi\rangle_0^{(1)}(t) &= \int_0^\infty \Lambda S_0^\tau \Psi(x) T(t-\tau) d\tau \bar{\mu}(dx); \\ &= \int_0^\infty G_{\Psi,X}^{(1)}(\tau) T(t-\tau) d\tau. \end{aligned} \quad (54) $$

The Green's function $G_{\Psi,X}^{(1)}(\tau)$ above is given by

$$ G_{\Psi,X}^{(1)}(\tau) = \int \Theta(\tau) \Lambda S_0^\tau \Psi(x) \bar{\mu}(dx), \quad (55) $$

where $\Lambda(\bullet) = X \cdot \nabla(\bullet)$ and $S_0^t(\bullet) = \exp(tF \cdot \nabla)(\bullet)$ is the semigroup of unperturbed Koopman operators$^8$, $(\cdot)$

denotes the inner product in $\mathcal{Y}$, and the Heaviside distribution $\Theta(\tau)$ enforces causality. If the unperturbed invariant measure $d\bar{\mu}(x)$ is smooth with respect to the standard Lebesgue measure, one has $\bar{\mu}(dx) = \tilde{\mu}(x)dx$, with $\tilde{\mu}(x)$ the density, and the Green's function can be written as follows:

$$ \begin{aligned} G_{\Psi,X}^{(1)}(\tau) &= \Theta(\tau) \int \frac{-\nabla \cdot (\tilde{\mu}(x)X)}{\tilde{\mu}(x)} S_0^\tau \Psi(x) \tilde{\mu}(x) dx \\ &= \Theta(\tau) C(\Phi, S_0^\tau \Psi); \end{aligned} \quad (56) $$

here $\Phi = -(\nabla \cdot (\tilde{\mu}(x)X))/\tilde{\mu}(x)$ and $C(A, S_0^\tau B)$ is the $\tau$-lagged correlation between the variables $A$ and $B$, and the average of $\Phi$ vanishes. Note that Eq. (56) is the appropriate generalization of the FDT for the nonautonomous, out-of-equilibrium system (52).

Given any specific choice of the forcing's time dependency $T(t)$ in Eqs. (52)-(54) and measuring the linear correction term $\langle\Psi\rangle_0^{(1)}(t)$ from a set of experiments, the same equations allow one to derive the appropriate Green's function. Therefore, using the output of a specific set of experiments or of GCM simulations, we achieve predictive power for any temporal pattern of the forcing $X(x)$.

Consider now the Fourier transform of Eq. (54):

$$ \langle\Psi\rangle_0^{(1)}(\omega) = \chi_{\Psi,X}^{(1)}(\omega) T(\omega), \quad (57) $$

where we have introduced the susceptibility $\chi_{\Psi,X}^{(1)}(\omega) = \mathcal{F}[G_{\Psi,X}^{(1)}]$, defined as the Fourier transform of the Green's function $G_{\Psi,X}^{(1)}(t)$. Under suitable integrability conditions, the fact that the Green's function $G_{\Psi,X}^{(t)}$ is causal is equivalent to saying that its susceptibility obeys the so-called Kramers-Kronig relations (Lucarini and Sarno, 2011; Ruelle, 2009); these provide integral constraints that link the real and imaginary parts, so that $\chi_{\Psi,X}^{(1)}(\omega) = i\mathcal{P}(1/\omega) \star \chi_{\Psi,X}^{(1)}(\omega)$, where $i = \sqrt{-1}$, $\star$ indicates the convolution product, and $\mathcal{P}$ stands for integration by parts. Extensions to the case of higher-order susceptibilities are also available (Lucarini, 2008a, 2009a; Lucarini and Colangeli, 2012; Lucarini et al., 2005b).

Instead of studying merely individual trajectories, one can study the evolution of ensembles of trajectories to obtain probabilistic estimates. The evolution of the measure $\rho$ driven by the dynamical system $\dot{x} = F(x)$ is described by the transfer or Perron-Frobenius operator $\mathcal{L}_0^t$ (Baladi, 2000; Chekroun et al., 2014; Villani, 2009), which is the adjoint of the Koopman or composition operator $S_0^t$; it is defined as follows:

$$ \int \mathcal{L}_0^t \rho(x) \Psi(x) dx = \int \rho(x) S_0^t \Psi(x) dx. \quad (58) $$

Assuming that no degenaries are present, the generator of the semigroup $\{\mathcal{L}_0^t\}$ is the so-called Liouville operator $L$ and it satisfies $\partial_t \rho = -\nabla(\rho F) = L\rho$.

<footnote>$^8$ In discrete-time dynamics given by $x_{k+1} = g(x_k)$, the Koopman operator is a linear operator acting on observables $h: \mathcal{Y} \to \mathbb{R}$ via $Uh(x) = h(f(x))$. In continuous-time dynamics, like Eq. (50) herein, this operator is replaced by a semigroup of operators, as in Eq. (55) above; see Budisić et al. (2012) for a good review.</footnote>
Let $\{(\sigma_k, \rho_k) : k = 1, \dots, \infty\}$ be the eigenpairs of $L$. Then the eigenvalues of $\mathcal{L}_0^t$ are given by $\{\exp(\sigma_k t)\}$, with the same eigenvectors $\{\rho_k\}$. Correspondingly, the invariant measure $\bar{\mu}$ is the eigenvector having the null eigenvalue $\sigma_1 = 0$ for $L$ or, $\forall t \ge 0$, unit eigenvalue of $\mathcal{L}_0^t$. Note that $\forall k \ge 2$, $\Re(\sigma_k) < 0$, which implies exponential decay of correlations for all the system's smooth observables, with an asymptotic rate that is given by the largest value of $\Re(\sigma_k)$, $k \ge 2$. Moreover, the presence of a small spectral gap between the unit eigenvalue of $\mathcal{L}_0^t$ and its other eigenvalues within the unit disk leads to a small radius of expansion for Ruelle's perturbative approach (Chekroun et al., 2014; Lucarini, 2016).

It turns out that, if one neglects the essential part of the spectrum, the susceptibility can be written as

$$ \chi_{\Psi,X}^{(1)}(\omega) = \sum_{k=1}^{\infty} \frac{\alpha_k \{\Psi, X\}}{\omega - \sigma_k}, \qquad (59) $$

where the factor $\alpha_k$ evaluates how the response of the system to the forcing $X$ for the observable $\Psi$ projects on the eigenvector $\{\rho_k\}$. The constants $\pi_k = i\sigma_k$ are usually referred to as Ruelle-Pollicott poles (Pollicott, 1985; Ruelle, 1986). Note that Eq. (59) implies that the Green's function corresponding to the susceptibility can be written as a weighted sum of exponential functions, if one neglects possible degeneracies.<fnref n="9" />

## 4. Climate Change Prediction via Ruelle Response Theory

A somewhat different approach to constructing the climate response to forcings focuses on computing it directly from Eq. (55), without relying on the applicability of the FDT, which fails in certain cases of geophysical relevance (e.g., Gritsun and Lucarini, 2017). The difficulty in applying this direct approach lies in the fact that the formula contains contributions from both stable and unstable directions in the tangent space (Ruelle, 2009).

Evaluating the contribution of the unstable directions is especially hard; hence, Abramov and Majda (2007) proposed a blended approach that also uses the FDT. Using adjoint methods in the direct approach has also yielded promising results (Wang, 2013). Faced with the so-called cold-start problem of climate simulations, Haselmann et al. (1993) suggested a heuristic approach to

computing a climate model's Green's function and applied it to study the relaxation to steady state of a coupled GCM's globally averaged surface temperatures.

Lucarini and associates (Lucarini et al., 2017; Lucarini and Sarno, 2011; Ragone et al., 2016) proposed to evaluate the Green's function using an experimental but rigorous approach, suggested by standard optics laboratory practice (Lucarini et al., 2005b). The idea is to use a set of carefully selected probe experiments — typically, step-like increases of the parameter of interest — to construct the Green's function and then, exploiting Eq. (55), use this operator to predict the response of the system to a temporal pattern of interest for the forcing.

Given a set of forced climate simulations and a background unperturbed one, such an approach allows one to construct the Green's functions response operators for as many observables as desired, global as well as local in space. Such a toolkit allows one to treat a continuum of scenarios of temporal patterns forcings, thus providing a general framework for improving climate change projections given in the form shown in Fig. 18.

The idea is to consider the set of equations describing an unperturbed climate evolution in the form $\dot{x} = F(x)$, with the vector field $X = X(x)$ in Eq. (52) as the 3-D radiative forcing associated with the increase of CO$_2$ concentration, and $\epsilon T(t)$ its time modulation. By plugging $T(t) = \epsilon \Theta(t)$ into Eq. (54), we have, for any climatic observable $\Psi$,

$$ \frac{d}{dt} \langle \Psi \rangle_0^{(1)}(t) = \epsilon G_{\Psi, [CO_2]}^{(1)}(t). \qquad (60) $$

We estimate $\langle \Psi \rangle_0^{(1)}(t)$ by taking the system's average of response over an ensemble of initial states and use the previous equation to derive our estimate of $G_{\Psi, [CO_2]}^{(1)}(t)$, by assuming linearity in the response.

Note that these Green's functions are specifically related to changes in the atmospheric CO$_2$ concentration, rather than to a generic radiative forcing. Indeed, for each climatic variable, the response to changes in the solar irradiance — e.g., via modulation of the parameter $\mu$ in Eqs. (13) — is different from the impact of changes in the CO$_2$ concentration, because the details of the radiative forcing are very different in the two cases. Thus, geoengineering proposals that aim to reduce solar irradiance by injecting sulphate aerosols in the stratosphere (e.g., Smith and Wagner, 2018) rely on flawed scientific reasoning (Bódai et al., 2020; Lucarini, 2013), as well as being of dubious practical help (Proctor et al., 2018) and hardly defensible from an ethical perspective (Lawrence et al., 2018).

Ragone et al. (2016) derived, moreover, a closed formula relating the TCR and the ECS. Indeed, with the definition of $TCR(\tau)$ given at the beginning of

<footnote n="9">⁹ Note that Eq. (59) mirrors the quantum-mechanical expressions for the electric susceptibility of atoms or molecules. In the latter case, the summation involves all the pairs of eigenstates of the system's unperturbed Hamiltonian operator; in each term, the poles' imaginary part corresponds to the energy difference between the pair of considered eigenstates; and the real part is the so-called line width of the transition, whose inverse is the life time; finally, the numerator is its so-called dipole strength (Cohen-Tannoudji et al., 2007).</footnote>
Sect. IV.C, one gets

$$ 
\begin{aligned}
ECS - TCR(\tau) &= INR(\tau) \\
&= \Delta \tilde{R}_{2 \times CO_2} \times P \int_{-\infty}^{\infty} \chi_{T_S, [CO_2]}^{(1)}(\omega') \\
&\qquad \frac{1 + \frac{\sin(\omega' \tau/2)}{\omega' \tau/2} \exp[-i\omega' \tau/2]}{2\pi i\omega} d\omega'.
\end{aligned} 
\quad (61) $$

The difference between ECS and TCR is a weighted integral of the susceptibility, accounting for the contribution of processes and feedbacks occurring at different time scales. The integral in Eq. (61) yields $\chi_{T_S, [CO_2]}^{(1)}(0)$ if $\tau \to 0$ (i.e., $TCR(0) = 0$), decreases with $\tau$, and vanishes in the limit $\tau \to \infty$. $INR(\tau)$ is a measure of the system's inertia at the time scale $\tau$, due to the overall contribution of the internal physical processes and characteristic time scales of the relevant climatic subsystems.

Figure 42 illustrates a climate prediction obtained by applying the response theory above to the open-source model PLASIM (see http://tiny.cc/zgk0bz), an atmospheric GCM of intermediate complexity coupled with a mixed-layer ocean model (Fraedrich, 2012; Fraedrich et al., 2005). The figure shows a good agreement between (a) the ensemble average of 200 simulations where the CO₂ concentration is increased at the rate of 1% until doubling (black curve), and then kept constant; and (b) the prediction done by convolving the Green's function of Eq. (60) with the time pattern of the forcing (blue curve). Note that the temperature increase predicted at year 70, when doubling of the CO₂ is reached, gives as estimate the $TCR \simeq 4.1$ K (blue curve). The true value of the $TCR \simeq 4.3$ K is given by the black line. The TCR is indeed smaller than the $ECS \simeq 4.8$ K, which corresponds to a good approximation to the temperature increase predicted at year 200.

The power of response theory lies in the fact that — once the Green's function of interest has been computed — one can predict the future evolution of bespoke observables, defined as needed. Furthermore, Lucarini (2018) has shown that it is possible to use certain classes of observables as surrogate forcings of other observables, in the spirit of the linear feedback analysis of Eq. (44). Given the notoriously more difficult prediction of precipitation, Lucarini et al. (2017) have found that predictions of changes in the globally averaged precipitation can be as good as those of the temperature.

Proceeding from global predictions to more localized ones, Fig. 43 shows the outcome of predicting the change in the zonal averages of the surface temperature. It is clear that response theory does a good job in reproducing the spatial patterns of temperature change, except for an underestimate of the temperature change in the high latitudes on the scale of few decades. Indeed, such mishap is due the strong polar amplification of the warming due to the ice-albedo feedback, which can only qual-

<figure><img src="image_51.png" /><figcaption>FIG. 42: Comparison between the climate model simulation (black) and response theory prediction (blue) for a TCR experiment using PLASIM, a GCM of intermediate complexity. The CO₂ concentration was ramped up by 1% per year to double its initial value. The upper and lower limit of the light-shaded bands are computed as two standard deviations of the ensemble distribution. Reproduced with permission from (Lucarini et al., 2017).</figcaption></figure>

itatively captured by a linear approach like the one used by Lucarini et al. (2017).

## 5. Slow Correlation Decay and Sensitive Parameter Dependence

Equation (59) implies that resonant amplification occurs if the forcing acts at a frequency $\omega$ that is close to the imaginary part of a subdominant Ruelle-Pollicott pole $\pi_k$, $k \ge 2$ that has a small real part, i.e., for some $k \ge 2$,

$$ |\omega - \Im(\pi_k)| \ll \omega, \quad |\Re(\pi_k)| \ll 1, \quad (62) $$

because the system's susceptibility is greatly enhanced at such an $\omega$. Conditions (62) are easily satisfied for a broadband forcing and a system that has a small spectral gap. Conversely, as discussed after Eq. (58), the presence of a small spectral gap is associated with a slow decay of correlation for smooth observables. We now provide two examples that show how linear response breaks down for forced systems possessing slow decay of correlations due to a small spectral gap in the unperturbed dynamics.

The first example is due to Chekroun et al. (2014), who investigated the response of the highly simplified ENSO model of Jin and Neelin (1993), to which Jin et al. (1994) added periodic forcing that led to chaotic behavior. Chekroun and coauthors focussed on how the response changes when one varies the model parameter $\delta$ that controls the travel time of the equatorially trapped waves, which play an essential role in the ENSO mechanism. Chekroun et al. (2014) found that when the spectral gap is small, the system exhibits rough dependence of its properties with respect to small modulations of the
<figure><img src="image_52.png" /><figcaption>FIG. 43: Patterns of climate response for the zonal average of the surface temperature. (a) Prediction via response theory (red and black lines not relevant in this context). b) Difference between the ensemble average of the direct numerical simulations and the predictions obtained using the response theory. Reproduced with permission from (Lucarini et al., 2017).</figcaption></figure>

parameters, since Ruelle's perturbative expansion breaks down even for a very small intensity of the forcing.

The second example shows the importance of correlation slowdown near saddle-node bifurcation points in the PLASIM model discussed already in the previous subsection. This intermediate-complexity model is multistable: its bifurcation diagram with respect to the solar insolation parameter $\mu$ is very similar to Fig. 21 of Sect. III.A see Lucarini et al. (2010a, 2013). PLASIM's bifurcation diagram is reproduced in Fig. 45a here, in which the present climate $S = S_0 = 1360 \text{ Wm}^{-2}$ is denoted by a bullet ● and the bifurcation points associated with the transition from the warm to the snowball state ($W-SB$) and from the snowball to the warm state ($SB-W$) are indicated by solid arrowheads ▲. For comparison with Fig. 21, note that $\mu = S/S_0$ there.

Tantet et al. (2018) investigated PLASIM's response to changes in the solar irradiance for $S < S_0$. They found that as the $W-SB$ transition nears, the lag correlation for a large-scale observable — namely the average

<figure><img src="image_53.png" /><figcaption>FIG. 44: Rough dependence (red dots) of the skewness of the equatorial ocean temperature's distribution with respect to changes in the parameter $\delta$ that controls the travel time of equatorially trapped waves in a simplified ENSO model; see text for details. Black points refer to nonchaotic windows where the parameter dependence is smooth; blue points indicate particularly sharp jumps in the parameter dependence. Reproduced with permission from Chekroun et al. (2014).</figcaption></figure>

equatorial near surface air temperature — decays more and more slowly, as a result of a narrowing spectral gap. Indeed, near the saddle-node bifurcations, the response of the system to perturbations is greatly amplified, and, by definition, becomes singular exactly at the bifurcation point.<fnref n="10" />

This second example leads us naturally to the investigation of the climate system's global stability properties and of its critical transitions, which are the subject of the next section.

## V. CRITICAL TRANSITIONS AND EDGE STATES

In Sect. III, we introduced several types of bifurcations that involve transitions between two or more regimes: saddle-node bifurcations whose pairing can lead to coexistence of two stable steady states (fixed points); pitchfork bifurcations that, in the presence of a mirror symmetry, can lead from one stable symmetric steady state

<footnote n="10">¹⁰ Scheffer et al. (2009) have proposed a set of “warning signals” when a dynamical system approaches a critical transition. Correlation slowdown is one of the main ones and it is verified in the present case for the back-to-back saddle-node bifurcations that make up the backbone of the PLASIM model’s hysteresis cycle. Colon et al. (2015), though, have found that — in a predator–prey system, modeled by an ODE system as well as by an agent-based mode — more complex critical transitions can behave rather differently than in the by-and-large saddle-node-bound early warning literature.</footnote>
<figure><img src="image_54.png" /><figcaption>FIG. 45: Slowdown of correlations near a saddle-node bifurcation point in the PLASIM model. (a) The model's bifurcation diagram as a function of the value of the solar irradiance $S$. The critical transitions Warm-to-Snowball and Snowball-to-Warm as well as the present-day conditions are indicated. Reproduced with permission from Lucarini et al. (2010a). (b) Decay of correlation of the equatorial mean surface temperature for various values of the solar irradiance $S$. The decay becomes considerably slower as $S$ approaches the critical value of $S \simeq 1265$ Wm$^{-2}$. Reproduced with permission from Tantet et al. (2018).</figcaption></figure>

to the coexistence of two steady states that are mirror images of each other; Hopf bifurcations that lead from a stable steady state to a stable periodic solution (limit cycle) and then again from the stable limit cycle to a stable torus on which an infinity of quasi-periodic solutions live; and nonlocal bifurcations associated with the existence of homoclinic and heteroclinic orbits that lead on to chaotic regimes.

Successive bifurcation scenarios (Dijkstra, 2013; Dijkstra and Ghil, 2005; Eckmann, 1981; Ghil and Chil-dress, 1987) that involve several of these bifurcations and additional ones lead from solutions with high symmetry in space and time to successively more complex and

chaotic ones. We saw in Sect. III that such bifurcation scenarios shed considerable light on the phenomenology of large-scale atmospheric, oceanic and coupled ocean-atmosphere flows. Finally, in Sects. III.G and IV.E, we introduced stochastic effects into the nonlinearly deterministic setting recalled above, and outlined the role of pullback attractors (PBAs) and of the invariant measures they support in the theory of nonautonomous and random dynamical systems (NDSs and RDSs).

Recently, the interest in bifurcations in the climate and environmental sciences has greatly increased due to the introduction of the concept of tipping points from the social sciences (Gladwell, 2000; Lenton et al., 2008). Clearly, a tipping point sounds a lot more threatening than a bifurcation point, especially when dealing with a hurricane or a dramatic and irreversible climate change.

Beyond the linguistic effectiveness of the term, it does also generalize the bifurcation concept in the context of open systems that are modeled mathematically by NDSs or RDSs. As we saw in Sects. III and II, the climate system — as well as its subsystems, namely the atmosphere, oceans, biosphere and cryosphere — are open and exchange time-dependent fluxes of mass, energy and momentum with each other and with outer space. Hence it is quite appropriate to consider this generalization. Following Ashwin et al. (2012), one distinguishes among three kinds of tipping points, as shown in Table II.

The investigation of systems possessing multiple attractors — which may include fixed points, limit cycles, strange attractors as well as PBAs or random attractors — is an active area of research, encompassing mathematics, as well as the natural and socio-economic sciences. Feudel et al. (2018), for instance, provide a brief review that introduces a special issue of the journal *Chaos* on the topic and sketch several interesting examples.

This section is devoted to providing a general framework for the study of multistability in the Earth system, to forced transitions between different regimes when multistability is present, and to the critical transitions taking place at classical bifurcation and more general tipping points. Recall, for motivation, two cases of multistability in the climate system.

First, planet Earth as a whole has boundary conditions that arise from its space environment in the solar system, which are compatible with at least two competing, stable regimes: today's relatively ice-free and warm climate and a deep freeze or snowball climate; see Sect. III.A. Second, the oceans' THC is — over a rather broad range of the parameters that control its heat and freshwater budget — bistable, with the current active state coexisting with a greatly reduced or even reversed circulation; see Sect. III.C.

In order to illustrate the main mathematical and physical aspects of the problem, we choose at first the specific example of coexistence between the warmer, Holocene-like, and the much colder, snowball state of planet Earth,
<figure><img src="image_55.png" /><figcaption>FIG. 46: Schematic diagram of B-Tipping: as a control parameter is being changed slowly, the number of stable equilibria is reduced from two to one. Reproduced with permission from Lenton et al. (2008).</figcaption></figure>

as discussed in Sects. III.A and IV.E.5; see, in particular, Figs. 21 and 45a.

## A. Bistability for Gradient Flows and EBMs

The theory of SDEs and of RDSs provides a comprehensive framework for deriving the probability of occurrence of coexisting regimes in multistable systems and for estimating the probability of stochastic-forcing-triggered transitions between them. A good starting point is gradient flow in the presence of additive white noise.

The governing SDE in this case is:

$$dx = -\nabla_x V(x)dt + \epsilon dW, \quad (63)$$

where $x \in \mathbb{R}^d$, the potential $V: \mathbb{R}^d \to \mathbb{R}$ is sufficiently smooth, and $dW$ is a vector of $d$ independent increments of Brownian motion, while $\epsilon$ determines the strength of the noise. The Fokker-Planck equation associated with Eq. (63) describes the evolution of the PDF $p_\epsilon(x, t)$ of an ensemble of trajectories obeying the SDE.

The stationary solution corresponding to the latter equation's invariant measure $\lim_{t \to \infty} p_\epsilon(x, t) = p_\epsilon(x)$ is given by a large-deviation law (Varadhan, 1966):

$$p_\epsilon(x) = \frac{1}{Z} \exp \left( - \frac{2V(x)}{\epsilon^2} \right) \quad (64)$$

where $Z$ is a normalization constant. The local minima of the potential $V$ are the stable fixed points of the deterministic dynamics. One obtains them by setting $\epsilon = 0$ in the SDE (63), and they correspond to the local maxima of $p$. Thus, for instance, the double-well potential of Fig. 20 corresponds to a bimodal PDF.

The mountain pass lemma (e.g., Bisgard, 2015) states that the two minima of $V$ that give rise to the two stable

fixed points have to be separated by an unstable one of the saddle type, which is a maximum of $V$ in one direction and a minimum in all the other ones. Such a saddle looks like a mountain pass on a topographic map, hence the name of the lemma. Ghil and Childress (1987, Sec. 10.4) discussed its application to 1-D EBMs.

In the limit of weak noise, trajectories starting near a local minimum of $V$ located at $x = x_1$ typically wait a long time before moving to the neighborhood of a different local minimum of $V$. The transitions are most likely to occur through the lowest energy saddle $x = x_s$ that links the initial basin of attraction to any other basin. Let us call the corresponding local minimum $x = x_2$.

The optimal path linking $x_1$ to $x_s$ minimizes the Freidlin and Wentzell (1984) action and is called an instanton. In the weak-noise limit $\epsilon \to 0$, the persistence time inside the basin of attraction of $x = x_1$ is

$$\bar{\tau}_{\epsilon} \propto \exp \left( \frac{2(V(x_s) - V(x_1))}{\epsilon^2} \right), \quad (65)$$

which is referred to as the Kramers (1940) formula.

Assume now, for simplicity, that we deal with a bistable system, as in Fig. 20 of Sect. III.A. If the potential $V$ depends on a slowly varying control parameter $\phi$, a B-tipping would involve, for instance, the number of local minima decreasing from 2 to 1, as a local minimum merges with a saddle point at $\phi = \phi_c$. This merging is illustrated in Fig. 46 for a simple 1-D case, like that of the 0-D EBM governed by Eqs. (13)$^{11}$.

When the system nears the tipping point, the persistence time in the shallow, metastable minimum is reduced, because $\Delta V \equiv V_s - V_1 \to 0$. Another flag anticipating the tipping point is the increase of the autocorrelation time $\tau$ in the metastable state, which is proportional to the inverse of the second derivative of the potential evaluated at the minimum, $\bar{\tau}_{\epsilon} \propto 1/V_{x_1}''$.

These two easily observable phenomena are the simplest manifestation of the slowing-down process (Scheffer et al., 2009), which we have discussed already in Sect. IV.E.5 when describing the findings of Tantet et al. (2018). As indicated in the footnote there, the warning signals of tipping-point approach have been mostly studied in the B-tipping case for saddle-node bifurcations. Much remains to be done in this context for more complex situations that may involve N-tipping, R-tipping or global bifurcations rather than local ones. For instance, Ditlevsen and Johnsen (2010) analyzed a high-resolution ice core record and excluded the possibility of interpreting Dansgaard-Oeschger events as B-tipping by noting the lack of the two early warning signals mentioned above. Their results leave us with the alternative

<footnote n="11">11 Note that the number of dimensions in phase space is 1, while that of the model in physical space is 0.</footnote>
<table><caption>TABLE II: Tipping points in open systems; see also Ashwin et al. (2012) and Fig. 46.</caption><thead><tr><td>Type</td><td>Cause</td><td>Mechanism</td></tr></thead><tbody><tr><td>B-Tipping</td><td>Bifurcation-due tipping</td><td>Slow change in a parameter leads to the system's passage through a classical bifurcation</td></tr><tr><td>N-Tipping</td><td>Noise-induced tipping</td><td>Random fluctuations lead to the system's crossing an attractor basin boundary</td></tr><tr><td>R-Tipping</td><td>Rate-induced tipping</td><td>Due to rapid changes in the forcing, the system loses track of a slow change in its attractors</td></tr></tbody></table>

of interpreting the associated irregular oscillations (e.g., Boers et al., 2018) either as the result of N-tipping for an S-shaped back-to-back saddle-node bifurcation or as B-tipping for a Hopf bifurcation; see Fig. 29 and its discussion in Sect. III.E.3 for similar ambiguities in interpreting changes in atmospheric LFV.

Returning to the right-hand side of the 0-D EBM of Eq. (13a), it can be written as minus the derivative of a potential $V(T)$, cf. Ghil (1976). The gain factor $\Lambda$ given in Eq. (47), as well as the ECS, is proportional to the inverse of the second derivative of the potential $V$ evaluated at the local minimum defining the reference climate $T_0$, ECS $\propto \Lambda \propto 1/V''_{T_0}$. Hence, near the critical transition, a small forcing in the right direction can lead the system to jump to the other basin of attraction, even in the absence of noise. Note, though, that, near such a transition, the linear stability analysis performed in Eqs. (46, 47) is no longer valid, and an accurate quantitative analysis of climate response requires taking into account the essential nonlinearity in the problem; see the Zaliapin and Ghil (2010) comments on Roe and Baker (2007).

Determining the saddle point that potentially connects, via instantons, two local minima is not at all trivial when looking at high- or even infinite-dimensional systems. Figure 47 reports the findings of Bódai et al. (2015) on the Ghil (1976) 1-D EBM. Panel (a) shows the bifurcation diagram of the model, in which the control variable is the globally averaged surface air temperature $\bar{T}$, and the tuning parameter is the insolation parameter $\mu = S/S_0$ introduced in Eq. (13b).

Panel (b) shows the zonally averaged temperatures of the 1-D model's steady states, for the reference conditions $T = T(x)$ at $\mu = 1$. The red and blue lines represent the two stable solutions — warm ($W$) and snowball ($SB$), respectively — while the green line is the unstable solution lying in-between. Note that, at each tipping point, the unstable state comes in contact with a stable one, and then both disappear, in accordance with the scenario of basin boundary crisis (Ott, 2002). In fact, the lack of stability of the steady-state solution described by the green line is apparent by observing that $\bar{T}$ in Fig. 47a decreases with increasing solar radiation — as opposed to the $W$ and $SB$ states — which is clearly not physical (Ghil, 1976; North et al., 1981) and can be loosely interpreted as the result of a negative heat capacity; see

a)

<figure><img src="image_56.png" /></figure>

<figure><img src="image_57.png" /><figcaption>FIG. 47: Bistability for the Ghil-Sellers 1-D EBM. (a) Simplified bifurcation diagram for the Ghil (1976) 1-D EBM; compare Figs. 21 and 45. The red, blue, and green segments correspond to the globally averaged surface air temperature for the warm, snowball, and Melancholia state, respectively, as a function of the relative solar irradiance $\mu$; see also Arcoya et al. (1998). (b) Meridional temperature profile for the warm (red), snowball (blue), and $M$-state (green) for this model under present-day conditions ($\mu = 1$), with the sine of latitude on the $x$-axis; compare Ghil (1976, Fig. 3a). The circles indicate the $M$-state estimate obtained using the edge tracking method. Adapted with permission from Bódai et al. (2015).</figcaption></figure>

discussion in Bódai et al. (2015).
<figure><img src="image_58.png" /><figcaption>FIG. 48: Schematic diagram of the edge tracking algorithm proposed by Skufca *et al.* (2006) and Schneider *et al.* (2007), as applied to the 1-D EBM of Fig. 47. Reproduced with permission from Bódai *et al.* (2015).</figcaption></figure>

## B. Finding the Edge States

In one phase space dimension, every first-order system $\dot{x} = f(x; \mu)$ is a gradient system, since one can always write $V(x; \mu) = -\int f(\xi, \mu)d\xi$. This is not so in two or more dimensions and one might wonder whether it is possible to compute unstable solutions for more general multistable systems, with no gradient property. Scott *et al.* (1999) presented an example of such an unstable stationary solution in a simple box model of the THC featuring competing stable solutions and no gradient structure. Moreover, the actual computation of the unstable steady state in Ghil (1976, Fig. 3a) did not rely on the system's potential $V(T(x))$ but on its 1-D character in physical space and on a shooting method to solve the Sturm-Liouville equation that results when setting $\partial T(x, t)/\partial t \equiv 0$ in the Ghil-Sellers model. Neither of these approaches (Ghil, 1976; Scott *et al.*, 1999) can easily be extended to more general systems.

In the remainder of this section, we discuss a more general paradigm of multistability for deterministic systems and, subsequently, how this paradigm may be useful for studying the problem of noise-induced transitions between competing states, without the simplifying assumption of a gradient, as in Eq. (63).

One can formally introduce general multistable systems as follows. We consider a smooth autonomous continuous-time dynamical system acting on a smooth finite-dimensional compact manifold $M$, described by a ODE system of the form $\dot{x} = F(x)$. The system is multistable if it possesses more than one asymptotically stable state, defined by the attractors $\{\Omega_j : j = 1, \dots, J\}$. The asymptotic state of an orbit is determined by its initial state, and the phase space is partitioned between the basins of attraction $\{B_j\}$ of the attractors $\{\Omega_j\}$ and the

boundaries $\{\partial B_\ell : \ell = 1, \dots, L\}$ separating these basins.

The basin boundaries can be strange geometrical objects with codimension smaller than one. Orbits initialized on the basin boundaries $\{\partial B_\ell\}$ are attracted towards invariant saddles $\{\Pi_{\ell k}: \ell = 1, \dots, L, k = 1, \dots, k_\ell\}$, where we allow for the existence of $k_\ell$ separate saddles in the basin boundary $\partial B_\ell$; these saddles are often referred to as edge states (cf. Schneider *et al.*, 2007; Skufca *et al.*, 2006).<fnref n="12" /> For nongrادي flows, the edge states, like the asymptotic states, can feature chaotic dynamics (Grebogi *et al.*, 1983; Ott, 2002; Robert *et al.*, 2000; Vollmer *et al.*, 2009). Lucarini and Bódai (2017, 2019a,b) refer to the chaotic edge states as Melancholia ($M$) states. Hence, we cannot expect to find easily edge states and, a fortiori, $M$-states for such general multistable systems.

The edge tracking algorithm of B. Eckhart, J. A. Yorke and associates (Schneider *et al.*, 2007; Skufca *et al.*, 2006) allows one to do so by constructing a shadowing trajectory that leads an orbit starting on the basin boundary toward the corresponding edge state, cf. Fig. 48; therein one uses a bisection method to control the instability associated with the trajectories' diverging away from the basin boundary. Bódai *et al.* (2015) first used edge tracking in a geophysical context to reproduce the unstable solution of Ghil (1976) with this more easily generalizable approach; compare green line and circles in Fig. 47b.

Lucarini and Bódai (2017) computed $M$-states for an intermediate complexity climate model with $O(10^5)$ degrees of freedom that couples PUMA, an atmospheric primitive equations model (Frisius *et al.*, 1998) with a modified version of the 1-D EBM of Ghil (1976) that acts as a surrogate ocean and contributes to meridional heat transport. Figure 49 shows the bifurcation diagram of this climate model, where the $W \to SB$ tipping point is located near $\mu = 0.98$, while the $SB \to W$ tipping point is located near $\mu = 1.06$. Just as in the 1-D EBM of Fig. 47, the tipping points are associated with basin boundary crises (Ott, 2002).

Over a wide range of $\mu$-values, the edge state features chaotic dynamics that arises from the atmospheric model's baroclinic instability; it leads to weather variability and to a limited predictability horizon. Since this instability is much faster than the climatic one due to the ice-albedo feedback, the basin boundary is a fractal set with near-zero codimension, in agreement with results obtained in low-dimensional cases (Grebogi *et al.*, 1983; Lai and Tél, 2011). In other words, the basin boundary has almost full Lebesgue measure. As a result, near the basin boundary there is virtually no predictability on the asymptotic state of the system, because infinitesimally

<footnote n="12">¹² Note that, strictly speaking, edge states refer to solutions separating long-lived turbulence and laminar flow, rather than describing the interface between truly separated attractor basins.</footnote>
<figure><img src="image_59.png" /><figcaption>FIG. 49: Bifurcation diagram for the intermediate-complexity climate model in Lucarini and Bódai (2017), drawn for the long-term, globally averaged ocean temperatures [⟨T_s⟩]. Bistability is found for a large range of values of the control parameter μ = S/S_0. The relevant solid lines are red for the Warm W-states, blue for the Snowball SB-states, and green for the Melancholia M-states, as in Fig. 47. Reproduced with permission from Lucarini and Bódai (2019a).</figcaption></figure>

close initial states have a high probability of belonging to separate basins of attraction.

### C. Invariant Measures and Noise-induced Transitions

As mentioned above, transitions across the basin boundaries of competing attractors are possible and quite likely in the presence of stochastic perturbations. Noise-induced escape from an attractor has, in fact, long been studied in the natural sciences (e.g., Grassberger, 1989; Hanggi, 1986; Kautz, 1987).

Let us generalize Eq. (63) as follows:

$$dx = F(x)dt + \epsilon s(x)dW; \quad (66)$$

here $x \in \mathbb{R}^d$, $F(x)dt$ is the drift term given by a vector flow field that admits multiple steady states, as discussed in the previous two subsections, $dW$ are the increments of a $d$-dimensional Brownian motion. The volatility matrix $s(x) \in \mathbb{R}^{d \times d}$ is such that $s(x)^T s(x)$ is the covariance matrix of the noise. Finally, the parameter $\epsilon \ge 0$ controls the noise intensity.

Recent extensions of the classical Freidlin-Wentzell (Freidlin and Wentzell, 1984) theory (Graham et al., 1991; Hamm et al., 1994; Lai and Tél, 2011) yield results that mirror closely those summarized before in the case of gradient flows with additive noise, cf. Eqs. (64) and (65), given suitable assumptions on the drift term and the volatility (e.g., Lucarini and Bódai, 2019a,b). In the weak-noise limit given by $\epsilon \to 0$, the invariant measure of the system can be written as a large-deviation law. Formally, one just has to replace

$V(x)$ by a general pseudo-potential $\Phi(x)$, which depends in a nontrivial way on $F(x)$ and $s(x)$. Moreover, the constant $Z$ in Eq. (64) is replaced by a function $\zeta(x)$, which is relatively unimportant because the behavior of the system depends mostly on the properties of $\Phi(x)$.

Certain general properties of $\Phi(x)$ apply to all choices of the noise law and, once a noise law is chosen, one can derive how the properties of the system change as a function of the parameter $\epsilon$. In general, regardless of the noise law, $\Phi(x)$ has local minima supported on the attractors $\{\Omega_j : j = 1, \dots, J\}$ of the deterministic dynamics. Correspondingly, the invariant measure has local maxima on the attractors. Moreover, $\Phi(x)$ has a constant value on the support of each edge state $\Pi_\ell$ and each attractor $\Omega_j$, respectively.

In the simplest case of just $J = 2$ attractors and $L = 1$ edge states then, in the weak-noise limit, transitions away from either attractor basin take place exponentially more likely along the instanton connecting the corresponding attractor with the edge state $\Pi_1$. Just as for the gradient flows in Sect. V.A, instantons can be calculated as minimizers of a suitably defined action (Beri et al., 2005; Grassberger, 1989; Kautz, 1987; Kraut and Feudel, 2002).

Consider now the case that more than one edge state separate a given $\Omega_j$ from the other attractors $\{\Omega_k : k \ne j\}$. Let us denote then by $\Pi_j$ the edge state for which the value of $\Phi$ is lowest. In this case, the most probable exit path connects $\Omega_j$ with $\Pi_j$, while the other escape channels are basically switched off in the weak-noise limit.

Note that, if the attractors and the edge states are more complex sets than isolated points, the instantons connecting them are not unique, because in principle any point of the attractor can be linked by an instanton to any point of the edge state, because the pseudo-potential is constant on $\Omega_j$ and $\Pi_j$. Indeed, the edge states are the gateways for the noise-induced escapes from the deterministic basins of attraction, as illustrated in the bistable example with nongrادي flow given below.

Lucarini and Bódai (2019a) introduced stochastic forcing into the climate model studied by Lucarini and Bódai (2017) as a fluctuating factor of the form $1 + \epsilon dW/dt$, where $dW$ is a Brownian motion that modulates the solar insolation parameter $\mu$. This forcing yields a multiplicative noise law because the energy input into the system depends on the product of $\mu$ times the co-albedo $(1 - \alpha)$, where the albedo $\alpha = \alpha(T)$ depends explicitly on the surface air temperature $T = T(x, t)$, which is a state variable of the model; compare with Eq. (13b).

In Fig. 50 we reproduce the estimate by Lucarini and Bódai (2019a) of a 2-D projection of the invariant measure for $\mu = 1$ from a climate simulation of $\approx 6.0 \times 10^4$ yr with fluctuations in the solar insolation of 1.5% on the scale of 100 yr. The simulation features $92 SB \to W$ and $W \to SB$ transitions, and the measure is projected onto the ([T_S], ΔT_S)-plane; here $[T_S]$ is the globally averaged surface air temperature, and ΔT_S measures the temper-
<figure><img src="image_60.png" /></figure>

FIG. 50: Invariant measure of the transitions between the two stable regimes of the stochastically perturbed climate model of Lucarini and Bódai (2017) via the *M*-state separating them, for $\mu = 1$. Main graph: projection of this invariant measure on the reduced space ($T_S$, $[\Delta T_S]$); $W$ attractor – red dot; $SB$ attractor – blue dot; and $M$-state – green dot. Red and blue dashed lines plot the $W \to SB$ and the $SB \to W$ instanton, respectively. Top left inset: marginal PDF with respect to $\Delta T_S$; bottom right inset: marginal PDF with respect to $[T_S]$; and center right inset: probability along the two instantons. Reproduced with permission from Lucarini and Bódai (2019a).

ature difference between the surface temperature in the low- and high-latitude regions.

The peaks of the PDFs are very close to the $W$ and $SB$ attractors, and the agreement improves even further when considering the two marginal PDFs (top left and bottom right insets). Both the $W \to SB$ and the $SB \to W$ instantons can be estimated: their starting and final points agree remarkably well with the attractors and the $M$-state. The instantons follow a path of monotonic descent that tracks closely the crests of the PDF, with the minimum occurring at the $M$-state.

These features are all in excellent agreement with the theoretical predictions for multistable systems with a generalized pseudo-potential.

## D. Nearing Critical Transitions

The general framework outlined in the previous subsection is also quite useful for studying the properties of such systems near a critical transition (e.g., Lucarini and Bódai, 2019a). In Fig. 49, it is clear that $\mu = 0.98$ is close to the $W \to SB$ tipping point. One finds that at low noise intensities — i.e., at a relative $\mu$-fluctuation smaller or equal to 1.4% on a centennial scale — it is extremely hard to escape from the $SB$ state. This finding happens to agree with snowball state simulations that used more detailed models as well (e.g., Crowley *et al.*,

<figure><img src="image_61.png" /></figure>

FIG. 51: Escape from the $W$ attractor near the $W \to SB$ critical transition of the stochastically perturbed climate model of Lucarini and Bódai (2017), for $\mu = 0.98$. (a) The expected value of the transition time obeys Eq. (65) in the weak-noise regime; the diagram is in log-linear coordinates and $\sigma$ on the $x$-axis is the relative fluctuation of the solar insolation $\mu$. (b) Estimate of the $W \to SB$ instanton and, in the inset, empirical density; both used 50 trajectories that escape to the $SB$ state. Reproduced with permission from Lucarini and Bódai (2019a)

2001; Ghil, 2019; Pierrehumbert, 2004)

Let us then focus on the escape from the $W$-state. Lucarini and Bódai (2019a) estimated the expected value of the transition time from the $W$ to the $SB$ state using 50 simulations per noise strength value. These values grow exponentially with the inverse of the square of the parameter $\epsilon$, as predicted by the theory; see Fig. 51a. Note that the difference between the pseudo-potential value $\Phi$ at the $M$-state and at the $W$-state is equal to half the slope of the straight line, in agreement with Eq. (65).

Finally, the escape transition paths from the $W$-state to the $M$-state are plotted in Fig. 51b. In the weak-noise case of relative $\mu$ fluctuation smaller than 1% on a centennial scale, the highest densities of these paths lie quite close to the instanton connecting the $W$ attractor to the $M$-state, and follow a path of decreasing probability.

These encouraging results suggest that the underlying
ing methodology of edge tracking and instanton estimation could be applied to observational and reanalysis data sets. As discussed in Sect. III.E.3 the episodic, or ‘particle’, approach to atmospheric LFV results in a Markov chain of regimes {$R_j : j = 1, \dots, J$} that involves preferential transition paths {$\Theta_{j,k} : j = 1, \dots, J, k = 1, \dots, J$} between them. Heretofore, these transition paths, as well as the regimes themselves, were estimated by purely statistical methods; see Ghil et al. (2018) and Fig. 29a herein. Deloncle et al. (2007) and Kondrashov et al. (2007) applied a random-forest algorithm (Breiman, 2001) to find the best real-time predictors of the next regime in a Markov chain, conditional upon the one currently occupied. The results were quite satisfying for the QG3 model of Marshall and Molteni (1993) and encouraging for an observational data set of 55 winters of Northern Hemisphere 700-hPa geopotential height anomalies (Deloncle et al., 2007; Kondrashov et al., 2007, respectively).

## E. Chaos-to-Chaos Transition

So far, we have seen that critical transitions associated with saddle-node bifurcations and mild generalizations thereof are fairly well understood by now. We conclude this section and the main part of the text with a somewhat more exotic example of chaos-to-chaos transition for a delay differential equation (DDE) system with and without stochastic perturbations.

We saw in Sect. II.B.1 and III.E.2 that ENSO plays a key role in the global climate on interannual-to-interdecadal scales. Hence, a large number of relatively simple models thereof exist to better understand its main features. Important mechanisms involved are air-sea interaction, equatorial wave dynamics, and radiative forcing by the seasonal cycle (Bjerknes, 1969). In particular, the role of the wave dynamics has been captured by introducing one or two delays into the governing equations of some of the simpler models (e.g., Ghil et al., 2008a; Tziperman et al., 1994b, and references therein).

Chekroun et al. (2018) studied the PBA of the seasonally forced Tziperman et al. (1994b) model both with and without stochastic perturbations. The model has two delays, associated with a positive and a negative feedback; these delays are based on the basin-crossing times of the eastward-traveling Kelvin waves and the westward-traveling Rossby waves. The control parameter *a* is the intensity of the positive feedback and the PBA undergoes a crisis that consists of a chaos-to-chaos transition; the authors refer to it as a strange PBA since model behavior is chaotic within it in the purely deterministic case.

The changes in the invariant, time-dependent measure $\mu_t$ supported on this ENSO model's PBA are illustrated in Figs. 52a–52d, as a function of the control parameter *a*. The PBA experiences a critical transition at a

value $a_*$, as illustrated in Fig. 52, where $h(t)$ is the thermocline depth anomaly from seasonal depth values at the domain's eastern boundary, with $t$ in years. Here $a = (1.12 + \delta)/180$ and $0.015700 < \delta_* < 0.015707$.

The transition in the Kolmogorov-Smirnov metric of the invariant measure's dependence on the parameter *a*, i.e. in $\mu_t = \mu_t(a)$, is quite sharp, according to Chekroun et al. (2018, Fig. 3; not shown here). The singular support of the measure is in full agreement with rigorous mathematical results, as well as with the numerical results (Chekroun et al., 2011; Ghil, 2017) on the random attractor of the stochastically perturbed Lorenz (1963) convection model.

The change in the PBA is clearly associated with the population lying towards the ends of the elongated filaments apparent in Figs. 52a–52d. This population is due to the occurrence of large warm, El Niño and cold, La Niña events. Thus, $\mu_t(a)$ encripts faithfully the disappearance of such extreme events as $a \nearrow a_*$.

Note that Mukhin et al. (2015b,a) have used neural-network methodology to predict such transitions between two types of chaotic regimes in inverse models learned from time series simulated by ENSO models of increasing complexity. The ENSO models were subjected to a linear change in a parameter meant to represent gradual global warming. Strikingly, transitions between regimes in the inverse model did occur out of sample, even when the training interval did not contain both types of behavior.

Returning to the results of Chekroun et al. (2018), they found that perturbing the Tziperman et al. (1994b) model by small additive noise eliminates the crisis. The explanation of this numerical observation is tied to the role played in ENSO dynamics by the interaction between the intrinsic frequency of the coupled ocean-atmosphere system (e.g., Neelin et al., 1998) and the seasonal forcing (Ghil et al., 2008a; Jin et al., 1994, 1996; Tziperman et al., 1994b). This interaction induces a Devil's staircase in model frequency, which has plausible counterparts in observations (Ghil and Robertson, 2000).

As shown by Ghil et al. (2008, Appendix B), a Devil's staircase step that corresponds to a rational rotation number can be smoothed out by a sufficiently intense noise. In fact, the narrower a Devil's staircase step is, the less robust is it to noise perturbations, while the wider ones are the most robust. The effect of noise on the paradigmatic example of such a staircase, the standard circle map, has been examined in greater depth by Marangio et al. (2019).

This example is just one step on the long road of using the tools of nonautonomous and random dynamical systems for a better understanding of major climatic phenomena and processes. Note, however, that the model, while relatively simple, is actually infinite-dimensional because of the dependence of a DDE solution on a full interval of initial values on the real axis.
<figure><img src="image_62.png" /><figcaption>FIG. 52: Embedding of the invariant, time-dependent measure $\mu_t$ supported on the PBA associated with the highly idealized ENSO DDE model of Tziperman et al. (1994b). The embedding is shown within the $(h(t), h(t+1))$-plane for $a = (1.12 + \delta)/180$, $t \simeq 147.64$ yr and, respectively, (a) $\delta = 0.0$; and (b) $\delta = 0.01500$; (c) $\delta = 0.0157$; and (d) $\delta = 0.015707$. The red curves in the four panels represent the singular support of the measure. Reproduced with permission from Chekroun et al. (2018).</figcaption></figure>

## VI. CONCLUDING REMARKS

The goal of this review was to highlight some of the key physical and mathematical ingredients that can help address the description, understanding and prediction of climate variability and climate change. Complementary aspects of observations, theory and numerics have been taken into consideration. We have leaned heavily on dynamical systems theory and nonequilibrium statistical mechanics, and have tried to present a coherent picture of

the time dependence of the climate system, its multiscale nature, and its multistability.

We have emphasized the complex interplay between intrinsic climate variability and the climate's response to perturbations. The topic is, in fact, relevant for three problems of great scientific relevance: (a) anthropogenic climate change; (b) coevolution of the Earth's climate and of the biosphere, and (c) the quest for life on other planets, along with the habitability of our own planet.

The presentation also aimed to show the extent to
which basic mathematical and physical tools can help solve the main challenges inherent to the climate sciences. These challenges cannot be overcome merely by increas- ing the resolution of numerical models and including in them more and more physical and biogeochemical pro- cesses. In addition, a balanced interplay of observations, modeling and theory is definitely needed to achieve the necessary progress.

This review is far from exhaustive: the authors had to
make hard and, obviously, personal choices on the topics
to be covered. We would like to briefly mention here
several additional approaches to the problems at hand
that have been developed in recent years.

* *Network Theory* (e.g. Barrat *et al.*, 2008; Newman, 2010) has provided a novel viewpoint for constructing a parsimonious yet efficient representation of many complex processes taking place in the Earth system (e.g., Donges *et al.*, 2009; Tsonis and Roebber, 2004; Tsonis *et al.*, 2006). Gozolchiani *et al.* (2011) and Wang *et al.* (2013) provide examples of use of networks for capturing specific climatic processes, while Boers *et al.* (2014) an example of applying network theory to climate prediction.

* *Extreme Value Theory (EVT)* (e.g., Coles, 2001; Embrechts *et al.*, 1999) has been used extensively in studying the statistical properties of rare hydrometeorological events, such as the occurrence of very intense rain or of extreme temperatures (Ghil *et al.*, 2011; Katz *et al.*, 2002, and references therein). Yet this classical methodology has been almost entirely neglected in IPCC reports dedicated to the study of extremes, cf. IPCC (2012). Recently, though, rapid advances in EVT theory for general observables of chaotic deterministic dynamical systems (Holland *et al.*, 2012; Lucarini *et al.*, 2016, 2014b) have led to the derivation of indicators of weather regimes based on extreme value statistics in the recurrence of atmospheric fields. In addition to the mere classification of such regimes, as mentioned at the end of Sect. V.D, this approach allows one to infer higher or lower instability of these regimes and, hence, their lower or higher predictability (Faranda *et al.*, 2017; Hochman *et al.*, 2019).

* *Large Deviation Theory* (e.g., Touchette, 2009; Varadhan, 1984) has been first used in the context of geophysical fluid dynamics for studying the self-organization and the multistability of turbulent flows and, specifically, of jet structures (Bouchet *et al.*, 2014; Bouchet and Venaille, 2012); see also Sects. V.A and V.C herein. More recently, it has been applied successfully in studying weather and climate extremes. First, it has helped formulate rare events algorithms able to nudge a climate model towards representing preferentially the class

of extreme events of interest (Ragone *et al.*, 2018). Second, it has provided a solid theoretical and numerical basis for the study of spatially extended or temporally persistent temperature extremes (Gálfi *et al.*, 2019). Large Deviation Theory has also been recently used to study multiscale and coupled atmosphere-ocean instabilities in a hierarchy of climate models (De Cruz *et al.*, 2018; Vannitsem and Lucarini, 2016).

* *Detection and Attribution Studies.* In a complex system like the climate, inferring causal relationships among events and phenomena is far from obvious. Nonetheless, doing so for well-defined weather and climate events is essential for building simplified models and, at a practical level, for causal attribution of weather- and climate-related events. This is quite important in the case of detection and attribution studies of anthropogenic climate change, especially when aiming to go beyond changes in mean climate properties, such as globally averaged temperatures, and on to determining to what extent the occurrence of an individual extreme event — e.g., of a hurricane or an extended drought — can be attributed to climate change (Adam, 2011; Allen, 2003).

Methodologies based on the causal counterfactual
theory of Pearl (2009) are being increasingly recog-
nised as a key instrument for providing a more rig-
orous basis for detection and attribution studies
(e.g., Hannart et al., 2016a). They also seem better
suited for defining reliably the link between anthro-
pogenic climate forcing and individual events (Han-
nart and Naveau, 2018). Combining these methods
with those of data assimilation, discussed herein
in Sect. II.A, appears to be well suited to refine
the distinction between the factual and counterfac-
tual world that separates causation from the lack
thereof (e.g., Hannart et al., 2016b).

* *Beyond Linear Response Theory.* As shown in Sect. IV.E, linear response theory can provide a systematic improvement upon the standard methodology of forward integration of model ensembles with perturbed initial states and parameters. While it does apply to systems out of thermodynamic equilibrium, it is still limited by its linearity to fairly small perturbations in parameters. In Sect. IV.D, we mentioned that one can use the Wasserstein distance to measure arbitrary changes between two invariant measures, whether supported on a time-independent, classical attractor or a time-dependent PBA. One might thus wish to apply the PBA-based methodologies herein directly to observational data sets or to the simulations of IPCC-class models, either instead of or in combination
with previously tested statistical methods.

In Sect. V.E, we showed that the PBA of an intermediate, but still infinite-dimensional ENSO model can undergo a chaos-to-chaos transition that involves major changes in its invariant measure. Moreover, these changes could be connected to a physically quite significant change in model behavior, namely in the number and size of extreme events, i.e. of the largest warm and cold events. Thus exploring similarly interesting changes in model PBAs and in the time-dependent invariant measures supported by them appears to be a very promising road toward a deeper understanding of climate variability and its interaction with both natural and anthropogenic forcing.

In the text above, we have pointed to several open scientific problems; see also Ghil (2019) for a historical perspective. An important matter that was left practically untouched, aside from the brief discussion of ENSO and its prediction in Sec. III.E.2, is the tropical atmosphere. Half the surface of the Earth lies in the Tropics, between 30 degrees N and 30 degrees S, and much of humanity lives there. As one gets closer to the Equator, the quasi-geostrophic approximation collapses and moist processes take on a much greater importance. Four of the classical references are Krishnamurti *et al.* (1979); Palmer (1951); Riehl (1954) and Yanai (1975). The current literature is too plethoric to Marshall here.

More broadly, we have already alluded at the end of Sec. I to the lack of a proper definition of climate and to the efforts to provide a mathematically rigorous definition that would clearly distinguish it from weather. We have also mentioned in Sec. III.E.1 the distinction made by Von Neumann (1960) between the short-, medium- and long-term prediction problem in the climate sciences.

At an even higher level of abstraction, one might wish to recall the parallel drawn by H. Poincaré between the deterministic unpredictability of weather and that of the position of the planets on the ecliptic (Poincaré, 1902, pp. 69-70). Harking back to these two giants of both physics and mathematics, we could say that: (i) the astronomical ephemerides are the analog of weather and of its short-term prediction; (ii) the problem of climate corresponds to that of the long-term stability of the solar system; and (iii) the secular variations of the orbital configurations in the solar system are the appropriate counterpart of the slow evolution of the climate system's low-frequency modes of variability (Chown, 2004).

Making these analogies stick at a level that might help determine more reliably climate sensitivity to human activities is clearly a worthwhile effort for the best mathematicians, physicists and climate scientists.

## ACKNOWLEDGMENTS

We have benefited immensely from discussions with many colleagues and friends in the last few years, but we would like to thank especially T. Bódai, N. Boers, F. Bouchet, M. D. Chekroun, J. I. Díaz, H. A. Dijkstra, D. Faranda, A. Feigin, U. Feudel, K. Fraedrich, G. Gallavotti, A. Gritsun, A. Groth, D. Holm, F.-F. Jin, D. Kondrashov, S. Kravtsov, J. Kurths, J. C. McWilliams, J. D. Neelin, S. Pierini, F. Ragone, D. Ruelle, A. Speranza, H. Swinney, O. Talagrand, A. Tantet, S. Vannitsem, J. Yorke, and J. Wouters. Two anonymous reviewers provided constructive comments. The authors acknowledge the financial support provided by Horizon 2020 through the project TiPES (grant No. 820970). This activity has been supported by the EIT Climate-KIC. EIT Climate-KIC is supported by EIT, a body of the European Union. VL also acknowledges the financial support provided by Horizon 2020 through the projects CRESCENDO (grant No. 641816) and Blue-ACTION (grant No. 727852) and by the Royal Society (grant No. IEC R2 170001).
## Appendix A: Acronyms

This Appendix contains two tables of acronyms used throughout the paper: Table III contains the scientific acronyms and Table IV the institutional ones.

<table><caption>TABLE III: Scientific acronyms</caption><thead><tr><td>Acronym</td><td>Meaning</td></tr></thead><tbody><tr><td>AS</td><td>Arakawa-Schubert (parametrization)</td></tr><tr><td>CNs</td><td>Complex networks</td></tr><tr><td>DDE</td><td>Delay differential equation</td></tr><tr><td>DNS</td><td>Direct numerical simulation</td></tr><tr><td>EBM</td><td>Energy Balance Model</td></tr><tr><td>ECS</td><td>Equilibrium Climate Sensitivity</td></tr><tr><td>ENSO</td><td>El Niño-Southern Oscillation</td></tr><tr><td>EMS</td><td>Empirical Mode Reduction</td></tr><tr><td>FDT</td><td>Fluctuation-Dissipation Theorem</td></tr><tr><td>GCM</td><td>General Circulation Model</td></tr><tr><td>GCM</td><td>Global Climate Model</td></tr><tr><td>GFD</td><td>Geophysical fluid dynamics</td></tr><tr><td>GHGs</td><td>Greenhouse Gases</td></tr><tr><td>GLE</td><td>Generalized Langevin equation</td></tr><tr><td>LFV</td><td>Low Frequency Variability</td></tr><tr><td>MJO</td><td>Madden-Julian Oscillation</td></tr><tr><td>MOC</td><td>Meridional overturning circulation</td></tr><tr><td>MZ</td><td>Mori-Zwanzig</td></tr><tr><td>NAO</td><td>North Atlantic Oscillation</td></tr><tr><td>NDS</td><td>Nonautonomous dynamical system</td></tr><tr><td>NWP</td><td>Numerical Weather Prediction</td></tr><tr><td>NSEs</td><td>Navier-Stokes Equations</td></tr><tr><td>ODE</td><td>Ordinary differential equation</td></tr><tr><td>PBA</td><td>Pullback attractor</td></tr><tr><td>PDE</td><td>Partial differential equation</td></tr><tr><td>PNA</td><td>Pacific North American (pattern)</td></tr><tr><td>PSA</td><td>Pacific South American (pattern)</td></tr><tr><td>QG</td><td>Quasi-geostrophic (flow, model)</td></tr><tr><td>RDS</td><td>Random dynamical system</td></tr><tr><td>SDE</td><td>Stochastic differential equation</td></tr><tr><td>TCR</td><td>Transient Climate Response</td></tr><tr><td>THC</td><td>Thermohaline Circulation</td></tr></tbody></table>

<table><caption>TABLE IV: Institutional acronyms</caption><thead><tr><td>Acronym</td><td>Meaning</td></tr></thead><tbody><tr><td>AR</td><td>Assessment Report</td></tr><tr><td>CMIP</td><td>Climate Model Intercomparison Project</td></tr><tr><td>ECMWF</td><td>European Centre for<br/>Mid-range Weather Forecast</td></tr><tr><td>ENSO</td><td>El-Niño-Southern Oscillation</td></tr><tr><td>IPCC</td><td>Intergovernmental Panel on Climate Change</td></tr><tr><td>NCAR</td><td>National Center for Atmospheric Research</td></tr><tr><td>NCEP</td><td>National Center for Environmental Prediction</td></tr><tr><td>PCMDI</td><td>Program for Climate Model<br/>Diagnostics and Intercomparison</td></tr><tr><td>SPM</td><td>Summary for Policy Makers</td></tr><tr><td>UNEP</td><td>United Nations Environmental Programme</td></tr><tr><td>WCRP</td><td>World Climate Research Programme</td></tr><tr><td>WMO</td><td>World Meteorological Organization</td></tr></tbody></table>
REFERENCES

Abramov, R V, and A. J. Majda (2007), "Blended response algorithms for linear fluctuation-dissipation for complex nonlinear dynamical systems," Nonlinearity 20 (12), 2793-2821.

Adam, D (2011), "Climate change in court," Nature Climate Change 1 (3), 127–130.

Allen, M (2003), "Liability for climate change," Nature 421 (6926), 891–892.

Andersen, K K, and Coauthors (2004), “High-resolution record of Northern Hemisphere climate extending into the last interglacial period,” Nature 431 (7005), 147–151.

Andronov, A A, and L. S. Pontryagin (1937), “Systèmes grossiers,” Dokl. Akad. Nauk SSSR 14 (5), 247–250.

Arakawa, A, and W. H. Schubert (1974), “Interaction of a cumulus cloud ensemble with the large-scale environment, Part I,” Journal of the Atmospheric Sciences 31 (3), 674–701.

Arcoya, D, J. I. Diaz, and L. Tello (1998), “S-shaped bifurcation branch in a quasilinear multivalued model arising in climatology,” Journal of Differential Equations 150 (1), 215–225.

Arnold, L (1998), *Random Dynamical Systems* (Springer-Verlag, New York/Berlin).

Arnold, Ludwig (1988), *Random Dynamical Systems* (Springer, New York).

Arnold, V (1983), *Geometric Methods in the Theory of Ordinary Differential Equations* (Springer, New York/Berlin).

Arnold, VI (2003), *Catastrophe Theory* (Springer Nature, Berlin).

Ashwin, P. S. Wieczorek, R. Vitolo, and P. Cox (2012), “Tipping points in open systems: bifurcation, noise-induced and rate-dependent examples in the climate system,” Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 370 (1962), 1166–1184.

Baladi, V (2000), *Positive Transfer Operators and Decay of Correlations* (World Scientific, Singapore).

Baladi, V (2008), “Linear response despite critical points,” Nonlinearity 21 (6), T81.

Balmaseda, MA, F. Hernandez, A. Storto, M.D. Palmer, O. Alves, L. Shi, G.C. Smith, T. Toyoda, M. Valdivieso, B. Barnier, D. Behringer, T. Boyer, Y-S. Chang, G.A. Chepurin, N. Ferry, G. Forget, Y. Fujii, S. Good, S. Guine�ut, K. Haines, Y. Ishikawa, S. Keeley, A. Khl, T. Lee, M.J. Martin, S. Masina, S. Masuda, B. Meyssignac, K. Mogensen, L. Parent, K.A. Peterson, Y.M. Tang, Y. Yin, G. Vernieres, X. Wang, J. Waters, R. Wedd, O. Wang, Y. Xue, M. Chevallier, J-F. Lemieux, F. Dupont, T. Kuragano, M. Kamachi, T. Awaji, A. Caltabiano, K. Wilmer-Becker, and F. Gaillard (2015), “The ocean re-analyses intercomparison project (ora-ip),” Journal of Operational Oceanography 8 (sup1), s80–s97.

Barenblatt, G L (1987), *Dimensional Analysis* (Gordon and Breach, New York).

Barnes, E A, and J. A. Screen (2015), “The impact of Arctic warming on the midlatitude jetstream: Can it? Has it? Will it?” WIREs Climate Change 6, 277–286.

Barnston, A G, M. K. Tippett, M. L. Heureux, S. Li, and D. G. DeWitt (2012), “Skill of real-time seasonal ENSO model predictions during 2002–2011 — is our capability improving?” Bull. Am. Meteorol. Soc. 93 (5), 631–651.

Barrat, A, M. Barthelemy, and A. Vespignany (2008), *Dynamical Processes in Complex Networks* (Cambridge University Press).

Batchelor, G (1974), *Introduction to Fluid Dynamics* (Cambridge University Press, Cambridge, UK).

Beck, C (1990), “Brownian motion from deterministic dynamics,” Physica A: Statistical Mechanics and its Applications 169 (2), 324–336.

Bell, T L (1980), “Climate sensitivity from fluctuation dissipation: Some simple model tests,” Journal of the Atmospheric Sciences 37 (8), 1700–1707.

Bellenger, H. E. Guilyardi, J. Leloup, M. Lengaigne, and J. Vialard (2014), “Enso representation in climate models: from cmip3 to cmip5,” Climate Dynamics 42 (7), 1999–2018.

Bellman, R E, and K. L. Cooke (1963), *Differential-Difference Equations* (Rand Corporation, Santa Monica, CA).

Bengtsson, L, M. Ghil, and E. Källén (1981), *Dynamic Meteorology: Data Assimilation Methods* (Springer).

Bensid, S, and J. I. Diaz (2019), “On the exact number of monotone solutions of a simplified Budyko climate model and their different stability,” Discrete and Continuous Dynamical Systems - B 24, 1033.

Benzi, R, P. Malguzzi, A. Speranza, and A. Sutera (1986), “The statistical properties of general atmospheric circulation: observational evidence and a minimal theory of bi-modality,” Q. J. R. Meteorol. Soc. 112 (473), 661–674.

Benzi, R, and A. Speranza (1989), “Statistical properties of low-frequency variability in the northern hemisphere,” Journal of Climate 2 (4), 367–379.

Berger, A, and S. Siegmund (2003), “On the gap between random dynamical systems and continuous skew products,” Journal of Dynamics and Differential Equations 15 (2-3), 237–279.

Beri, S, R. Mannella, D. G. Luchinsky, A. N. Silchenko, and P. V. E. McClintock (2005), “Solution of the boundary value problem for optimal escape in continuous stochastic systems and maps,” Phys. Rev. E 72, 036131.

Berloff, P. A. Hogg, and W. K. Dewar (2007), “The turbulent oscillator: A mechanism of low-frequency variability of the wind-driven ocean gyres,” J. Phys. Oceanogr. 37, 2363–2386.

Berner, J, U. Achatz, L. Batté, L. Bengtsson, A. de la Cámara, H. M. Christensen, M. Colangeli, D. R. B. Coleman, D. Crommelin, S. I. Dolaptchiev, and C. L. Franzke (2017), “Stochastic parameterization: Toward a new view of weather and climate models,” Bulletin of the American Meteorological Society 98 (3), 565–588.

Bhattacharya, K, M. Ghil, and I. L. Vulis (1982), “Internal variability of an energy-balance model with delayed albedo effects,” Journal of the Atmospheric Sciences 39, 1747–1773.

Bisgard, J (2015), “Mountain passes and saddle points,” SIAM Review 57 (2), 275–292.

Bjerknes, J P (1969), “Atmospheric teleconnections from the equatorial Pacific,” Mon. Weather Rev. 97, 163–172.

Bjerknes, V (1904), “Das Problem der Wettervorhersage, betrachtet vom Standpunkte der Mechanik und der Physik,” Meteorologische Zeitschrift 21, 1–7.

Bódai, T, G. Károlyi, and T. Tél (2013), “Driving a conceptual model climate by different processes: Snapshot attractors and extreme events,” Phys. Rev. E 87, 022822.

Bódai, T, Gy. Károlyi, and T. Tél (2011), “A chaotically driven model climate: extreme events and snapshot attrac-
tors,” Nonlinear Processes in Geophysics **18** (5), 573–580.

Bódai, T, V. Lucarini, and F. Lunkeit (2020), “Can we use linear response theory to assess geoengineering strategies?” *Chaos: An Interdisciplinary Journal of Nonlinear Science* **30** (2), 023124, https://doi.org/10.1063/1.5122255.

Bódai, T, V. Lucarini, F. Lunkeit, and R. Boschi (2015), “Global instability in the Ghil–Sellers model,” *Clim. Dyn.* **44** (11), 3361–3381.

Bódai, T, and T. Tél (2012), “Annual variability in a con- ceptual climate model: Snapshot attractors, hysteresis in extreme events, and climate sensitivity,” *Chaos: An Inter- disciplinary Journal of Nonlinear Science* **22** (2), 023110.

Boers, N, B. Bookhagen, H. M. J. Barbosa, N. Marwan, J. Kurths, and J. A. Marengo (2014), “Prediction of ex- treme floods in the eastern central andes based on a com- plex networks approach,” *Nature Communications* **5**, 5199 EP –.

Boers, N, M. Ghil, and D.-D. Rousseau (2018), “Ocean circulation, ice shelf, and sea ice interactions explain Dansgaard–Oeschger cycles,” *Proceedings of the National Academy of Sciences* **115** (47), E11005–E11014.

Boiseau, M, M. Ghil, and A. Juillet-Leclerc (1999), “Trends and interdecadal variability from South-Central Pacific coral records,” *Geophys. Res. Lett.* **26**, 2881–2884.

Bond, G, W. Showers, M. Cheseby, R. Lotti, P. Almasi, P. de- Menocal, P. Priori, H. Cullen, I. Hajdas, and G. Bonani (1995), “A pervasive millennial-scale cycle in the North At- lantic Holocene and glacial climates,” *Science* **278**, 1257– 1265.

Bony, S, B. Stevens, D. M. W. Frierson, C. Jakob, M. Kageyama, R. Pincus, T. G. Shepherd, S. C. Sher- wood, A. P. Siebesma, A. H. Sobel, M. Watanabe, and M. J. Webb (2015), “Clouds, circulation and climate sen- sitivity,” *Nat Geosci* **8**, 10.1038/ngeo2398.

Boos, W R, and J. V. Hurley (2013), “Thermodynamic bias in the multimodel mean boreal summer monsoon,” *J. Climate* **26**.

Bouchet, F, J. Laurie, and O. Zaboronski (2014), “Langevin dynamics, large deviations and instantons for the quasi- geostrophic model and two-dimensional euler equations,” *Journal of Statistical Physics* **156** (6), 1066–1092.

Bouchet, F, and J. Sommeria (2002), “Emergence of intense jets and Jupiter’s Great Red Spot as maximum-entropy structures,” *J. Fluid Mech.* **464**, 165–207.

Bouchet, F, and A. Venaille (2012), “Statistical mechanics of two-dimensional and geophysical flows,” *Physics Reports* **515** (5), 227–295.

Branstator, G (1987), “A striking example of the atmo- sphere’s leading traveling pattern,” *J. Atmos. Sci.* **44** (16), 2310–2323.

Breiman, L (2001), “Random forests,” *Machine Learning* **45** (1), 5–32.

Broecker, W S (1991), “The great ocean conveyor,” *Oceanog- raphy* **4**, 79–89.

Bryan, F O (1986), “High-latitude salinity effects and inter- hemispheric thermohaline circulations,” *Nature* **323**, 301– 304.

Budivsić, M, R. Mohr, and I. Mezić (2012), “Applied Koop- manism,” *Chaos: An Interdisciplinary Journal of Nonlinear Science* **22** (4), 047510.

Budyko, M I (1969), “The effect of solar radiation variations on the climate of the Earth,” *Tellus* **21**, 611–619.

Cane, M A, and S. E. Zebiak (1985), “A theory for El Niño and the Southern Oscillation,” *Science* **228**, 1084–1087.

Caraballo, T, and X. Han (2017), *Applied Nonautonomous and Random Dynamical Systems: Applied Dynamical Systems* (Springer Science + Business Media).

Carrassi, A. M. Bocquet, L. Bertino, and G. Evensen (2018), “Data assimilation in the geosciences: An overview of methods, issues, and perspec- tives,” *WIREs Climate Change* **9** (5), e535, https://onlinelibrary.wiley.com/doi/pdf/10.1002/wcc.535.

Carton, J A, and B. S. Giese (2008), “A reanalysis of ocean climate using Simple Ocean Data Assimilation (SODA),” *Mon. Wea. Rev.* **136** (8), 2999–3017.

Carvalho, A N, J. Langa, and J. C. Robinson (2013), “The pullback attractor,” in *Attractors for infinite-dimensional non-autonomous dynamical systems*, Applied Mathemati- cal Sciences, Vol. 182 (Springer New York) pp. 3–22.

Cessi, P (2019), “The global overturning circulation,” *Annual Review of Marine Science* **11**, 249–270.

Cessi, P, and G. R. Ierley (1995), “Symmetry-breaking mul- tiple equilibria in quasi-geostrophic, wind-driven flows,” *J. Phys. Oceanogr.* **25**, 1196–1205.

Cessi, P, and W. R. Young (1992), “Multiple equilibria in two-dimensional thermohaline circulation,” *J. Fluid Mech.* **241**, 291–309.

Cvitanović, Predrag (1988), “Invariant measurement of strange sets in terms of cycles,” *Phys. Rev. Lett.* **61**, 2729– 2732.

Chang, C P, M. Ghil, M. Latif, and J. M. Wallace, Eds. (2015), *Climate Change: Multidecadal and Beyond* (World Scientific Publishing Co./Imperial College Press).

Chang, K-Il, M. Ghil, K. Ide, and C-C. A. Lai (2001), “Tran- sition to aperiodic variability in a wind-driven double-gyre circulation model,” *J. Phys. Oceanogr.* **31**, 1260–1286.

Charney, J G (1947), “The dynamics of long waves in a baro- clinic westerly current,” *J. Meteorol.* **4** (5), 136–162.

Charney, J G (1971), “Geostrophic turbulence,” *J. Atmos. Sci.* **28** (6), 1087–1095.

Charney, J G, A. Arakawa, D. J. Baker, B. Bolin, R. E. Dick- inson, R. M. Goody, C. E. Leith, H. M. Stommel, and C. I. Wunsch (1979), *Carbon Dioxide and Climate: A Scientific Assessment* (National Academy of Sciences, Washington, DC).

Charney, J G, and J. G. DeVore (1979), “Multiple flow equi- libria in the atmosphere and blocking,” *J. Atmos. Sci.* **36**, 1205–1216.

Charney, J G, R. Fjørtoft, and J. von Neumann (1950a), “Nu- merical integration of the barotropic vorticity equation,” *Tellus* **2**, 237–254.

Charney, J G, R. Fjørtoft, and J. Von Neumann (1950b), “Numerical integration of the barotropic vorticity equa- tion,” *Tellus* **2**, 38–54.

Charney, J G, M. Halem, and R. Jastrow (1969), “Use of incomplete historical data to infer the present state of the atmosphere,” *J. Atmos. Sci.* **26**, 1160–1163.

Charney, J G, J. Shukla, and K. C. Mo (1981), “Compari- son of a barotropic blocking theory with observation,” *J. Atmos. Sci.* **38** (4), 762–779.

Chavanis, P H, and J. Sommeria (1996), “Classification of self-organized vortices in two-dimensional turbulence: the case of a bounded domain,” *Journal of Fluid Mechanics* **314**, 267–297.

Chekroun, M D, M. Ghil, and J. D. Neelin (2018), “Pull- back attractor crisis in a delay differential ENSO model,” in *Advances in Nonlinear Geosciences*, edited by A. Tsonis (Springer) pp. 1–33.
Chekroun, M D, J. D. Neelin, D. Kondrashov, J. C. McWilliams, and M. Ghil (2014), “Rough parameter dependence in climate models and the role of Ruelle-Pollicott resonances,” Proceedings of the National Academy of Sciences 111 (5), 1684–1690.

Chekroun, M D, E. Simonnet, and M. Ghil (2011), “Stochastic climate dynamics: Random attractors and time-dependent invariant measures,” Physica D: Nonlinear Phenomena 240 (21), 1685–1700.

Chen, F, and M. Ghil (1996), “Interdecadal variability in a hybrid coupled ocean-atmosphere model,” J. Phys. Oceanogr. 26, 1561–1578.

Cheng, X, and J. M. Wallace (1993), “Cluster analysis of the Northern Hemisphere wintertime 500 pHa heigh field: Spatial patterns,” J. Atmos. Sci. 50, 2674–2696.

Chorin, A, and P. Stinis (2007), “Problem reduction, renormalization, and memory,” Communications in Applied Mathematics and Computational Science 1 (1), 1–27.

Chorin, A J, O. H. Hald, and R. Kupferman (2002), “Optimal prediction with memory,” Physica D: Nonlinear Phenomena 166, 239–257.

Chown, M (2004), “Chaotic heavens,” New Scientist 181 (2436), 32–35.

Cionni, I, G. Visconti, and F. Sassi (2004), “Fluctuation dissipation theorem in a general circulation model,” Geophysical Research Letters 31 (9), 10.1029/2004GL019739.

Cohen-Tannoudji, C, J. Dupont-Roc, and G. Grynberg (2007), *Photons and Atoms: Introduction to Quantum Electrodynamics* (Wiley, New York).

Coles, S (2001), *An Introduction to Statistical Modeling of Extreme Values* (Springer Science & Business Media).

Colon, C, D. Claessen, and M. Ghil (2015), “Bifurcation analysis of an agent-based model for predator–prey interactions,” *Ecological Modelling* 317, 93–106.

Compo, G P, J. S. Whitaker, P. D. Sardeshmukh, N. Matsui, R. J. Allan, X. Yin, B. E. Gleason, R. S. Vose, G. Rutledge, P. Bessemoulin, S. Brönnimann, M. Brunet, R. I. Crouthamel, A. N. Grant, P. Y. Groisman, P. D. Jones, M. C. Kruk, A. C. Kruger, G. J. Marshall, M. Maugeri, H. Y. Mok, O. Nordli, T. F. Ross, R. M. Trigo, X. L. Wang, S. D. Woodruff, and S. J. Worley (2011), “The twentieth century reanalysis project,” *Quarterly Journal of the Royal Meteorological Society* 137 (654), 1–28.

Cooper, F C, and P. H. Haynes (2011), “Climate sensitivity via a nonparametric fluctuation-dissipation theorem,” *Journal of the Atmospheric Sciences* 68 (5), 937–953.

Cox, P M, C. Huntingford, and M. S. Williamson (2018), “Emergent constraint on equilibrium climate sensitivity from global temperature variability,” *Nature* 553, 319–322.

Cronin, T N (2010), *Paleoclimates: Understanding Climate Change Past and Present* (Columbia University Press, New York).

Crowley, T J, W. T. Hyde, and W. R. Peltier (2001), “CO₂ levels required for deglaciation of a ”near-snowball” Earth,” *Geophysical Research Letters* 28 (2), 283–286.

Cushman-Roisin, B, and J.-M. Beckers (2011a), *Introduction to Geophysical Fluid Dynamics: Physical and Numerical Aspects* (Academic Press).

Cushman-Roisin, B, and J-M. Beckers (2011b), *Introduction to Geophysical Fluid Dynamics: Physical and Numerical Aspects* (Academic Press, 2nd Edition) 875 pp.

Cvitanović, P, and B Eckhardt (1991), “Periodic orbit expansions for classical smooth flows,” J. Phys. Math. Gen. 24.

Da Costa, E D, and A. C. Colin de Verdiére (2004), “The 7.7 year North Atlantic oscillation,” Q. J. R. Meteorol. Soc. 128A, 797–817.

Dansgaard, W, S. J. Johnsen, H. B. Clausen, D. Dahl-Jensen, N. S. Gundestrup, C. U. Hammer, C. S. Hvidberg, J. P. Steffensen, A. E. Sveinbjörnsdottir, J. Jouzel, and G. Bond (1993), “Evidence for general instability of past climate from a 250-kyr ice-core record,” Nature 364 (6434), 218–220.

Davini, P, and F. D'Andrea (2016), “Northern Hemisphere atmospheric blocking representation in global climate models: Twenty years of improvements?” Journal of Climate 29 (24), 8823–8840.

De Cruz, L, S. Schubert, J. Demaeyer, V. Lucarini, and S. Vannitsem (2018), “Exploring the Lyapunov instability properties of high-dimensional atmospheric and climate models,” *Nonlinear Processes in Geophysics* 25 (2), 387–412.

Dee, D P, and Coauthors (2011), “The ERA-Interim reanalysis: configuration and performance of the data assimilation system,” Q. J. R. Meteorol. Soc. 137 (656), 553–597.

Dell'Aquila, A, V. Lucarini, P. M. Ruti, and S. Calmanti (2005), “Hayashi spectra of the Northern Hemisphere mid-latitude atmospheric variability in the NCEP and ERA-40 reanalyses,” *Climate Dynamics* 25, 639–652.

Dell'Aquila, A, P. M. Ruti, S. Calmanti, and V. Lucarini (2007), “Southern Hemisphere Mid-latitude Atmospheric Variability of the NCEP-NCAR and ECMWF Reanalyses,” J. Geophys. Res. 112, D08106.

Deloncle, A, R. Berk, F. D'Andrea, and M. Ghil (2007), “Weather regime prediction using statistical learning,” *Journal of the Atmospheric Sciences* 64 (5), 1619–1635.

Delworth, T, S. Manabe, and R. J. Stouffer (1993), “Interdecadal variations of the thermohaline circulation in a coupled ocean-atmosphere model,” J. Climate 6, 1993–2011.

Delworth, T L, and M. E. Mann (2000), “Observed and simulated multidecadal variability in the Northern Hemisphere,” *Clim. Dyn.* 16, 661–676.

Demaeyer, J, and S. Vannitsem (2017), “Stochastic parametrization of subgrid-scale processes in coupled ocean-atmosphere systems: benefits and limitations of response theory,” *Quarterly Journal of the Royal Meteorological Society* 143 (703), 881–896.

Dickey, J O, M. Ghil, and S. L. Marcus (1991), “Extratropical aspects of the 40–50 day oscillation in length-of-day and atmospheric angular momentum,” J. Geophys. Res.: Atmospheres 96 (D12), 22643–22658.

Dijkstra, H A (2005), *Nonlinear Physical Oceanography: A Dynamical Systems Approach to the Large Scale Ocean Circulation and El Niño*, 2nd ed. (Springer Science+Business Media., Berlin-Heidelberg, Germany).

Dijkstra, H A (2007), “Characterization of the multiple equilibria regime in a global ocean model,” Tellus A 59, 695–705.

Dijkstra, H A (2013), *Nonlinear Climate Dynamics* (Cambridge University Press, Cambridge, UK).

Dijkstra, H A, and G. Burgers (2002), “Fluid mechanics of El Niño variability,” Annu. Rev. Fluid Mech. 34, 531–558.

Dijkstra, H A, and M. Ghil (2005), “Low-frequency variability of the large-scale ocean circulation: A dynamical systems approach,” *Reviews of Geophysics* 43 (3), RG3002.

Dijkstra, H A, and C. A. Katsman (1997), “Temporal variability of the wind-driven quasi-geostrophic double gyre ocean circulation: Basic bifurcation diagrams,” Geophys.
Astrophys. Fluid Dyn. 85, 195-232.

Ditlevsen, P D, and S. J. Johnsen (2010), “Tipping points: Early warning and wishful thinking,” Geophysical Research Letters **37** (19), 10.1029/2010GL044486.

Doblas-Reyes, Francisco J, Javier Garcia-Serrano, Fabian Lienert, Aida Pinto Biescas, and Luis R. L. Rodrigues (2013), “Seasonal climate predictability and forecasting: status and prospects,” Wiley Interdisciplinary Reviews: Climate Change **4** (4), 245–268.

Dobrushin, R (1970), “Prescribing a system of random variables by conditional distributions,” Theory of Probability & Its Applications **15** (3), 458–486.

Donges, J F, Y. Zou, N. Marwan, and J. Kurths (2009), “Complex networks in climate dynamics,” The European Physical Journal Special Topics **174** (1), 157–179.

Driver, R D (1977), *Ordinary and Delay Differential Equations* (Springer-Verlag, New York) re-issued in 2012 by Springer Science & Business Media.

Drótos, G, T. Bódai, and T. Tél (2015), “Probabilistic concepts in a changing climate: a snapshot attractor picture,” *Journal of Climate*, *Journal of Climate* 10.1175/JCLI-D-14-00459.1.

Duane, G S, C. Grabow, F. Selten, and M. Ghil (2017), “Introduction to focus issue: Synchronization in large networks and continuous media-data, models, and supermodels,” *Chaos* **27**, 126601.

Duplessy, J-C, and N. J. Shackleton (1985), “Response of global deep-water circulation to Earth’s climatic change 135,000-107,000 years ago,” *Nature* **316**, 500–507.

E, W, and J. Lu (2011), “Multiscale modeling,” *Scholarpedia* **6** (10), 11527, revision #91540.

Eade, R, D. Smith, A. Scaife, E. Wallace, N. Dunstone, L. Hermanson, and N. Robinson (2014), “Do seasonal-to-decadal climate predictions underestimate the predictability of the real world?” *Geophysical Research Letters* **41** (15), 5620–5628.

Eady, E T (1949), “Long waves and cyclone waves,” *Tellus* **1** (3), 33–52.

Eckmann, J P (1981), “Roads to turbulence in dissipative dynamical systems,” *Rev. Mod. Physics* **53**, 643–654.

Eckmann, J P, and D. Ruelle (1985), “Ergodic theory of chaos and strange attractors,” *Rev. Mod. Phys.* **57**, 617–656.

Edwards, P N (2010), *A Vast Machine: Computer Models, Climate Data, and the Politics of Global Warming* (MIT Press, Cambridge).

Egger, J (1978), “Dynamics of blocking highs,” *Journal of the Atmospheric Sciences* **35** (10), 1788–1801.

Einstein, A (1905), “Über die von der molekularinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen,” *Annalen der Physik* **322** (8), 549–560; reprinted in *Investigations on the Theory of the Brownian Movement, five articles by A. Einstein*, R. Furth (ed.) and A. D. Cowper (transl.), 1956, Dover Publ., New York, 122 pp.

Emanuel, K A (1994), *Atmospheric Convection* (Oxford University Press on Demand).

Embrechts, P, C. Klüppelberg, and T. Mikosch (1999), *Modelling Extremal Events for Insurance and Finance* (Springer-Verlag, New York).

Epstein, E S (1988), “Long-range weather prediction: Limits of predictability and beyond,” *Weather and Forecasting* **3**, 69–75.

Esper, J, E. R. Cook, and F. H. Schweingruber (2002), “Low-frequency signals in long tree-ring chronologies and the reconstruction of past temperature variability,” *Science* **295** (5563), 2250–2253.

Eyring, V, L. Bock, A. Lauer, M. Righi, M. Schlund, B. Andela, E. Arnone, O. Bellprat, B. Brötz, L.-P. Caron, N. Carvalhais, I. Cionni, N. Cortesi, B. Crezee, E. Davin, P. Davini, K. Debeire, L. de Mora, C. Deser, D. Docquier, P. Earnshaw, C. Ehbrecht, B. K. Gier, N. Gonzalez-Reviriego, P. Goodman, S. Hagemann, S. Hardiman, B. Hassler, A. Hunter, C. Kadow, S. Kindermann, S. Koirala, N. V. Koldunov, Q. Lejeune, V. Lembo, T. Lovato, V. Lucarini, F. Massonnet, B. Müller, A. Pandde, N. Pérez-Zanón, A. Phillips, V. Predoi, J. Russell, A. Sellar, F. Serva, T. Stacke, R. Swaminathan, V. Torralba, J. Vegas-Regidor, J. von Hardenberg, K. Weigel, and K. Zimmermann (2019), “Esmvaltool v2.0 – extended set of large-scale diagnostics for quasi-operational and comprehensive evaluation of earth system models in cmip,” *Geoscientific Model Development Discussions* **2019**, 1–81.

Eyring, V, S. Bony, G. A. Meehl, C. A. Senior, B. Stevens, R. J. Stouffer, and K. E. Taylor (2016a), “Overview of the Coupled Model Intercomparison Project Phase 6 (CMIP6) experimental design and organization,” *Geoscientific Model Development* **9** (5), 1937–1958.

Eyring, V, M. Righi, A. Lauer, M. Evaldsson, S. Wenzel, C. Jones, A. Anav, O. Andrews, I. Cionni, E. L. Davin, C. Deser, C. Ehbrecht, P. Friedlingstein, P. Gleckler, K.-D. Gottschaldt, S. Hagemann, M. Juckes, S. Kindermann, J. Krasting, D. Kunert, R. Levine, A. Loew, J. Mäkelä, G. Martin, E. Mason, A. S. Phillips, S. Read, C. Rio, R. Roehrig, D. Senftleben, A. Sterl, L. H. van Ulft, J. Walton, S. Wang, and K. D. Williams (2016b), “ESMValTool (v1.0) – a community diagnostic and performance metrics tool for routine evaluation of Earth system models in CMIP,” *Geoscientific Model Development* **9** (5), 1747–1802.

Faranda, D, G. Messori, and P. Yiou (2017), “Dynamical proxies of North Atlantic predictability and extremes,” *Scientific Reports* **7**, 41278 EP –.

Farrell, B F, and P. J. Ioannou (1996), “Generalized stability theory. I: Autonomous operators,” *J. Atmos. Sci.* **53**, 2025–2040.

Feliks, Y, M. Ghil, and A. W. Robertson (2010), “Oscillatory climate modes in the Eastern Mediterranean and their synchronization with the North Atlantic Oscillation,” *J. Climate* **23**, 4060–4079.

Feliks, Y, M. Ghil, and A. W. Robertson (2011), “The atmospheric circulation over the North Atlantic as induced by the SST field,” *J. Climate* **24**, 522–542.

Feliks, Y, M. Ghil, and E. Simonnet (2004), “Low-frequency variability in the mid-latitude atmosphere induced by an oceanic thermal front,” *J. Atmos. Sci.* **61**, 961–981.

Feliks, Y, M. Ghil, and E. Simonnet (2007), “Low-frequency variability in the mid-latitude baroclinic atmosphere induced by an oceanic thermal front,” *J. Atmos. Sci.* **64**, 97–116.

Fenichel, N (1979), “Geometric singular perturbation theory for ordinary differential equations,” *Journal of Differential Equations* **31** (1), 53–98.

Ferranti, L, S. Corti, and M. Janousek (2015), “Flow-dependent verification of the ECMWF ensemble over the Euro-Atlantic sector,” *Quarterly Journal of the Royal Meteorological Society* **141** (688), 916–924.
Feudel, U, A. N. Pisarchik, and K. Showalter (2018), "Multistability and tipping: From mathematics and physics to climate and brain—minireview and preface to the focus issue," Chaos: An Interdisciplinary Journal of Nonlinear Science 28 (3), 033501.

Fraedrich, K (2012), "A suite of user-friendly global climate models: hysteresis experiments," The European Physical Journal Plus 127 (5), 53.

Fraedrich, K, and H. Bottger (1978), "A wavenumber frequency analysis of the 500 mb geopotential at 50° n," J. Atmos. Sci. 35, 745-750.

Fraedrich, K, H. Jansen, E. Kirk, U. Luksch, and F. Lunkeit (2005), "The planet simulator: Towards a user friendly mode," Meteorologische Zeitschrift 14 (3), 299-304.

Francis, J A, and S. J. Vavrus (2012), "Evidence linking Arctic amplification to extreme weather in mid-latitudes," Geophysical Research Letters 39 (6).

Franzke, C L E, T. J. O'Kane, J. Berner, P. D. Williams, and V. Lucarini (2015), "Stochastic climate theory and modeling," Wiley Interdisciplinary Reviews: Climate Change 6 (1), 63-78.

Freidlin, M I, and A.D. Wentzell (1984), *Random Perturbations of Dynamical Systems* (Springer, New York).

Frisius, T, F. Lunkeit, K. Fraedrich, and I. N. James (1998), "Storm-track organization and variability in a simplified atmospheric global circulation model," Quarterly Journal of the Royal Meteorological Society 124 (548), 1019-1043.

Gálfi, V M, V. Lucarini, and J. Wouters (2019), "A large deviation theory-based analysis of heat waves and cold spells in a simplified model of the general circulation of the atmosphere," Journal of Statistical Mechanics: Theory and Experiment 2019 (3), 033404.

Gallavotti, G, and E. G. D. Cohen (1995), "Dynamical ensembles in stationary states," Journal of Statistical Physics 80 (5-6), 931-970.

Ghil, M (1976), "Climate stability for a Sellers-type model," J. Atmos. Sci. 33, 3-20.

Ghil, M (1989), "Meteorological data assimilation for oceanographers. Part I: Description and theoretical framework," Dyn. Atmos. Oceans 13, 171-218.

Ghil, M (1994), "Cryothermodynamics: the chaotic dynamics of paleoclimate," Physica D 77, 130-159.

Ghil, M (2001), "Hilbert problems for the geosciences in the 21st century," Nonlin. Processes Geophys. 8, 211-222.

Ghil, M (2002), "Natural climate variability," in *Encyclopedia of Global Environmental Change, Vol. 1*, edited by T. E. Munn, M. MacCracken, and J. Perry (J. Wiley and Sons, Chichester/New York) pp. 544-549.

Ghil, M (2015), "A mathematical theory of climate sensitivity or, How to deal with both anthropogenic forcing and natural variability?" in *Climate Change : Multidecadal and Beyond*, edited by C. P. Chang, M. Ghil, M. Latif, and J. M. Wallace (World Scientific Publishing Co./Imperial College Press) pp. 31-51.

Ghil, M (2017), "The wind-driven ocean circulation: Applying dynamical systems theory to a climate problem," Discr. Cont. Dyn. Syst.-A 37 (1), 189-228.

Ghil, M (2019), "A century of nonlinearity in the geosciences," Earth and Space Science 6, 1007-1042.

Ghil, M, M. R. Allen, M. D. Dettinger, K. Ide, D. Kondrashov, M. E. Mann, A. W. Robertson, A. Saunders, Y. Tian, F. Varadi, and P. Yiou (2002), "Advanced spectral methods for climatic time series," Rev. Geophys. 40 (1), 31-341.

Ghil, M, M. D. Chekroun, and E. Simonnet (2008), "Climate dynamics and fluid mechanics: Natural variability and related uncertainties," Physica D: Nonlinear Phenomena 237 (14-17), 2111-2126.

Ghil, M, M. D. Chekroun, and G. Stepan (2015), "A collection on 'Climate dynamics: multiple scales and memory effects, Introduction," R. Soc. Proc. A 471, 20150097.

Ghil, M, and S. Childress (1987), *Topics in Geophysical Fluid Dynamics: Atmospheric Dynamics, Dynamo Theory, and Climate Dynamics* (Springer-Verlag, Berlin/Heidelberg/New York).

Ghil, M, A. Groth, D. Kondrashov, and A. W. Robertson (2018), "Extratropical sub-seasonal-to-seasonal oscillations and multiple regimes: The dynamical systems view," in *The Gap Between Weather and Climate Forecasting: Sub-Seasonal to Seasonal Prediction*, edited by A. W. Robertson and F. Vitart, Chap. 6 (Elsevier) in press.

Ghil, M, M. Halem, and R. Atlas (1979), "Time-continuous assimilation of remote-sounding data and its effect on weather forecasting," Mon. Wea. Rev. 107, 140-171.

Ghil, M, M. Kimoto, and J. D. Neelin (1991), "Nonlinear dynamics and predictability in the atmospheric sciences," Rev. Geophys., Supplement (U.S. Nat'l Rept. to Int'l Union of Geodesy and Geophys. 1987-1990) 29, 46-55.

Ghil, M, and P. Malanotte-Rizzoli (1991), "Data assimilation in meteorology and oceanography," in *Advances in Geophysics*, Vol. 33 (Academic Press, New York) pp. 141-266.

Ghil, M, and K. C. Mo (1991), "Intraseasonal oscillations in the global atmosphere. Part 1: Northern Hemisphere and tropics," J. Atmos. Sci. 48, 752-779.

Ghil, M, A Mulhaupt, and P Pestiaux (1987), "Deep water formation and Quaternary glaciations," Climate Dynamics 2, 1-10.

Ghil, M, and A. W. Robertson (2000), "Solving problems with GCMs: General circulation models and their role in the climate modeling hierarchy," in *General Circulation Model Development: Past, Present and Future*, edited by D. A. Randall (Academic Press, New York) pp. 285-325.

Ghil, M, and A. W. Robertson (2002), "'Waves' vs. 'particles'" in the atmosphere's phase space: A pathway to long-range forecasting?" Proc. Natl. Acad. Sci. USA 99, 2493-2500.

Ghil, M, P. Yiou, S. Hallegatte, B. D. Malamud, P. Naveau, A. Soloviev, P. Friederichs, V. Keilis-Borok, D. Kondrashov, V. Kossobokov, O. Mestre, C. Nicolis, H. W. Rust, P. Shebalin, M. Vrac, A. Witt, and I. Zaliapin (2011), "Extreme events: dynamics, statistics and prediction," Nonlinear Processes in Geophysics 18 (3), 295-350.

Ghil, M, I. Zaliapin, and S. Thompson (2008a), "A delay differential model of ENSO variability: parametric instability and the distribution of extremes," Nonlinear Processes in Geophysics 15 (3), 417-433.

Gildor, H, and E. Tziperman (2001), "A sea ice climate switch mechanism for the 100-kyr glacial cycles," J. Geophys. Res. 106, 9117-9133.

Gill, A E (1982), *Atmosphere-Ocean Dynamics* (Academic Press, New York, U.S.A.).

Ginelli, F, P. Poggi, A. Turchi, H. Chaté, R. Livi, and A. Politi (2007), "Characterizing dynamics with covariant lyapunov vectors," Phys. Rev. Lett. 99, 130601.

Gladwell, M (2000), *The Tipping Point: How Little Things Can Make a Big Difference* (Little Brown, Boston).

Gleckler, P J, C. Doutriaux, P. J. Durack, K. E. Taylor, Y. Zhang, D. N. Williams, E. Mason, and J. Servonnat
(2016), "A more powerful reality test for climate models," Eos Trans. AGU **97**.

Gómez-Leal, I., L. Kaltenegger, V. Lucarini, and F. Lunkeit (2018), "Climate sensitivity to carbon dioxide and the moist greenhouse threshold of earth-like planets under an increasing solar forcing," The Astrophysical Journal **869** (2), 129.

Goody, R (2000), "Sources and sinks of climate entropy," Q. J. R. Meteorol. Soc. **126**, 1953-1970.

Gozolchiani, A., S. Havlin, and K. Yamasaki (2011), "Emergence of El Niño as an autonomous component in the climate network," Physical Review Letters **107** (14), 148501, arXiv:1010.2605 [physics.ao-ph].

Graham, R., A. Hamm, and T. Tél (1991), "Nonequilibrium potentials for dynamical systems with fractal attractors or repellers," Phys. Rev. Lett. **66**, 3089-3092.

Grasman, J (1987), *Asymptotic Methods for Relaxation Oscillations and Applications* (Springer Science & Business Media).

Grassberger, P (1989), "Noise-induced escape from attractors," Journal of Physics A: Mathematical and General **22** (16), 3283.

Grebogi, C., E. Ott, and J. A. Yorke (1983), "Fractal basin boundaries, long-lived chaotic transients, and unstable-unstable pair bifurcation," Physical Review Letters **50**, 935-938.

Gritsun, A., G. Branstator, and A. J. Majda (2008), "Climate response of linear and quadratic functionals using the fluctuation-dissipation theorem," J. Atmos. Sci. **65**.

Gritsun, Andrey, and Grant Branstator (2007), "Climate response using a three-dimensional operator based on the fluctuation-dissipation theorem," Journal of the Atmospheric Sciences **64** (7), 2558-2575.

Gritsun, Andrey, and Valerio Lucarini (2017), "Fluctuations, response, and resonances in a simple atmospheric model," Physica D: Nonlinear Phenomena **349**, 62-76.

de Groot, S R, and P Mazur (1984), *Non-Equilibrium Thermodynamics* (Dover).

Guckenheimer, J., and P. Holmes (1983), *Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields*, 2nd ed. (Springer-Verlag, Berlin/Heidelberg).

Guckenheimer, John, Kathleen Hoffman, and Warren Weckesser (2003), "The forced van der Pol equation I: The slow flow and its bifurcations," SIAM Journal on Applied Dynamical Systems **2** (1), 1-35.

Haines, K (1994), "Low-frequency variability in atmospheric middle latitudes," Surveys in Geophysics **15**, 1-61.

Hale, J K (1977), *Theory of Functional Differential Equations* (Springer-Verlag, New York).

Hamm, A., T. Tél, and R. Graham (1994), "Noise-induced attractor explosions near tangent bifurcations," Physics Letters A **185** (3), 313-320.

Hanggi, P (1986), "Escape from a metastable state," Journal of Statistical Physics **42** (1), 105-148.

Hannart, A., A. Carrassi, M. Bocquet, M. Ghil, P. Naveau, M. Pulido, J. Ruiz, and P. Tandeo (2016b), "DADA: data assimilation for the detection and attribution of weather and climate-related events," Climatic Change **136** (2), 155-174.

Hannart, A., and P. Naveau (2018), "Probabilities of causation of climate changes," Journal of Climate **31** (14), 5507-5524, https://doi.org/10.1175/JCLI-D-17-0304.1.

Hannart, A., J. Pearl, F. E. L. Otto, P. Naveau, and M. Ghil (2016a), "Causal counterfactual theory for the attribution

of weather and climate-related events," Bull. Amer. Meteorol. Soc. **97** (1), 99-110.

Hassanzadeh, P., Z. Kuang, and B. F. Farrell (2014), "Responses of midlatitude blocks and wave amplitude to changes in the meridional temperature gradient in an idealized dry GCM," Geophysical Research Letters **41** (14), 5223-5232.

Hasselmann, K (1976), "Stochastic climate models. I: Theory," Tellus **28**, 473-485.

Hasselmann, K., R. Sausen, E. Maier-Reimer, and R. Voss (1993), "On the cold start problem in transient simulations with coupled atmosphere-ocean models," Climate Dynamics **9** (2), 53-61.

Hasson, S., V. Lucarini, and S. Pascale (2013), "Hydrological cycle over South and Southeast Asian river basins as simulated by PCMDI/CMIP3 experiments," Earth System Dynamics **4** (2), 199-217.

Hawkins, E., R. S. Smith, L. C. Allison, J. M. Gregory, T. J. Woollings, H. Pohlmann, and B. Cuevas (2011), "Bistability of the atlantic overturning circulation in a global climate model and links to ocean freshwater transport," Geophysical Research Letters **38** (1), copyright - Copyright 2011 by American Geophysical Union; Last updated - 2014-08-09.

Hayashi, Y (1971), "A generalized method for resolving disturbances into progressive and retrogressive waves by space Fourier and time cross-spectral analysis," J. Meteorol. Soc. Jpn. **49**, 125-128.

Heinrich, H (1988), "Origin and Consequences of Cyclic Ice Rafting in the Northeast Atlantic Ocean during the Past 130,000 Years," Quaternary Research **29**, 142-152.

Held, I M (2001), "The partitioning of the poleward energy transport between the tropical ocean and atmosphere," J. Atmos. Sci. **58** (8), 943-948.

Held, I M (2005), "The gap between simulation and understanding in climate modeling," Bull. Am. Meteorol. Soc. **86**, 1609-1614.

Held, I M, and M. J. Suarez (1974), "Simple albedo feedback models of the ice caps," Tellus **26**, 613-629.

von der Heydt, Anna S, Henk A. Dijkstra, Roderik S. W. van de Wal, Rodrigo Caballero, Michel Crucifix, Gavin L. Foster, Matthew Huber, Peter Köhler, Eelco Rohling, Paul J. Valdes, Peter Ashwin, Sebastian Bathiany, Tijn Berends, Loes G. J. van Bree, Peter Ditlevsen, Michael Ghil, Alan M. Haywood, Joel Katzav, Gerrit Lohmann, Johannes Lohmann, Valerio Lucarini, Alice Marzocchi, Heiko Pälike, Itzel Ruvalcaba Baroni, Dirk Simon, Appy Sluijs, Lennert B. Stap, Alexis Tantet, Jan Viebahn, and Martin Ziegler (2016), "Lessons on climate sensitivity from past climate changes," Current Climate Change Reports **2** (4), 148-158.

Hochman, A., P. Alpert, T. Harpaz, H. Saaroni, and G. Messori (2019), "A new dynamical systems perspective on atmospheric predictability: Eastern Mediterranean weather regimes as a case study," Science Advances **5** (6), 10.1126/sciadv.aau0936.

Hoffman, P F., A. J. Kaufman, G. P. Halverson, and D. P. Schrag (2002), "On the initiation of a snowball earth," Science **281**, 1342.

Hoffman, P F, and D. P. Schrag (2002), "The snowball earth hypothesis: testing the limits of global change," Terra Nova **14**, 129.

Holland, Mark P, Renato Vitolo, Pau Rabassa, Alef E. Sterk, and Henk W. Broer (2012), "Extreme value laws in dynamical systems under physical observables," Physica D:
Nonlinear Phenomena **241** (5), 497–513.

Holton, J R, and G. J. Hakim (2013), *An Introduction to Dynamic Meteorology, 4th ed.* (Academic Press, San Diego, CA).

Hoskins, B J, M. E. McIntyre, and A. W. Robertson (1985), “On the use and significance of isentropic potential vorticity maps,” Quarterly Journal of the Royal Meteorological Society **111**, 877–946.

Hu, G, T. Bódai, and V. Lucarini (2019), “Effects of stochastic parametrization on extreme value statistics,” Chaos: An Interdisciplinary Journal of Nonlinear Science **29** (8), 083102, https://doi.org/10.1063/1.5095756.

Huybers, P (2005), “Comment on ”Hockey sticks, principal components, and spurious significance” by S. McIntyre and R. McKittrick,” Geophys. Res. Lett. **32** (20), L20705.

Imbrie, J, and K. P. Imbrie (1986), *Ice Ages: Solving the Mystery* (2nd Edn., Harvard Univ. Press, Cambridge, Mass.).

IPCC, (2001), *Climate Change 2001: The Scientific Basis. Contribution of Working Group I to the Third Assessment Report of the Intergovernmental Panel on Climate Change*, edited by J. T. Houghton et al. (Cambridge University Press, Cambridge, UK).

IPCC, (2001), *Third Assessment Report of the Intergovernmental Panel on Climate Change*, edited by Houghton et al. (Cambridge University Press).

IPCC, (2007), *Climate Change 2007 - The Physical Science Basis: Working Group I Contribution to the Fourth Assessment Report of the IPCC*, edited by S. Solomon et al. (Cambridge University Press, Cambridge, UK and New York, NY, USA).

IPCC, (2012), *Managing the Risks of Extreme Events and Disasters to Advance Climate Change Adaptation. A Special Report of Working Groups I and II of the Intergovernmental Panel on Climate Change*, edited by C.B. Field et al. (Cambridge University Press, Cambridge, UK, and New York, USA).

IPCC, (2014a), *Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change*, edited by T. Stocker et al. (Cambridge University Press, Cambridge).

IPCC, (2014b), *Climate Change 2014. Impacts, Adaptation and Vulnerability. Contribution of Working Group II to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change*, edited by C.B. Field et al. (Cambridge University Press, Cambridge).

IPCC, (2014c), *Climate Change 2014. Mitigation of Climate Change. Contribution of Working Group III to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change*, edited by O. Edenhofer et al. (Cambridge University Press, Cambridge).

Itoh, H, and M. Kimoto (1996), “Multiple attractors and chaotic itinerancy in a quasigeostrophic model with realistic topography: Implications for weather regimes and low-frequency variability,” J. Atmos. Sci. **53** (15), 2217–2231.

Itoh, H, and M. Kimoto (1997), “Chaotic itinerancy with preferred transition routes appearing in an atmospheric model,” Physica D **109** (3-4), 274–292.

Jiang, N, J. D. Neelin, and M. Ghil (1995a), “Quasiquadrennial and quasi-biennial variability in the equatorial Pacific.,” Clim. Dyn. **12**, 101–112.

Jiang, S, F.-F. Jin, and M. Ghil (1995b), “Multiple equilibria and aperiodic solutions in a wind-driven double-gyre, shallow-water model,” J. Phys. Oceanogr. **25**, 764–786.

Jin, F-F, and M. Ghil (1990), “Intraseasonal oscillations in the extratropics: Hopf bifurcation and topographic instabilities,” J. Atmos. Sci. **47** (24), 3007–3022.

Jin, F-F, and J. D. Neelin (1993), “Modes of interannual tropical ocean-atmosphere interaction - a unified view. I: Numerical results.” J. Atmos. Sci. **50**, 3477–3503.

Jin, F-F, J. D. Neelin, and M. Ghil (1994), “El Niño on the devil’s staircase: Annual subharmonic steps to chaos,” Science **264**, 70–72.

Jin, F-F, J. D. Neelin, and M. Ghil (1996), “El Niño/Southern Oscillation and the annual cycle: Subharmonic frequency-locking and aperiodicity,” Physica D **98**, 442–465.

Jones, C K R T (1995), “Geometric singular perturbation theory,” in *Dynamical Systems*, Lecture Notes in Mathematics, Vol. 1609, edited by R. Johnson (Springer) pp. 44–118.

Jouzel, J, N. I. Barkov, J. M. Barnola, M. Bender, J. Chappelaz, C. Genthon, V. M. Kotlyakov, V. Lipenkov, C. Lorius, J. R. Petit, D. Raynaud, G. Raisbeck, C. Ritz, T. Sowers, M. Stievenard, F. Yiou, and P. Yiou (1991), “Extending the Vostok ice-core record of paleoclimate to the penultimate glacial period.” Nature **364**, 407–412.

Just, W, H. Kantz, C. Rödenbeck, and M. Helm (2001), “Stochastic modelling: replacing fast degrees of freedom by noise,” Journal of Physics A: Mathematical and General **34** (15), 3199–3213.

Kalmay, E (2003), *Atmospheric Modeling, Data Assimilation and Predictability* (Cambridge University Press, Cambridge, UK).

Karamperidou, C, and Coauthors (2015), “The response of ENSO flavors to mid-Holocene climate: Implications for proxy interpretation,” Palaeoceanography **30**, 527–547.

Karspeck, Alicia R, Gokhan Danabasoglu, Jeffrey Anderson, Svetlana Karol, Nancy Collins, Mariana Vertenstein, Kevin Raeder, Tim Hoar, Richard Neale, Jim Edwards, and Anthony Craig (2018), “A global coupled ensemble data assimilation system using the community earth system model and the data assimilation research testbed,” Quarterly Journal of the Royal Meteorological Society **144** (717), 2404–2430.

Katsman, C A, H. A. Dijkstra, and S. S. Drijfhout (1998), “The rectification of the wind-driven circulation due to its instabilities,” J. Mar. Res. **56**, 559–587.

Katz, RW, M.B. Parlange, and P. Naveau (2002), “Statistics of extremes in hydrology,” Advances in Water Resources **25**, 1287–1304.

Kautz, RL (1987), “Activation energy for thermally induced escape from a basin of attraction,” Physics Letters A **125** (6), 315–319.

Kennett, J P, and L. D. Stott (1991), “Abrupt deep-sea warming, palaeoceanographic changes and benthic extinctions at the end of the Palaeocene,” Nature **353**, 225–229.

Kharin, VV, F.W. Zwiers, and X. Zhang (2005), “Intercomparison of near surface temperature and precipitation extremes in AMIP-2 simulations, reanalyses and observations,” Journal of Climate **18**, 5201–5223.

Kim, Y-H, and M-H Kim (2013), “Examination of the global Lorenz energy cycle using MERRA and NCEP-reanalysis 2,” Clim. Dyn. **40**, 1499–1513.

Kimoto, M, and M. Ghil (1993a), “Multiple flow regimes in the Northern Hemisphere winter. Part I: Methodology and hemispheric regimes,” J. Atmos. Sci. **50**, 2625–2643.

Kimoto, M, and M. Ghil (1993b), “Multiple flow regimes in the Northern Hemisphere winter. Part II: Sectorial regimes
and preferred transitions,” J. Atmos. Sci. 50, 2645–2673.

Kistler, R, and Coauthors (2001), “The NCEP-NCAR 50-year reanalysis: Monthly means CD-ROM and documentation,” Bull. Amer. Meteorol. Soc. 82, 247–267.

Kleidon, A (2009), “Non-equilibrium thermodynamics and maximum entropy production,” *Naturwissenschaften* **96**, 653–677.

Kleidon, A (2010), “Life, hierarchy, and the thermodynamic machinery of planet Earth,” Physics Life Rev. **7** (4), 424–460.

Kleidon, A, and R. Lorenz, Eds. (2005), *Non-equilibrium Thermodynamics and the Production of Entropy* (Springer, Berlin).

Klein, R (2010), “Scale-dependent models for atmospheric flows,” Annu. Rev. Fluid Mech. **42** (1), 249–274.

Kloeden, P E, and M. Rasmussen (2011), *Nonautonomous Dynamical Systems* (American Mathematical Society).

Kondrashov, D, M. D. Chekroun, and M. Ghil (2015), “Data-driven non-Markovian closure models,” Physica D **297**, 33–55.

Kondrashov, D, M. D. Chekroun, A. W. Robertson, and M. Ghil (2013), “Low-order stochastic model and ”past-noise forecasting” of the Madden-Julian Oscillation,” Geophys. Res. Lett. **40**, 5305–5310.

Kondrashov, D, Y. Feliks, and M. Ghil (2005a), “Oscillatory modes of extended Nile River records (A.D. 622–1922),” *Geophysical Research Letters* **32**, L10702.

Kondrashov, D, K. Ide, and M. Ghil (2004), “Weather regimes and preferred transition paths in a three-level quasigeostrophic model,” J. Atmos. Sci. **61** (5), 568–587.

Kondrashov, D, S. Kravtsov, and M. Ghil (2006), “Empirical mode reduction in a model of extratropical low-frequency variability,” J. Atmos. Sci. **63** (7), 1859–1877.

Kondrashov, D, S. Kravtsov, A. W. Robertson, and M. Ghil (2005b), “A hierarchy of data-based ENSO models,” J. Climate **18** (21), 4425–4444.

Kondrashov, D, J. Shen, R. Berk, F. DAndrea, and M. Ghil (2007), “Predicting weather regime transitions in northern hemisphere datasets,” Climate dynamics **29** (5), 535–551.

Kramers, HA (1940), “Brownian motion in a field of force and the diffusion model of chemical reactions,” Physica **7** (4), 284–304.

Kraut, S, and U. Feudel (2002), “Multistability, noise, and attractor hopping: The crucial role of chaotic saddles,” Phys. Rev. E **66**, 015207.

Kravtsov, S, P. Berloff, W. K. Dewar, M. Ghil, and J. C. McWilliams (2007), “Dynamical origin of low-frequency variability in a highly nonlinear mid-latitude coupled model,” J. Climate **19**, 6391–6408.

Kravtsov, S, D. Kondrashov, and M. Ghil (2005), “Multi-level regression modeling of nonlinear processes: Derivation and applications to climatic variability,” J. Climate **18** (21), 4404–4424.

Kravtsov, S, D. Kondrashov, and M. Ghil (2009), “Empirical model reduction and the modeling hierarchy in climate dynamics and the geosciences,” in *Stochastic Physics and Climate Modeling*, edited by T. N. Palmer and P. Williams (Cambridge University Press) pp. 35–72.

Kravtsov, S, D. Kondrashov, Igor Kamenkovich, and M. Ghil (2011), “An empirical stochastic model of sea-surface temperatures and surface winds over the Southern Ocean,” Ocean Science **7** (6), 755–770.

Krishnamurti, T N, L. Stefanova, and V. Misra (1979), *Tropical Meteorology, An Introduction* (Springer).

Kubo, R (1957), “Statistical-mechanical theory of irreversible processes. i. general theory and simple applications to magnetic and conduction problems,” Journal of the Physical Society of Japan **12** (6), 570–586.

Kubo, R (1966), “The fluctuation-dissipation theorem,” Reports on Progress in Physics **29** (1), 255–284.

Kuhlbrodt, T, A. Griesel, M. Montoya, A. Levermann, M. Hofmann, and S. Rahmstorf (2007), “On the driving processes of the atlantic meridional overturning circulation,” Reviews of Geophysics **45** (2), 10.1029/2004RG000166.

Kushnir, Y (1987), “Retrograding wintertime low-frequency disturbances over the North Pacific ocean,” J. Atmos. Sci. **44** (19), 2727–2742.

Kushnir, Y (1994), “Interdecadal variations in North Atlantic sea surface temperature and associated atmospheric conditions,” J. Phys. Oceanogr. **7**, 141–157.

Lagerstrom, P A (1988), *Matched Asymptotic Expansions: Ideas and Techniques* (Springer Science & Business Media).

Lai, Y-C, and T. Tél (2011), *Transient Chaos* (Springer, New York).

Laloyaux, P, E. de Boisseson, M. Balmaseda, J.-R. Bidlot, S. Broennimann, R. Buizza, P. Dalhgren, D. Dee, L. Haimberger, H. Hersbach, Y. Kosaka, M. Martin, P. Poli, N. Rayner, E. Rustemeier, and D. Schepers (2018), “CERA-20C: A coupled reanalysis of the twentieth century,” Journal of Advances in Modeling Earth Systems **10** (5), 1172–1195.

Lamb, H H (1972), *Climate: Present, Past and Future* (Methuen, London).

Langen, P L, and V. A. Alexeev (2005), “Estimating 2 x CO2 warming in an aquaplanet GCM using the fluctuation-dissipation theorem,” Geophysical Research Letters **32** (23), L23708.

Lawrence, M G, S. Schäfer, H. Muri, V. Scott, A. Oschlies, N. E. Vaughan, O. Boucher, H. Schmidt, J. Haywood, and J. Scheffran (2018), “Evaluating climate geoengineering proposals in the context of the Paris Agreement temperature goals,” Nature Communications **9** (1), 3734.

Ledrappier, F, and L.-S. Young (1988), “Entropy formula for random transformations,” Prob. Theory Related Fields **80**, 217–240.

Lee, T, T. Awaji, M. A. Balmaseda, E. Greiner, and D. Stammer (2009), “Ocean state estimation for climate research,” Oceanography **22**, 10.5670/oceanog.2009.74.

Legras, B, and M. Ghil (1985), “Persistent anomalies, blocking, and variations in atmospheric predictability,” J. Atmos. Sci. **42**, 433–471.

Leith, C E (1975), “Climate response and fluctuation dissipation,” *Journal of the Atmospheric Sciences*, Journal of the Atmospheric Sciences **32** (10), 2022–2026.

Leith, C E, and R. H. Kraichnan (1972), “Predictability of turbulent flows,” J. Atmos. Sci. **29**, 1041–1058.

Lembo, V, F. Lunkeit, and V. Lucarini (2019), “Thediato (v1.0) – a new diagnostic tool for water, energy and entropy budgets in climate models,” Geoscientific Model Development **12** (8), 3805–3834.

Lenton, T M, H. Held, E. Kriegler, J. W. Hall, W. Lucht, S. Rahmstorf, and H. J. Schellnhuber (2008), “Tipping elements in the Earth’s climate system,” Proc. Natl. Acad. Sci. USA **105**, 1786–1793.

Li, Z-X, K. Ide, H. Le Treut, and M. Ghil (1997), “Atmospheric radiative equilibria in a simple column model,” Climate Dyn. **13**, 429–440.
Lindzen, R S (1986), "Stationary planetary waves, blocking, and interannual variability," in *Advances in Geophysics*, Vol. 29 (Elsevier) pp. 251-273.

Liu, J., J. A. Curry, H. Wang, M. Song, and R. M. Hor-
ton (2012), "Impact of declining arctic sea ice on winter
snowfall," *Proceedings of the National Academy of Sciences*
**109** (11), 4074–4079.

Lorenz, E N (1963), "Deterministic nonperiodic flow," J. At-
mos. Sci. **20**, 130–141.

Lorenz, E N (1967), *The Nature and Theory of the General
Circulation of the Atmosphere* (World Meteorological Or-
ganization, Geneva, Switzerland).

Lorenz, E N (1969a), "The predictability of a flow which pos-
sesses many scales of motion," Tellus **21**, 289-307.

Lorenz, E N (1969b), "Three approaches to atmospheric pre-
dictability," Bull. Am. Meteorol. Soc. **50**, 345-349.

Lorenz, E N (1976), "Nondeterministic theories of climatic
change," Quatern. Res. **6** (4), 495-506.

Lorenz, E N (1979), "Forced and free variations of weather and climate," J. Atmos. Sci **36**, 1367-1376.

Lorenz, E N (1984), "Irregularity: a fundamental property of the atmosphere," Tellus A **36** (2), 98-110.

Lorenz, E N (1996), "Predictability: a problem partly solved,"
in *Seminar on Predictability, 4-8 September 1995* (Euro-
pean Centre for Medium-range Weather Forecasts, Read-
ing, UK) pp. 1-18.

Lorenz, EN (1955), "Available potential energy and the main-
tenance of the general circulation," Tellus **7**, 157-167.

Lu, Z., Z. Fu, L. Hua, N. Yuan, and L. Chen (2018), “Evaluation of ENSO simulations in CMIP5 models: A new perspective based on correlation phase transition in complex networks,” Scientific Reports **8** (1), 14912.

Lucarini, V (2002), "Towards a definition of climate science,"
Intl. J. Environ. Pollution **18** (5), 413–422.

Lucarini, V (2008a), "Response theory for equilibrium and non-equilibrium statistical mechanics: Causality and generalized kramers-kronig relations," Journal of Statistical Physics **131**, 543–558.

Lucarini, V (2008b), "Validation of climate models," in *Encyclopaedia of Global Warming and Climate Change*, edited by G. Philander (SAGE, Thousand Oaks, California) pp. 1053–1057.

Lucarini, V (2009a), "Evidence of dispersion relations for the nonlinear response of the Lorenz 63 system," Journal of Statistical Physics **134**, 381–400.

Lucarini, V (2009b), "Thermodynamic efficiency and en-
tropy production in the climate system," Phys. Rev. E **80**,
021118.

Lucarini, V (2012), “Stochastic perturbations to dynamical
systems: A response theory approach,” Journal of Statisti-
cal Physics **146** (4), 774–786.

Lucarini, V (2013), “Modeling complexity: the case of cli-
mate science,” in *Models, Simulations, and the Reduction of
Complexity*, edited by U Gähde, S Hartmann, and JH Wolf
(De Gruyter) pp. 229–254.

Lucarini, V (2016), "Response operators for Markov processes in a finite state space: Radius of convergence and link to the response theory for Axiom A systems," Journal of Statistical Physics **162** (2), 312–333.

Lucarini, V (2018), “Revising and extending the linear re-
sponse theory for statistical mechanical systems: Evaluat-
ing observables as predictors and predictands,” Journal of
Statistical Physics **173** (6), 1698–1721.

Lucarini, V, R. Blender, C. Herbert, F. Ragone, S. Pascale, and J. Wouters (2014), “Mathematical and physical ideas for climate science,” Rev. Geophys. **52** (4), 809–859.

Lucarini, V, and T. Bódai (2017), “Edge states in the climate system: exploring global instabilities and critical transitions,” *Nonlinearity* **30** (7), R32.

Lucarini, V, and T. Bódai (2019a), “Transitions across
melancholia states in a climate model: Reconciling the de-
terministic and stochastic points of view,” Phys. Rev. Lett.
**122**, 158701.

Lucarini, V, and T. Bódai (2019b), “Global Stability Properties of the Climate: Melancholia States, Invariant Measures, and Phase Transitions,” arXiv e-prints , arXiv:1903.08348arXiv:1903.08348 [physics.ao-ph].

Lucarini, V, S. Calmanti, and V. Artale (2005a), “Destabilization of the thermohaline circulation by transient changes in the hydrological cycle,” Climate Dynamics **24** (2), 253–262.

Lucarini, V, S. Calmanti, and V. Artale (2007), “Experimental mathematics: Dependence of the stability properties of a two-dimensional model of the atlantic ocean circulation on the boundary conditions,” Russian Journal of Mathematical Physics **14** (2), 224–231.

Lucarini, V, and M. Colangeli (2012), “Beyond the linear fluctuation-dissipation theorem: the role of causality,” Journal of Statistical Mechanics: Theory and Experiment **2012** (05), P05013.

Lucarini, V, D. Faranda, A. C. G. M. M. de Freitas, J. M. M. de Freitas, M. Holland, T. Kuna, M. Nicol, M. Todd, and S. Vaienti (2016), *Extremes and Recurrence in Dynamical Systems* (Wiley).

Lucarini, V, K. Fraedrich, and F. Lunkeit (2010a), “Thermodynamic analysis of snowball earth hysteresis experiment: Efficiency, entropy production, and irreversibility,” Q. J. Royal Met. Soc. **136**, 2–11.

Lucarini, V, K. Fraedrich, and F. Lunkeit (2010b), “Thermodynamics of climate change: generalized sensitivities,” Atmospheric Chemistry and Physics **10** (20), 9729–9737.

Lucarini, V, and A. Gritsun (2020), “A new mathematical framework for atmospheric blocking events,” Climate Dynamics **54** (1), 575–598.

Lucarini, V, T. Kuna, D. Faranda, and J. Wouters (2014b), “Towards a general theory of extremes for observables of chaotic dynamical systems,” Journal of Statistical Physics **154**.

Lucarini, V, S. Pascale, V. Boschi, E. Kirk, and N. Iro (2013), “Habitability and multistability in earth-like planets,” Astronomische Nachrichten **334** (6), 576–588.

Lucarini, V, and F Ragone (2011), “Energetics of climate
models: Net energy balance and meridional enthalpy trans-
port,” Rev. Geophys. **49**, RG1001.

Lucarini, V, F. Ragone, and F. Lunkeit (2017), “Predicting climate change using response theory: Global averages and spatial patterns,” Journal of Statistical Physics **166** (3), 1036–1064.

Lucarini, V, J. J. Saarinen, K.-E. Peiponen, and E. M. Vartiainen (2005b), *Kramers-Kronig relations in Optical Materials Research* (Springer, New York).

Lucarini, V, and S. Sarno (2011), “A statistical mechanical approach for the computation of the climatic response to general forcings,” Nonlin. Processes Geophys **18**, 7–28.

Lucarini, Valerio, and Peter H. Stone (2005), “Thermohaline circulation stability: A box model study. part i: Uncoupled model,” Journal of Climate **18** (4), 501–513,
https://doi.org/10.1175/JCLI-3278.1.

Lynch, P (2008), "The ENIAC Forecasts: A Re-creation," Bull. Am. Meteorol. Soc. **89** (1), 45–55.

Madden, RA, and P.R. Julian (1971), "Detection of a 40-50 day oscillation in the zonal wind in the Tropical Pacific," J. Atmos. Sci. **28**, 702–708.

Majda, A J, B. Gershgorin, and Y. Yuan (2009), “Low-frequency climate response and fluctuation–dissipation the-orems: Theory and practice,” *Journal of the Atmospheric Sciences*, *Journal of the Atmospheric Sciences*, **67**(4), 1186–1201.

Majda, A J, and X. Wang (2006), *Nonlinear Dynamics and Statistical Theories for Basic Geophysical Flows* (Cambridge University Press, Cambridge).

Maloney, E D, and D. L. Hartmann (2000), “Modulation of hurricane activity in the Gulf of Mexico by the Madden-Julian oscillation,” *Science* **287**, 2002–2004.

Manabe, Syukuro, and Robert F. Strickler (1964), “Thermal equilibrium of the atmosphere with a convective adjust-ment,” *Journal of the Atmospheric Sciences* **21** (4), 361–385.

Mann, M E, R. S. Bradley, and M. K. Hughes (1998), “Global-scale temperature patterns and climate forcing over the past six centuries,” *Nature* **392**, 779–787.

Mann, M E, R. S. Bradley, and M. K. Hughes (1999), “Northern Hemisphere temperatures during the past mil-lennium: Inferences, uncertainties, and limitations,” *Geo-phys. Res. Lett.* **26** (6), 759–762.

Mann, M E, Z. Zhang, M. K. Hughes, R. S. Bradley, S. K. Miller, and F. Rutherford, S. and Ni (2008), “Proxy-based reconstructions of hemispheric and global surface temper-ature variations over the past two millennia,” *Proc. Natl. Acad. Sci. USA* **105** (36), 13252–13257.

Marangio, L, J. Sedro, S. Galatolo, A. Di Garbo, and M. Ghil (2019), “Arnold maps with noise: Differentiability and non-monotonicity of the rotation number,” *arXiv preprint arXiv:1904.11744*.

Marconi, U Marini Bettolo, A. Puglisi, L. Rondoni, and A. Vulpiani (2008), “Fluctuation-dissipation: Response theory in statistical physics,” *Phys. Rep.* **461**, 111.

Marotzke, J (2000), “Abrupt climate change and thermo-haline circulation: Mechanisms and predictability,” *Proc. Natl. Acad. Sci. USA* **97**, 1347–1350.

Marques, CAF, A Rocha, and J Corte-Real (2010), “Com-parative energetic of ERA-40, JRA-25 and NCEP-R2 re-analysis, in the wave number domai,” *Dyn Atmos Oceans* **50**, 375–399.

Marques, CAF, A Rocha, J Corte-Real, JM Castanheira, J Ferreira, and P Melo-Goncalves (2009), “Global atmo-spheric energetic from NCEP-reanalysis 2 and ECMWF-ERA40 reanalysis,” *Int. J. Climatol.* **29**, 159–174.

Marshall, J, C. Hill, L. Perelman, and A. Ad-croft (1997), “Hydrostatic, quasi-hydrostatic, and non-hydrostatic ocean modelling,” *J. Geophys. Res.* **102**, 5733–5752.

Marshall, J, and F. Molteni (1993), “Toward a dynamical understanding of atmospheric weather regimes,” *J. Atmos. Sci.* **50**, 1993–2014.

Martinson, D G, K. Bryan, M. Ghil, M. M. Hall, T. R. Karl, E. S. Sarachik, S. Sorooshian, and L. D. Talley, Eds. (1995), *Natural Climate Variability on Decade-to-Century Time Scales* (National Academies Press, Washington, DC, U.S.A.).

McIntyre, S, and R. McKittrick (2005), “Hockey sticks, principal components, and spurious significance,” *Geophys. Res. Lett.* **32** (3), n/a–n/a.

McWilliams, J C (2006), *Fundamentals of Geophysical Flus-ids* (Cambridge University Press, Cambridge, UK and New York, NY, USA).

McWilliams, J C (2019), “A perspective on the legacy of Edward Lorenz,” *Earth and Space Sciences* in press.

Meacham, S P (2000), “Low frequency variability of the wind-driven circulation,” *J. Phys. Oceanogr.* **30**, 269–293.

Meehl, G A, L. Goddard, G. Boer, R. Burgman, G. Branstator, C. Cassou, S. Corti, G. Danabasoglu, F. Doblas-Reyes, E. Hawkins, A. Karspeck, M. Kimoto, A. Kumar, D. Matei, J. Mignot, R. Msadek, A. Navarra, H. Pohlmann, M. Riebecker, T. Rosati, E. Schneider, D. Smith, R. Sutton, H. Teng, G. J. van Oldenborgh, G. Vecchi, and S. Yeager (2014), “Decadal climate prediction: An update from the trenches,” *Bulletin of the American Meteorological Society* **95** (2), 243–267.

Melbourne, I, and A. M. Stuart (2011), “A note on diffusion limits of chaotic skew-product flows,” *Nonlinearity* **24** (4), 1361–1367.

Members, PALAEOSENS Project, E. J. Rohling, E. J. Rohling, A. Sluijs, H. A. Dijkstra, P. Köhler, R. S. W. van de Wal, A. S. von der Heydt, D. J. Beerling, A. Berger, P. K. Bijl, M. Crucifix, R. DeConto, S. S. Drijfhout, A. Fedorov, G. L. Foster, A. Ganopolski, J. Hansen, B. Hönisch, H. Hooghiemstra, M. Huber, P. Huybers, R. Knutti, D. W. Lea, L. J. Lourens, D. Lunt, V. Masson-Delmotte, M. Medina-Elizalde, B. Otto-Bliesner, M. Pagani, H. Pälike, H. Renssen, D. L. Royer, M. Siddall, P. Valdes, J. C. Zachos, and R. E. Zeebe (2012), “Making sense of palaeoclimate sensitivity,” *Nature* **491**, 683 EP –.

Merkin, V G, D. Kondrashov, M. Ghil, and B. J. Anderson (2016), “Data assimilation of low-altitude magnetic perturbations into a global magnetosphere model,” *Space Weather* **14**, 165–184.

Merryfield, William J, Johanna Baehr, Lauriane Batté, Emily J. Becker, Amy H. Butler, Caio A. S. Coelho, Gokhan Danabasoglu, Paul A. Dirmeyer, Francisco J. Doblas-Reyes, Daniela I. V. Domeisen, et al. (2020), “Current and emerging developments in subseasonal to decadal prediction,” *Bulletin of the American Meteorological Society* (2020).

Mitchell, J M (1976), “An overview of climate variability and its causal mechanisms,” *Quatern. Res.* **6**, 481–493.

Mo, KC, and M. Ghil (1987), “Statistics and dynamics of persistent anomalies,” *J. Atmos. Sci.* **44** (5), 877–902.

Molteni, F M (2002), “Weather regimes and multiple equilibria,” in *Encyclopedia of Atmospheric Sciences*, edited by J. R. Holton (Academic Press) pp. 2577–2586.

Monge, G (1781), “Mémoire sur la théorie des déblais et des remblais,” *Histoire de l’Académie Royale des Sciences* **1**, 666–704.

Mori, H (1965), “Transport, collective motion, and Brownian motion,” *Progress of Theoretical Physics* **33** (3), 423–455.

Moron, V, R. Vautard, and M. Ghil (1998), “Trends, interdecadal and interannual oscillations in global sea-surface temperature,” *Clim. Dyn.* **14**, 545–569.

Mukhin, D, D. Kondrashov, E. Loskutov, A. Gavrilov, A. Feigin, and M. Ghil (2015b), “Predicting critical transitions in ENSO models. Part II: Spatially dependent models,” *Journal of Climate* **28** (5), 1962–1976, https://doi.org/10.1175/JCLI-D-14-00240.1.