# The physics of climate variability and climate change

Michael Ghil, Valerio Lucarini

► To cite this version:

Michael Ghil, Valerio Lucarini. The physics of climate variability and climate change. Reviews of Modern Physics, 2020, 92, 10.1103/RevModPhys.92.035002. insu-03726977

HAL Id: insu-03726977

https://insu.hal.science/insu-03726977v1

Submitted on 28 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# The Physics of Climate Variability and Climate Change

Michael Ghil

Geosciences Department and Laboratoire de Météorologie Dynamique (CNRS and IPSL),
Ecole Normale Supérieure and PSL University,
F-75231 Paris Cedex 05,
France
Department of Atmospheric and Oceanic Sciences,
University of California,
Los Angeles, CA 90095-1565,
USA

Valerio Lucarini

Department of Mathematics and Statistics,
University of Reading,
Reading, RG66AX,
UK
Centre for the Mathematics of Planet Earth,
University of Reading,
Reading, RG66AX,
UK
CEN - Institute of Meteorology,
University of Hamburg,
Hamburg, 20144,
Germany

(Dated: February 27, 2020)

The climate system is a forced, dissipative, nonlinear, complex and heterogeneous system that is out of thermodynamic equilibrium. The system exhibits natural variability on many scales of motion, in time as well as space, and it is subject to various external forcings, natural as well as anthropogenic. This paper reviews the observational evidence on climate phenomena and the governing equations of planetary-scale flow, as well as presenting the key concept of a hierarchy of models as used in the climate sciences. Recent advances in the application of dynamical systems theory, on the one hand, and of nonequilibrium statistical physics, on the other, are brought together for the first time and shown to complement each other in helping understand and predict the system's behavior. These complementary points of view permit a self-consistent handling of subgrid-scale phenomena as stochastic processes, as well as a unified handling of natural climate variability and forced climate change, along with a treatment of the crucial issues of climate sensitivity, response, and predictability.

PACS numbers: 92.05.Df, 92.30.Bc, 92.60.Ry, 92.70.Np, 05.45.-a, 05.70.Ln

## CONTENTS

I. Introduction and Motivation 2

A. Basic Facts of the Climate Sciences 2

B. More Than “Just” Science 4

1. The Intergovernmental Panel on Climate Change 4

2. Hockey Stick Controversy and Climate Blogs 5

C. This Review 6

II. A Brief Introduction to Climate Dynamics 7

A. Climate Observations: Direct and Indirect 7

1. Instrumental Data and Reanalyses 7

2. Proxy Data 11

B. Climate Variability on Multiple Time Scales 12

1. A Survey of Climatic Time Scales 13

2. Atmospheric Variability in Mid-latitudes 15

C. Basic Properties and Fundamental Equations 17

1. Governing Equations 17

2. Approximate Balances and Filtering 18

3. Quasi-Geostrophy and Weather Forecasting 21

D. Climate Prediction and Climate Model Performance 22

1. Predicting the State of the System 24

2. Predicting the System's Statistical Properties 24

3. Metrics for Model Validation 26

III. Climate Variability and the Modeling Hierarchy 27

A. Radiation Balance and Energy-Balance Models 27

B. Other Atmospheric Processes and Models 30

C. Oscillations in the Oceans' Thermohaline Circulation 31

1. Theory and Simple Models 31

2. Bifurcation Diagrams for GCMs 33

D. Bistability, Oscillations and Bifurcations 34

1. Bistability and Steady-State Bifurcations 34

2. Oscillatory Instabilities and Hopf Bifurcations 35

E. Main Modes of Variability 36

1. Modes of Variability and Extended-Range Prediction 36

2. Coupled Atmosphere-Ocean Modes of Variability 37

3. Atmospheric Low-Frequency Variability 38
<table><tr><td>F. Internal Variability and Routes to Chaos</td><td>40</td></tr><tr><td>1. A simple model of the double-gyre circulation</td><td>41</td></tr><tr><td>2. Bifurcations in the double-gyre problem</td><td>41</td></tr><tr><td>G. Multiple Scales: Stochastic and Memory Effects</td><td>44</td></tr><tr><td>1. Fast Scales and Their Deterministic Parametrization</td><td>44</td></tr><tr><td>2. An Example: Convective Parametrization</td><td>45</td></tr><tr><td>3. Stochastic Parametrizations</td><td>45</td></tr><tr><td>4. Modeling Memory Effects</td><td>47</td></tr></table>

<table><tr><td>IV. Climate Sensitivity and Response</td><td>50</td></tr><tr><td>A. A Simple Framework for Climate Sensitivity</td><td>50</td></tr><tr><td>B. Climate Sensitivity: Uncertainties and Ambiguities</td><td>51</td></tr><tr><td>C. Transient Climate Response (TCR)</td><td>53</td></tr><tr><td>D. Beyond Climate Sensitivity</td><td>53</td></tr><tr><td>E. A General Framework for Climate Response</td><td>54</td></tr><tr><td>1. Pullback Attractors (PBAs)</td><td>54</td></tr><tr><td>2. Fluctuation Dissipation and Climate Change</td><td>56</td></tr><tr><td>3. Ruelle Response Theory</td><td>56</td></tr><tr><td>4. Climate Change Prediction via Ruelle Response Theory</td><td>58</td></tr><tr><td>5. Slow Correlation Decay and Sensitive Parameter Dependence</td><td>59</td></tr></table>

<table><tr><td>V. Critical Transitions and Edge States</td><td>60</td></tr><tr><td>A. Bistability for Gradient Flows and EBMs</td><td>62</td></tr><tr><td>B. Finding the Edge States</td><td>64</td></tr><tr><td>C. Invariant Measures and Noise-induced Transitions</td><td>65</td></tr><tr><td>D. Nearing Critical Transitions</td><td>66</td></tr><tr><td>E. Chaos-to-Chaos Transition</td><td>67</td></tr></table>

<table><tr><td>VI. Concluding Remarks</td><td>68</td></tr><tr><td>Acknowledgments</td><td>70</td></tr><tr><td>A. Acronyms</td><td>71</td></tr><tr><td>References</td><td>72</td></tr></table>

# I. INTRODUCTION AND MOTIVATION

## A. Basic Facts of the Climate Sciences

The climate system is a forced, dissipative, chaotic system that is out of equilibrium and whose complex natural variability arises from the interplay of positive and negative feedbacks, instabilities and saturation mechanisms. These processes span a broad range of spatial and temporal scales and include many chemical species and all physical phases. The system's heterogeneous phenomenology includes the mycrophysics of clouds, cloud-radiation interactions, atmospheric and oceanic boundary layers, and several scales of turbulence (Ghil, 2019); it evolves, furthermore, under the action of large-scale agents that drive and modulate its evolution, mainly differential solar heating and the Earth's rotation and gravitation.

As is often the case, the complexity of the physics is interwoven with the chaotic character of the dynamics. Moreover, the climate system's large natural variability on different time scales is strongly affected by relatively small changes in the forcing, anthropogenic as well as natural (Ghil and Childress, 1987; Lucarini et al., 2014; Peixoto and Oort, 1992)

<figure><img src="image_2.png" /><figcaption>FIG. 1: Globally averaged energy fluxes in the Earth system [W m⁻²]. The fluxes on the left represent solar radiation in the visible and the ultraviolet, those on the right terrestrial radiation in the infrared, and those in the middle nonradiative fluxes. Reproduced from Trenberth et al. (2009) ©American Meteorological Society; used with permission.</figcaption></figure>

On the macroscopic level, climate is driven by differences in the absorption of solar radiation throughout the depth of the atmosphere, as well as in a narrow surface layer of the ocean and of the soil; the system's actual governing equations are given in Sect. II.C below. The prevalence of absorption at the surface and in the atmosphere's lower levels leads, through several processes, to compensating vertical energy fluxes — most notably, fluxes of infrared radiation throughout the atmosphere and convective motions in the troposphere; see Fig. 1.

More solar radiation is absorbed in the low latitudes, leading to horizontal energy fluxes as well. The atmosphere's large-scale circulation is, to first order, a result of these horizontal and vertical fluxes arising from the gradients in solar radiation absorption, in which the hydrological cycle plays a key role as well. The ocean circulation, in turn, is set into motion by surface or near-surface exchanges of mass, momentum and energy with the atmosphere: the so-called wind-driven component of the circulation is due mainly to the wind stress and the thermohaline one is due mainly to buoyancy fluxes (Dijkstra, 2005; Dijkstra and Ghil, 2005). The coupled atmospheric and oceanic circulation reduces the temperature differences between tropics and polar regions with respect to that on an otherwise similar planet with no horizontal energy transfers (Held, 2001; Lorenz, 1967; Lucarini and Ragone, 2011; Peixoto and Oort, 1992). At steady state, the convergence of enthalpy transported by the atmosphere and by the oceans compensates the radiative imbalance at the top of the atmosphere; see Fig. 2.

The classical theory of the general circulation of the atmosphere (Lorenz, 1967) describes in further detail how the mechanisms of energy generation, conversion and dissipation produce the observed circulation, which deviates
substantially from the highly idealized, zonally symmet-
ric picture sketched so far. According to Lorenz (1955),
atmospheric large-scale flows result from the conversion
of available potential energy — which is produced by
the atmosphere’s differential heating — into kinetic en-
ergy, and the Lorenz (1967) energy cycle is completed
by energy cascading to smaller scales to be eventually
dissipated. McWilliams (2019) provides an up-to-date
criticism of and further perspective on this theory.

Overall, the climate system can be seen as a ther-
mal engine capable of transforming radiative heat into
mechanical energy with a given, highly suboptimal
efficiency, given the many irreversible processes that
make it less than ideal (Kleidon and Lorenz, 2005; Lu-
carini, 2009b; Lucarini et al., 2014; Pauluis and Held,
2002). This conversion occurs through genuinely three-
dimensional (3-D) baroclinic instabilities (Charney, 1947;
Eady, 1949) that are triggered by large temperature gra-
dients and would break zonal symmetry even on a so-
called aqua-planet, with no topographic or thermal asym-
metries at its surface. These instabilities give rise to a
negative feedback, as they tend to reduce the temper-
ature gradients they feed upon by favoring the mixing
between masses of fluids at different temperatures.

Note that, while these baroclinic and other large-scale
instabilities do act as negative feedbacks, they cannot be
treated as diffusive, Onsager (1931)-like processes. Faced
with the Earth system's complexity discussed herein and
illustrated by Fig. 3, the closure of the coupled thermody-
namical equations governing the general circulation of the
atmosphere and oceans would provide a self-consistent
theory of climate. Such a theory should able to con-
nect instabilities and large-scale stabilizing processes on
longer spatial and temporal scales, and to predict its re-
sponse to a variety of forcings, both natural and anthro-
pogenic (Ghil, 2015; Ghil and Childress, 1987; Lucarini,
2009b; Lucarini et al., 2014). This goal is being actively
pursued but is still out of reach at the time of this writing
(e.g., Ghil, 2019, and references); see also Sects. IV and
V herein. The observed persistence of spatial gradients in
chemical concentrations and temperatures, as well as the
associated mass and energy fluxes, are a basic signature
of the climate system's intrinsic disequilibrium.

Figure 3 emphasizes, moreover, that the fluid and the
solid parts of the Earth system are coupled on even longer
time scales, on which geochemical processes become of
paramount importance (Kleidon, 2009; Rothman et al.,
2003). In contrast, closed, isolated systems cannot main-
tain disequilibrium and have to evolve towards homoge-
nous thermodynamical equilibrium, as a result of the sec-
ond law of thermodynamics (Prigogine, 1961).

Studying the climate system’s entropy budget provides
a good global perspective on this system. The Earth
as a whole absorbs shortwave radiation carried by low-
entropy solar photons at T<sub>Sun</sub> ≈ 6000 K and emits in-
frared radiation to space via high-entropy thermal pho-

<figure><img src="image_3.png" /><figcaption>FIG. 2: Meridional distribution of net radiative fluxes and of horizontal enthalpy fluxes. (a) Observed zonally averaged radiative imbalance at the top of the atmosphere from the ERBE experiment (1985–1989). (b) Inferred meridional enthalpy transport from ERBE observations (solid line) and estimate of the atmospheric enthalpy transport from two reanalysis datasets (ECMWF and NCEP). Reproduced from Trenberth and Caron (2001) c ○ American Meteorological Society; used with permission.</figcaption></figure>

tons at T<sub>Earth</sub> ≈ 255 K (Lucarini et al., 2014; Peixoto and Oort, 1992). Besides the viscous dissipation of kinetic energy, many other irreversible processes—such as turbulent diffusion of heat and chemical species, irreversible phase transitions associated with various hydrological processes, and chemical reactions involved in the biogeochemistry of the planet—contribute to the total material entropy production (Goody, 2000; Kleidon, 2009).

These and other important processes appear in the
schematic diagram of Fig. 3. In general, in a forced dissi-
pative system, entropy is continuously produced by irre-
versible processes and at steady state, this production is
balanced by a net outgoing flux of entropy at the system’s
boundaries (de Groot and Mazur, 1984; Prigogine, 1961);
in the case at hand, this flux leaves mainly through the
top of the atmosphere (Goody, 2000; Lucarini, 2009b).
Thus, on average, the climate system’s entropy budget is
balanced, just like its energy budget.
The phenomenology of the climate system is commonly approached by focusing on distinct and complementary aspects that include:

* wave-like features such as Rossby waves or equatorially trapped waves (e.g., Gill, 1982), which play a key role in the transport of energy, momentum, and water vapor, as well as in the study of atmospheric, oceanic and coupled-system predictability;

* particle-like features such as hurricanes, extratropical cyclones or oceanic vortices (e.g., McWilliams, 2019; Salmon, 1998), which strongly affect the local properties of the climate system and of its subsystems and subdomains;

* turbulent cascades, which are of crucial importance in the development of large eddies through the mechanism of geostrophic turbulence (Charney, 1971), as well as in mixing and dissipation within the planetary boundary layer (Zilitinkevich, 1975).

Each of these points of view is useful and they do overlap and complement each other (Ghil and Robertson, 2002; Lucarini et al., 2014) but neither by itself provides a comprehensive understanding of the properties of the climate system. It is a key objective of this review paper to provide the interested and motivated reader with the tools for achieving such a comprehensive understanding with predictive potential.

<figure><img src="image_4.png" /><figcaption>FIG. 3: Schematic diagram representing forcings, dissipative and mixing processes, gradients of temperature and chemical species, and coupling mechanisms across the Earth system. Bluish (reddish) colors refer to the fluid (solid) Earth. Reproduced from Kleidon (2010) with permission.</figcaption></figure>

While much progress has been achieved (Ghil, 2019), understanding and predicting the dynamics of the cli-

mate system faces — on top of all the difficulties that are intrinsic to any nonlinear, complex system out of equilibrium — the following additional obstacles that make it especially hard to grasp fully:

* the presence of well-defined subsystems — the atmosphere, the oceans, the cryosphere — characterized by distinct physical and chemical properties and widely differing time and space scales;

* the complex processes coupling these subsystems;

* the continuously varying set of forcings that result from fluctuations in the incoming solar radiation and in the processes, both natural and anthropogenic, that alter the atmospheric composition;

* the lack of scale separation between different processes, which requires a profound revision of the standard methods for model reduction, and calls for unavoidably complex parametrization of subgrid-scale processes in numerical models;

* the lack of detailed, homogeneous, high-resolution and long-lasting observations of climatic fields that leads to the need for combining direct and indirect measurements when trying to reconstruct past climate states preceding the industrial era; and, last but not least,

* the fact that we only have one realization of the processes that give rise to climate evolution in time.

For all these reasons, it is far from trivial to separate the climate system's response to various forcings from its natural variability in the absence of time-dependent forcings. In simpler words, and as noted already by Lorenz (1979), it is hard to separate forced and free climatic fluctuations (Lucarini et al., 2014, 2017; Lucarini and Sarno, 2011). This difficulty is a major stumbling block on the road to a unified theory of climate evolution (Ghil, 2015, 2017) but some promising ideas for overcoming it are emerging and will be addressed in Sects. IV and V herein; see also Ghil (2019).

## B. More Than “Just” Science

### 1. The Intergovernmental Panel on Climate Change

Besides the strictly scientific aspects of climate research, much of the recent interest in it has been driven by the accumulated observational and modeling evidence on the ways humans influence the climate system. In order to review and coordinate the research activities carried out by the scientific community in this respect, the United Nations Environment Programme (UNEP) and the World Meteorological Organization (WMO) established in 1988 the Intergovernmental Panel on Climate
Change (IPCC); its assessments reports (ARs) are issued every 4–6 years. By compiling systematic reviews of the scientific literature relevant to climate change, the ARs summarize the scientific progress, the open questions, and the bottlenecks regarding our ability to observe, model, understand and predict the climate system’s evolution.

More specifically, it is the IPCC Working Group I that focuses on the physical basis of climate change; see IPCC (2001, 2007, 2014a) for the three latest reports — AR3, AR4 and AR5 — in this area. Working Groups II and III are responsible for the reports that cover the advances in the interdisciplinary areas of adapting to climate change and of mitigating its impacts; see IPCC (2014b,c) for the contributions of Working Groups II and III, respectively, to AR5. AR6 is currently in preparation, cf. https://www.ipcc.ch/assessment-report/ar6/.

Moreover, the IPCC supports the preparation of special reports on themes which are of interest across two of the working groups, e.g. climatic extremes (IPCC, 2012), or across all three of them. The IPCC experience and working group structure is being replicated for addressing climate change at regional level, as in the case for the Hindu Kush Himalayan region, called sometimes the “third pole” (Wester et al., 2019).

The IPCC reports are based on the best science available and are policy-relevant but not policy-prescriptive. Their multi-stage review is supposed to guarantee neutrality but they are still inherently official, UN-sanctioned documents and have to bear the imprimatur of the IPCC's 195 member countries. Their release thus leads to considerable and often adversarial debates involving a variety of stakeholders from science, politics, civil society, and business; they also affect media production, cinema, videogames, and art at large, and are reflected more and more by them.

Climate change has thus become an increasingly central topic of discussion in the public arena, involving all levels of decision making from local through regional and on to global. In recent years, climate services have emerged as a new area at the intersection of science, technology, policy making and business. They emphasize tools to enable climate change adaptation and mitigation strategies, and have benefited from large public investments like the European Union's Copernicus Programme, https://climate.copernicus.eu.

The lack of substantial progress in national governments' and international bodies' tackling climate change has recently led to the rapid growth of global, young-people-driven grassroots movements like Extinction Rebellion, https://rebellion.earth/, and FridaysForFuture, https://www.fridaysforfuture.org. Some countries, like the United Kingdom, have declared a state of climate emergency, cf. https://www.bbc.co.uk/news/uk-politics-48126677, and some influential media outlets have started to use the expression climate

crisis instead of climate change, https://tinyurl.com/y2v2jwzy. While of great consequence, the present paper does not further dwell on such socio-economic and political issues.

## 2. Hockey Stick Controversy and Climate Blogs

Mann et al. (1999, Fig. 3(a)) produced a temperature reconstruction from proxy data, cf. Sect. II.A, for the last 1000 yr, shown as part of the blue curve in Fig. 4. This curve was arguably the most striking, and hence controversial, scientific result contained in AR3 (IPCC, 2001) and it was dubbed, for obvious reasons, the hockey stick. The AR3 report combined into one figure — Fig. 1(b) of the Summary for Policy Makers (SPM) — this blue curve and the red curve shown in Fig. 4 here, which was based on instrumental data over the last century-and-one-half; see Fig. 1(a) of the SPM. This superposition purported to demonstrate that recent temperature increase was unprecedented over the last two millennia, in both values attained and rate of change.

Figure 1(b) of the AR3's SPM received an enormous deal of attention from the social and political forces wishing to underscore the urgency of tackling anthropogenic climate change. For the opposite reasons, the paper and its authors were the subject of intense political and judicial scrutiny and attack by other actors in the controversy, claiming that the paper was both politically motivated and scientifically unsound.

McIntyre and McKitrick (2005), among others, strongly criticized the results of Mann et al. (1999), claiming that the statistical procedures used for smoothly combining the diverse proxy records used — including tree rings, coral records, ice cores, and long historical records, with their diverse sources and ranges of uncertainty — into a single multiproxy record, and the latter with instrumental records, were marred by biases and underestimation of the actual statistical uncertainty. Later papers criticized in turn the statistical methods of McIntyre and McKitrick (2005) and confirmed the overall correctness of the hockey stick reconstruction (e.g., Huybers, 2005; Mann et al., 2008; PAGES, 2013; Taricco et al., 2009). NRC (2006) provided an excellent review of the state of our knowledge concerning the last two millennia of climate change and variability.

This controversy included the notorious “Climategate” incident, in which data hacked from the computer of a well-known UK scientist were used to support the thesis that scientific misconduct and data falsification had been routinely used to support the hockey stick reconstruction. These claims were later dismissed but they did lead to an important change in the relationship between the climate sciences, society, and politics, and in the way climate scientists interact among themselves and with the public. In certain countries — e.g., the UK and
<figure><img src="image_5.png" /></figure>

FIG. 4: Surface air temperature record for the last two millennia. Green dots show the 30-year average of the latest PAGES 2k reconstruction (PAGES, 2013), while the red curve shows the global mean temperature, according to HadCRUT4 data from 1850 onwards; the original “hockey stick” of Mann *et al.* (1999) is plotted in dark blue and its uncertainty range in light blue. Graph by Klaus Bitterman.

the USA — stringent rules have been imposed to ascertain that scientists working in governmental institutions have to publicly reveal the data they use in the preparation of scientific work, if formally requested to do so.

Faced with the confusion generated by the polemics, several leading scientists started blogs — e.g., www.realclimate.org, www.ClimateAudit.org, www.climate-lab-book.ac.uk, and judithcurry.com — in which scientific literature and key ideas are presented for a broader audience and debated outside the traditional media of peer-reviewed journals or public events, such as conferences and workshops. Most contributions are of high quality, but sometimes arguments appear to sink to the level of a bitter strife between those in favor and those against the reality of climate change and of the anthropogenic contribution to its causes.

### C. This Review

The main purpose of this review paper is to bring together a substantial body of literature published over the last few decades in the geosciences, as well as in mathematical and physical journals, and provide a comprehensive picture of climate dynamics. Moreover, this picture should appeal to a readership of physicists, and help stimulate interdisciplinary research activities.

For decades meteorology and oceanography, on the one side, and physics, on the other, have had a relatively low level of interaction, with by-and-large separate scientific gatherings and scholarly journals. Recent developments in dynamical systems theory, both finite- and infinite-

dimensional, as well as in random processes and statistical mechanics, have created a common language that makes it possible, at this time, to achieve a higher level of communication and mutual stimulation.

The key aspects of the field that we want to tackle here are the natural variability of the climate system, the deterministic and random processes that contribute to this variability, its response to perturbations, and the relations between internal and external causes of observed changes in the system. Moreover, we will present tools for the study of critical transitions in the climate system, which can help understand and possibly predict the potential for catastrophic climate change.

In Sect. II, we provide an overview for non-specialists of the way climate researchers collect and process information on the state of the atmosphere, the land surface, and the oceans. Next, the conservation laws and the equations that govern climatic processes are introduced.

An important characteristic of the climate system is the already mentioned coexistence and nonlinear interaction of multiple subsystems, processes and scales of motion. This state of affairs entails two important consequences that are also addressed in Sect. II. First is the need for scale-dependent filtering: on the positive side, this filtering leads to simplified equations; on the negative one, it calls for so-called parametrization of unresolved processes, i.e. for the representation of subgrid-scale processes in terms of the resolved, larger-scale ones. Second is the fact that no single model can encompass all the subsystems, processes and scales, hence the need for resorting to a hierarchy of models. The section ends by discussing present-day standard protocols for climate modeling and the associated problem of evaluating the models’ performance in a coherent way.

Section III treats in greater depth climate variability. We describe the most important modes of climate variability and provide an overview of the coexistence of several equilibria in the climate system, and of their dependence on parameter values. While the study of bifurcations and exchange of stability in the climate system goes back to the work of E. N. Lorenz, H. M. Stommel and G. Veronis in the 1960s (e.g., Dijkstra, 2013; Ghil, 2019; Ghil and Childress, 1987), a strikingly broadened interest in these matters has been stimulated by the borrowing from the social sciences of the term “tipping points” (Gladwell, 2000; Lenton *et al.*, 2008).

Proceeding beyond multiple equilibria, we show next how complex processes give rise to the system’s internal variability by successive instabilities setting in, competing and eventually leading to the quintessentially chaotic nature of the evolution of climate. This section concludes by addressing the need for using random processes to model the faster and smaller scales of motion in multiscale systems, and by discussing markovian and non-markovian approximations for the representation of the neglected degrees of freedom. We will also discuss top-
down vs. data-driven approaches.

Section IV delves into the analysis of climate response. The response to external forcing of a physico-chemical system out of equilibrium is the overarching concept we use in clarifying the mathematical and physical bases of climate change. We critically appraise climate models as numerical laboratories and review ways to test their skill at simulating past and present changes, as well as at predicting future ones. The classical concept of equilibrium climate sensitivity is critically presented first, and we discuss its merits and limitations.

We present next the key concepts and methods of nonautonomous and random dynamical systems, as a framework for the unified understanding of intrinsic climate variability and forced climate change, and emphasize the key role of pullback attractors in this framework. These concepts have been introduced only quite recently into the climate sciences, and we show how pullback attractors and the associated dynamical-systems machinery provide an excellent setting for studying the statistical mechanics of the climate system as an open system.

This system is subject to variations in the forcing and in its boundary conditions on all time scales. Such variations include, on different time scales, the incoming solar radiation, the position of the continents, and the sources of aerosols and greenhouse gases. We further introduce time-dependent invariant measures on a parameter-dependent pullback attractor, and the Wasserstein distance between such measures, as the main ingredients for a more geometrical treatment of climate sensitivity in the presence of large and sudden changes in the forcings.

We then outline, in the context of non-equilibrium statistical mechanics, Ruelle's response theory as an efficient and flexible tool for calculating climate response to small and moderate natural and anthropogenic forcings, and reconstruct the properties of the pullback attractor from a suitably defined reference background state. The response of a system near a tipping point is studied, and we emphasize the link between properties of the autocorrelation of the unperturbed system and its vicinity to the critical transition, along with their implications in terms of tell-tale properties of associated time series.

Section V is devoted to discussing multistability in the climate system, and the critical transitions that occur in the vicinity of tipping points in systems possessing multiple steady states. The corresponding methodology is then applied to the transitions between a fully frozen, so-called snowball state of our planet and its warmer states. These transitions have played a crucial role in modulating the appearance of complex life forms. We introduce the concept of an edge state, a dynamical object that has helped explain bistability in fluid mechanical systems, and argue that such states will also yield a more complete picture of tipping points in the climatic context. Finally, we present an example of a more exotic chaos-to-chaos critical transition that occurs in a delay

differential equation model for the Tropical Pacific.

In Section VI, we briefly summarize this review paper's main ideas and introduce complementary research lines that have not been discussed herein, as well as a couple of the many still open questions. Appendix A contains a table of scientific acronyms and one of institutional acronyms that are used throughout the paper.

We have started this section by characterizing the climate system and giving a broad-brush description of its behavior. But we have not defined the concept of climate as such, since we do not have as yet a consensual definition of what the climate, as opposed to weather, really is. An old adage states that “Climate is what you expect, weather is what you get.”

Clearly this implies that stochastic and ergodic approaches must play a role in disentangling the proper types of averaging on the multiple time and space scales involved. A fuller understanding of the climate system's behavior should eventually lead to a proper definition of climate. Mathematically rigorous work aimed at such a definition is being undertaken but is far from complete (e.g., F. Flandoli, Lectures at the Institut Henri Poincaré, Paris, October 2019; cf. http://users.dma.unipi.it/flandoli/IHP_contribution_Flandoli_Tonello_v3.pdf).

## II. A BRIEF INTRODUCTION TO CLIMATE DYNAMICS

### A. Climate Observations: Direct and Indirect

A fundamental difficulty in the climate sciences arises from humanity's insufficient ability to collect data of standardized quality, with sufficient spatial detail and of sufficient temporal coverage. Instrumental data sets have substantial issues of both synchronic and diachronic coherence. Moreover, such data sets only extend, at best, for about one-to-two centuries into the past. In this section, we cover first instrumental data sets and then so-called historical and proxy data sets, which use indirect evidence on the value of meteorological observables before the industrial era.

#### 1. Instrumental Data and Reanalyses

Since the establishment of the first meteorological stations in Europe and in North America in the 19th century, the extent and quality of the network of observations and the technology supporting the collection and storage of data have rapidly evolved. Still, at any given time, the spatial density of data changes dramatically across the globe, with much sparser observations over the oceans and over land areas characterized by low population density or a low degree of technological development (e.g., Ghil and Malanotte-Rizzoli, 1991, Fig. 1).
<figure><img src="image_6.png" /><figcaption>FIG. 5: Schematic diagram representing the evolution of the observing network for weather and climate data. The dotted vertical line corresponds to the International Geophysical Year (IGY). Courtesy of Dick Dee.</figcaption></figure>

Starting in the late 1960s, polar orbiting and geostationary satellites have led to a revolution in collecting weather, land surface and ocean surface data. Space-borne instruments are now remotely sensing many climatic variables from the most remote areas of the globe; for instance, they measure the overall intensity and spectral features of emitted infrared and reflected visible and ultraviolet radiation, and complex algorithms relate their raw measurements to the actual properties of the atmosphere, such as temperature and cloud cover.

Figure 5 represents schematically the evolution of the observational network for climatic data, while Fig. 6 portrays the instruments that today comprise the Global Observing System of the World Meteorological Organization (WMO), the United Nations agency that coordinates the collection and quality check of weather and climate data over the entire globe.

Since the early 20th century, the daily measurements have grown in number by many orders of magnitude and cover now more regularly the entire globe, even though large swaths of the Earth still feature relatively sparse observations. Figures 7 and 8 illustrate the coverage and variety of the observing systems available at present to individual researchers and practitioners, as well as to environmental and civil-protection agencies. Note that the so-called conventional network of ground-based weather stations and related observations has evolved since the Global Weather Experiment in the late 1970s, cf. Fig. 1 in Bengtsson et al. (1981) but only marginally so: it is the remote-sensing observations that have increased tremendously in number, variety and quality.

The number and quality of oceanographic observations was several orders of magnitude smaller than that of meteorological ones in the 1980s (Ghil and Malanotte-Rizzoli, 1991; Munk and Wunsch, 1982). Here also the advent of space-borne altimetry for sea surface heights (SSHs), scatterometry for surface winds, and other

<figure><img src="image_7.png" /><figcaption>FIG. 6: An illustration of the instruments and platforms that comprise the World Meteorological Organization's (WMO's) Global Observing System (GOS). From the COMET© website at http://meted.ucar.edu/ of the University Corporation for Atmospheric Research (UCAR), sponsored in part through a cooperative agreements with the National Oceanic and Atmospheric Administration (NOAA), U.S. Department of Commerce (DOC). ©1997–2016 UCAR; all rights reserved.</figcaption></figure>

remote-sensing methods has revolutionized the field (e.g., Robinson, 2010). This number, however, is still smaller by at least one order of magnitude than that of atmospheric observations since, as pointed out by Munk and Wunsch (1982), the oceans' interior is not permeable to exploration by electromagnetic waves. This is a fundamental barrier hindering our ability to directly observe the deep ocean.

Observational data for the atmosphere and oceans are at any rate sparse, irregular, and of different degrees of accuracy, while in many applications one has to obtain the best estimate, with known error bars, of the state of the atmosphere or oceans at a given time and with a given, uniform spatial resolution. More often than not, this estimate also needs to include meteorological, oceanographic or coupled-system variables, such as vertical wind velocity or surface heat fluxes, that can only be observed very poorly or not at all.

The very active field of data assimilation has developed to bridge the gap between the observations that are, typically, discrete in both time and space, and the continuum of the atmospheric and oceanic fields. Data assimilation — as distinct from polynomial interpolation, statistical regression or the inverse methods used in solid-earth geophysics — first arose in the late 1960s from the needs of numerical weather prediction (NWP), on the one hand, and the appearance of time-continuous data streams from satellites, on the other (Charney et al., 1969; Ghil et al., 1979). NWP is, essentially, an initial-value problem for the partial differential equations (PDEs) governing large-
Data Coverage, 10 April 2009

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

FIG. 7: Maps of point observations from the World Meteorological Organization's (WMO's) Global Observing System (GOS) on 10 April 2009: (a) synoptic weather station and ship reports; (b) upper-air station reports; (c) buoy observations; (d) aircraft wind and temperature; (e) wind profiler reports; (f) temperature and humidity profiles from Global Positioning System (GPS) radio occultation; and (g) observations from citizen weather observers. The tropics are the bright areas bordered by ±30° latitude. From the COMET© website at http://meted.ucar.edu/ of the University Corporation for Atmospheric Research (UCAR), sponsored in part through cooperative agreements with the National Oceanic and Atmospheric Administration (NOAA), U.S. Department of Commerce (DOC). ©1997-2016 UCAR; all rights reserved.

scale atmospheric flows (Richardson, 1922) that needed a
complete and accurate initial state every 12 or 24 hours.

Data assimilation combines partial and inaccurate ob-
servational data with a dynamic model, based on phys-
ical laws, that governs the evolution of the continuous

medium under study in order to provide the best esti-
mates of the state of the medium. This model also is
subject to errors, due to incomplete knowledge of the
smaller-scale processes, numerical discretization errors
and other factors. Given these two sources of infor-
<figure><img src="image_9.png" /><figcaption>FIG. 8: Geographic distribution of observing systems: (a) geostationary satellite observations; (b) and (c) polar-orbiting satellite soundings; (d) ocean surface scatterometer-derived winds; and (e, f) TRMM Microwave Imager (TMI) orbits. Each color represents the coverage of a single satellite. Observations in (b) and (c) represent vertical layers and area-averaged values. The tropics are marked by the lighter areas bordered by ±30° latitude. Same source as Fig. 7. ©1997-2016 UCAR; all rights reserved.</figcaption></figure>

mation, observational and physico-mathematical, there are three types of problems that can be formulated and solved, given measurements over a time interval {$t_0 \le t \le t_1$}: filtering, smoothing and prediction; see Fig. 9.

Filtering involves obtaining a best-possible estimate of the state $\mathbf{X}(t)$ at $t = t_1$, smoothing at all times $t_0 \le t \le t_1$, and prediction at times $t > t_1$. Filtering and prediction are typically used in NWP, and can be considered as the generation of short “video loops,” while smoothing is typically used in climate studies, and resembles the generation of long “feature movies” (Ghil and Malanotte-Rizzoli, 1991).

Figure 10 illustrates a so-called *forecast-assimilation cycle*, as used originally in NWP: at evenly spaced, pre-selected times {$t_k : k = 1, 2, \dots, K$}, one obtains an *analysis* of the state $\mathbf{X}(t_k)$ by combining the observations over some interval preceding the time $t_k$ with the forecast from the previous state $\mathbf{X}(t_{k-1})$ (Bengtsson *et al.*, 1981; Kalnay, 2003). Many variations on this relatively simple

scheme have been introduced in adapting it to oceano-graphic data (Ghil and Malanotte-Rizzoli, 1991, and references therein) and to space plasmas (e.g., Merkin *et al.*, 2016) or to actually using the time-continuous stream of remote-sensing data (Ghil *et al.*, 1979). Carrassi *et al.* (2018) have recently provided a comprehensive review of data assimilation methodology and applications.

Analyses are routinely used for numerical weather forecasts and take advantage of the continuous improvements of models and observations. But climate studies require data of consistent spatial resolution and accuracy over long time intervals, over which an operational NWP center might have changed its numerical model or its data assimilation scheme, as well as its raw data sources. To satisfy this need, several NWP centers have started in the 1990s to produce so-called *reanalyses* that use the archived data over multidecadal time intervals, typically since World War II, as well as the best model and data assimilation method available at the time of the reanal-
<figure><img src="image_10.png" /><figcaption>FIG. 9: Schematic diagram of filtering (F), smoothing (S) and prediction (P); green filled circles are observations. Based on Wiener (1949).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIG. 10: Schematic diagram of a forecast-assimilation cycle that is used for constructing the best estimates of the state of the atmosphere, oceans or both through the procedure of data assimilation: observational data are dynamically interpolated using the a meteorological, oceanographic or coupled model to yield the analysis products. The red arrow corresponds to a longer forecast, made only from time to time. Greater detail for the case of operational weather prediction appears in Ghil (1989, Fig. 1).</figcaption></figure>

ysis project. For obvious reasons of computational cost, reanalyses are often run at a lower spatial resolution than the latest version in operational use at that time.

Some leading examples of such diachronically coherent reanalyses for the atmosphere are those produced by the European Centre for Medium-range Weather Forecasts (ECMWF: Dee and Coauthors, 2011), the NCEP-NCAR Reanalysis produced in collaboration by the U.S. National Centers for Environmental Prediction (NCEP) and the National Center for Atmospheric Research (NCAR: Kistler and Coauthors, 2001), and the JRA-25 reanalysis produced by the Japan Meteorological Agency (JMA: Onogi and Coauthors, 2007). While these reanalyses agree pretty well for fields that are relatively well observed, such as the geopotential field (see Sect. II.C herein) over the continents of the Northern Hemisphere, substantial differences persist in their fields over the Southern Hemisphere or those that are poorly observed or not at all (Dell’Aquila et al., 2005, 2007; Kharin et al., 2005; Kim and Kim, 2013; Marques et al., 2010, 2009).

Recently, Compo et al. (2011) produced a centennial reanalysis from 1871 to the present by assimilating only surface pressure reports and using observed monthly sea-surface temperature and sea-ice distributions as boundary conditions, while Poli et al. (2016) provided a similar product for the time interval 1899–2010, where the surface pressure and the surface winds were assimilated. These enterprises are motivated by the need to provide a benchmark for testing the performance of climate models for the late 19th and the 20th century.

A similar need has arisen for the oceans: on the one hand, several much more sizeable data sources have become available through remote sensing, and led to detailed ocean modeling; on the other, the study of the coupled climate system requires a more uniform data set, albeit less accurate than for the atmosphere alone. Thus

the equivalent of a reanalysis for the oceans had to be produced, in spite of the fact that the equivalent of NWP for the oceans did not exist. A good example of a diachronically coherent data set for the global ocean is the Simple Ocean Data Assimilation (SODA: Carton and Giese, 2008). More recently, the community of ocean modelers and observationalists has delivered several ocean reanalyses able to provide a robust estimate of the state of the ocean (Balmaseda et al., 2015; Lee et al., 2009).

Finally, by relying on recent advances in numerical methods and in the increased availability of observational data, as well as of increased performance of computing and storage capabilities, coupled atmosphere-ocean data assimilation systems have been constructed (e.g., Penny and Hamill, 2017). Vannitsem and Lucarini (2016) provide a theoretical rationale for the need of coupled data assimilation schemes to be able to deal effectively with the climate system’s multiscale instabilities. These coupled systems play a key role in efforts to produce seamless weather, subseasonal-to-seasonal (S2S), seasonal, and interannual climate predictions (Palmer et al., 2008; Robertson and Vitart, 2018); they have already been used for constructing climate reanalyses (e.g., Karspeck et al., 2018; Laloyaux et al., 2018).

## 2. Proxy Data

As mentioned already repeatedly, and discussed in greater detail in Sect. II.B below, climate variability covers a vast range of time scales, and the information we can garner from the instrumental record is limited to the last century or two. Even so-called historical records only extend to the few millennia of a literate humanity (Lamb, 1972). In order to extend our reach beyond this eyeblink of the planet’s life, it is necessary to resort to indirect measures of past climatic conditions able to inform us on its state thousands or even millions of years ago.
Climate proxies are physical, chemical or biological characteristics of the past that have been preserved in various natural repositories and that can be correlated with the local or global state of the atmosphere, oceans or cryosphere at that time. Paleoclimatologists and geochemists currently take into consideration multiple proxy records, including coral records (Boiseau *et al.*, 1999; Karamperidou and Coauthors, 2015) and tree rings (Esper *et al.*, 2002) for the last few millennia, as well as marine-sediment (Duplessy and Shackleton, 1985; Taricco *et al.*, 2009) and ice-core (Andersen and Coauthors, 2004; Jouzel *et al.*, 1991) records for the last two million years of Earth history, the Quaternary. Glaciation cycles, i.e. an alternation of warmer and colder climatic episodes, dominated the latter era. The chemical and physical characteristics, along with the accumulation rate of the samples, require suitable calibration and dating, and are used thereafter to infer some of the properties of the climate of the past (e.g., Cronin, 2010; Ghil, 1994, and references therein).

Proxies differ enormously in terms of precision, uncertainties in the values and dating, and spatio-temporal extent, nor do they cover homogeneously the Earth. It is common practice to combine and cross-check many different sources of data to have a more precise picture of the past (Cronin, 2010; Imbrie and Imbrie, 1986). Recently, data assimilation methods have started to be applied to this problem as well, using simple models and addressing the dating uncertainties in particular (e.g., Roques *et al.*, 2014). Combining the instrumental and proxy data — with their very different characteristics of resolution and accuracy — is a complex, and sometimes controversial, exercise in applied statistics. An important example is that of estimating the globally averaged surface air temperature record well before the industrial era, cf. Fig. 4 and previous discussion on the hockey stick controversy in Sect. I.B.2.

## B. Climate Variability on Multiple Time Scales

The presence of multiple scales of motions — in space and in time — in the climate system can be summarized through so-called Stommel diagrams. Figure 11a presents the original Stommel (1963) diagram, in which a somewhat idealized spectral density associated with the oceans' variability was plotted in logarithmic space and time scales, while identifying characteristic oceanic phenomena whose variance exceeds the background level. Stommel diagrams describe the spatial-temporal variability in a climatic subdomain by associating different, phenomenologically well-defined dynamical features — such as cyclones and long waves in the atmosphere or meanders and eddies in the ocean — with specific ranges of scales; they emphasize relationships between spatial and temporal scales of motion. Usually, specific dynamical

features are associated with specific approximate balances governing the properties of the evolution equations of the geophysical fluids, cf. Sect. II.C.2.

In Figs. 11b and 11c, a qualitative Stommel diagram portrays today's estimates of the main range of spatial and temporal scales in which variability is observed for the oceans and the atmosphere, respectively. One immediately notices that larger spatial scales are typically associated with longer temporal scales, in both the atmosphere and oceans. The two plots clearly show that, for both geophysical fluids, a “diagonal” of high spectral density in the wavelength–frequency plane predominates. As the diagonal reaches the size of the planet in space, the variability can no longer maintain this proportionality of scales and it keeps increasing in time scales, which are not bounded, while the spatial ones are.

Both extratropical cyclones in the atmosphere and eddies in the ocean are due to baroclinic instability, but their characteristic spatial extent in the ocean is ten times smaller and their characteristic duration is one hundred times longer than in the atmosphere. In Fig. 11c, three important meteorological scales are explicitly mentioned: the microscale (small-scale turbulence), the mesoscale (e.g. thunderstorms and frontal structures), and the synoptic scale (e.g. extratropical cyclones).

Given the different dynamical variability ranges in space and time, different classes of numerical models, based on different dynamical balances, can simulate explicitly only one or a few such dynamical ranges. The standard way of modeling processes associated with a particular range of scales is to “freeze” processes on slower time scales or to prescribe their slow, quasi-adiabatic effect on the variability being modeled, in order to handle the processes that are too large or too slow in scale to be included in the model.

As for the faster processes, these are “parametrized,” i.e. one attempts to model their net effect on the variability of interest. Such parametrizations have been, until fairly recently, purely deterministic but have started over the last decade or so to be increasingly stochastic (Palmer and Williams, 2009, and references therein). We will discuss the mathematics behind parametrizations and provide a few examples in Sect. III.G.

To summarize, there are about 15 orders of magnitude in space and in time that are active in the climate system, from continental drift at millions of years to cloud processes at hours and shorter. The presence of such a wide range of scales in the system provides a formidable challenge for its direct numerical simulation. There is no numerical model that can include all the processes that are active on the various spatial and temporal scales and can run for $10^7$ simulated years. Ockham's razor and its successors, including Poincaré's parsimony principle (Poincaré, 1902), suggest that, if we had one, it wouldn't necessarily be such a good tool for developing scientific insight, rather than just a gigantic simulator not helping
scientists to distinguish the forest from the trees.

Still, there are increasing efforts for achieving “seem-
less prediction” across time and space scales (e.g., Palmer
et al., 2008; Robertson and Vitart, 2018). Merryfield
et al. (2020) provide a comprehensive review of the most
recent efforts for bridging the gap between S2S and
seasonal-to-decadal prediction.

1. A Survey of Climatic Time Scales

Combining proxy and instrumental data allows one to gather information not only on the mean state of the climate system, but also on its variability on many different scales. An artist’s rendering of climate variability on all time scales is provided in Fig. 12a. The first version of this figure was produced by Mitchell (1976) and many versions thereof have circulated since. The figure is meant to provide semi-quantitative information on the spectral power $S = S(\omega)$, where the angular frequency $\omega$ is $2\pi$ times the inverse of the oscillation period; $S(\omega)$ is supposed to give the amount of variability in a given frequency band for a generic climatic variable, although one has typically in mind the globally averaged surface air temperature. Unlike in the Stommel diagrams of Fig. 11, there is no information on the spatial scales of interest.

This power spectrum is not computed directly by spec-
tral analysis from a time series of a given climatic quan-
tity, such as (local or global) temperature; indeed, there
is no single time series that is 10⁷ years long and has a
sampling interval of hours, as the figure would suggest.
Figure 12a includes, instead, information obtained by an-
alyzing the spectral content of many different time series,
for example, the spectrum of the 335-year long record of
Central England temperatures in Fig. 12b. This time
series is the longest instrumentally measured record of
temperatures; see, though, Kondrashov et al. (2005a) for
Nile River water levels. Given the lack of earlier instru-
mental records, one can imagine, but cannot easily con-
firm, that the higher-frequency spectral features might
have changed, in amplitude, frequency or both, over the
course of climatic history.

With all due caution in its interpretation, Fig. 12a
reflects three types of variability: (i) sharp lines that
correspond to periodically forced variations, at one day
and one year; (ii) broader peaks that arise from internal
modes of variability; and (iii) a continuous portion of the
spectrum that reflects stochastically forced variations, as
well as deterministic chaos (Ghil, 2002).

Between the two sharp lines at 1 day and 1 year lies
the synoptic variability of mid-latitude weather systems,
concentrated at 3–7 days, as well as intraseasonal vari-
ability, i.e. variability that occurs on the time scale of
1–3 months. The latter is also called low-frequency at-
mospheric variability, a name that refers to the fact that
this variability has lower frequency, or longer periods,

<figure><img src="image_12.png" /><figcaption>(a) Stommel diagram of ocean variability</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(b) Chelton diagram of ocean variability</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(c) COMET diagram of atmospheric variability</figcaption></figure>

FIG. 11: Idealized wavelength-and-frequency power spectra for the climate system. (a) The original Stommel (1963) diagram representing the spectral density (vertical coordinate) of the oceans' variability as a function of the spatial and temporal scale. Reproduced with permission. (b) Diagram representing qualitatively the main features of the ocean variability; courtesy of D. Chelton. (c) Same as b), describing here the variability of the atmosphere, ©The COMET program.
<figure><img src="image_15.png" /><figcaption>(a) Composite power spectrum of climate</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(b) Spectrum of Central England Temperatures</figcaption></figure>

FIG. 12: Power spectra of climate variability across time scales. (a) An artist's rendering of the composite power spectrum of climate variability for a generic climatic variable, from hours to millions of years; it shows the amount of variance in each frequency range. (b) Spectrum of the Central England Temperature time series from 1650 to the present. Each peak in the spectrum is tentatively attributed to a physical mechanism; see Plaut et al. (1995) for details. From Ghil (2002). ©John Wiley and Sons, Ltd. Reproduced with permission.

than the life cycle of weather systems. Intraseasonal vari-
ability comprises phenomena such as the Madden–Julian
oscillation of winds and cloudiness in the tropics or the
alternation between episodes of zonal and blocked flow
in mid-latitudes (Ghil and Childress, 1987; Ghil et al.,
1991; Ghil and Mo, 1991; Haines, 1994; Molteni, 2002).

Immediately to the left of the seasonal cycle in Fig. 12a
lies interannual, i.e. year-to-year, variability. An impor-
tant component of this variability is the El Niño phe-
nomenon in the Tropical Pacific: once about every 2–7
years, the sea-surface temperatures (SSTs) in the East-
ern Tropical Pacific increase by one or more degrees
over a time interval of about one year. This SST varia-
tion is associated with changes in the trade winds over
the Tropical Pacific and in sea level pressures (Bjerknes,
1969; Philander, 1990); an East–West seesaw in the lat-

ter is called the Southern Oscillation. The combined El
Niño–Southern Oscillation (ENSO) phenomenon arises
through large-scale interaction between the equatorial
Pacific and the atmosphere above. Equatorial wave dy-
namics in the ocean plays a key role in setting ENSO's
time scale (Cane and Zebiak, 1985; Dijkstra and Burgers,
2002; Neelin et al., 1998, 1994).

The greatest excitement among climate scientists, as
well as the public, is more recently being generated by in-
terdecadal variability, i.e. climate variability on the time
scale of a few decades, the time scale of an individual
human's life cycle (Martinson et al., 1995). Figure 12b
represents a "blow-up" of the interannual-to-interdecadal
portion of Fig. 12a. The broad peaks are due to the cli-
mate system's internal processes: each spectral compo-
nent can be associated, at least tentatively, with a mode
of interannual or interdecadal variability (Plaut et al.,
1995). Thus the rightmost peak, with a period of 5.5
years, can be attributed to the remote effect of ENSO's
low-frequency mode (Ghil et al., 2002; Ghil and Robert-
son, 2000), while the 7.7-year peak captures a North
Atlantic mode of variability that arises from the Gulf
Stream's interannual cycle of meandering and intensifi-
cation (Dijkstra and Ghil, 2005, and references therein).
The two interdecadal peaks, near 14 and 25 years, are
also present in global records, instrumental as well as pa-
leoclimatic (Delworth and Mann, 2000; Ghil et al., 2002;
Kushnir, 1994; Mann et al., 1998; Moron et al., 1998).

Finally, the leftmost part of Fig. 12a represents pale-
oclimatic variability. The information summarized here
comes exclusively from proxy indicators of climate; see
Sect. II.A.2.

The presence of near-cyclicity is manifest, in this range,
in the broad peaks present in Fig. 12a between roughly
1 kyr and 1 Myr. The two peaks at about 20 kyr and
40 kyr reflect variations in Earth's orbit, while the dom-
inant peak at 100 kyr remains to be convincingly ex-
plained (Ghil, 1994; Gildor and Tziperman, 2001; Imbrie
and Imbrie, 1986). Quaternary glaciation cycles provide
a fertile testing ground for theories of climate variabil-
ity for two reasons: (i) they represent a wide range of
climatic conditions; and (ii) they are much better docu-
mented than earlier parts of paleoclimatic history.

Within these glaciation cycles, there is higher-
frequency variability prominent in North Atlantic paleo-
climatic records. These are the Heinrich (1988) events,
marked by a sediment layer that is rich in ice-rafted
debris, whose near-periodicity is of 6–7 kyr, and the
Dansgaard-Oeschger cycles (Dansgaard et al., 1993) that
provide the peak at around 1–2.5 kyr in Fig. 12a. Rapid
changes in temperature, of up to one half of the am-
plitude of a typical glacial–interglacial temperature dif-
ference, occurred during Heinrich events and somewhat
smaller ones over a Dansgaard-Oeschger cycle. Progres-
sive cooling through several of the latter cycles, followed
by an abrupt warming, defines a Bond cycle (Bond et al.,
1995). None of these higher-frequency phenomena can be
directly connected to orbital or other cyclic forcings.

In summary, climate variations range from the large-
amplitude climate excursions of past millennia to smaller-
amplitude fluctuations on shorter time scales. Several
spectral peaks of variability can be clearly related to forc-
ing mechanisms; others cannot. In fact, even if the exter-
nal forcing were constant in time — i.e., if no systematic
changes in insolation or atmospheric composition, such
as trace gas or aerosol concentrations, would occur —
the climate system would still display variability on many
time scales. This statement is clearly true for interannual
ENSO variability in the equatorial Pacific, as discussed
above. We will understand multiscale climate variability
better in Sects. III and IV, where we will look in greater
detail at natural climate variability and climate response
to forcings, respectively.

## 2. Atmospheric Variability in Mid-latitudes

Mid-latitude atmospheric variability during boreal
winter — when the winds are stronger and the variability
is larger in the Northern Hemisphere — has long been a
major focus of dynamic meteorology. The intent of this
section is to motivate the reader to appreciate the com-
plexity of large-scale atmospheric dynamics by focusing
on a relatively well understood aspect thereof. We shall
see that fairly diverse processes contribute to the spectral
features discussed in connection with Fig. 12.

The synoptic disturbances that are most closely associ-
ated with mid-latitude weather have characteristic time
scales of the order of 3–10 days, with a corresponding
spatial scale of the order of 1000–2000 km (Holton and
Hakim, 2013). They roughly correspond to the famil-
iar eastward-propagating cyclones and anticyclones, and
emerge as a result of the process of baroclinic instabil-
ity, which converts available potential energy of the zonal
flow into eddy kinetic energy.

This conversion is a crucial part of the Lorenz (1967,
1955) energy cycle and it occurs through the lowering of
the center of mass of the atmospheric system undergoing
an unstable development. Baroclinic instability (Char-
ney, 1947; Holton and Hakim, 2013; Vallis, 2006) is ac-
tive when the meridional temperature gradient or, equiv-
alently (see below), the vertical wind shear are strong
enough. These conditions are more readily verified in
the winter season, which features a large equator-to-
pole temperature difference and a strong mid-latitude jet
(Holton and Hakim, 2013; Speranza, 1983).

The space-time spectral analysis introduced by
Hayashi (1971) and refined by Pratt (1976) and by
Fraedrich and Bottger (1978) builds upon the idea of the
Stommel diagrams in Fig. 11. In addition, it provides
information about the direction and speed at which the
atmospheric eddies move and associates to each range

of spatial and temporal scales a corresponding weight in
terms of spectral power. This information may be ob-
tained, in the first instance, by Fourier analysis of a one-
dimensional spatial field and it allows one to reconstruct
the propagation of atmospheric waves. This analysis is
usually carried out in the so-called zonal, i.e. west-to-east
direction; see Sect. II.C below.

Next, one can compute the power spectrum in the fre-
quency domain for each spatial Fourier component, and
then average the results across consecutive winters to
derive a climatology of winter atmospheric waves. The
difficulty here lies in the fact that straightforward space-
time decomposition will not distinguish between standing
and traveling waves: a standing wave will give two spec-
tral peaks corresponding to waves that travel east- and
westward at the same speed and with the same phase.
This problem can only be circumvented by making as-
sumptions regarding a given wave's nature. For instance,
we may assume complete coherence between the east-
ward and westward components of standing waves and
attribute the incoherent part of the spectrum to actual
traveling waves.

Figure 13 shows the spectral properties of the win-
ter 500-hPa geopotential height field meridionally aver-
aged across the mid-latitudes of the northern hemisphere
(specifically, between 30° and 75° N) for the time interval
1957–2002. The properties of all the waves, as well as the
standing, eastward- and westward-traveling waves appear
in panels (a)–(d), respectively. As discussed later, the
500-hPa height field provides a synthetic yet comprehen-
sive picture of the atmosphere’s synoptic to large-scale
dynamics.

Figure 13c clearly shows that the eastward-
propagating waves are dominated by synoptic vari-
ability, concentrated between 3–12-day periods and
zonal wavenumbers 5–8; note that a single cyclone or
anticyclone counts for half-a-wavelength. In addition,
the slanting high-variability ridge in $\overline{H}_E$ indicates the
existence of a statistically defined dispersion relation
that relates frequency and wavenumber, in agreement
with the basic tenets of baroclinic instability theory
(Holton and Hakim, 2013).

When looking at the westward propagating variance in
Fig. 13d, one finds that the dominant portion of the vari-
ability is associated with low-frequency, planetary-scale
Rossby waves. Finally, Fig. 13b shows the contribution to
the variance given by standing waves, which correspond
to large-scale, geographically locked and persistent phe-
nomena like blocking events. Note that westward prop-
agating and stationary waves provide the lion's share of
the overall variability of the atmospheric field; see also
Kimoto and Ghil (1993a, Fig. 7).

The dynamics and energetics of planetary waves are
still under intensive scrutiny. Descriptions and expla-
nations of several highly nonlinear aspects thereof are
closely interwoven with those of blocking events, identi-
<figure><img src="image_17.png" /><figcaption>FIG. 13: Variance $\bar{H}$ of the winter (December-January-February) atmospheric fields in the mid-latitudes of the Northern Hemisphere. (a) Total variance $\bar{H}_T$; (b) variance associated with standing waves $\bar{H}_S$; (c) variance associated with eastward-propagating waves $\bar{H}_E$; and (d) variance associated with westward-propagating waves $\bar{H}_W$. Based on NCEP-NCAR reanalysis data (Kistler and Coauthors, 2001). See text for details. From Dell'Aquila et al. (2005), with permission from the American Geophysical Union.</figcaption></figure>

fied as persistent, large-scale deviations from the zonally symmetric general circulation (Benzi et al., 1986; Benzi and Speranza, 1989; Charney and DeVore, 1979; Ghil and Childress, 1987; Kimoto and Ghil, 1993a; Legras and Ghil, 1985). Weeks et al. (1997, Fig. 1) provide a fine example of contrast between a blocking event and the climatologically more prevalent zonal flow.

Persistent blocking events affect strongly the weather for up to a month over continental-size areas. Such persistence offers some hope for extended predictability of large-scale flows, and of the associated synoptic-scale weather, beyond the usual predictability of the latter, which is believed not to exceed 10–15 days (Lorenz, 1996, and references therein).

Today, the most highly resolved and physically detailed NWP models are reasonably good at predicting the persistence of a blocking event, once the model's initial state is within that event, but not at predicting the onset of such an event or its collapse (Ferranti et al., 2015). Likewise, the capability of necessarily lower-resolution climate models to simulate the spatio-temporal statistics of such events is far from perfect; in fact, relatively limited

improvement has been realized in the last two decades (Davini and D'Andrea, 2016). Weeks et al. (1997) reproduced successfully in the laboratory key features of the dynamics and statistics of blocking events.

Ghil and Robertson (2002) reviewed several schematic descriptions of the mid-latitude atmosphere's low-frequency variability (LFV) as jumping between a zonal regime and a blocked one or, more generally, a small number of such regimes. This coarse-graining of the LFV's phase space and Markov chain representation of the dynamics continues to inform current efforts at understanding what atmospheric phenomena can be predicted beyond 10–15 days and how. Recently, analyses based on dynamical systems theory have associated blocked flow configurations with higher instability of the atmosphere as a whole (Faranda et al., 2017; Lucarini and Gritsun, 2020; Schubert and Lucarini, 2016), as predicted by Legras and Ghil (1985).

The effect of global warming on the statistics of blocking events has become recently a matter of considerable controversy. The sharper increase of near-surface temperatures in the Arctic than near the Equator is
fairly well understood (e.g., Ghil, 1976, Fig. 7) and has been abundantly documented in recent observations (e.g., Walsh, 2014, Fig. 8). This decrease of pole-to-equator temperature difference $\Delta T$ is referred to as polar amplification of global warming.

Francis and Vavrus (2012) and Liu et al. (2012) have suggested that reduced $\Delta T$ slows down the prevailing westerlies and increases the north-south meandering of the subtropical jet stream, resulting in more frequent blocking events. This suggestion seems to agree with fairly straightforward arguments of several authors on the nature of blocking (Charney and DeVore, 1979; Ghil and Childress, 1987; Legras and Ghil, 1985), as illustrated, for instance, in Ghil and Robertson (2002, Fig. 2): The strength $\psi_A^*$ of the driving jet in the figure is proportional to $\Delta T$, in accordance with standard quasi-geostrophic flow theory,<fnref n="1" /> and lower jet speeds $\psi_A$ favor blocking. Ruti et al. (2006) have also presented observational evidence for a nonlinear relation between the strength of the subtropical jet and the probability of occurrence of blocking events, in agreement with the bent-resonance theory of mid-latitude LFV proposed by Benzi et al. (1986).

Considerable evidence against the apparently straightforward argument for global warming as the cause of an increase in midlatitude blocking has accumulated, too, from both observational and modeling studies (Barnes and Screen, 2015; Hassanzadeh et al., 2014, and references therein). The issue is far from settled, as are many questions about other regional effects of global warming.

## C. Basic Properties and Fundamental Equations

### 1. Governing Equations

The evolution of the atmosphere, the oceans, the soil, and the ice masses can be described by using the continuum approximation, in which these subsystems are represented by field variables that depend on three spatial dimensions and time. For each of the climatic subdomains, we consider the following field variables: the density $\rho$ and the heat capacity at constant volume $C$, with the specific expression for $\rho$ and $C$ defining the thermodynamics of the medium; the concentration of the chemical species $\{\xi_k: 1 \le k \le K\}$ contained in the medium, and present in different phases, e.g. the salt dissolved in the oceans or the water vapor in the atmosphere; the three components $\{v_i: 1 \le i \le 3\}$ of the velocity vector; the temperature $T$; the pressure $p$; the heating rate $J$; and the gravitational potential $\Phi$.

Note that, by making the thin-shell approximation $H/R \ll 1$, where $H$ is the vertical extent of the geo-

physical fluid and $R$ is the radius of the Earth, we can assume that the gravitational potential at the local sea level is zero and we can thus safely use the approximation $\Phi = gz$, where $\Phi$ is then called the geopotential, $g$ is the gravity at the surface of the Earth, and $z$ is the geometric height above sea level. Moreover, one has to take into account that the climate system is embedded in a noninertial frame of reference that rotates with an angular velocity $\Omega$ with components $\{\Omega_i: 1 \le i \le 3\}$.

The PDEs that govern the evolution of the field variables are based on the budget of mass, momentum and energy. When the fluid contains several chemical species, their separate budgets also have to be accounted for (Vallis, 2006). In order to have a complete picture of the Earth system, one should in principle also consider the evolution of biological species. Doing so, however, is well beyond our scope here, even though present Earth system models do attempt to represent biological processes, albeit in a simplified way; see later discussion in Sect. II.D.

The mass budget equations for the constituent species can be written as follows:

$$ \partial_t(\rho\xi_k) = -\partial_i(\rho\xi_k v_i) + D_{\xi_k} + L_{\xi_k} + S_{\xi_k}. \quad (1) $$

Here $\partial_t$ is the partial derivative in time and $\partial_i$ in the $x_i$-direction; $D_{\xi_k}$, $L_{\xi_k}$ and $S_{\xi_k}$ are the diffusion operator, phase changes, and local mass budget due to other chemical reactions that are associated with $k$.

The momentum budget's $i$th component is written:

$$ \partial_t(\rho v_i) = -\partial_j(\rho v_j v_i) - \partial_i p + \rho \partial_i \Phi - 2\rho \epsilon_{ijk} \Omega_j v_k + T_i + F_i. \quad (2) $$

Here the Levi-Civita antisymmetric tensor $\epsilon_{ijk}$ has been used to write the Coriolis force; $T_i$ indicates direct mechanical forcings, e.g. those resulting from luni-solar tidal effects; $F_i = -\partial_j\tau_{ij}$ corresponds to friction, with $\{\tau_{ij}\}$ the stress tensor; and summation over equal indices is used. Equation (2) is just a forced version of the momentum equation in the Navier-Stokes equations (NSEs), written in a rotating frame of reference.

A general state equation valid for both fluid envelopes of the Earth, i.e., the atmosphere and the oceans, is

$$ \rho = g(T, p, \xi_1, \dots, \xi_K). \quad (3) $$

As a first approximation, one can take $K=1$, where $\xi_1 = \xi$ is, respectively, moisture for the atmosphere and salinity for the oceans. In the interests of brevity, we will restrict ourselves here to the atmosphere.

In general, one can write the specific energy of the climate system as the sum of the specific internal, potential, and kinetic energies, taking into account also the contributions coming from chemical species in different phases. In order to get manageable formulas, some approximations are necessary (e.g., Peixoto and Oort, 1992).

<footnote n="1">¹ See Sects. II.C.2 and II.C.3 for geostrophic balance, quasi-geostrophy and their consequences.</footnote>
Neglecting reactions other than the phase changes between the liquid and gas phases of water, the expression of the specific energy in the atmosphere is

$$e = c_v T + \Phi + v_j v_j / 2 + Lq,$$

where $c_v$ is the specific heat at constant volume for the gaseous atmospheric mixture, $L$ is the latent heat of evaporation, and $q = \rho \xi$ is the specific humidity. In this formula, we neglect the heat content of the liquid and solid water and the heat associated with the phase transition between solid and liquid water. Instead, the approximate expression for the specific energy in the oceans is

$$e = c_W T + \Phi + v_j v_j / 2,$$

where $c_W$ is the specific heat at constant volume of water, while neglecting the effects of salinity and of pressure. Finally, for the specific energy of soil or ice, we can take $e = c_{\{S,I\}} T + \Phi$, respectively.

After some nontrivial calculations, one derives the following general equation for the local energy budget:

$$\partial_t(\rho e) = -\partial_j(\rho \varepsilon v_j) - \partial_j Q_j^{\mathrm{SW}} - \partial_j Q_j^{\mathrm{LW}} \\ - \partial_j J_j^{\mathrm{SH}} - \partial_j J_j^{\mathrm{LH}} - \partial_j(v_i \tau_{ij}) + v_i T_i, \quad (4)$$

where $e$ is the energy per unit mass and $\varepsilon = e + p/\rho$ is the enthalpy per unit mass. The energy source and sink terms can be written as the sum of the work done by the mechanical forcing $v_i T_i$ and of the respective divergences of the shortwave (solar) and longwave (terrestrial) components of the Poynting vector, $Q_j^{\mathrm{SW}}$ and $Q_j^{\mathrm{LW}}$; of the turbulent sensible and latent heat fluxes, $J_j^{\mathrm{SH}}$ and $J_j^{\mathrm{LH}}$; and of the scalar product of the velocity field with the stress tensor $v_i \tau_{ij}$.

Equation (4) is written in a conservative form, with the right-hand side containing only the sum of flux divergences, except for the last term, which is negligible. By taking suitable volume integrals of Eq. (4) and assuming steady-state conditions, one derives meridional enthalpy transports from the zonal budgets of energy fluxes (Ghil and Childress, 1987; Lucarini et al., 2014; Lucarini and Ragone, 2011; Peixoto and Oort, 1992); recall Figs. 2(a,b). The presence of inhomogeneous absorption of shortwave radiation due to the geometry of the Sun-Earth system and of the physico-chemical properties of the climatic subdomains determines the presence of non-equilibrium conditions for the climate system, as already discussed in Sect. I.A.

Various versions of Eqs. (1)–(4) have been studied for over a century, and well-established thermodynamical and chemical laws describe accurately the phase transitions and reactions of the climate system’s constituents. Finally, quantum mechanics allows one to calculate in detail the interaction between matter and radiation.

Still, despite the fact that climate dynamics is governed by well-known evolution equations, it is way beyond our

scientific abilities to gain a complete picture of the mathematical properties of the solutions. In fact, many fundamental questions are still open regarding the basic NSEs in a homogeneous fluid, without phase transitions and rotation (Temam, 1984, 1997). In particular, even for the basic NSEs, providing analytical, closed-form solutions is only possible in some highly simplified cases that are either linear or otherwise separable (Batchelor, 1974).

As already discussed above, Eqs. (1)–(4) can simulate a range of phenomena that spans many orders of magnitude in spatial and temporal scales. Hence, it is virtually impossible to construct numerical codes able to represents explicitly all the ongoing processes, at the needed resolution in space and time. It is thus necessary to parametrize the processes that occur at subgrid scales and cannot, therewith, be directly represented. Among the most important processes of this type are cloud-radiation interactions and turbulent diffusion. We will discuss the theoretical framework behind the formulation of parametrizations in Sect. III.G.

## 2. Approximate Balances and Filtering

Equations (1)–(4) are too general and contain too many wave propagation speeds and instabilities in order to properly focus on certain classes of phenomena of interest. For instance, the fluid’s compressibility only plays a role in the propagation of very fast acoustic waves, whose energy is negligible, compared to that of the winds and currents we are interested in here.

Therefore, starting with Charney (1947), more or less systematic approximations have been introduced to filter out undesirable waves from the equations of motion and to study particular classes of phenomena and processes. Depending on the climate subsystem and the range of scales of interest, different approximations can be adopted. For example, if one considers ice sheets and bedrock on the time scale of millennia, it is reasonable to assume that $v_j \simeq 0$, i.e., to remove the flow field from the evolution equation (Ghil, 1994; Saltzman, 2001, and references therein). Obviously, this approximation is not valid if one wants to describe explicitly the motion of the ice sheets on shorter time scales (Paterson, 1981).

More precisely, the filtering process consists in applying a set of mathematical approximations into the fundamental governing equations — usually considered to be a suitable generalization of the classical NSEs, like Eqs. (1)–(4). The purpose of this filtering is to exclude from the filtered system certain physical processes that are heuristically assumed to play only a minor role in its dynamics at the time and space scale under investigation. The magnitudes of all the terms in the governing equations for a particular type of motion are estimated by dimensional analysis (Barenblatt, 1987), and various classes of simplified equations have been derived
by considering distinct asymptotic regimes in the associated scales (Ghil and Childress, 1987; Klein, 2010; McWilliams, 2006; Pedlosky, 1987).

The approximations adopted rely on assuming that the continuous medium, whether fluid or solid, obeys suitably defined time-independent relations — or undergoes only small departures from such relations — and that these relations result from the balance of two or more dominating forces. Imposing such balances leads to reducing the number of independent field variables of the system that obey a set of evolution equations.

In meteorological terminology, a prognostic variable whose tendency, i.e. time derivative, appears in the full equations, may thus become a diagnostic variable, which only appears in a nondifferentiated form in the filtered equations. Additionally, the imposition of dynamical balances leads to removing specific wave motions from the range of allowed dynamical processes. Below we shall see two of the most important examples of filtering, which, among other things, are essential for the practical implementation of numerical models of the atmosphere and oceans (Cushman-Roisin and Beckers, 2011a; Holton and Hakim, 2013; Vallis, 2006).

a. Hydrostatic Balance. A classical example of filtering is the hydrostatic approximation. In a local Cartesian coordinate system (x, y, z), we define by $\hat{z}=\nabla\Phi/|\nabla\Phi|$ the direction perpendicular to a geopotential isosurface $\Phi \equiv \text{const}$. An obvious stationary solution of the NSEs is given in these coordinates by the time-independent hydrostatic balance equation:

$$ \rho_h g = -\partial_z p_h, \quad (5) $$

where the subscript 'h' denotes this particular solution.

On sufficiently large spatial and temporal scales, many geophysical flows — e.g., atmospheric and oceanic, as well as continental surface and ground water — are near hydrostatic equilibrium. In general, stable hydrostatic equilibrium is achieved when fluid with lower specific entropy lies below fluid with higher specific entropy; in this case, Eq. (5) is obeyed for all (x, y, z) within the domain occupied by the fluid.

When this condition is broken because of external forcing or internal heating, say, the stratification is not stable and the fluid readjusts so as to reestablish the hydrostatic equilibrium. This readjustment is achieved by vertical convective motions, in which available potential energy is transformed into kinetic energy responsible for vertical motions that can be locally much faster than the large-scale flow. Violations of hydrostatic balance thus exist only on short time and space scales.

A large class of models used in studying, simulating and attempting to predict climate variability and climate change is based on a particular simplification of the full set of Eqs. (1)–(4). This simplification leads to

the so-called primitive equations, which filter out nonhydrostatic motions. As a result, sound waves are excluded from the solutions, a fact that greatly facilitates the numerical implementation of such models (Cushman-Roisin and Beckers, 2011a; Washington and Parkinson, 2005).

The primitive equations are derived by assuming that the time-independent hydrostatic balance given in Eq. (5) does apply at all times, even when the flow is time dependent. One thus assumes that the vertical acceleration of the fluid is everywhere much smaller than gravity. Neglecting the vertical acceleration altogether leads to a monotonic relationship between the vertical coordinate $z$ and pressure $p$.<fnref n="2" />

One can hence replace $z$ by $p$ as a vertical coordinate, and rewrite the modified Eqs. (1)–(4) using a new coordinate system in which one also replaces, typically, the horizontal coordinates $(x, y)$ by the longitude $\lambda$ and the latitude $\phi$, respectively. In such a coordinate system, $v_3 \equiv \omega \equiv \partial_t p$ is the change of pressure with time, the three-dimensional (3-D) velocity field is nondivergent, $\partial_j v_j = 0$, and the density is automatically constant and set to 1 (Holton and Hakim, 2013; Vallis, 2006).

Despite the use of primitive equations, climate models aim at simulating a system that nonetheless features nonhydrostatic motions. Additionally, geophysical flows obeying primitive equations can lead to unstable vertical configurations of the fluids. Typical climate models are formulated so as to quickly eliminate local non-hydrostatic conditions and to parametrize the balanced condition that is thus recovered in terms of variables that are explicitly represented in the model, following the idea of convective adjustment originally proposed by Manabe and Strickler (1964). Finding optimal ways to parametrize the effect of convection on hydrostatic climate models has become a key research area in climate dynamics; see Emanuel (1994) for a classical treatise on the problem of atmospheric convection and Plant and Yano (2016) for a recent overview.

The primitive equations were at the core of both climate models and weather prediction models for three decades, from about the 1970s through the 1990s. Non-hydrostatic models — constructed through discretization of the full NSEs with forcing and rotation — have become increasingly available over the last decade, first as limited-area models and, more recently, even as global ones (e.g., Marshall et al., 1997).

These models require extremely high resolution in space and time and are computationally quite expensive; as they start to enter operational use, they require state-of-the-art hard- and software, in terms of both computing

<footnote n="2">² The validity of this approximation, even locally in space and at most times, explains Pascal’s experimental observation during the famed 1648 Puy-de-Dôme ascension that a barometer’s reading decreases monotonically with altitude.</footnote>
<figure><img src="image_18.png" /><figcaption>(a) Geostrophic balance in Northern Hemisphere</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(b) Synoptic conditions at 500 hPa over U.S.</figcaption></figure>

FIG. 14: Geostrophic balance. (a) Schematics of geostrophically balanced flow in the Northern Hemisphere: at a given z-level, the pressure gradient force (upward-pointing arrow), and Coriolis force (downward-pointing arrow) cancel out and the flow (horizontal arrow) is parallel to the isobars. (b) Synoptic conditions for the 500-hPa level over the United States at 0000 Greenwich mean time (GMT) on 12 July 2019. The dark gray lines indicate the isolines of geopotential height $z = \Phi/g$ (in units of 10 m), and the barbed blue arrows indicate the direction of the wind; the barbs indicate the wind speed, each short (long) barb corresponding to a speed of 5 (10) knots, with 2 knots = 1.03 $ms^{-1}$. Reproduced with permission from the NOAA-National Weather Service https://www.spc.noaa.gov/obswx/maps/.

and storage, as well as sophisticated post-processing of the output. Moreover, for any practical purpose, sufficiently detailed initial and boundary data are not available and methodological problems reappear in connection with formulating such data. Still, the use of nonhydrostatic models allows, in principle, for bypassing the key problem of parametrizing convection.

b. Geostrophic Balance. Another important example of filtering is the one associated with time-independent purely horizontal balanced flows, in which the horizon-

tal components of the pressure gradient force and the Coriolis force cancel out in the NSEs. Such flows are termed geostrophic; etymologically, this means Earth-turning. The nondimensional parameter that determines the validity of this approximation is the Rossby number $Ro = U/f_0L$, where U is a characteristic horizontal velocity, L a characteristic horizontal extent, and $f_0 = 2\Omega \sin \phi_0$; here $\Omega$ is the modulus of Earth's angular frequency of rotation around its axis, and $\phi_0$ is a reference latitude.

The geostrophic approximation provides a rather good diagnostic description of the flow when $Ro \ll 1$. Such small $Ro$-values prevail for the atmosphere on synoptic and planetary scales, such as those of extratropical weather systems — say poleward of 30°, where $f_0$ is large enough — and in the free atmosphere, i.e. above the planetary boundary layer, where frictional forces can be comparable to the Coriolis force. In the oceans, this approximation is extremely accurate everywhere except very close to the equator. The smallness of Rossby number in the oceans arises from oceanic currents being orders of magnitude slower than atmospheric winds, so that $Ro$ nears unity only where $f_0$ is extremely small.

Purely geostrophic flow fields are obtained via a zeroth-order expansion in $Ro$ of the NSEs. Fluid motion is introduced by considering small perturbations $\rho_g$ and $p_g$ that break the translation symmetry of the basic, purely hydrostatic density and pressure fields, $\rho_h$ and $p_h$, but preserve the geostrophic balance, i.e. one sets $\rho = \rho_h + \rho_g$ and $p = p_h + p_g$, respectively.

Letting $\hat{x}$ be locally the zonal direction and $\hat{y}$ be locally the meridional direction, one can write (Holton and Hakim, 2013) geostrophic balance at each height $z$ as:

$$ \rho_h u f_0 = -\partial_y p_g, \quad (6a) $$

$$ \rho_h v f_0 = \partial_x p_g, \quad (6b) $$

where the derivatives are taken at constant $z$, $(u, v) = (v_1, v_2)$ and, furthermore, the geostrophic perturbation itself is in hydrostatic equilibrium, i.e. $g\rho_g = -\partial_z p_g$. The combined hydrostatic and geostrophic balance constrain atmospheric and oceanic motions so that the velocity field is uniquely determined by the pressure field. In such doubly balanced flows, the winds or currents are parallel to the isobars at a given geopotential height $z$, as shown in Fig. 14a, rather than perpendicular to the isobars, as in nonrotating fluids. Moreover, the vertical component of the velocity field vanishes, $w = v_3 \equiv 0$.

Using the pressure coordinate system $(x, y, p)$ described above, where $p$ plays the role of the vertical coordinate, it is possible to express the geostrophic balance as follows (Holton and Hakim, 2013):

$$ f_0 u_g = -\partial_y \Phi, \quad (7a) $$

$$ f_0 v_g = \partial_x \Phi, \quad (7b) $$

where, this time, derivatives are taken at constant $p$. The vertical derivative $\partial_p \Phi$ can, furthermore, be expressed
in terms of $\rho_g$. For the atmosphere, one can make the simplifying assumption of dry conditions, for which the equation of state is described by the gas constant $R$ and the heat capacity at constant pressure $C_p$.

Assuming, moreover, a horizontally homogeneous background state with $T_h = T_h(p)$, introducing the potential temperature $\Theta_h = T_h(p/p_s)^{R/C_p}$, with $p_s = \text{const.}$ the reference surface pressure, and letting $\sigma := -RT_h d \ln \Theta_h/dp$, which defines the stratification of the background state, one obtains:

$$ \partial_p \Phi = -\frac{Rp}{\sigma} T_g. \qquad (8) $$

Thus, in the geostrophic approximation, the scalar field $\Phi$ provides complete diagnostic information on the state of the fluid: its horizontal derivatives give us the velocity field, while the vertical derivative gives us the geopotential, or mass field, via perturbations with respect to the background temperature field (Holton and Hakim, 2013).

Essentially, geophysical fluid dynamics (GFD) is the study of large-scale flows dominated by the combined hydrostatic and geostrophic balances. The former is a consequence of the shallowness of planetary flows, the latter of their rotation (Cushman-Roisin and Beckers, 2011a; Ghil and Childress, 1987; McWilliams, 2006; Pedlosky, 1987). The next subsection will introduce the more complex situation of finite-but-small perturbations from this combined balance.

Although modern computers allow a fully ageostrophic (beyond geostrophic approximation) description of most geophysical fluids, geostrophic balance remains a fundamental tool of theoretical research, as well as being used for practical applications in everyday weather prediction. Figure 14b illustrates a typical mid-latitude synoptic situation: at a given pressure level, the winds blow, to a very good approximation, parallel to the isolines of geopotential height, and the speed of the wind is higher where the gradient of geopotential height is stronger. An appendix in Ghil et al. (1979) describes the extent to which diagnostic relations based on the geostrophic approximation are still used in the analysis and prediction of mid-latitude weather, as simulated by advanced numerical models. Indeed, geostrophic approximation is implicit in the day-to-day reading of weather maps on synoptic-to-planetary scales.

### 3. Quasi-Geostrophy and Weather Forecasting

The diagnostic, i.e. time-independent nature of the geostrophic balance implies that the ageostrophic terms, although relatively small, are important for the time evolution of the flow. A planetary flow that evolves slowly in time compared to $1/f_0$ can be described using quasi-geostrophic theory, namely a perturbative theory that expands the NSEs in $R_0$ and truncates at first order.

The use of the quasi-geostrophic approximation effectively filters out solutions that correspond to higher-speed atmospheric or oceanic inertia-gravity waves; the latter, also called Poincaré waves, are gravity waves modified by the presence of rotation (Cushman-Roisin and Beckers, 2011a; Ghil and Childress, 1987; Holton and Hakim, 2013; McWilliams, 2006; Pedlosky, 1987; Vallis, 2006). This approximation breaks down, though, near frontal discontinuities and in other situations in which the ageostrophic component of the velocity field plays an important advective role, and it has to be improved upon by retaining higher-order terms.

In quasi-geostrophic theory, a fundamental role is played by the quasi-geostrophic potential vorticity<fnref n="3" />:

$$ q = \frac{1}{f_0} \Delta_2 \Phi + f + \partial_p \left( \frac{f_0}{\sigma} \partial_p \Phi \right). \qquad (9) $$

Here $\Delta_2$ is the horizontal Laplacian in $p$-coordinates, the first term $\zeta_g = \Delta_2 \Psi/f_0$ is the relative vorticity of the geostrophic flow, $f$ is the planetary vorticity, and the last term is the stretching vorticity; the static stability $\sigma$, which measures the stratification of the fluid was introduced in Eq. (8).

The potential vorticity $q$ thus combines dynamical and thermodynamical information on the fluid flow. In more abstract terms, one can write:

$$ q - f = \mathcal{L}\Phi, \qquad (10) $$

where $\mathcal{L}$ is a modified three-dimensional Laplacian operator. Using suitable boundary conditions, one can invert the expression above and derive the geopotential field from the vorticity field,

$$ \Phi = \mathcal{L}^{-1}(q - f). \qquad (11) $$

In the absence of forcing and dissipation, $q$ is conserved by geostrophic motions, so that:

$$ \partial_t q + J \left( \frac{\Phi}{f_0}, q \right) = 0, \qquad (12) $$

where $J(A, B) = \partial_x A \partial_y B - \partial_y A \partial_x B$ is the Jacobian operator, so that $J(\Phi/f_0, q) = (1/f_0)(\partial_x \Phi / \partial_y q - \partial_y \Phi / \partial_x q)$ describes vorticity advection by the geostrophic velocity field $(u_g, v_g)$, cf. Eq. (6). Note that, in the limit $\sigma \to \infty$, i.e. if one assumes infinitely stable stratification, the third term on the right-hand side of Eq. (9) drops out, and Eq. (12) describes the conservation of the absolute vorticity in a barotropic, two-dimensional flow.

Quasi-geostrophic theory, as introduced by Charney (1947), arguably provided a crucial advance for the understanding of the dynamics of planetary flows and provided the foundation for the successful start of NWP

<footnote n="3">³ More general formulations of potential vorticity are also quite important in GFD (e.g., Hoskins et al., 1985).</footnote>
(Charney *et al.*, 1950a). The filtering associated with this theory was, in particular, instrumental in eliminating the numerical instability that marred the pioneering weather forecast experiment of Richardson (1922).

In fact, L. F. Richardson had used a more accurate set of equations of motion for the atmosphere than the set of Eqs. (9)–(12). Somewhat counterintuitively, the quasi-geostrophic model eventually provided a more accurate prediction tool, because it represented more robustly the dynamical processes at the scales of interest. The inconclusive result of Richardson’s one-step, 6-hour numerical experiment was received by the meteorological community at the time as proof that NWP, as proposed by Bjerknes (1904), was not possible.

We outline herewith how the quasi-geostrophic potential-vorticity equations can be used to perform stable numerical weather forecasts. Assume that, at time $t = t_0$, we have information on the field $\Phi_0 = \Phi(t_0)$.

* Step 0: Using Eq. (10), one computes the corresponding potential vorticity field $q_0 = q(t_0)$;

* Step A: Equation (12) is used next to predict the potential vorticity field at time $t_1$ as follows: $q_1 = q(t_1) = q(t_0) + J(\Phi_0/f_0, q_0)\Delta t$, where $\Delta t = t_1 - t_0$;

* Step B: Equation (11) is used now to infer $\Phi_1$;

* Step C: Go to Step A and predict $q_2$, and so on.

Three decades after Richardson (1922), J. G. Charney, R. Fjørtoft and J. von Neumann performed the first successful numerical weather forecasts on the ENIAC computer in Princeton by using essentially the procedure detailed above for the simpler case of barotropic, two-dimensional flow (Charney *et al.*, 1950a). In this case, the standard vorticity is given by the sum of relative and planetary vorticity and it plays a role analogous to potential vorticity in quasi-geostrophic three-dimensional flow. The procedure was adapted to describe the evolution of the 500-hPa field, which deviates the least from the behavior of an idealized atmosphere of homogeneous density with no vertical velocities. The simulation was performed on a limited domain that covered North America and the adjacent ocean areas.

The authors investigated the issues associated with horizontal boundary conditions, as well as with the numerical stability of the integration. Their success (Charney *et al.*, 1950a) paved the way for the theory’s successive applications to physical oceanography and climate dynamics as a whole (Dijkstra, 2005; Ghil and Childress, 1987; Pedlosky, 1987). The saga of this scientific and technological breakthrough was told many times; see Lynch (2008) for a good recent account.

## D. Climate Prediction and Climate Model Performance

A key area of interest in the climate sciences is the development and testing of numerical models used for simulating the past, present, and future of the climate system. As discussed later in the paper, climate models differ enormously in terms of scientific scope, computational cost, and flexibility, so that one has to consider systematically a hierarchy of climate models, rather than one model that could incorporate all subsystems, processes, and scales of motion (Ghil, 2001; Ghil and Robertson, 2000; Held, 2005; Lucarini, 2002).

Figure 15a shows the so-called Bretherton horrendogram (NAC, 1986), which displays the full range of subsystems one needs to deal with when addressing the Earth system as a whole, along with some of the many interactions among these subsystems. The climate modeling community has slowly caught up with the complexities illustrated by F. Bretherton and his colleagues in this diagram, three decades ago.

Models of different levels of complexity and detail are suited for addressing different kinds of questions, according, in particular, to the main spatial and temporal scales of interest. At the top of the hierarchy, one finds global climate models, known in the 1970s and ‘80s as general circulation models. Both these phrases share the acronym GCMs, with the change in name reflecting a change in emphasis from understanding the planetary-scale circulation — first of the atmosphere, then of the oceans — to simulating and predicting the global climate. GCMs aim to represent, at the highest computationally achievable resolution, the largest number of physical, chemical, and biological processes of the Earth system; see an early overview in Randall (2000).

The atmosphere and oceans, whose modeling relies on the equations discussed in Sect. II.C.1 above, are still at the core of the Earth system models being developed today. Following improvements in basic scientific knowledge, as well as in computing and data storage capabilities, the latter models include an increasing number of physical, chemical and biological processes. They also rely on a much higher spatial and temporal resolution of the fields of interest.

Figure 15b illustrates schematically the process of model development across the last three decades, from the first to the fourth assessment report: FAR–SAR–TAR–AR4. The graph does not cover the development of the last 5–10 years, which has largely dealt with the inclusion of eco-biological modules and the above-mentioned nonhydrostatic effects. Currently available state-of-the-art models result from the cumulative efforts of generations of climate scientists and coding experts, with hundreds of individual researchers contributing to different parts of the code. As in the early stages of post-WWII NWP, climate modelling exercises now are some of the heaviest users of civilian High Performance Computing.
<figure><img src="image_20.png" /><figcaption>(a) Bretherton horrendogram</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(b) Climate model evolution</figcaption></figure>

FIG. 15: The Earth system, its components and its modelling. (a) The NAC (1986) horrendogram that illustrates the main components of the Earth system and the interactions amongst them. (b) Evolution of climate models across the first four IPCC assessment reports, ranging from the early 1990s to the mid-2000s (IPCC, 2007).
# 1. Predicting the State of the System

In spite of its remarkable progress, climate modeling and prediction face several kinds of uncertainties: First, uncertainties in predicting the state of the system at a certain lead time, given the uncertain knowledge of its state at the present time. Geophysical flows are typically chaotic, as are other processes in the system. Hence, as discussed in the next two sections, the climate system depends sensitively on its initial data, as suggested already by Poincaré (1908) and recognized more fully by Lorenz (1963). Following up on terminology introduced by Lorenz (1976), these are *uncertainties of the first kind*, namely small errors in the initial data that can lead at later times to large errors in the flow pattern.

Assume that two model runs start from nearby initial states. In a nearly linear regime, in which the phase space distance between the system's orbits is small, their divergence rate can be studied by considering the spectrum of Lyapunov exponents and, in particular, the algebraically largest ones. A chaotic system has at least one positive Lyapunov exponent, and physical instabilities that act on distinct spatial and temporal scales are related to distinct positive exponents (Eckmann and Ruelle, 1985). These instabilities can be described by using the formalism of covariant Lyapunov vectors (Ginelli et al., 2007).

Addressing the uncertainties of the first kind and providing good estimates for the future state of the system is the classical goal of NWP. After the initially exponential increase of errors, nonlinear effects kick in, and the orbits in the ensemble populate the attractor of the system, so that any predictive skill is lost. In such an ensemble, several simulations using the same model are started with slightly perturbed initial states, and the ensemble of orbits produced is used to provide a probabilistic estimate of how the system will actually evolve (Palmer, 2017).

Figure 16a illustrates the main phases of error growth from an ensemble of initial states in the Lorenz (1963) model and Fig. 16b outlines how an ensemble forecast system actually works. Obvious limitations are related to the computational difficulties of running a sufficient number of ensemble members in order to obtain an accurate estimate, as well as with the fact that an initial ellipsoid of states tends to become flattened out in time, making the error estimate more laborious. Thus, care has to be exercised in the choice of the initial states, which is done in practice by taking advantage of the reanalyses described earlier in Sect. II.A.1.

Recently, similar methods are being used on longer time scales to perform experimental climate predictions, seen as an initial value problem (Eade et al., 2014). These predictions are performed on seasonal (Doblas-Reyes et al., 2013) to decadal (Meehl et al., 2014) time scales; they aim not only to reach a better understanding of the multiscale dynamics of the climate system, but also to help achieve the many potential benefits of

skillful medium-to-long-term predictions. One interesting example of such an application is the wind energy market (Torralba et al., 2017).

The World Climate Research Programme (WCRP) has created an encompassing project focusing on the challenges of S2S prediction; see https://www.wcrp-climate.org/s2s-overview and Robertson and Vitart (2018). We will delve more specifically into extended-range forecasting in Sect. III.E.1 and see how its performance depends on a model's ability to capture the natural modes of variability of the climate system on different time scales.

# 2. Predicting the System's Statistical Properties

A key goal of climate modeling is to capture the system's statistical properties — i.e., its mean state and its variability — and its response to forcings of different nature. These problems will be treated thoroughly in the next two sections but are rendered quite difficult by a second set of uncertainties. Uncertainties in model formulation, as well as an unavoidably limited knowledge of the external forcings are referred to as *uncertainties of the second kind* (Lorenz, 1967, 1976; Peixoto and Oort, 1992) and limit intrinsically the possibility of providing realistic simulations of the statistical properties of the climate system; they affect, in particular, severely the modeling of abrupt climate changes and of the processes that may lead to such changes.

These deficiencies are related to uncertainties in many key parameters of the climate system, as well as to the fact that each model may represent incorrectly certain processes that are relevant on the temporal and spatial scales of interest or that it may miss them altogether. These types of uncertainty are termed parametric and structural uncertainty, respectively (Lucarini, 2013).

A growing number of comprehensive climate models are available to the international scientific community for studying the properties of the climate system and predicting climate change; currently there are about 50 such models (IPCC, 2014a). Still, many of these models have in common a substantial part of their numerical code, as they originate from a relatively small number of models, atmospheric and oceanic, that were originally developed in the 1960s and 1970s. Figure 17 illustrates this family tree for the case of atmospheric GCMs. Hence the widely noted fact that the climate simulations produced by this large number of models often fall into classes that bear marked similarities and do not necessarily resemble to a desired extent the observed climate evolution over the last century or so.

IPCC's evolution over three decades has led to a coordination and restructuring of modeling activities around the world. In order to improve comparisons among distinct models and the replicability of investigations aimed
<figure><img src="image_22.png" /><figcaption>(a) Evolution in phase space</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(b) Evolution in time</figcaption></figure>

FIG. 16: Schematic diagram of the evolution of an ensemble of initial states in a chaotic system. (a) A small ball of initial states in the Lorenz (1963) convection model evolves initially in phase space according to the stretching and contracting directions associated with positive and negative Lyapunov exponents, until nonlinear effects become important and the set of initial states propagated by the system's vector field populates its attractor. Reproduced with permission from Slingo and Palmer (2011). (b) Evolution in time of such an initial ball, in a generic chaotic system. Reproduced with permission from Swinbank et al. (2016).

at climate change, the Program for Climate Model Diagnostics and Intercomparison (PCMDI), through its climate model intercomparison projects (CMIPs), defines standards for the modeling exercises to be performed by research groups that wish to participate in a given assessment report (AR) and provide projections of future climate change; see http://cmip-pcmdi.llnl.

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

FIG. 17: Family tree of the main atmospheric general circulation models. The tree clearly shows that many of today's state-of-the-art models share a — smaller or larger — portion of their genes, i.e., of the basic ideas and parameter values that went into the coding. Reproduced with permission from Edwards (2010).

gov/cmip5/data_portal.html for CMIP5 and https://www.wcrp-climate.org/wgcm-cmip/wgcm-cmip6 for CMIP6, respectively.

The PCMDI's CMIPs have also supported a single website that gathers climate model outputs contributing to the IPCC-initiated activities, thereby providing a unique opportunity to evaluate the state-of-the-art capabilities of climate models in simulating the climate system's past, present and future behavior. A typical IPCC-style package of experiments includes simulating the climate system under various conditions, such as:

* a reference state, e.g., a statistically stationary preindustrial state with fixed parameters;

* industrial era and present-day conditions, including known natural and anthropogenic forcings;

* future climate projections, performed by using a set of future scenarios of greenhouse gas and aerosol emissions and land-use change with some degree of realism, as well as idealized ones (e.g. instantaneous doubling of CO₂ concentration).

The latter changes of greenhouse gas, aerosol concentrations and land use follow a prescribed evolution over a given time window, and are then fixed at a certain value
to observe the relaxation of the system to a new stationary state. Each such evolution of the forcing is defined a “scenario” in ARs 1–4 (IPCC, 2001, 2007), and a Representative Concentration Pathway (RCP) in AR5 (IPCC, 2014a). Note that each AR has involved an increasing number of models, people — reaching many hundreds by the latest AR — and hence acronyms.

Each scenario or RCP is a representation of the expected greenhouse gas and aerosol concentrations resulting from a specific path of industrialization and change in land use, as provided by Working Groups II and III to Working Group I; see Fig. 18 for AR5’s simulations of the 20th and 21st century with two extreme RCPs. At the same time, the attribution of unusual climatic conditions to specific climate forcing is far from being a trivial matter (Allen, 2003; Hannart et al., 2016a).

In the most recent CMIP exercise, CMIP6, the scope of model intercomparisons has grown beyond just testing future climate response to specific forcing scenarios or RCPs. In order to test more comprehensively the performance of climate models, various standardized intercomparison projects focus on the analysis of specific subsystems, processes or time scales; see https://search.es-doc.org. Eyring et al. (2016a) provide a useful summary of the strategy for PCMDI’s CMIP6 project and of the scientific questions it addresses.

### 3. Metrics for Model Validation

The standardization of climate model outputs promoted by the CMIPs has also helped address a third major uncertainty in climate modeling: what are the best metrics for analyzing a model’s outputs and evaluating its skill? Note that, in this context, the term metric does not have its usual mathematical meaning of a distance in function space but refers instead to a statistical estimator of model performance, whether quadratic or not.

The validation, or auditing — i.e., the overall evaluation of the faithfulness — of a set of climate models, is a delicate operation, which can be decomposed into two distinct but related procedures. The first one is model intercomparison, which aims to assess the consistency of different models in the simulation of certain physical phenomena in a certain time frame. The second procedure is model verification, whose goal is to compare a model’s outputs with corresponding observed or reconstructed scalar quantities or fields (Lucarini, 2008b).

In principle, there are numerous ways to construct a metric by taking any reasonable function of a climate model’s variables. Nevertheless, even if several observables are mathematically well defined, their physical relevance and robustness can differ widely. There is no a priori valid criterion for selecting a good climatic observable, even though taking into account basic physical properties of the climate system can provide useful guidance.

There is thus no unique recipe for testing climate models, a situation that is in stark contrast with more traditional areas of physics. For instance, in high-energy physics, the variables mass, transition probability or cross-section are suggested by the very equations that one tries to solve or to study experimentally.

In the absence of dissipation and of sources and sinks, certain scalars—such as the atmosphere’s or the oceans’ total mass, energy and momentum—are also conserved by the fundamental equations of Sect. II.C.1. Local quantities, though, like total rain over India during the summer monsoon or the intensity of the subtropical jet over North America during a given winter month, often matter in evaluating a climate model’s skill: current models still have substantial difficulties in simulating the statistics of major regional processes, such as ENSO in the Tropical Pacific (Bellenger et al., 2014; Lu et al., 2018), the Indian Monsoon (Boos and Hurley, 2013; Hasson et al., 2013; Turner and Annamalai, 2012), and mid-latitude LFV associated with blocking (Davini and D’Andrea, 2016; Woollings et al., 2018).

From the end user’s point of view, it is important to check how realistic the modeled fields of practical interest are. But if the aim is to define strategies for the radical improvement of model performance — beyond incremental advances often obtained at the price of large increases in computer power — it is crucial to fully understand the differences among models in their representation of the climate system, and to decide whether specific physical processes are correctly simulated by a given model.

Additional issues, practical as well as epistemological, emerge when we consider the actual process of comparing theoretical and numerical investigations with observational or reanalysis data. Model results and approximate theories can often be tested only against observational data from a sufficiently long past history, which may pose problems of both accuracy and coverage, as mentioned in Sect. II.A.1.

To summarize, difficulties emerge in evaluating climate model performance because (i) we always have to deal with three different kinds of attractors: the attractor of the real climate system, its reconstruction from observations, and the attractors of the climate models; and (ii) because of the high dimensionality of both the phase space and the parameter space of these attractors.

In order to address these issues, multi-variable metrics are currently used to try to assess the skill of available climate models. Thus, Fig. 19 shows a diagram describing the performance of the 42 models participating in PCMDI’s CMIP5. Even superficial analysis of the diagram indicates that no model is the best for all the variables under consideration. An improved assessment package is given by PCMDI’s Metrics Package (Gleckler et al., 2016) and by the ESMValTool package (Eyring et al., 2016b). Recently, Lembo et al. (2019) released TheDiato, a flexible diagnostic tool able to evaluate com-
<figure><img src="image_25.png" /></figure>

FIG. 18: State-of-the-art climate model outputs for various climate change scenarios. Left panel: Change in the globally averaged surface temperature as simulated by climate models included in IPCC (2014a). Vertical bands indicate the range of model outputs, and the colors correspond to different Representative Concentration Pathways (RCPs). Right panel: Spatial patterns of temperature change — i.e., 2081–2100 average with respect to the present — for the two most extreme RCPs. Reproduced with permission from (IPCC, 2014a)

previously the energy, entropy and water budgets, and their transports for climate models. This package will become part of the second generation of the ESMValTool package (Eyring et al., 2019).

Finally, in order to describe synthetically and comprehensively the outputs of a growing number of climate models, it has become common to consider multi-model ensembles and focus the attention on the ensemble mean and the ensemble spread. Mean and spread have been taken as the, possibly weighted, first two moments of the models' outputs for the metric under study (Tebaldi and Knutti, 2007); see Fig. 18.

This approach merges information from different attractors and the resulting statistical estimators cannot be interpreted in the standard way, with the mean approximating the true field and the standard deviation describing its uncertainty. Such a naive interpretation relies on an implicit assumption that the set is a probabilistic ensemble formed by equivalent realizations of a given process, and that the underlying probability distribution is unimodal; see Parker (2010) for a broader epistemological discussion of the issues.

While the models in such an “ensemble of opportunity” may be related to each other, as shown in Fig. 17, they are by no means drawn from the same distribution. A number of alternative approaches for uncertainty quantification in climate modeling have been proposed but they go beyond the scope of the present review.

### III. CLIMATE VARIABILITY AND THE MODELING HIERARCHY

#### A. Radiation Balance and Energy-Balance Models

The concepts and methods of the theory of deterministic dynamical systems (Andronov and Pontryagin, 1937; Arnold, 1983; Guckenheimer and Holmes, 1983) have been applied first to simple models of atmospheric and oceanic flows, starting about half-a-century ago (Lorenz, 1963; Stommel, 1961; Veronis, 1963). More powerful computers now allow their application to fairly realistic and detailed models of the atmosphere, ocean, and the coupled atmosphere-ocean system. We start herewith by presenting such a hierarchy of models.

This presentation is interwoven with that of the successive bifurcations that lead from simple to more complex solution behavior for each climate model. Useful tools for comparing model behavior across the hierarchy and with observations are provided by ergodic theory (Eckmann and Ruelle, 1985; Ghil et al., 2008). Among these tools, advanced methods for the analysis and prediction of uni- and multivariate time series play an important role (Ghil et al., 2002, and references therein).

The concept of a modeling hierarchy in climate dynamics was introduced by Schneider and Dickinson (1974). Several authors have discussed lately in greater detail the role of such a hierarchy in the understanding and prediction of climate variability (Dijkstra and Ghil, 2005; Ghil
<figure><img src="image_26.png" /><figcaption>FIG. 19: Example of model performance evaluation. This diagram shows how 42 models participating in the Fifth Climate Model Intercomparison Project (CMIP5) fare in terms of representation of the seasonal cycle during the time interval 1980–2005 for 13 climate variables. Values are normalized and perfect agreement with the observations is given by 0. See text for details. Reproduced with permission from IPCC (2014a).</figcaption></figure>

and Robertson, 2000; Held, 2005; Lucarini, 2002; Saltzman, 2001). At present, the best-developed hierarchy is for atmospheric models. These models were originally developed for weather simulation and prediction on the time scale of hours to days. Currently they serve — in a stand-alone mode or coupled to oceanic and other models — to address climate variability on all time scales.

The first rung of the modeling hierarchy for the atmosphere is formed by zero-dimensional (0-D) models, where the number of dimensions, from zero to three, refers to the number of independent space variables used to describe the model domain, i.e. to physical-space dimensions. Such 0-D models essentially attempt to follow the evolution of the globally averaged air temperature at surface as a result of changes in global radiative balance:

$$c \frac{d\bar{T}}{dt} = R_i - R_o, \qquad (13a)$$

$$R_i = \mu Q_0 \{1 - \alpha(\bar{T})\}, \qquad (13b)$$

$$R_o = \sigma m(\bar{T})(\bar{T})^4. \qquad (13c)$$

Here $R_i$ and $R_o$ are incoming solar radiation and outgoing terrestrial radiation. The heat capacity $c$ is that of the global atmosphere, plus that of the global ocean or some fraction thereof, depending on the time scale of interest: one might only include in $c$ the ocean mixed layer

when interested in subannual time scales but the entire ocean when studying paleoclimate (e.g., Saltzman, 2001).

The rate of change of $\bar{T}$ with time $t$ is given by $d\bar{T}/dt$, while $Q_0$ is the solar radiation received at the top of the atmosphere, also called the “solar constant,” $\sigma$ is the Stefan-Boltzmann constant, and $\mu$ is an insolation parameter equal to unity for present-day conditions. To have a closed, self-consistent model, the planetary reflectivity or albedo $\alpha$ and grayness factor $m$ have to be expressed as functions of $\bar{T}$; $m = 1$ for a perfectly black body and $0 < m < 1$ for a grey body like planet Earth.

There are two kinds of one-dimensional (1-D) atmospheric models, for which the single spatial variable is latitude or height, respectively. The former are so-called energy-balance models (EBMs), which consider the generalization of the model (2.1) for the evolution of surface-air temperature $T = T(x, t)$, say,

$$c(x) \frac{\partial T}{\partial t} = R_i - R_o + D. \qquad (14)$$

The terms $R_i$ and $R_o$ are similar to those given in Eqs. (13b) and (13c) for the 0-D case above, but can now be functions of the meridional coordinate $x$ — latitude, co-latitude, or sine of latitude — as well as of time $t$ and temperature $T$. The horizontal heat-flux term $D$
describes the convergence of the heat transport across latitude belts; it typically contains first and second partial derivatives of T with respect to x, while c(x) represents the system's space-dependent heat capacity.

Thus Eq. (14) corresponds physically to a nonlinear heat or reaction-diffusion equation, and mathematically to a nonlinear parabolic PDE. Hence the rate of change of local temperature T with respect to time also becomes a partial derivative, $\partial T(x, t)/\partial t$. Two such models were introduced independently by Budyko (1969) in the then Soviet Union and by Sellers (1969) in the United States.

<figure><img src="image_27.png" /><figcaption>FIG. 20: Scalar double-well potential function $V(T)$; the warm and the “deep-freeze” or snowball states correspond to the system's two stable fixed points, separated by an unstable one.</figcaption></figure>

The first striking results of theoretical climate dynamics were obtained in showing that Eq. (14) could have two stable steady-state solutions, depending on the value of the insolation parameter $\mu$, cf. Eq. (13b)$^4$.

This multiplicity of stable steady states, or physically possible stationary climates of our planet, can be explained, in its simplest form, already by the 0-D model of Eq. (13). Note that the time derivative of the global temperature $\bar{T}$ in Eq. (13a) can be written as minus the derivative of a potential $V(\bar{T}) = -\int \{R_i(\bar{T}) - R_o(\bar{T})\}d\bar{T}$, viz. $d\bar{T}/dt = -dV(\bar{T})/d\bar{T}$. In the case of bistability, the two local minima of V correspond to the stable solutions, or fixed points, and the local maximum of V corresponds to the unstable solution; see the sketch in Fig. 20.

The physical explanation lies in the fact that — for a fairly broad range of $\mu$-values around $\mu = 1.0$ — the curves for $R_i$ and $R_o$ as a function of $\bar{T}$ intersect in three

points. One of these points corresponds to the present climate (highest $\bar{T}$-value), and another one to an ice-covered planet (lowest $\bar{T}$-value); both of these are stable, while the third one, i.e. the intermediate $\bar{T}$-value, is unstable.

To obtain this result, it suffices to make two assumptions: (i) that $\alpha = \alpha(\bar{T})$ in Eq. (13b) is a piecewise-linear function of $\bar{T}$ — or, more generally, a monotonically increasing one with a single inflection point — with high albedo at low temperature, due to the presence of snow and ice, and low albedo at high $\bar{T}$, due to their absence; and (ii) that $m = m(\bar{T})$ in Eq. (13c) is a smooth, increasing function of $\bar{T}$ that captures in its simplest from the “greenhouse effect” of trace gases and water vapor.

The EBM modelers (Ghil, 1976; Held and Suarez, 1974; North, 1975) called the ice-covered state a “deep freeze.” The possibility of such a state in Earth history was met with considerable incredulity by much of the climate 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. \qquad (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).

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. \qquad (16) $$

<footnote>$^4$After the publication of the Budyko (1969) and Sellers (1969) models, it became apparent that massive nuclear explosions could induce a nuclear winter: namely, reducing the incoming solar radiation by a dramatic increase in atmospheric particulate matter could potentially trigger an even greater disaster for life on Earth than nuclear war itself (e.g., Turco et al., 1983). Studies to this effect have been influential in reducing the size of the nuclear arsenals at the end of the Cold War. Strikingly, the two contributions came almost simultaneously from scientists belonging to the Cold War's two opposing geopolitical blocks.</footnote>
<figure><img src="image_28.png" /><figcaption>FIG. 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.]</figcaption></figure>

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_c, T_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_d, T_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.

## B. 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 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.C.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.D, 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 pre-
cisely, 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_c, T_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.E.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.

### C. Oscillations in the Oceans' Thermohaline Circulation

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

Basically, the THC is due to denser water sinking,

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

FIG. 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_{in} - R_{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_{rain} - P_{evap}$ for the ocean and $P_{snow} - P_{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.]

lighter water rising, and water mass continuity closing the circuit through near-horizontal flow between the areas of rising and sinking<fnref n="5" /> 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

<footnote n="5">⁵ A 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).</footnote>
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 paleo-circulations; the association of the latter with salinity effects that exceed the thermal ones (e.g., during the Palaeocene, ≈ 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., Lucarini and Stone, 2005; Rooth, 1982; Scott et al., 1999; Titz et al., 2002).

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 non-linear PDEs with two independent space variables, say latitude and depth (e.g., Cessi and Young, 1992; Lucarini et al., 2005a, 2007; Quon and Ghil, 1992). 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 anti-symmetric 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

<figure><img src="image_30.png" /><figcaption>FIG. 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.]</figcaption></figure>

and Ghil (2005):

$$ \dot{X} = f(X; \mu) = \mu X - X^3. \qquad (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.

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

## 2. Bifurcation Diagrams for GCMs

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.B 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 I<fnref n="6" />. Such oscillatory behavior seems to match more closely the instrumentally recorded THC variability, as

<figure><img src="image_31.png" /><figcaption>(a) Intermediate-complexity models</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>FIG. 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.</figcaption></figure>

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 $λ_{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.

The self-sustained THC oscillations in question are

<footnote n="6">⁶ Anomalies in the table are defined as the difference between the monthly mean value of a variable and its climatological mean.</footnote>
<table><caption>TABLE I: Oscillations in the oceans' thermohaline circulation. See also Ghil (1994).</caption><thead><tr><th>Time scale</th><th>Phenomena</th><th>Mechanism</th></tr></thead><tbody><tr><td>Decadal</td><td>· Local migration of surface anomalies. in the NW corner of the ocean basin<br>· Gyre advection in mid-latitudes</td><td>· Localized surface density anomalies due to surface coupling<br>· Gyre advection</td></tr><tr><td>Centennial</td><td>Loop-type, meridional circulation</td><td>Conveyor-belt advection of density anomalies</td></tr><tr><td>Millennial</td><td>Relaxation oscillation, with “flushes” and superimposed decadal fluctuations</td><td>Bottom water warming, due to strong braking effect of salinity forcing</td></tr></tbody></table>

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.D.2 below.

## D. Bistability, Oscillations and Bifurcations

In Sects. III.A, III.B and III.C.1, we have introduced bistability of steady-state solutions via saddle-node and pitchfork bifurcations, while in Sect. III.C.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.

### 1. Bistability and Steady-State Bifurcations

Section III.B 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 and 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:

$$ \dot{X} = \mu - X^2, \quad (18a) $$

$$ \dot{X} = (\mu - \mu_0) + (X - X_0)^2, \quad (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; Arnold, 2003; Guckenheimer and Holmes, 1983; 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.B 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.C.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, \quad (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 \to -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
<figure><img src="image_33.png" /></figure>

<figure><img src="image_34.png" /><figcaption>FIG. 25: Dependence of THC solutions on two parameters in a hybrid coupled model; the two parameters are the atmosphere-ocean coupling coefficient $\lambda_{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_{ao}$; this curve was obtained by fitting a parabola to the model's numerical-simulation results, shown as full and open circles.</figcaption></figure>

[From Chen and Ghil (1996) with permission from the American Meteorological Society.]

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., Ghil, 2017; Jiang et al., 1995b; Speich et al., 1995b).

## 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, \quad (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.

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} \ge 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{r} \ge 0$, according to

$$ \dot{\rho} = 2\rho(\mu + c\rho), \quad (21a) $$

$$ \dot{\theta} = \omega. \quad (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$.

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
<figure><img src="image_35.png" /><figcaption>FIG. 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.]</figcaption></figure>

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.

## E. 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.A-II.B-II.C. We review here some of the most important ones, and comment on the general features of such modes.

### 1. Modes of Variability and Extended-Range Prediction

Several large-scale spatial patterns of atmospheric co-variability 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.D, 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; Charney *et al.*, 1950b; Richardson, 1922) 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 (Leith and Kraichnan, 1972; Lorenz, 1969a,b; Thompson, 1957). 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
<figure><img src="image_36.png" /><figcaption>FIG. 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 \le \theta < 2\pi$. [After Ghil and Childress (1987) with permission from Springer Science+Business Media.]</figcaption></figure>

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; Ghil and Robertson, 2002; Tribbia and Anthes, 1987, 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.

## 2. Coupled Atmosphere-Ocean Modes of Variability

The best known of all these modes is the El Niño-Southern Oscillation (ENSO: Dijkstra, 2005; Philander, 1990, and references therein), mentioned earlier in Sect. II.B. 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.

Such strong episodes recurred every 2-7 years during the instrumental record of roughly 150 years, and the ENSO phenomenon, with its 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., Ghil *et al.*, 2002; Rasmuson *et al.*, 1990).

Even larger variance, accompanied by lesser regularity, is associated with a quasi-quadrennial mode (Ghil *et al.*, 2002; Jiang *et al.*, 1995a), 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
<figure><img src="image_37.png" /><figcaption>FIG. 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 https://www.meted.ucar.edu/ams/wim_2014/9b.html.]</figcaption></figure>

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

### 3. Atmospheric Low-Frequency Variability

As discussed in Sect. II.A, 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.B, 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’ descrip-

tion, 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.B. 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.B. 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; Jin and Ghil, 1990; Legras and Ghil, 1985). The potential role of these results in the controversy surrounding the effect of anthropogenic polar amplification on blocking frequency was mentioned in Sect. II.B.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 op-
posed 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⁻¹, three distinct oscillatory modes that do not seem
to synchronize: in the Tropics the MJO, with a period
≈ 50 days, and in mid-latitudes the Branstator-Kushnir
wave, with a period ≈ 30 days, along with the topo-
graphic oscillatory mode, with a period ≈ 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 re-
spect to the radius of the Earth, and several telecon-
nections mentioned so far extend across continents and
oceans, cf. Fig. 28 herein.

Be that as it may, let us now concentrate on the topo-
graphic ‘wave’ mode and its relationship with the block-
ing and zonal ‘particles’. Ghil *et al.* (2018) recently re-
viewed the evidence provided by a hierarchy of models for
the presence of a Hopf bifurcation that arises from the in-
teraction 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 Mar-
shall 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 conclu-
sions of the Ghil *et al.* (2018) review paper on the obser-
vational, theoretical and modeling literature of multiple
regimes and oscillatory modes of subseasonal variability.
Four complementary approaches to explaining this vari-
ability are illustrated in Fig. 29.

One approach to persistent anomalies in mid-latitude
atmospheric flows on subseasonal time scales is to con-
sider 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 = π/2 or t = 3π/2 for a sine wave sin(t); or else they

<figure><img src="image_38.png" /><figcaption>FIG. 29: Schematic overview of atmospheric low-frequency variability (LFV) mechanisms; see text for details. Reproduced from Ghil *et al.* (2018), with permission from Elsevier.</figcaption></figure>

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 tri-
ads between a topographic Rossby wave and two free
Rossby waves (Egger, 1978; Ghil and Childress, 1987;
Trevisan and Buzzi, 1980). Neither version of this line of
thought, though, explains the organization of the persis-
tent anomalies into distinct flow regimes.

Rossby *et al.* (1939) initiated a different, genuinely
nonlinear approach by suggesting that multiple equi-
libria 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 at-
mospheric flows could occur. This multiple-equilibria ap-
proach was then pursued vigorously in the 1980s (Char-
ney and DeVore, 1979; Charney *et al.*, 1981; Ghil and
Childress, 1987; Legras and Ghil, 1985) and it is illus-
trated in Fig. 29 by the sketch labeled (a): one version
of the sketch illustrates models that concentrated on the
B–Z dichotomy (Benzi *et al.*, 1986; Charney and De-
Vore, 1979; Charney *et al.*, 1981), the other on models
(e.g., Legras and Ghil, 1985) that allowed for the pres-
ence of additional clusters, found by Kimoto and Ghil
(1993a) and Smyth *et al.* (1999), among others, in obser-
vations. 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 oscil-
latory 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 aris-
ing 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; Cvitanović 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 in 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 (Kondrashov et al., 2015, 2013, 2006; Kravtsov et al., 2005, 2009) 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.

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

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: $O(100)$ km in the ocean vs. $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.E.3, does correspond in the ocean to a time scale of years to decades.
<figure><img src="image_39.png" /><figcaption>FIG. 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.</figcaption></figure>

## 1. A simple model of the double-gyre circulation

The simplest model that includes several of the most pertinent mechanisms described in Sect. II.B 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:

$$ \partial_t q + J(\psi, q) - \nu \Delta^2 \psi + \mu \Delta \psi = -\tau \sin\left(\frac{2\pi y}{L_y}\right), (22a) $$

$$ q = \Delta \psi - \lambda_R^{-2} \psi + \beta y. \qquad (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.C.3, so that $J(\psi, q) = \psi_x q_y - \psi_y q_x = \mathbf{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_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 (Dijkstra and Ghil, 2005; Jiang et al., 1995b).

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 $S$,

$$ S[\psi(x, y)] = -\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.

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

a. 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><img src="image_40.png" /><figcaption>FIG. 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 $τ$. 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).</figcaption></figure>

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 $τ$ (Cessi and Ierley, 1995; Jiang et al., 1995b). As $τ$ 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, $τ = τ_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 voracity 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.

b. *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 (Dijkstra and Katsman, 1997; Sheremet et al., 1997; Speich et al., 1995a). 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 **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 **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 **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 **Pitchfork** bifurcation in the figure, and arises in much more complex situations and models (e.g., Simonnet, 2005; Simonnet et al., 2003b).

In fact, such merging is quite common in small-
dimensional dynamical systems with symmetry, as ex-
emplified 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 pe-
riod — 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
(Dijkstra, 2005; Dijkstra and Ghil, 2005; Ghil, 2015,
2017; Ghil and Childress, 1987). It is precisely this kind
of near-periodic recurrence that is identified in the cli-
mate sciences as LFV.

c. *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 con-
firmed 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 (Chang *et al.*, 2001; Meacham, 2000; Nadiga and
Luce, 2001; Simonnet *et al.*, 2005, 2003a,b, 1995). 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 (Ghil and Childress, 1987; Guckenheimer
and Holmes, 1983). In effect, the quasi-geostrophic model
reviewed herein undergoes a genuine homoclinic bifurca-
tion 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 bi-
furcations. 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 Simon-
net 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 or-
bits 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 ana-
lyzed so far — it is reinforced by results with much more
detailed models in the hierarchy (e.g., Dijkstra, 2005; Di-
jkstra and Ghil, 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 climat-
ically important 7-year and 14-year oscillations of the
oceanic circulation and the related atmospheric phenom-
ena in and around the North Atlantic basin (Da Costa
and Colin de Verdiére, 2004; Dijkstra, 2005; Dijkstra and
Ghil, 2005; Feliks et al., 2010, 2011, 2004, 2007; Kon-
drashov et al., 2005a; Moron et al., 1998; Plaut et al.,
1995; Simonnet et al., 2005, 2003b; Wunsch, 1999). 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 (Berloff et al.
2007; Dijkstra and Ghil, 2005; Feliks et al., 2007; Ghil
and Childress, 1987; Gill, 1982; Kravtsov et al., 2007;
Pedlosky, 1987, 1996; Simonnet et al., 2003b; Stommel,
1965) surely play a fundamental role, as they do in the ob-
served 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 under-
standing 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 situ-
ation, dynamical systems theory by itself cannot be used
as a full explanation of complex, observed behavior re-
sulting from successive bifurcations that are rooted in
simple stationary or periodic solutions.