community, as incompatible with existing paleoclimatic evidence at the time.
Geochemical evidence led in the early 1990s to the discovery of a snowball or,
at least, slushball Earth prior to the emergence of multicellular life, 600
Myr B.P. Hoffman  _et al._ (2002); Hoffman and Schrag (2002). This discovery
did not lead, however, to more enthusiasm for the theoretical prediction of
such a state, almost two decades earlier.

The bifurcation diagram of an 1-D EBM, like the one of Eq. (14), is shown in
Fig. 21. It displays the model’s mean temperature $\bar{T}$ as a function of
the fractional change $\mu$ in the insolation $Q_{0}=Q_{0}(x)$ at the top of
the atmosphere. The ‘$S$’-shaped curve in the figure arises from two back-to-
back saddle-node bifurcations. Bensid and Diaz (2019) have recently provided a
mathematically rigorous treatment of the bifurcation diagram of a Budyko-type
1-D EBM.

The normal form of the first saddle-node bifurcation is

$\dot{X}=\mu-X^{2}.$ (15)

Here $X$ stands for a suitably normalized form of $\bar{T}$ and $\dot{X}$ is
the rate of change of $X$, while $\mu$ is a parameter that measures the stress
on the system, in particular a normalized form of the insolation parameter in
Eq. (13b).

The uppermost branch corresponds to the steady-state solution $X=+\mu^{1/2}$
of Eq. (15) and it is stable. This branch matches rather well Earth’s present-
day climate for $\mu=1.0$; more precisely the steady-state solution
$T=T(x;\mu)$ of the full 1-D EBM (not shown) matches closely the annual mean
temperature profile from instrumental data over the last century (Ghil, 1976).

Figure 21: Bifurcation diagram for the solutions of an energy-balance model
(EBM), showing the global-mean temperature $\bar{T}$ vs. the fractional change
$\mu$ of insolation at the top of the atmosphere. The arrows pointing up and
down at about $\mu=1.4$ indicate the stability of the branches: towards a
given branch if it is stable and away if it is unstable. The other arrows show
the hysteresis cycle that global temperatures would have to undergo for
transition from the upper stable branch to the lower one and back. The angle
$\gamma$ gives the measure of the present climate’s sensitivity to changes in
insolation. [After Ghil and Childress (1987) with permission from Springer
Science+Business Media.]

The intermediate branch starts out at the left as the second solution,
$X=-\mu^{1/2}$, of Eq. (15) and it is unstable. It blends smoothly into the
upper branch of a coordinate-shifted and mirror-reflected version of Eq. (15),
say

$\dot{X}=(\mu-\mu_{0})+(X-X_{0})^{2}.$ (16)

This branch, $X=X_{0}+(\mu_{0}-\mu)^{1/2}$, is also unstable. Finally, the
lowermost branch in Fig. 21 is the second steady-state solution of Eq. (16),
$X=X_{0}-(\mu_{0}-\mu)^{1/2}$, and it is stable, like the uppermost branch.
The lowermost branch corresponds to an ice-covered planet at the same distance
from the Sun as Earth.

The fact that the upper-left bifurcation point $(\mu_{\rm c},T_{\rm c})$ in
Fig. 21 is so close to present-day insolation values created great concern in
the climate dynamics community in the mid-1970s, when these results were
obtained. Indeed, much more detailed computations (see below) confirmed that a
reduction of about 2–5% of insolation values would suffice to precipitate
Earth into a “deep freeze.” The great distance of the lower-right bifurcation
point $(\mu_{\rm d},T_{\rm d})$ from present-day insolation values, on the
other hand, suggests that one would have to nearly double atmospheric opacity,
say, for the Earth’s climate to jump back to more comfortable temperatures.

The results above follow Ghil (1976). Held and Suarez (1974) and North (1975)
obtained similar results, and a detailed comparison between EBMs appears in
Chapter 10 of Ghil and Childress (1987). Ghil (1976), rigorously, and then
North  _et al._ (1979), numerically, pointed out that a double-well potential,
like the one sketched in Fig. 20, does exist for higher-dimensional (North
_et al._ , 1979) and even infinite-dimensional (Ghil, 1976) versions of an
EBM. In higher dimensional cases, the maximum of the potential shown in Fig.
20 is replaced by a saddle, or “mountain pass” (e.g., Ghil and Childress,
1987, Sec. 10.4). In this case, dimensionality refers to phase space, rather
than physical space; see a more detailed discussion in Sect. V.

###  III.2 Other Atmospheric Processes and Models

The 1-D atmospheric models in which the details of radiative equilibrium are
investigated with respect to a height coordinate $z$ (geometric height,
pressure, etc.) are often called radiative-convective models (Ramanathan and
Coakley, 1978). This name emphasizes the key role that convection plays in
vertical heat transfer. While these models preceded historically EBMs as rungs
on the modeling hierarchy, it was only recently shown that they, too, could
exhibit multiple equilibria (Li  _et al._ , 1997). The word “equilibrium,”
here and in the rest of this article, refers simply to a steady state of the
model, rather than a to true thermodynamic equilibrium.

Two-dimensional (2-D) atmospheric models are also of two kinds, according to
the third space coordinate that is not explicitly included. Models that
resolve explicitly two horizontal coordinates, on the sphere or on a plane
tangent to it, tend to emphasize the study of the dynamics of large-scale
atmospheric motions. They often have a single layer or two. Those that resolve
explicitly a meridional coordinate and height are essentially combinations of
EBMs and radiative-convective models and emphasize therewith the thermodynamic
state of the system, rather than its dynamics.

Yet another class of “horizontal” 2-D models is the extension of EBMs to
resolve zonal, as well as meridional surface features, in particular land-sea
contrasts. We shall see in Sec. III.3.2 how such a 2-D EBM is used, when
coupled to an oceanic model.

Schneider and Dickinson (1974) and Ghil and Robertson (2000) discuss
additional types of 1-D and 2-D atmospheric models, along with some of their
main applications. Finally, to encompass and resolve the main atmospheric
phenomena with respect to all three spatial coordinates, as discussed in Sect.
II.4, GCMs occupy the pinnacle of the modeling hierarchy (e.g., Randall,
2000).

The mean zonal temperature’s dependence on the insolation parameter $\mu$ — as
obtained for 1-D EBMs and shown in Fig. 21 here — was confirmed, to the extent
possible, by using a simplified GCM, coupled to a “swamp” ocean model
(Wetherald and Manabe, 1975). More precisely, forward integrations with a GCM
cannot confirm the presence of the intermediate, unstable branch. Nor was it
possible in the mid-70s, when this numerical experiment was carried out, to
reach the deep-freeze stable branch, as it was called at the time, because of
the GCM’s computational limitations.

Still, the normal form of the saddle-node bifurcation, given by Eq. (15),
suggests a parabolic shape of the upper, present-day–like branch near the
upper-left bifurcation point in our figure, namely $(\mu_{\rm c},T_{\rm c})$.
This parabolic shape is characteristic of the dependence of a variable that
represents the model solution on a parameter that represents the intensity of
the forcing in several types of bifurcations; moreover, this shape is not
limited to the bifurcation’s normal form, like Eqs. (15) or (16), but is much
more general. The GCM simulations supported quite well a similar shape for the
globally averaged temperature profiles of the GCM’s five vertical layers, cf.
Wetherald and Manabe (1975, Fig. 5). See discussion below in Sects. IV.5.5 and
V.

Ghil and Robertson (2000) also describe the separate hierarchies that have
grown over the last quarter-century in modeling the ocean and the coupled
ocean–atmosphere system. More recently, an overarching hierarchy of Earth
system models — that encompass all the subsystems of interest, atmosphere,
biosphere, cryosphere, hydrosphere and lithosphere — has been developing.
Eventually, the partial results about each subsystem’s variability, outlined
in this section and the next one, will have to be verified from one rung to
the next of the full Earth system modeling hierarchy.

###  III.3 Oscillations in the Oceans’ Thermohaline Circulation

####  III.3.1 Theory and Simple Models

Historically, the thermohaline circulation (THC) (e.g., Dijkstra and Ghil,
2005; Kuhlbrodt  _et al._ , 2007) was first among the climate system’s major
processes to be studied using a very simple mathematical model. Stommel (1961)
formulated a two-box model and showed that it possesses multiple equilibria.

A sketch of the Atlantic Ocean’s THC and its interactions with the atmosphere
and cryosphere on long time scales is shown in Fig. 22. These interactions can
lead to climate oscillations with multi-millennial periods, such as the
Heinrich events (e.g., Ghil, 1994, and references therein), and are summarized
in the figure’s caption. An equally schematic view of the global THC is
provided by the widely known “conveyor belt” (Broecker, 1991) diagram. The
latter diagram captures greater horizontal, 2-D detail but it does not
commonly include the THC’s interactions with water in both its gaseous and
solid phases, which our Fig. 22 here does include.

Figure 22: Schematic diagram of an Atlantic meridional cross section from
North Pole (NP) to South Pole (SP), showing mechanisms likely to affect the
thermohaline circulation (THC) on various time-scales. Changes in the
radiation balance $R_{\rm in}-R_{\rm out}$ are due, at least in part, to
changes in extent of Northern Hemisphere (NH) snow and ice cover $V$, and to
how these changes affect the global temperature $T$; the extent of Southern
Hemisphere ice is assumed constant, to a first approximation. The change in
hydrologic cycle expressed in the terms $P_{\rm rain}-P_{\rm evap}$ for the
ocean and $P_{\rm snow}-P_{\rm abl}$ for the snow and ice is due to changes in
ocean temperature. Deep-water formation in the North Atlantic Subpolar Sea
(North Atlantic Deep Water: NADW) is affected by changes in ice volume and
extent, and regulates the intensity $C$ of the THC; changes in Antarctic
Bottom Water (AABW) formation are neglected in this approximation. The THC
intensity $C$ in turn affects the system’s temperature, and is also affected
by it. [After Ghil  _et al._ (1987) with permission from Springer Nature.]

Basically, the THC is due to denser water sinking, lighter water rising, and
water mass continuity closing the circuit through near-horizontal flow between
the areas of rising and sinking555A complementary point of view suggests,
instead, that surface winds and tides play a major role in the driving and
maintenance of the large-scale ocean circulation Wunsch (2002, 2013). effects
of temperature $T$ and salinity $S$ on the ocean water’s density,
$\rho=\rho(T,S)$, oppose each other: the density $\rho$ decreases as $T$
increases and it increases as $S$ increases. It is these two effects that give
the thermohaline circulation its name, from the Greek words for $T$ and $S$.
In high latitudes, $\rho$ increases as the water loses heat to the air above
and, if sea ice is formed, as the water underneath is enriched in brine. In
low latitudes, $\rho$ increases due to evaporation but decreases due to
sensible heat flux into the ocean.

For the present climate, the temperature effect is stronger than the salinity
effect, and ocean water is observed to sink in certain areas of the high-
latitude North Atlantic and Southern Ocean — with very few and limited areas
of deep-water formation elsewhere — and to rise everywhere else. Thus, in a
thermohaline regime, $T$ is more important than $S$ and hence comes before it.

During certain remote geological times, deep water may have formed in the
global ocean near the equator; such an overturning circulation of opposite
sign to that prevailing today has been dubbed halothermal, $S$ before $T$. The
quantification of the relative effects of $T$ and $S$ on the oceanic water
masses’ buoyancy in high and low latitudes is far from complete, especially
for paleocirculations; the association of the latter with salinity effects
that exceed the thermal ones (e.g., during the Palaeocene, $\simeq 57$ Myr
ago, cf. Kennett and Stott, 1991) is thus rather tentative.

To study the reversal of the abyssal circulation, due to the opposite effects
of $T$ and $S$, Stommel (1961) considered a model with two pipes connecting
two boxes. He showed that the system of two nonlinear, coupled ordinary
differential equations that govern the temperature and salinity differences
between the two well-mixed boxes has two stable steady-state solutions; these
two steady states are distinguished by the direction of flow in the upper and
the lower pipe.

Stommel’s paper was primarily concerned with distinct local convection
regimes, and hence vertical stratifications, in the North Atlantic and the
Mediterranean or the Red Sea, say. Today, we mainly think of one box as
representing the low latitudes and the other one the high latitudes in the
global THC (Marotzke, 2000, and references therein). Subsequently, many other
THC models that used more complex box-and-pipe geometries have been proposed
and studied (e.g., Rooth, 1982; Scott  _et al._ , 1999; Titz  _et al._ , 2002;
Lucarini and Stone, 2005).

The next step in the hierarchical modeling of the THC is that of 2-D
meridional-plane models, in which the temperature and salinity fields are
governed by coupled nonlinear PDEs with two independent space variables, say
latitude and depth (e.g., Quon and Ghil, 1992; Cessi and Young, 1992; Lucarini
_et al._ , 2005a, 2007). Given boundary conditions for such a model that are
symmetric about the Equator, like the equations themselves, one expects a
symmetric solution, in which water either sinks near the poles and rises
everywhere else (thermohaline) or sinks near the Equator and rises everywhere
else (halothermal). These two symmetric solutions would correspond to the two
equilibria of the Stommel (1961) two-box model; see Thual and McWilliams
(1992) for a discussion of the relationship between 2-D and box models.

In fact, symmetry breaking can occur, leading gradually from a symmetric two-
cell circulation to an antisymmetric one-cell circulation. In between, all
degrees of dominance of one cell over the other are possible. A situation
lying somewhere between the two seems to resemble most closely the meridional
overturning diagram of the Atlantic Ocean in Fig. 22.

The gradual transition is illustrated by Fig. 23 and it can be described by a
pitchfork bifurcation, cf. Dijkstra and Ghil (2005):

$\dot{X}=f(X;\mu)=\mu X-X^{3}.$ (17)

Here $X$ measures how asymmetric the solution is, so that $X=0$ is the
symmetric branch, and $\mu$ measures the stress on the system, in particular a
normalized form of the buoyancy flux at the surface. For $\mu<0$ the symmetric
branch is stable, while for $\mu>0$ the two branches $X=\pm\mu^{1/2}$ inherit
its stability.

Figure 23: Transition from a symmetric to an increasingly asymmetric
meridional ocean circulation. The streamfunction plots represent results from
an idealized 2-D model of thermosolutal convection in a rectangular domain,
for a prescribed value of the Rayleigh number and increasing values of the
nondimensional intensity $\gamma$ of the salinity flux at the surface. (a)
$\gamma=0.40$; (b) $\gamma=0.50$: and (c) $\gamma=0.55$. [After Quon and Ghil
(1992) with permission from Cambridge University Press.]

In the 2-D THC problem, the left cell dominates on one of the two branches,
while the right cell dominates on the other: for a given value of $\mu$, the
two stable steady-state solutions — on the $\\{X=+\mu^{1/2}\\}$ branch and on
the $\\{X=-\mu^{1/2}\\}$ branch, respectively — are mirror images of each
other. The idealized THC in Fig. 22, with the North Atlantic Deep Water
extending to the Southern Ocean’s polar front, corresponds to one of these two
branches. In theory, therefore, a mirror-image circulation, with the Antarctic
Bottom Water extending to the North Atlantic’s polar front, is equally
possible.

More recently, Cessi (2019, and references therein) have argued that the
meridional overturning is actually powered by momentum fluxes and not by
buoyancy fluxes. Because of the ongoing arguments on this topic, the THC is
termed increasingly the meridional overturning circulation (MOC).

####  III.3.2 Bifurcation Diagrams for GCMs

Time scale |  Phenomena |  Mechanism
---|---|---
Decadal |  $\cdot$ Local migration of surface anomalies. in the NW corner of the ocean basin |  $\cdot$ Localized surface density anomalies due to surface coupling
|  $\cdot$ Gyre advection in mid-latitudes |  $\cdot$ Gyre advection
Centennial |  Loop-type, meridional circulation |  Conveyor-belt advection of density anomalies
Millennial |  Relaxation oscillation, with “flushes” and superimposed decadal fluctuations |  Bottom water warming, due to strong braking effect of salinity forcing
Table 1: Oscillations in the oceans’ thermohaline circulation. See also Ghil
(1994).

Bryan (1986) was the first to document transition from a two-cell to a one-
cell circulation in a simplified ocean GCM with idealized, symmetric forcing.
In Sect. III.2 atmospheric GCMs confirmed essentially the EBM results. Results
of coupled ocean–atmosphere GCMs, however, have led to questions about whether
the presence of more than one stable THC equilibrium is actually realistic.
The situation with respect to the THC’s pitchfork bifurcation (17) is thus
subtler than it was with respect to Fig. 21 for radiative equilibria.

Dijkstra (2007) showed, however, by comparing observational and reanalysis
data with high-end ocean GCMs, that the current Atlantic Ocean’s THC is
actually in its multiple-equilibria regime. Climate models of intermediate
complexity did indeed find the Atlantic Ocean to be multistable (e.g.,
Rahmstorf  _et al._ , 2005), while several of the GCMs used in recent IPCC ARs
did not show such a behavior, with an exception coming from the study by
Hawkins  _et al._ (2011); see Fig. 24. Thus a failure of high-end models in
the hierarchy to confirm results obtained on the hierarchy’s lower rungs does
not necessarily imply it is the simpler models that are wrong. To the
contrary, such a failure might well indicate that the high-end models, no
matter how detailed, may still be rather imperfect; see also Ghil (2015) for a
summary.

Internal variability of the THC — with smaller and more regular excursions
than the huge and totally irregular jumps associated with bistability — was
studied intensively in the late 1980s and the 1990s. These studies placed
themselves on various rungs of the modeling hierarchy, from box models through
2-D models and all the way to ocean GCMs. A summary of the different kinds of
oscillatory variability found in the latter appears in Table 1666Anomalies in
the table are defined as the difference between the monthly mean value of a
variable and its climatological mean.. Such oscillatory behavior seems to
match more closely the instrumentally recorded THC variability, as well as the
paleoclimatic records for the recent geological past, than bistability.

The (multi)millennial oscillations interact with variability in the surface
features and processes shown in Fig. 22. Chen and Ghil (1996), in particular,
studied some of the interactions between atmospheric processes and the THC.
They used a so-called hybrid coupled model, namely a 2-D EBM, coupled to a
rectangular-box version of the North Atlantic rendered by a low-resolution
ocean GCM. This hybrid model’s regime diagram is shown in Fig. 25(a). A steady
state is stable for higher values of the coupling parameter $\lambda_{\rm ao}$
or of the EBM’s diffusion parameter $d$. Interdecadal oscillations with a
period of 40–50 years are self-sustained and stable for lower values of these
parameters.

(a) Intermediate-complexity models

(b) IPCC-class climate model

Figure 24: Multistability of the THC for (a) Earth models of intermediate
complexity, reproduced with permission from Rahmstorf  _et al._ (2005); and
(b) for an IPCC-class climate model, reproduced with permission from Hawkins
_et al._ (2011). Parameter controlling the freshwater input in the North
Atlantic basin on the abscissa and the THC’s intensity on the ordinate.

The self-sustained THC oscillations in question are characterized by a pair of
vortices of opposite sign that grow and decay in quadrature with each other in
the ocean’s upper layers. Their centers follow each other anti-clockwise
through the northwestern quadrant of the model’s rectangular domain. Both the
period and the spatio-temporal characteristics of the oscillation are thus
rather similar to those seen in a fully coupled GCM with realistic geometry
(Delworth  _et al._ , 1993). The transition from a stable equilibrium to a
stable limit cycle, via Hopf bifurcation, in this hybrid coupled model, is
shown in Fig. 25(b), and we will further elaborate upon it in Sect. III.4.2
below.

(a) Regime diagram (b) Bifurcation diagram

Figure 25: Dependence of THC solutions on two parameters in a hybrid coupled
model; the two parameters are the atmosphere–ocean coupling coefficient
$\lambda_{\rm ao}$ and the atmospheric thermal diffusion coefficient $d$. (a)
Schematic regime diagram. The full circles stand for the model’s stable steady
states, the open circles for stable limit cycles, and the solid curve is the
estimated neutral stability curve between the former and the latter. (b) Hopf
bifurcation curve at fixed $d$ = 1.0 and varying $\lambda_{\rm ao}$; this
curve was obtained by fitting a parabola to the model’s numerical-simulation
results, shown as full and open circles. [From Chen and Ghil (1996) with
permission from the American Meteorological Society.]

###  III.4 Bistability, Oscillations and Bifurcations

In Sects. III.1, III.2 and III.3.1, we have introduced bistability of steady-
state solutions via saddle-node and pitchfork bifurcations, while in Sect.
III.3.2, we mentioned oscillatory solutions and Hopf bifurcation as the
typical way the latter are reached as a model parameter changes. We start this
subsection by summarizing the bifurcations of steady-state solutions, often
referred to as equilibria, and will then introduce and discuss the normal
forms of Hopf bifurcation, both supercritical and subcritical.

####  III.4.1 Bistability and Steady-State Bifurcations

Section III.2 introduced EBMs and explained how the present climate and a
totally ice-covered Earth can result as coexisting stable steady states over a
certain parameter range. The normal forms of a supercritical and a subcritical
saddle-node bifurcation were given as Eqs. (15) and (16), respectively. Here
the criticality is defined as the existence of the two equilibria, stable an
unstable, to the right or the left of the critical or bifurcation point.

An $S$-shaped bifurcation curve, such as the one that appears in Fig. 21 for
the 1-D EBM of Eq. (14), can be obtained easily by soldering smoothly together
the back-to-back saddle-node bifurcations of Eqs. (15) and (16). While there
is no generic normal form for such a curve, here is a simple example:

$\displaystyle\dot{X}$ $\displaystyle=$ $\displaystyle\mu-X^{2},$ (18a)
$\displaystyle\dot{X}$ $\displaystyle=$
$\displaystyle(\mu-\mu_{0})+(X-X_{0})^{2},$ (18b)

with $\mu_{0}=1$ and $X_{0}=-1/2$.

Note that both the sub- and supercritical saddle-node bifurcations are
structurally stable, i.e. they persist as the system of evolution equations,
be it a system of ordinary or partial differential equations, is smoothly
perturbed (Andronov and Pontryagin, 1937; Guckenheimer and Holmes, 1983;
Arnold, 2003; Temam, 1997) by a small amount. This robustness is the reason
why saddle-node bifurcations — and other elementary bifurcations called of co-
dimension one, i.e., depending on a single parameter, like the Hopf
bifurcation — are so important and can, in practice, be a guide through the
hierarchy of models, in the climate sciences and elsewhere. A striking example
was provided in Sect. III.2 for the supercritical saddle-node bifurcation that
can be found in the 3-D GCM of Wetherald and Manabe (1975, Fig. 5), as well as
the 1-D EBM reproduced in Fig. 21 herein.

The next kind of bifurcation that leads to bistability of stationary states is
the pitchfork bifurcation introduced in Sect. III.3.1, in connection with the
THC, and whose normal form is given by Eq. (17). This bifurcation, however,
only arises in the presence of a mirror symmetry in the model under study,
such as seen in Fig. 23a. One suspects, therefore, that it is nor structurally
stable with respect to perturbations of the dynamics that do not preserve this
symmetry.

A simple example is given by the following perturbed pitchfork bifurcation:

$\dot{X}=f(X;\mu,\epsilon)=X(\mu-X^{2})+\epsilon,$ (19)

where $\epsilon$ is a small parameter. The bifurcation diagrams for
$\epsilon=0$ and $\epsilon=+0.1$ are given in Figs. 26(a,b), respectively.
Clearly, a nonzero value of $\epsilon$ breaks the $X\rightarrow-X$ symmetry of
Eq. (17), since it is no longer the case that $f(X)=f(-X)$. Hence, the
pitchfork breaks up into a continuous uppermost branch that is stable for all
$X$-values, and a supercritical saddle-node bifurcation, whose lower branch is
stable. If $\epsilon<0$, it will be the lowermost branch that is stable for
all $X$-values (not shown), and the upper branch of the saddle-node
bifurcation above it that is stable. In both cases, above a critical value
associated with the saddle-node bifurcation, three solution branches exist,
with the middle one that is unstable and the other two that are stable.

An obvious way that the symmetric pitchfork bifurcation present in 2-D models
of the THC can be broken, as illustrated in Fig. 26, is simply by considering
3-D models with zonally asymmetric basin boundaries. Perturbed pitchfork
bifurcations were also found in shallow-water models of the wind-driven ocean
circulation (e.g., Jiang  _et al._ , 1995b; Speich  _et al._ , 1995b; Ghil,
2017).

####  III.4.2 Oscillatory Instabilities and Hopf Bifurcations

In spite of the considerable detail and 3-D character of the ocean GCM
involved in the hybrid coupled model of Chen and Ghil (1996), it is clear that
the numerically obtained bifurcation diagram in Fig. 25b is of a very simple,
fundamental type. The normal form of such a Hopf bifurcation is given by

$\dot{z}=(\mu+i\omega)z+c(z\bar{z})z,$ (20)

where $z=x+iy$ is a complex variable, while $c$ and $\omega$ are nonnegative
and $\mu$ is a real parameter. Clearly, for small $z$, this equation describes
a rotation in the complex plane with an increasing radius $|z|$ when $\mu>0$,
i.e., it contains the possibility of an oscillatory instability, while the
cubic term corresponds to an increasing modification of this simple rotation
away from the origin.

Figure 26: Bifurcation diagram of the pitchfork bifurcation in Eq. (19), (a)
for $\epsilon=0$; and (b) for $\epsilon=0.1.$ Solid lines indicate stable
solutions, and dashed lines indicate unstable ones. [After Dijkstra and Ghil
(2005) with permission from Elsevier.]

The parsimonious complex notation above follows Arnold (1983) and Ghil and
Childress (1987) and is, we believe, more suggestive and transparent than the
more common one that uses separately the two real variables $x$ and $y$ (e.g.,
Guckenheimer and Holmes, 1983). The advantages of the former are apparent when
introducing polar variables via $z=\rho^{1/2}e^{i\theta},$ with
$\rho=z\bar{z}\geq 0.$ One can then separate the flow induced by Eq. (20) into
a constant rotation with angular velocity $\omega$ and a radius $r=\rho^{1/2}$
that either increases or decreases as $\dot{\rho}\gtrless 0$, according to

$\displaystyle\dot{\rho}$ $\displaystyle=$ $\displaystyle 2\rho(\mu+c\rho),$
(21a) $\displaystyle\dot{\theta}$ $\displaystyle=$ $\displaystyle\omega.$
(21b)

Equation (21a) is quadratic in $\rho$ and has the two roots $\rho=0$ and
$\rho=-\mu/c.$ The former corresponds to a fixed point at the origin $z=0$ in
Eq. (20), while the latter only exists when $c\mu<0$.

(a) supercritical Hopf (b) subcritical Hopf

Figure 27: Bifurcation diagram of the Hopf bifurcation in Eq. (20): (a)
supercritical Hopf bifurcation for $c=-1$; and (b) subcritical Hopf
bifurcation for $c=+1.$ Solid lines indicate stable solutions (s), and dashed
lines indicate unstable ones (u); it is common to plot $|z|=-\mu^{1/2}$ along
with $|z|=\mu^{1/2},$ in order to emphasize that one is actually projecting
onto the $(|z|,\mu)$-plane a paraboloid-shaped, one-parameter family of limit
cycles with given $r=|z|=\mu^{1/2}$ and $0\leq\theta<2\pi.$ [After Ghil and
Childress (1987) with permission from Springer Science+Business Media.]

We thus anticipate a solution that is a circle with radius $r=(-\mu/c)^{1/2}$
when $c\mu\neq 0$ and the two parameters have opposite signs. The simpler case
is that of the nonzero radii given by Eq. (21a) for positive stability
parameter $\mu$ and negative saturation parameter $c$: it corresponds to the
supercritical Hopf bifurcation sketched in Fig. 27a. In this case, the stable
fixed point at the origin loses its stability as $\mu$ changes sign, and
transmits it to a limit cycle of parabolically increasing radius.

The opposite case, of $c>0$, is plotted in Fig. 27b: it corresponds to
subcritical Hopf bifurcation. Here, in fact, the stability of the fixed point
is also lost at $\mu=0$ but there is no oscillatory solution for $\mu>0$ at
all, since the higher-order terms in either Eq. (20) or (21a) cannot stabilize
the linearly unstable oscillatory solution, which spirals out to infinity. It
is not just that the parabolas in Figs. 27a and 27b point in opposite
directions, as for the saddle-node bifurcations in Eqs. (18a) and (18b) or
Eqs. (15) and (16), but the two figures are topologically distinct.

###  III.5 Main Modes of Variability

The atmosphere, ocean, and the coupled ocean–atmosphere climate system have
many modes of variability, as initially discussed in Sects. II.1-II.2-II.3. We
review here some of the most important ones, and comment on the general
features of such modes.

####  III.5.1 Modes of Variability and Extended-Range Prediction

Several large-scale spatial patterns of atmospheric covariability have been
studied, starting in the second half of the 19th century. Lorenz (1967) and
Wallace and Gutzler (1981) provided good reviews of the earlier studies. The
earliest work tended to emphasize “centers of action,” where the variability
is strongest Teisserenc de Bort (1881), while more recently, it is so-called
teleconnections between such centers of action that have been emphasized. J.
Namias (1910–1997) played a key role in developing the interest in such
teleconnections, by applying systematically the use of their spatio-temporal
properties to the development of operational extended-range weather
forecasting (e.g., Namias, 1968).

Returning to the discussion of prediction in Sect. II.4, it is important to
recall John von Neumann’s (1903–1957) important distinction (Von Neumann,
1960) between weather and climate prediction. To wit, short-term NWP is the
easiest — i.e., it is a pure initial-value problem; long-term climate
prediction is next easiest — it corresponds to studying the system’s
asymptotic behavior; while intermediate-term prediction is hardest — both
initial and boundary values are important. Von Neumann’s role in solving the
NWP problem by integrating the discretized equations that govern large-scale
atmospheric motions (Bjerknes, 1904; Richardson, 1922; Charney  _et al._ ,
1950b) is well known. In fact, he also played a key role in convening the
conference on the “Dynamics of Climate” that was held at Princeton’s Institute
for Advanced Studies in October 1955, and whose proceedings Pfeffer (1960)
were finally published three years after Von Neumann’s untimely death.

Today, routine NWP is quite skillful for several days, but we also know that
detailed prediction of the weather is limited in theory by the exponential
growth of small errors (Lorenz, 1963) and by their turbulent cascading from
small to large scales (Thompson, 1957; Lorenz, 1969a, b; Leith and Kraichnan,
1972). The theoretical limit of detailed prediction — in the sense of
predicting future values of temperature, wind and precipitation at a certain
point, or within a small volume, in time and space — is of the order of 10–15
days (Epstein, 1988; Tribbia and Anthes, 1987).

In the sense of Von Neumann (1960), short-term prediction is being improved by
meteorologists and oceanographers through more-and-more accurate
discretization of the governing equations, increased horizontal and vertical
resolution of the numerical models, improved observations and data
assimilation methodologies, and improved parametrization of subgrid-scale
processes.

Important strides in solving the theoretical problems of the climate system’s
asymptotic behavior are being taken by the use of idealized models — either by
simplifying the governing equations in terms of the number of subsystems and
of active physical processes, e.g., by eliminating phase transitions or
chemical processes in the atmosphere — or by systematic model reduction to
small or intermediate-size sets of ordinary differential or stochastic
differential equations (e.g., Chang  _et al._ , 2015; Palmer and Williams,
2009). We shall return to the latter in Sects. IV and V.

What, then, if anything, and how, and how accurately can climate-related
scalars or fields be predicted beyond the limits of NWP? In other words, what
can be done about the gap between short-term and asymptotic prediction of
climate? These issues have been actively pursued for the last three decades
(e.g., Epstein, 1988; Tribbia and Anthes, 1987; Ghil and Robertson, 2002, and
many others). Atmospheric, oceanic and coupled modes of variability play an
important role in extended- and long-range forecasting.

The key idea is that a mode that is stationary or oscillatory can, by its
persistence or periodicity, contribute a fraction of variance that is
predictable for the mode’s characteristic time or, at least, for a substantial
fraction thereof. We start therefore with a brief review of some of the most
promising modes of variability.

####  III.5.2 Coupled Atmosphere–Ocean Modes of Variability

The best known of all these modes is the El Niño–Southern Oscillation (ENSO:
Philander, 1990; Dijkstra, 2005, and references therein), mentioned earlier in
Sect. II.2. The ENSO phenomenon is particularly strong over the Tropical
Pacific, but it affects temperatures and precipitations far away and over a
large area of the globe. Some of the best documented and statistically most
significant ones of these teleconnections are illustrated in Fig. 28. Just as
an example, destructive droughts over Northeast Brazil, Southeast Africa or
cold spells over Florida are often associated with particularly strong warm
ENSO episodes.

Figure 28: Teleconnection pattern for a warm ENSO episode, during the boreal
winter months December-January-February. Both colors and labels indicate warm-
vs.-cold and wet-vs.-dry anomalies, with anomalies being defined as the
difference between a monthly mean value of a variable and its climatological
mean. [Reproduced, with permission, from
%%␣here␣and␣elsewhere␣we␣use␣the␣climatological␣meaning␣of␣the␣word,%%␣instantaneous␣or,␣according␣to␣context,https://www.meted.ucar.edu/ams/wim_2014/9b.html.]

Such strong episodes recurred every 2–7 years during the instrumental record
of roughly 150 years, and the ENSO phenomenon, with it alternation of warm
episodes (El Niños) and cold ones (La Niñas), is quite irregular. Still, there
is a marked tendency for year-to-year alternation of, not necessarily strong,
El Niños and La Niñas; this alternation is associated with a well-known quasi-
biennial near-periodicity (e.g., Rasmusson  _et al._ , 1990; Ghil  _et al._ ,
2002).

Even larger variance, accompanied by lesser regularity, is associated with a
quasi-quadrennial mode Jiang  _et al._ (1995a); Ghil  _et al._ (2002),
sometimes just called the low-frequency ENSO mode. The positive interference
of these two modes generates large ENSO events that visually coincide with the
instrumentally recorded ones (Jiang  _et al._ , 1995a, Fig. 9).

These ENSO features have been used since 1992 for real-time forecasting that
essentially relies on predicting the oscillatory modes of two scalar indices
that capture much of the ENSO variability, namely the Southern Oscillation
Index (SOI) and the so-called Niño-3 index, obtained by averaging mean-monthly
sea surface temperatures (SSTs) over an area of the Eastern Tropical Pacific.
For the time being, such a data-driven forecast appears to still be quite
competitive with those made by high-end, detailed GCMs Barnston  _et al._
(2012).

####  III.5.3 Atmospheric Low-Frequency Variability

As discussed in Sect. II.1, datasets for the atmosphere are both longer and
more plentiful than for the oceans. Thus, there are quite a few modes of
variability that have been detected and described, especially for the Northern
Hemisphere, where both the human population and the data sets are denser.

As discussed in Sect. II.2, in connection with Fig. 12, atmospheric phenomena
are designated as having low, or intraseasonal, frequency if their
characteristic time is longer than the life cycle of a mid-latitude weather
system but still shorter than a season. We are thus talking here of intrinsic
variability, as opposed to the externally forced seasonal cycle, the latter
being easier to understand. More recently, one is also referring to this
variability as _subseasonal_ (e.g., Robertson and Vitart, 2018).

There are essentially two approaches for describing this intrinsic,
subseasonal variability: (i) as being episodic or intermittent, and (ii) as
being oscillatory. The two approaches are complementary, as we shall see, and
they have been dubbed the ‘particle’ vs. the ‘wave’ description, by a crude
analogy with quantum mechanics Ghil and Robertson (2002).

The key ingredient of the particle approach is provided by so-called
_persistent anomalies_ or regimes. The best known among these are blocking vs.
zonal flow, which were already mentioned in Sect. II.2. Other well-known
persistent features, especially during boreal winter, are the positive and
negative phase of the North Atlantic Oscillation (NAO) and the Pacific North
American (PNA) pattern. A rich literature exists on the reliable
identification, description and modeling of these patterns; see Ghil  _et al._
(2018, and references therein). In the Southern Hemisphere, there has also
been some interest in an approximate counterpart of the PNA, called the
Pacific South American (PSA) pattern Mo and Ghil (1987).

The key ingredient of the wave approach is provided by oscillatory modes that
do not necessarily possess exact periodicities, but rather the broad spectral
peaks that were illustrated in Fig. 12 and discussed already to some extent in
Sect. II.2. Probably the best known example of this type in the subseasonal
band is the Madden-Julian Oscillation (MJO: Madden and Julian, 1971). It has a
near-periodicity of roughly 50 days and affects winds and precipitation in the
Tropics, being strongest in the Indo-Pacific sector. Like the much lower-
frequency ENSO phenomenon, important extratropical effects have been
documented (e.g., Maloney and Hartmann, 2000).

In spite of the considerable amount of observational, theoretical and modeling
work dedicated to the MJO, it is still incompletely understood and not very
well simulated or predicted. Some of the reasons for these difficulties
include the key role played in its mechanism by transitions among liquid and
gaseous phases of water in tropical clouds, its multiscale character, and the
pronounced interactions with the oceans Zhang (2005).

Two extratropical modes in the subseasonal band are the Branstator
(1987)–Kushnir (1987) wave, and the 40-day mode associated with the
topographic instability first described by Charney and DeVore (1979) in a low-
order model. The Branstator-Kushnir wave is, like the MJO, an eastward
traveling wave, while the 40-day mode is a standing wave anchored by the
topography.

Charney and DeVore (1979) emphasized the bimodality of the solutions of a
model version with only three Fourier modes — one stable steady state being
zonal and the other being blocked — and did not pursue further the fact that,
in a more highly resolved version of the model, oscillatory solutions did
appear. M. Ghil and associates did clarify the role of the higher meridional
modes in the Hopf bifurcation that gives rise to the oscillatory solutions
Ghil and Childress (1987); Legras and Ghil (1985); Jin and Ghil (1990). The
potential role of these results in the controversy surrounding the effect of
anthropogenic polar amplification on blocking frequency was mentioned in Sect.
II.2.2.

While clearly mid-latitude weather systems, like the tropical ones, also
involve precipitation, their large-scale properties seem to be much less
affected by wet processes; the former arise essentially from purely dynamical,
as opposed to largely thermodynamical mechanisms. For the sake of simplicity,
we will thus try to illustrate in further detail the complementarity of the
wave and the particle approach for the extratropical topographic oscillation.

Given the recent interest in the physical literature for synchronization in
continuous media, cf. Duane  _et al._ (2017) and references therein, it might
be challenging to this readership that one can accomodate in a fairly narrow
frequency band, between roughly $0.1$ and $0.01$ day-1, three distinct
oscillatory modes that do not seem to synchronize: in the Tropics the MJO,
with a period $\simeq 50$ days, and in mid-latitudes the Branstator-Kushnir
wave, with a period $\simeq 30$ days, along with the topographic oscillatory
mode, with a period $\simeq 40$ days; see, for instance, Dickey  _et al._
(1991, Fig. 11). Not only are these three frequencies fairly close. i.e., the
detuning fairly small, but the characteristic wave lengths of all three of
these oscillatory modes are quite large with respect to the radius of the
Earth, and several teleconnections mentioned so far extend across continents
and oceans, cf. Fig. 28 herein.

Be that as it may, let us now concentrate on the topographic ‘wave’ mode and
its relationship with the blocking and zonal ‘particles’. Ghil  _et al._
(2018) recently reviewed the evidence provided by a hierarchy of models for
the presence of a Hopf bifurcation that arises from the interaction of the
large-scale westerly flow in mid-latitudes with the topography of the Northern
Hemisphere.

Certain spatial features of the phases of this mode do present striking
similarities to the blocked and zonal flows that appear not just as the two
stable equilibria in the highly idealized model of Charney and DeVore (1979),
but also as unstable, though long-lived, patterns in much more realistic
models, like the three-level, quasi-geostrophic (QG3) model originally
formulated by Marshall and Molteni (1993). The latter model is still widely
used to study the nonlinear dynamics of large-scale, mid-latitude atmospheric
flows (e.g., Kondrashov  _et al._ , 2004; Lucarini and Gritsun, 2020).

More broadly, we summarize here the relevant conclusions of the Ghil  _et al._
(2018) review paper on the observational, theoretical and modeling literature
of multiple regimes and oscillatory modes of subseasonal variability. Four
complementary approaches to explaining this variability are illustrated in
Fig. 29.

Figure 29: Schematic overview of atmospheric low-frequency variability (LFV)
mechanisms; see text for details. Reproduced from Ghil  _et al._ (2018), with
permission from Elsevier.

One approach to persistent anomalies in mid-latitude atmospheric flows on
subseasonal time scales is to consider them simply as due to slowing down of
Rossby waves or to their linear interference (Lindzen, 1986). This approach is
illustrated in the sketch labeled (c) within the figure: zonal flow $Z$ and
blocked flow $B$ are simply slow phases of an harmonic oscillation, like the
neighborhood of $t=\pi/2$ or $t=3\pi/2$ for a sine wave $\sin(t)$; or else
they are due to an interference of two or more linear waves, like the one
occurring for a sum $A\sin(t)+B\sin(3t)$ near $t=(2k+1)\pi/2$. A more
ambitious, quasi-linear version of this approach is to study long-lived
resonant wave triads between a topographic Rossby wave and two free Rossby
waves (Egger, 1978; Trevisan and Buzzi, 1980; Ghil and Childress, 1987).
Neither version of this line of thought, though, explains the organization of
the persistent anomalies into distinct flow regimes.

Rossby  _et al._ (1939) initiated a different, genuinely nonlinear approach by
suggesting that multiple equilibria may explain preferred atmospheric flow
patterns. These authors drew an analogy between such equilibria and hydraulic
jumps, and formulated simple models in which similar transitions between
faster and slower atmospheric flows could occur. This multiple-equilibria
approach was then pursued vigorously in the 1980s (Charney and DeVore, 1979;
Charney  _et al._ , 1981; Legras and Ghil, 1985; Ghil and Childress, 1987) and
it is illustrated in Fig. 29 by the sketch labeled (a): one version of the
sketch illustrates models that concentrated on the $B$–$Z$ dichotomy (Charney
and DeVore, 1979; Charney  _et al._ , 1981; Benzi  _et al._ , 1986), the other
on models (e.g., Legras and Ghil, 1985) that allowed for the presence of
additional clusters, found by Kimoto and Ghil (1993a) and Smyth  _et al._
(1999), among others, in observations. The latter include opposite phases of
the NAO and PNA anomalies ($PNA,RNA$ and $BNAO$ in sketch (a) of Fig. 29). The
LFV dynamics in this approach is given by the preferred transition paths
between the two or more regimes; see again Table 1 in Ghil  _et al._ (2018)
and references therein.

A third approach is associated with the idea of oscillatory instabilities of
one or more of the multiple fixed points that can play the role of regime
centroids. Thus, Legras and Ghil (1985) found a 40-day oscillation arising by
Hopf bifurcation off their blocked regime $B$, as illustrated in sketch (b) of
the figure.

An ambiguity arises, though, between this point of view and the complementary
possibility that the regimes are just slow phases of such an oscillation,
caused itself by the interaction of the mid-latitude jet with topography that
gives rise to a supercritical Hopf bifurcation. Thus, Kimoto and Ghil (1993a,
b) found, in their observational data, closed paths within a Markov chain
whose states resemble well-known phases of an intraseasonal oscillation. Such
a possibility was confirmed in the QG3 model by Kondrashov  _et al._ (2004).
Furthermore, multiple regimes and intraseasonal oscillations can coexist in a
two-layer model on the sphere within the scenario of “chaotic itinerancy”
(Itoh and Kimoto, 1996, 1997).

Lucarini and Gritsun (2020) observed that blockings occur when the system’s
trajectory is in the neighborhood of a specific class of unstable periodic
orbits (UPOs). UPOs in general are natural modes of variability that cover a
chaotic system’s attractor Cvitanović (1988); Cvitanovíc and Eckhardt (1991).
Here, the UPOs that correspond to blockings have a higher degree of
instability compared to UPOs associated with zonal flow; thus blockings are
associated with anomalously unstable atmospheric states, as suggested
theoretically by Legras and Ghil (1985) and confirmed experimentally by Weeks
_et al._ (1997). Different regimes may be associated with different bundles of
UPOs, a conjecture that could also explain the efficacy of Markov chains fin
describing the transitions between qualitatively different regimes.

Sketch (d) in the figure refers to the role of stochastic processes in S2S
variability and prediction, whether it be noise that is white in time — as in
Hasselmann (1976) or in linear inverse models (LIMs: Penland, 1989, 1996;
Penland and Ghil, 1993; Penland and Sardeshmukh, 1995) — or red in time, as in
certain nonlinear data-driven models (Kravtsov  _et al._ , 2005, 2009;
Kondrashov  _et al._ , 2006, 2013, 2015) or even non-Gaussian (Sardeshmukh and
Penland, 2015). Stochastic processes may enter into models situated on various
rungs of the modeling hierarchy, from the simplest conceptual models to high-
resolution GCMs. In the former, they may enter via stochastic forcing, whether
additive or multiplicative, Gaussian or not (e.g., Kondrashov  _et al._ ,
2015), while in the latter, they may enter via stochastic parametrizations of
subgrid-scale processes (e.g., Palmer and Williams, 2009, and references
therein).

Figure 29 summarizes some of the key dynamical mechanisms of mid-latitude
subseasonal variability, as discussed in this section and in Ghil  _et al._
(2018), without attempting to provide a definitive answer as to which approach
to modeling and prediction of this variability is or will be most productive
in the near future.

###  III.6 Internal Variability and Routes to Chaos

In the present section, we illustrate on hand of the wind-driven ocean
circulation a sequence of successive bifurcations that lead from a highly
symmetric, steady-state circulation to much more finely structured, irregular,
possibly chaotic oscillations. Various issues arise in pursuing such a
bifurcation sequence across a hierarchy of models and on to the observational
data.

Figure 30: A map of the main oceanic currents: warm currents in red and cold
ones in blue. Reproduced from Ghil  _et al._ (2008), with permission from
Elsevier.

Mid-latitude oceanic gyres appear clearly in Fig. 30 below, in the four major
extratropical ocean basins, namely the North and South Atlantic, and the North
and South Pacific. The large, subtropical ocean gyres are formed by a
poleward-flowing western boundary current, an equatorward-flowing eastern
boundary current, and the roughly zonally flowing currents that connect these
two coastal currents off the equator and on the poleward basin side,
respectively. These gyres are characterized by so-called anticyclonic
rotation, clockwise in the Northern and anticlockwise in the Southern
Hemisphere. In the North Pacific and the North Atlantic, they are accompanied
by smaller, cyclonically rotating gyres, while in the Southern Hemisphere,
such subpolar gyres are missing and replaced by the Antarctic Circumpolar
Current.

The basic phenomenology of these gyres and the detailed physical mechanisms
that give rise to it are described in several books (e.g., Cushman-Roisin and
Beckers, 2011b; Ghil and Childress, 1987; Gill, 1982; Pedlosky, 1996; Sverdrup
_et al._ , 1946; Vallis, 2006) and review papers Dijkstra and Ghil (2005);
Ghil (2017). Clearly the sharp western boundary currents — like the Gulf
Stream in the North Atlantic. the Kuroshio and its cross-basin extension in
the North Pacific, and the Brazil Current in the South Atlantic — as well as
the more diffuse eastern boundary currents — like the Canaries Current in the
North Atlantic and the California and Peru Currents in the North and South
Pacific, respectively — play a major role in carrying heat poleward and colder
waters equatorward. Hence these gyres’ interannual and interdecadal
variability is a major contributor to climate variability.

To illustrate the bifurcation sequence that might lead to this oceanic LFV, we
use a highly idealized model of the wind-driven double-gyre circulation in a
rectangular geometry. Note that the counterparts of synoptic weather systems
in the ocean are eddies and meanders that have much shorter spatial scales
than in the atmosphere, but considerably longer time scales:
$\mathcal{O}(100)$ km in the ocean vs. $\mathcal{O}(1000)$ km in the
atmosphere, but order of several months in time vs. order of merely several
days. Thus the definition we used for LFV in the atmosphere, cf. Sect.
III.5.3, does correspond in the ocean to a time scale of years to decades.

####  III.6.1 A simple model of the double-gyre circulation

The simplest model that includes several of the most pertinent mechanisms
described in Sect. II.2 is governed by the barotropic quasi-geostrophic
equations. We consider an idealized, rectangular basin geometry and simplified
forcing that mimics the distribution of vorticity due to the wind stress, as
sketched in Simonnet  _et al._ (2005, Fig. 2). In this idealized model, the
amounts of subpolar and subtropical vorticity injected into the basin are
equal and the rectangular domain $\Omega=(0,L_{x})\times(0,L_{y})$ is
symmetric about the axis of zero wind stress curl $y=L_{y}/2$.

The barotropic $2$-D quasi-geostrophic equations in this idealized setting
are:

$\displaystyle\partial_{t}q$ $\displaystyle+$ $\displaystyle
J(\psi,q)-\nu\Delta^{2}\psi+\mu\Delta\psi=-\tau\sin\big{(}\frac{2\pi
y}{L_{y}}\big{)},$ (22a) $\displaystyle q$ $\displaystyle=$
$\displaystyle\Delta\psi-\lambda_{\rm R}^{-2}\psi+\beta y.$ (22b)

Here $x$ points east and $y$ points north, while $q$ and $\psi$ are the
potential vorticity and the streamfunction, respectively, and the Jacobian $J$
gives the advection of potential vorticity by the flow, as already discussed
in Sect. II.3.3 , so that $J(\psi,q)=\psi_{x}q_{y}-\psi_{y}q_{x}={\bf
u}\cdot\nabla q$.

The physical parameters are the strength of the planetary vorticity gradient
$\beta=\partial f/\partial y$, the Rossby radius of deformation
$\lambda_{\mathrm{R}}^{-2}$, the eddy-viscosity coefficient $\nu$, the bottom
friction coefficient $\mu$, and the wind-stress intensity $\tau$. One
considers here free-slip boundary conditions $\psi=\Delta^{2}\psi=0$; the
qualitative results described below do not depend on the choice of homogeneous
boundary conditions Jiang  _et al._ (1995b); Dijkstra and Ghil (2005).

The nonlinear system of PDEs (22) is an infinite-dimensional dynamical system
and one can thus study its bifurcations as the parameters change. Two key
parameters are the wind stress intensity $\tau$ and the eddy viscosity $\nu$:
as $\tau$ increases the solutions become rougher, while an increase in $\nu$
renders them smoother.

An important property of (22) is its mirror symmetry in the $y=L_{y}/2$ axis.
This symmetry can be expressed as invariance with respect to the discrete
$\mathbb{Z}_{2}$ group ${\mathcal{S}}$,

${\mathcal{S}}\left[\psi(x,y)\right]=-\psi(x,L_{y}-y);$

any solution of (22) is thus accompanied by its mirror-conjugated solution.
Hence the prevailing bifurcations are of either the symmetry-breaking or the
Hopf type.

####  III.6.2 Bifurcations in the double-gyre problem

The development of a comprehensive nonlinear theory of the double-gyre
circulation over the last two decades has gone through four main steps. These
four steps can be followed through the bifurcation tree in Fig. 31.

##### Symmetry-breaking bifurcation.

The “trunk” of the bifurcation tree is plotted as the solid blue line in the
lower part of the figure. When the forcing $\tau$ is weak or the dissipation
$\nu$ is large, there is only one steady solution, which is antisymmetric with
respect to the mid-axis of the basin. This solution exhibits two large gyres,
along with their $\beta$-induced western boundary currents. Away from the
western boundary, such a near-linear solution (not shown) is dominated by so-
called Sverdrup balance between wind stress curl and the meridional mass
transport Gill (1982); Sverdrup (1947).

Figure 31: Generic bifurcation diagram for the barotropic quasi-geostrophic
model of the double-gyre problem: the asymmetry of the solution is plotted
versus the intensity of the wind stress $\tau$. The streamfunction field is
plotted for a steady-state solution associated with each of the three
branches; positive values in red and negative ones in blue. After Simonnet
_et al._ (2005).

The first generic bifurcation of this quasi-geostrophic model was found to be
a genuine pitchfork bifurcation that breaks the system’s symmetry as the
nonlinearity becomes large enough with increasing wind stress intensity $\tau$
Jiang  _et al._ (1995b); Cessi and Ierley (1995). As $\tau$ increases, the
near-linear Sverdrup solution that lies along the solid blue line in the
figure develops an eastward jet along the mid-axis, which penetrates farther
into the domain and also forms two intense recirculation vortices, on either
side of the jet and near the western boundary of the domain.

The resulting more intense, and hence more nonlinear solution is still
antisymmetric about the mid-axis, but loses its stability for some critical
value of the wind-stress intensity, $\tau=\tau_{\rm P}$. This value is
indicated by the filled square on the symmetry axis of Fig. 31 and is labeled
“Pitchfork” in the figure.

A pair of mirror-symmetric solutions emerges and it is plotted as the two red
solid lines in the figure’s middle part. The streamfunction fields associated
with the two stable steady-state branches have a rather different vorticity
distribution and they are plotted in the two small panels to the upper-left
and upper-right of Fig. 31. In particular, the jet in such a solution exhibits
a large, stationary meander, reminiscent of the semi-permanent one that occurs
in the Gulf Stream, just downstream of Cape Hatteras. These asymmetric flows
are characterized by one recirculation vortex being stronger in intensity than
the other; accordingly the eastward jet is deflected either to the southeast,
as is the case in the observations for the North Atlantic, or to the
northeast.

##### Gyre modes.

The next step in the theoretical treatment of the problem was taken in part
concurrently with the first one above Jiang  _et al._ (1995b) and in part
shortly thereafter Speich  _et al._ (1995a); Dijkstra and Katsman (1997);
Sheremet  _et al._ (1997). It involved the study of time-periodic
instabilities that arise through Hopf bifurcation from either an antisymmetric
or an asymmetric steady flow. Some of these studies treated wind-driven
circulation models limited to a stand-alone, single gyre Pedlosky (1996);
Sheremet  _et al._ (1997); such a model concentrates on the larger,
subtropical gyre while neglecting the smaller, subpolar one.

The overall idea was to develop a full, generic picture of the time-dependent
behavior of the solutions in more turbulent regimes, by classifying the
various instabilities in a comprehensive way. However, it quickly appeared
that a particular kind of instability leads to so-called gyre modes Jiang  _et
al._ (1995b); Speich  _et al._ (1995a), and was prevalent across the full
hierarchy of models of the double-gyre circulation; furthermore, this
instability triggers the lowest nonzero frequency present in all such models
Dijkstra (2005); Dijkstra and Ghil (2005).

These gyre modes always appear after the first pitchfork bifurcation, and it
took several years to understand their genesis: gyre modes arise as two
eigenvalues merge — one of the two is associated with a symmetric
eigenfunction and responsible for the pitchfork bifurcation, the other one
with an antisymmetric eigenfunction Simonnet and Dijkstra (2002). This merging
is marked by a filled circle on the left branch of antisymmetric stationary
solutions and is labeled as $\bf M$ in Fig. 31.

Such a merging phenomenon is not a bifurcation in the term’s usual meaning:
Although it corresponds to a topological change in phase space, the
oscillatory behavior at and near $\bf M$ is damped. Nevertheless, this
oscillatory eigenmode is eventually destabilized through a Hopf bifurcation,
which is indicated in Fig. 31 by a heavy dot marked “Hopf,” from which a
stylized limit cycle emerges. A mirror-symmetric $\bf M$ and Hopf bifurcation
also occur on the right branch of stationary solutions, but have been omitted
for visual clarity. This merging is as generic as the $\bf Pitchfork$
bifurcation in the figure, and arises in much more complex situations and
models (e.g., Simonnet  _et al._ , 2003b; Simonnet, 2005).

In fact, such merging is quite common in small-dimensional dynamical systems
with symmetry, as exemplified by the unfolding of codimension-2 bifurcations
of Bogdanov-Takens type Guckenheimer and Holmes (1983). In particular, the
fact that gyre modes trigger the longest, multi-annual periodicity of the
model is due to the frequency of this mode growing quadratically with the
control parameter from zero — i.e., from infinite period — until nonlinear
saturation sets in (e.g., Simonnet and Dijkstra, 2002; Simonnet  _et al._ ,
2009).

More generally, Hopf bifurcations give rise to features that recur more-or-
less periodically in fully turbulent planetary-scale flows, atmospheric,
oceanic and coupled Ghil and Childress (1987); Dijkstra (2005); Dijkstra and
Ghil (2005); Ghil (2015, 2017). It is precisely this kind of near-periodic
recurrence that is identified in the climate sciences as LFV.

##### Global bifurcations.

The bifurcations studied so far — in this subsection, as well as in the
preceding ones — are collectively known as local bifurcations: they result
from an instability of a specific solution that arises at a particular value
of a control parameter. This term is meant to distinguish them from the global
bifurcations that will be studied forthwith.

The importance of the gyre modes was further confirmed through an even more
puzzling discovery. Several authors realized, independently of each other,
that the low-frequency dynamics of their respective double-gyre models was
driven by intense relaxation oscillations of the jet Simonnet  _et al._
(1995); Meacham (2000); Chang  _et al._ (2001); Nadiga and Luce (2001);
Simonnet  _et al._ (2003a, b, 2005). These relaxation oscillations, already
described by Jiang  _et al._ (1995b) and Speich  _et al._ (1995b), were now
attributed to a homoclinic bifurcation, which is no longer due to a linear
instability of an existing solution but to a so-called homoclinic
reconnection, whose character is global in phase space Guckenheimer and Holmes
(1983); Ghil and Childress (1987). In effect, the quasi-geostrophic model
reviewed herein undergoes a genuine homoclinic bifurcation that is generic
across the full hierarchy of double-gyre models.

This bifurcation is due to the growth and eventual merging of the two limit
cycles, each of which arises from either one of the two mutually symmetric
Hopf bifurcations. The corresponding bifurcation is marked in Fig. 31 by a
filled circle and labeled “Homoclinic”; the reconnecting orbit itself is
illustrated in the figure by a stylized lemniscate, and plotted accurately in
Simonnet  _et al._ (2005, Fig. 2). This global bifurcation is associated with
chaotic behavior of the flow due to the Shilnikov phenomenon Nadiga and Luce
(2001); Simonnet  _et al._ (2005), which induces Smale horseshoes in phase
space.

The connection between such homoclinic bifurcations and gyre modes was not
immediately obvious, but Simonnet  _et al._ (2005) emphasized that the two
were part of a single, global dynamical phenomenon. The homoclinic bifurcation
indeed results from the unfolding of the gyre modes’ limit cycles. This
familiar dynamical scenario is again well illustrated by the unfolding of a
codimension-$2$ Bogdanov-Takens bifurcation, where the homoclinic orbits
emerge naturally.

Since homoclinic orbits have an infinite period, it was natural to hypothesize
that the gyre-mode mechanism, in this broader, global-bifurcation context,
gave rise to the observed $7$-yr and $14$-yr North Atlantic oscillations.
Although this hypothesis may appear a little farfetched — given the simplicity
of the double-gyre models analyzed so far — it is reinforced by results with
much more detailed models in the hierarchy (e.g., Dijkstra and Ghil, 2005;
Dijkstra, 2005; Vannitsem  _et al._ , 2015).

The successive-bifurcation theory appears therewith to be fairly complete for
barotropic, single-layer models of the double-gyre circulation. This theory
also provides a self-consistent, plausible explanation for the climatically
important 7-year and 14-year oscillations of the oceanic circulation and the
related atmospheric phenomena in and around the North Atlantic basin Dijkstra
and Ghil (2005); Dijkstra (2005); Moron  _et al._ (1998); Wunsch (1999); Da
Costa and Colin de Verdiére (2004); Plaut  _et al._ (1995); Feliks  _et al._
(2004, 2007, 2010, 2011); Kondrashov  _et al._ (2005a); Simonnet  _et al._
(2003b, 2005). The dominant 7- and 14-year modes of this theory survive,
moreover, perturbation by seasonal-cycle changes in the intensity and
meridional position of the westerly winds Sushama  _et al._ (2007).

In baroclinic models, with two or more active layers of different density,
baroclinic instabilities Dijkstra and Ghil (2005); Feliks  _et al._ (2007);
Gill (1982); Ghil and Childress (1987); Pedlosky (1987); Stommel (1965);
Pedlosky (1996); Simonnet  _et al._ (2003b); Berloff  _et al._ (2007);
Kravtsov  _et al._ (2007) surely play a fundamental role, as they do in the
observed dynamics of the oceans. However, it is not known to what extent
baroclinic instabilities can destroy gyre-mode dynamics. The difficulty lies
in a deeper understanding of the so-called rectification process Katsman  _et
al._ (1998), which arises from the nonzero mean effect of the baroclinic
eddying and meandering of the flow on its barotropic component.

Roughly speaking, rectification drives the dynamics farther away from any
stationary solutions. In this situation, dynamical systems theory by itself
cannot be used as a full explanation of complex, observed behavior resulting
from successive bifurcations that are rooted in simple stationary or periodic
solutions.

Other tools from statistical mechanics and nonequilibrium thermodynamics have,
therefore, to be considered Bouchet and Sommeria (2002); Farrell and Ioannou
(1996); Chavanis and Sommeria (1996); Lucarini and Sarno (2011); Lucarini  _et
al._ (2014); 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.

###  III.7 Multiple Scales: Stochastic and Memory Effects

In Sect. II.2, 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.

####  III.7.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}$ yr${}^{-1}=f_{1}\leq f\leq
f_{2}=10^{0}$ yr-1, and write a model of this variability as

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

We have seen that oscillatory modes of both the THC, cf. Sect. III.4 , and of
the wind-driven circulation, cf. Sec. III.6, 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\leq f_{0}<10^{-2}$ 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,

$\bm{\mu}=\bm{\mu}(\epsilon t),\quad 0<\epsilon\ll 1,$ (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

$\displaystyle x^{\prime}$ $\displaystyle=$ $\displaystyle F(x,y;\epsilon),$
(25a) $\displaystyle y^{\prime}$ $\displaystyle=$ $\displaystyle\epsilon
G(x,y;\epsilon).$ (25b)

Here $\mathbf{z}=(x,y)^{\rm T}$, with $x$ the fast and $y$ the slow variable,
and $\mathbf{H}=(F,G)^{\rm T}$, while $(\cdot)^{\rm T}$ designates the
transpose and $(\cdot)^{\prime}=\mathrm{d}(\cdot)/\mathrm{d}t$ 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

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

in which the dot stands for differentiation $\mathrm{d}(\cdot)/\mathrm{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

$\displaystyle x^{\prime}$ $\displaystyle=$ $\displaystyle F(x,y;0),$ (27a)
$\displaystyle y^{\prime}$ $\displaystyle=$ $\displaystyle 0,$ (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

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

The algebraic equation $F(x,y;0)=0$ defines the critical manifold
$\cal{M}_{\rm 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.3.2) or the wind-
driven circulation of the oceans (Sect. III.6.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}^{\rm
T},\mathbf{Y}^{\rm T})^{\rm 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.2 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.

####  III.7.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.3.2 and III.7.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_{\rm B}(\lambda)+F(\lambda);$ (29)

here $M_{\rm B}(\lambda)$ is the nonnegative mass flux through the cloud base,
and $J\otimes M_{\rm 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_{\rm B}(\lambda)+F(\lambda).$ (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_{\rm{adj}}$ of a cloud ensemble to the mean state, $\epsilon\sim
1/\tau_{\mathrm{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;$ (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_{\rm B}^{2}/K,\quad\tau_{\rm D}=K/D.$ (32)

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

####  III.7.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

${\rm d}u=-\lambda u{\rm d}t+\eta(t).$ (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{\rm d}W$, 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}^{\prime}$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

$\displaystyle\dot{x}$ $\displaystyle=$
$\displaystyle\epsilon^{-1}f_{0}(y^{(\epsilon)})+f_{1}(x^{(\epsilon)},y^{(\epsilon)}),\;x^{(\epsilon)}(0)=x_{0},$
(34a) $\displaystyle\dot{y}^{(\epsilon)}$ $\displaystyle=$
$\displaystyle\epsilon^{-2}g(y^{(\epsilon)}),\;y^{(\epsilon)}(0)=y_{0},$ (34b)

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

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
the chaotic and fast $y^{(\epsilon)}$ induces, as $\epsilon\to 0$, a white-
noise driving of the slow $x$. Note that one needs $d^{\prime}\geq 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^{\prime}}$, 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)\rightarrow_{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}{\overline{F}}(X(s))\mathrm{d}s+\sigma W(t).$ (35)

Here $W$ is the Brownian motion with variance $\sigma^{2}$, such that the
white noise in Eq. (33) be given by $\eta=\sigma\mathrm{d}W$. Moreover,
${\overline{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\mathrm{with}\;x\in\mathbb{R}^{d},\;y\in\mathbb{R}^{d^{\prime}}$
(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.5.

####  III.7.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., Tziperman  _et al._ , 1994a; Ghil  _et al._ , 2015, 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.

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

$\displaystyle\dot{x}$ $\displaystyle=$ $\displaystyle a_{11}x+a_{12}y,$ (37a)
$\displaystyle\dot{y}$ $\displaystyle=$ $\displaystyle a_{21}x+a_{22}y.$ (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)\mathrm{d}s,$

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

$\dot{x}=a_{11}x+\int_{0}^{t}K(t-s)x(s)\mathrm{d}s+f(t).$ (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.7.1 and III.7.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).

##### 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{\mathbf{x}}=-\mathbf{A}\mathbf{x}+\mathbf{B}(\mathbf{x},\mathbf{x})+\mathbf{L}(\mathbf{x},\mathbf{r}_{t}^{l},\xi_{t},t),\;0\leq
l\leq\mathbf{L}-1.$ (39)

Here $\mathbf{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 $-\mathbf{A}\mathbf{x}$ and $\mathbf{B}(\mathbf{x},\mathbf{x})$
represent, respectively, the linear dissipation and the quadratic self-
interactions of these modes, while $\mathbf{L}$ is a time-dependent operator
that is bilinear in the resolved variables $\mathbf{x}$ and the unresolved
ones $\mathbf{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

$\mathrm{d}\mathbf{r}_{t}^{l}=\mathbf{M}_{l}(\mathbf{x},\mathbf{r}_{t}^{0},\ldots,\mathbf{r}_{t}^{l})\mathrm{d}+\mathbf{r}_{t}^{l+1}\mathrm{d}t.$

At each level $l$, the coupling between the variable $\mathbf{r}_{t}^{l}$ and
the previous-level variables
$(\mathbf{x},\mathbf{r}_{t}^{0},\ldots,\mathbf{r}_{t}^{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}=\xi_{t}$ has a
lag$-1$ vanishing autocorrelation. The stochastic residuals
$\mathbf{r}_{t}^{l}$, obtained in this recursive minimization procedure, are
ordered in decreasing order of decorrelation time, from $\mathbf{r}_{t}^{0}$
to $\xi_{t}$.

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

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

$\displaystyle\mathrm{d}x$ $\displaystyle=$
$\displaystyle(f(x)+r)\mathrm{d}t,$ (40a) $\displaystyle\mathrm{d}r$
$\displaystyle=$ $\displaystyle(\gamma x-\alpha
r)\mathrm{d}t+\mathrm{d}W_{t},$ (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)\mathrm{d}s+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)\mathrm{d}s+W_{t}(\omega),$ (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).

##### 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\times 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 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).

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

Kondrashov  _et al._ (2005b) fitted the global SST field between
$30^{\circ}$S–$60^{\circ}$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 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^{\circ}$S–$5^{\circ}$N, $150^{\circ}$–$90^{\circ}$W) over
which SST anomalies are averaged to obtain the Niño-3 index, already mentioned
in Sect. III.5.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.

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

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.

##### 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 described by the following evolution equations:

$\displaystyle\dot{\mathbf{X}}$ $\displaystyle=$ $\displaystyle
f_{\mathbf{X}}(\mathbf{X})+\epsilon\Psi_{\mathbf{X}}(\mathbf{X},\mathbf{Y}),$
(42a) $\displaystyle\dot{\mathbf{Y}}$ $\displaystyle=$ $\displaystyle
f_{\mathbf{Y}}(\mathbf{Y})+\epsilon\Psi_{\mathbf{Y}}(\mathbf{X},\mathbf{Y});$
(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)\mathrm{d}s.$ (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 $\mathbf{Y}$ variables, at $\epsilon=0$. In the limit of
infinite time scale separation between the $\mathbf{X}$ and $\mathbf{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_{\mathbf{Y}}=0$ in the master–slave system (42), the non-Markovian term
is identically zero, as expected.

Figure 34: Schematic representation of the three terms contributing to the
parametrization of the fast variables $\mathbf{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).

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.5.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 (2018b, a) 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.

###  IV.1 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.1 where the net
radiation $R=R_{\rm i}-R_{\rm 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 $c{\mathrm{d}T}/{\mathrm{d}t}=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,\ldots,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\left(T,\alpha_{1}(T),\ldots,\alpha_{N}(T)\right)$. 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 CO2 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

$\displaystyle\Delta\tilde{R}$
$\displaystyle=-{\displaystyle\frac{\mathrm{d}R}{\mathrm{d}T}}\bigg{|}_{T=T_{0}}\Delta
T+\mathcal{O}\left((\Delta\,T)^{2}\right)$ (44)
$\displaystyle={\displaystyle\frac{\partial R}{\partial
T}}\bigg{|}_{T=T_{0}}\Delta\,T-\sum_{k=1}^{N}{\displaystyle\frac{\partial
R}{\partial\alpha_{k}}}{\displaystyle\frac{\partial\alpha_{k}}{\partial
T}\bigg{|}_{T=T_{0}}\Delta T}$
$\displaystyle+\mathcal{O}\left((\Delta\,T)^{2}\right).$

Here, $\mathcal{O}(x)$ is a function such that $\mathcal{O}(x)\leq C\,x$ 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

$\displaystyle\frac{1}{\lambda_{0}(T_{0})}=$
$\displaystyle{\displaystyle-\frac{\partial R}{\partial
T}}\bigg{|}_{T=T_{0}},$ (45a) $\displaystyle f_{k}(T_{0})=$
$\displaystyle-\lambda_{0}(T_{0}){\displaystyle\frac{\partial
R}{\partial\alpha_{k}}\frac{\partial\alpha_{k}}{\partial
T}}\bigg{|}_{T=T_{0}},$ (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}\left((\Delta\,T)^{2}\right),$
(46)

which readily leads to

$\Delta
T=\Lambda(T_{0})\Delta\tilde{R}+\mathcal{O}\left((\Delta\,T)^{2}\right).$ (47)

Here

$\Lambda(T_{0})=-\frac{\mathrm{d}R}{\mathrm{d}T}\bigg{|}_{T=T_{0}}=\frac{\lambda_{0}((T_{0}))}{1-\sum_{k=1}^{N}f_{k}(T_{0})}$
(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.1, 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., Watson and Lovelock, 1983; Zeng and
Neelin, 2000; Rombouts and Ghil, 2015, 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:

$\displaystyle ECS$
$\displaystyle=\Lambda(T_{0})\Delta\tilde{R}_{2\times{\rm{CO}_{2}}}$ (49)
$\displaystyle={\displaystyle\frac{\lambda_{0}((T_{0}))}{1-\sum_{k=1}^{N}f_{k}(T_{0})}}\Delta\tilde{R}_{2\times{\rm
CO}_{2}}.$

Note that the radiative forcing is, to a good approximation, proportional to
the logarithm of the CO2 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 CO2
concentration is twice as large as the ECS; see Fig. 35.

Figure 35: Test of the linear scaling of the long-term climate response with
respect to the CO2 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 CO2
concentration. Reproduced with permission from Pfister and Stocker (2017).
Figure 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 CO2 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-CO2 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).

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 CO2 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).

###  IV.2 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 IPCC
(2001, 2007, 2014a). In fact, Charney  _et al._ (1979) estimated the ECS
uncertainty as $1.5$–$4.5$ K for CO2 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); 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.

Figure 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_{[{\rm CO}_{2},{\rm
LI}]}$ on the $y$-axis refers to sensitivity with respect to CO2 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).

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
CO2 concentration, so that $\Delta\tilde{R}_{2\times{\rm CO}_{2}}$ depends on
the concentration’s reference value. In fact, this dependence is weak across a
large range of CO2 concentrations, but it is greatly strengthened by optical
saturation effects in the CO2 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 CO2 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 irrandiance requires
exceedingly high CO2 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 CO2 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).

Figure 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/.

###  IV.3 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
CO2 concentration. TCR is defined as the change in the globally averaged
surface air temperature recorded at the time at which CO2 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.5 herein.

Figure 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 CO2 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 CO2
level, i.e., at $2\times{\rm CO}_{2}$ or at $4\times{\rm CO}_{2}$. Reproduced
with permission from IPCC (2001, Fig. 9.1).

###  IV.4 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,
$\mathbf{X}_{0}=\mathbf{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, 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_{\rm W}$ between the two invariant measures. The Wasserstein
distance or “earth mover’s distance” $d_{\rm 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_{\rm 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.

Figure 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 CO2 concentration, say, blue dash-dotted line) changes suddenly,
global temperature (red heavy solid) undergoes a transition. (a) Only the mean
temperature $\bar{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.

###  IV.5 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.

####  IV.5.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 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)$ (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\geq t_{0}$, with $\phi(s,s)=\mathbf{1}$, and
it thus generates a two-parameter semi-group. In the autonomous case, time-
translational invariance reduces the latter to a one-parameter semigroup
since, $\forall t\geq 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(\textrm{d}x)$ Eckmann and Ruelle (1985); Ledrappier and Young
(1988).777Among 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), 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.

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

$\mathrm{d}{x}=F(x,t)\mathrm{d}t+g(x)\mathrm{d}\eta,$ (51)

where $\eta=\eta(t;\omega)$ is a Wiener process, while $\omega$ labels the
particular realization of this random process, and $\mathrm{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 Chekroun  _et al._ (2011); Carvalho  _et al._ (2013)
and Ghil  _et al._ (2008, Appendix A). Careful numerical applications of PBAs
to explain the wind-driven circulation and the THC are now available Sevellec
and Fedorov (2015); Pierini  _et al._ (2016, 2018).

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._ (2011);
Bódai and Tél (2012); Bódai  _et al._ (2013); Drótos  _et al._ (2015).

Figure 41: A cartoon depicting the pullback attractor (PBA) for a
nonautonomous system as in Eq. (50). The set $\mathbf{B}$ of finite Lebesgue
measure, initialized at time $t=t_{0}-\tau$, evolves towards the set
$\mathbf{A}(t_{0})$. The construction is exact in the limit
$\tau\rightarrow\infty$; for a random attractor, see Ghil  _et al._ (2008,
Figs. A.1 and A.2). Adapted from Sevellec and Fedorov (2015), with permission.

Note that, in the most recent IPCC reports IPCC (2001, 2007, 2014a) and
according to the standard protocols described in Sect. II.4.2 , future climate
projections 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.

####  IV.5.2 Fluctuation Dissipation and Climate Change