2.1 Model hierarchy

In Dijkstra (2013), I suggested to classify climate models according to the two traits “scales” and “processes” (Fig. 3). Here the trait “scales” refers to both spatial and temporal scales as there exists a relation between both: on smaller spatial scales usually faster processes take place. “Processes” refers to either physical, chemical or biological processes taking place in the different climate subsystems. Models with a limited number of processes and scales are usually referred to as conceptual climate models. Examples are box models of the ocean circulation (Stommel, 1961) and models of glacial–interglacial cycles (Crucifix, 2012), all formulated by small-dimensional systems of ordinary differential equations. Limiting the number of processes, scales can be added by discretizing the governing partial differential equations spatially up to three dimensions. A higher spatial resolution and inclusion of more processes will give models located in the upper right part of the diagram, so-called Earth system models (ESMs). Between the conceptual models and ESMs are so-called intermediate complexity models which are spatially extended (described by partial differential equations) but with fewer scales and/or processes (Fig. 3).

Figure 3Organization of climate models according to the two model traits: number of processes and number of scales (Dijkstra, 2013).

Download

Any spatially extended climate model consists of a set of conservation laws, which are formulated as a set of coupled partial differential equations, that can be written in general form as (Griffies, 2004)

where ℒ and ℳ are linear operators, 𝒩 is a nonlinear operator, u is the state vector, ℱ contains the forcing of the system, and λ indicates the dependence of the operators on parameters. Appropriate boundary and initial conditions have to be added to this set of equations for a well-posed problem.

When Eq. (1) is discretized, eventually a set of differential equations with algebraic constraints arises, which can be written as

where x∈ℝⁿ is the state vector, n its dimension, M_λ is a (often singular) matrix of which every zero row is associated with an algebraic constraint, L_λ is the discretized version of ℒ_λ, and F_λ and N_λ are the finite-dimensional versions of the forcing and the nonlinear operator, respectively.

When noise is added, the evolution of the flow can generally be described by a stochastic differential-algebraic equation of the form

where X_t is now the stochastic state vector, f_λ(X_t) contains the linear and nonlinear processes, and the forcing. The quantity ${\mathit{W}}_t\in R^{n_{\mathrm{w}}}$ is a vector of n_w-independent standard Brownian motions (Gardiner, 2009) and ${\mathit{g}}_{\mathit{\lambda }}\left({\mathit{X}}_t\right)\in R^{n\times n_{\mathrm{w}}}$. Here, we use the Itô interpretation to represent unresolved processes but the Stratonovich interpretation can also be used (depending on the nature of processes represented).

2.2 Continuation methods

These methods form part of the numerical bifurcation analysis toolbox; here we are restricted to a single parameter λ. Finding steady states of the system (Eq. 2) versus λ amounts to solving

The idea of pseudo-arc-length continuation (Keller, 1977; Seydel, 1994) is to parametrize branches of solutions $\mathrm{\Gamma }\left(s\right)\equiv \left(\mathit{x}\right(s),\mathit{\lambda }(s\left)\right)$ with an arc-length parameter s (or an approximation of it, thus the term “pseudo”) and choose s as the continuation parameter. An additional equation is obtained by approximating the normalization condition of the tangent $\dot{\mathrm{\Gamma }}\left(s\right)=\left(\dot{\mathit{x}}\left(s\right),\dot{\mathit{\lambda }}\left(s\right)\right)$ to the branch Γ(s), where the dot refers to the derivative with respect to s, with ${\left|\dot{\mathrm{\Gamma }}\right|}^{\mathrm{2}}=\mathrm{1}$. More precisely, for a given solution (x₀,λ₀), the next solution (x,λ) is required to satisfy the constraint

where ${\dot{\mathrm{\Gamma }}}_{\mathrm{0}}=\left({\dot{\mathit{x}}}_{\mathrm{0}},{\dot{\mathit{\lambda }}}_{\mathrm{0}}\right)$ is the normalized direction vector of the solution family Γ(s) at (x₀,λ₀) and Δs is an appropriately small step size. Equation (5) stipulates that the projection of $(\mathit{x},\mathit{\lambda })-({\mathit{x}}_{\mathrm{0}},{\mathit{\lambda }}_{\mathrm{0}})$ onto $({\dot{\mathit{x}}}_{\mathrm{0}},{\dot{\mathit{\lambda }}}_{\mathrm{0}})$ has the value Δs. Efficient solution methods for high-dimensional systems of the form (Eqs. 4–5) are presented in De Niet et al. (2007) and Thies et al. (2009).

Suppose that the deterministic part of Eq. (3) has a stable fixed point ${\mathit{x}}_{\mathit{\lambda }}^{*}$ for a given range of parameter values. Then linearization of Eq. (3) around the deterministic steady state yields, with ${\mathit{Y}}_t={\mathit{X}}_t-{\mathit{x}}_{\mathit{\lambda }}^{*}$ (Kuehn, 2012),

where A_λ(x)≡(D_xf_λ)(x) is the Jacobian matrix and B_λ(x)=g_λ(x).

In the special case that M_λ is a non-singular matrix, the Eq. (6) can be rewritten (dropping the arguments and subscripts on the matrices) in Itô form as

which represents an n-dimensional Ornstein–Uhlenbeck process. The corresponding stationary covariance matrix C is then determined from the generalized Lyapunov equation:

The Gaussian probability density function at the fixed point can then be computed directly from C. When M is singular, special methods have been devised to cope with the singular part (Baars et al., 2017); also in this case, a generalized Lyapunov equation determines the covariance matrix C. Efficient solution methods for high-dimensional versions of Eq. (8) are presented in Baars et al. (2017).

3.1 Tipping points

An overview of possible tipping elements in the Earth's system is given in Lenton et al. (2008) and Steffen et al. (2018) and includes the Arctic winter sea ice, the marine ice sheets (MISs), the Amazon rainforest and the Atlantic Meridional Overturning Circulation (AMOC). Several of the associated transitions are thought to be associated with the existence of a multiple equilibrium regime associated with a back-to-back saddle-node structure, in particular the collapse of the AMOC and that of the MIS.

For a back-to-back saddle-node bifurcation there are two transition scenarios possible, called (i) bifurcation tipping and (ii) noise-induced tipping (Ditlevsen and Johnsen, 2010). In case (i) the parameter crosses a value at one of the saddle-node bifurcations, and in case (ii) a finite amplitude perturbation in the state vector causes a transition (even for a fixed value of the parameter). In the non-autonomous case also rate-induced tipping (Ashwin et al., 2012) is possible. For both cases (i) and (ii), it is crucial to determine the extent of the multiple equilibrium regime (Fig 4): this has been investigated in detail in spatially extended models for the following problems.

• AMOC. In a substantial number of papers, the bifurcation diagrams for both spatially two- and three-dimensional (intermediate complexity) ocean-only models of the AMOC have been determined (see chap. 6 in Dijkstra, 2005). The most advanced result is for a global ocean model coupled to an energy balance atmospheric model (den Toom et al., 2012), capturing ocean–atmosphere feedbacks, where also a back-to-back saddle-node structure was found.

  Numerical bifurcation analyses provided the basis for the stability indicator $\mathrm{\Sigma }=M_{\mathrm{ov}}^{\mathrm{s}}-M_{\mathrm{ov}}^{\mathrm{n}}$ of the multiple equilibrium regime of the AMOC (Huisman et al., 2010). Here, M_ov is the AMOC-induced freshwater transport, and the superscripts n and s indicate the northern and southern boundaries of the Atlantic, respectively. Although the central idea was already formulated in Rahmstorf (1996) and de Vries and Weber (2005), only bifurcation analysis of high-dimensional discretized ocean models provided more rigorous support for the use of this indicator. Although this stability indicator has its problems (Gent, 2018) and needs to be extended in a coupled ocean–atmosphere context, it is often used to investigate the stability of the AMOC (Hawkins et al., 2011).

  For a spatially two-dimensional ocean-only model, the covariance matrices C were determined from solving the Lyapunov equation in Eq. (8) in Baars et al. (2017) for the case of noise in the freshwater forcing. While here it served only to test the new Lyapunov equation solver (RAILS), the methodology was extended recently to compute (noise-induced) transition probabilities of the AMOC and to relate that probability to the stability indicator Σ (Castellana et al., 2019). Such transitions are thought to be involved in the Dansgaard–Oeschger (DO) events (Ditlevsen and Johnsen, 2010).

• MIS. The explicit computation of the bifurcation structure of a spatially one-dimensional marine ice-sheet model (with a moving grounding line) has been carried out only recently (Mulder et al., 2018). Here also a back-to-back saddle-node structure is found, of course compatible with the many initial value problem studies of this model (Schoof, 2007). The gravitational effect of a marine ice sheet on sea level has a stabilizing influence on the ice sheet as seen through a shift in the bifurcation diagram. While this was found in many different model studies (Gomez et al., 2010), the precise mechanism could be deduced from the bifurcation diagram shift (Mulder et al., 2018). For a stochastic MIS model, the covariance matrices C were determined for each stable state in Mulder et al. (2018). Typical noise in the accumulation leads to grounding line variations on the order of 1000 m, while for sea level noise this is about 100 m. In the multiple equilibrium regime, the study of the transition probabilities indicated that, for both noise types, it is more likely to jump from a large ice sheet state to a small ice sheet state than vice versa.

Figure 4The canonical bifurcation diagram with the back-to-back saddle node indicating two stable states (a, c) and an unstable state (b). Bifurcation tipping occurs when the parameter λ crosses the value at L₁ or L₂. Noise-induced tipping (e.g., from state a to state c) can occur through a perturbation in the state vector (for fixed λ).

Download

In high-dimensional climate models, also so-called edge states or Melancholy states have been computed, for example in a coupled atmospheric sea-ice model investigating ice-covered/ice-free multi-stability (Lucarini and Bódai, 2017). This edge state is a saddle embedded in the boundary between the two basins of attraction of the stable climate states.

3.2 Patterns of sea surface temperature (SST) variability

It is remarkable that on interannual-to-multidecadal timescales the variability in sea surface temperature is organized in large-scale patterns (Fig. 5). These patterns have been detected by using multi-variate analysis on long data sets, such as the HadISST, for example through principal component analysis where the patterns are then contained in the empirical orthogonal functions (EOFs). Well-known and much-studied patterns are those of the El Niño–Southern Oscillation (ENSO), the Atlantic Multidecadal Oscillation (Enfield et al., 2001) and the Pacific Decadal Oscillation (PDO) (Mantua et al., 1997), as shown in Fig. 5.

Figure 5Overview of patterns of climate variability (AMO, PDO and ENSO) as determined in Deser et al. (2010) with accompanying time series.

Numerical bifurcation analysis has been applied to several spatially extended models, in particular ENSO, PDO and the AMO.

• ENSO. The cornerstone intermediate complexity model is the Zebiak–Cane (ZC) model of which the behavior has been extensively analyzed (Zebiak and Cane, 1987). Numerical bifurcation analysis of different versions of the ZC model were performed (see chap. 7 in Dijkstra, 2005), and in each of them a Hopf bifurcation occurs once the coupling strength between the equatorial ocean and atmosphere crosses a critical value. The period of the leading mode is in the interannual range and determined by basin modes, just as in the recharge–discharge oscillator model (Jin, 1997). The spatial pattern of the leading mode is localized into the cold tongue region of the mean (steady) state and shares many similarities with the first EOF from observations.

• AMO. One of the intermediate complexity models which has been used is a spatially three-dimensional model of the North Atlantic (see chap. 8 in Dijkstra, 2013). The mean (steady) state is generated by a horizontal atmospheric surface buoyancy field which drives a meridional overturning circulation in the ocean model. Numerical bifurcation analysis of this model has shown that the background state destabilizes through a multi-decadal leading mode, and hence a Hopf bifurcation occurs. The timescale of the leading mode can be linked to the basin crossing time of temperature anomalies, and its spatial pattern shares many features with the observed pattern (Kushnir, 1994) when a representation of the continents is considered. The variability can be easily excited through noise in the heat flux, even when the leading mode is decaying in the deterministic case (see chap. 8 in Dijkstra, 2013).

• PDO. The same bifurcation analysis as for the AMO was applied to a model of two ocean basins which are connected by the Southern Ocean (von der Heydt and Dijkstra, 2007). It was found that the PDO cannot be related to a single mode of variability which arises through a Hopf bifurcation (as for ENSO and the AMO). The key here is that destabilization of the mean (steady) flow can only occur when there is sinking (by the AMOC) in the northern part of the basin. Such sinking is absent in the North Pacific for the present-day climate. Indeed, modern views of the PDO indicate that several different mechanisms are likely important for the existence of PDO variability (Newman et al., 2016).

My interpretation of these results is that several of these SST patterns (but not all) appear through a normal mode which destabilizes the mean state through positive feedbacks; the presence of negative feedbacks causes the oscillatory behavior. In this case, the associated Hopf bifurcation (of a spatially extended model) provides both the dominant timescale of variability and its spatial pattern. The elegant structure of leading modes in ocean models and the ZC model was presented in Dijkstra (2016). The key to why normal modes can be dominant in this variability may be that the nonlinearity in these models is rather weak (it involves advection of heat/salt and not of momentum). Hence the mean state is not modified (rectified) much due to the nonlinear interactions (contrary to the processes shown in the next section).

3.3 Collective interactions and emergent behavior

In the previous two subsections climate system variability phenomena were attributed to low-order dynamics. However, there are many phenomena which are intrinsically caused by the collective interaction of multiple instabilities. Clearly, the role of numerical bifurcation theory becomes quite limited in determining the behavior of these (in general) chaotic dynamical systems; I briefly describe below two examples.

• Ocean western boundary current variability (WBC). In pure wind-driven barotropic shallow-water models, the bifurcation structure consists of an imperfect pitchfork bifurcation followed by several Hopf bifurcations containing two types of modes: Rossby modes and so-called gyre modes. These gyre modes are eventually responsible for homoclinic orbits which lead to ultra-low-frequency behavior (see chap. 5 in Dijkstra, 2005). When another layer is added, baroclinic instabilities lead to a range of normal modes (Simonnet et al., 2003) which all destabilize. Hence the eventual emergent behavior is induced by their collective interactions. Such behavior can lead to low-frequency behavior in the statistical steady state which has been referred to as the turbulent oscillator (Berloff et al., 2007).

• Midlatitude atmospheric variability (MAV). Analysis of barotropic spatially extended models of midlatitude atmospheric flows has shown that there are many (highly) unstable equilibria (Legras and Ghil, 1983; Crommelin, 2003). Also here the resulting variability occurs through a collective interaction and low-frequency variability emerges in the statistical equilibrium state (Crommelin et al., 2004). Approaches to understand the dynamics on the attractor have been proposed through transfer operator methods (Tantet et al., 2015).

The cases briefly described above are examples of strongly nonlinear systems, where the nonlinearities occur in the momentum advection and where the mean state is strongly modified through rectification. Of course, there are many more examples of such geophysical systems, in particular on timescales up to interannual both in the ocean (internal waves) and the atmosphere (weather).