# Orbital insolation variations, intrinsic climate variability, and Quaternary glaciations

Keno Riechers, Takahito Mitsui, Niklas Boers, Michael Ghil

► To cite this version:

Keno Riechers, Takahito Mitsui, Niklas Boers, Michael Ghil. Orbital insolation variations, intrinsic climate variability, and Quaternary glaciations. Climate of the Past, 2022, 18, pp.863-893. 10.5194/cp-18-863-2022. insu-03726911

HAL Id: insu-03726911

https://insu.hal.science/insu-03726911

Submitted on 22 Jul 2022

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

# Orbital insolation variations, intrinsic climate variability, and Quaternary glaciations

Keno Riechers<fnref n="1" />,<fnref n="2" />, Takahito Mitsui<fnref n="2" />,<fnref n="3" />, Niklas Boers<fnref n="2" />,<fnref n="3" />,<fnref n="4" />, and Michael Ghil<fnref n="5" />,<fnref n="6" />

**Correspondence:** Keno Riechers (riechers@pik-potsdam.de)

Received: 8 October 2021 – Discussion started: 14 October 2021

Revised: 1 March 2022 – Accepted: 3 March 2022 – Published: 21 April 2022

**Abstract.** The relative role of external forcing and of intrinsic variability is a key question of climate variability in general and of our planet's paleoclimatic past in particular. Over the last 100 years since Milankovic's contributions, the importance of orbital forcing has been established for the period covering the last 2.6 Myr and the Quaternary glaciation cycles that took place during that time. A convincing case has also been made for the role of several internal mechanisms that are active on timescales both shorter and longer than the orbital ones. Such mechanisms clearly have a causal role in Dansgaard-Oeschger and Heinrich events, as well as in the mid-Pleistocene transition. We introduce herein a unified framework for the understanding of the orbital forcing's effects on the climate system's internal variability on timescales from thousands to millions of years. This framework relies on the fairly recent theory of non-autonomous and random dynamical systems, and it has so far been successfully applied in the climate sciences for problems like the El Niño-Southern Oscillation, the oceans' wind-driven circulation, and other problems on interannual to interdecadal timescales. Finally, we provide further examples of climate applications and present preliminary results of interest for the Quaternary glaciation cycles in general and the mid-Pleistocene transition in particular.

## 1 Introduction and motivation

In the early 20th century, Milutin Milankovic presented his theory of ice ages (Milankovitch, 1920). Based on his own calculations and on insightful suggestions from Wladimir Köppen and Alfred Wegener (Imbrie and Imbrie, 1986), he proposed that the transitions between glacial and interglacial climate conditions were primarily caused by variations of incoming solar radiation, which by that time was known to vary in a quasi-periodic manner on slow timescales of tens of thousands to hundreds of thousands of years (Poincaré, 1892–1899). These variations of insolation, which arise as a consequence of the gravitational interaction of the Earth with the other planets and with its own Moon, are typically referred to as orbital forcing.

The orbital forcing comprises variations in (i) the eccentricity of the Earth's orbit around the sun with dominant spectral peaks around 400 and 100 kyr; (ii) the obliquity, or axial tilt, i.e., the angle between the Earth's rotational and its orbital axis, with dominant periodicity around 41 kyr; and (iii) the climatic precession, which determines the phase of the summer solstice along the Earth's orbit and has its most pronounced spectral power around 23 and 19 kyr (Berger, 1978).

For 2 centuries or so of modern geology, records of our planet's physical and biological past were merely discrete sequences of strata with specific properties, like coloration

<footnote n="1">¹Department of Mathematics and Computer Science, Freie Universität Berlin, Berlin, Germany</footnote>

<footnote n="2">²Potsdam Institute for Climate Impact Research, Potsdam, Germany</footnote>

<footnote n="3">³Earth System Modelling, School of Engineering & Design, Technical University of Munich, Munich, Germany</footnote>

<footnote n="4">⁴Department of Mathematics and Global Systems Institute, University of Exeter, Exeter, UK</footnote>

<footnote n="5">⁵Geosciences Department and Laboratoire de Météorologie Dynamique (CNRS and IPSL), Ecole Normale Supérieure and PSL University, Paris, France</footnote>

<footnote n="6">⁶Department of Atmospheric and Oceanic Science, University of California at Los Angeles, Los Angeles, USA</footnote>
and composition (Imbrie and Imbrie, 1986). This state of af-
fairs led, after the initial success of the Milankovitch (1920)
theory of the ice ages, to severe criticism of the temporal
mismatch between insolation minima and glaciation maxima
(e.g., Flint, 1971).

The advent of marine-sediment cores after World War II
led, for the first time, to the availability of records that were
more or less continuous in time. Like all climate records,
these cores covered limited time intervals and did so with
limited resolution and with inaccuracies in absolute dating,
as well as in the quantities being measured. Moreover, they
posed the problem of inverting proxy records of isotopic and
microbiotic counts to physical quantities like temperature
and precipitation.

In spite of these limitations, the spectral analysis of deep-sea records allowed Hays et al. (1976) to overcome the difficulties previously encountered by the orbital theory of Quaternary glaciations, in particular the absence of the imprint of precessional and obliquity peaks in glaciation proxy records. Specifically, Hays et al. (1976) were able to create a composite record – back to over 400 kyr b2k, i.e., over 400 000 years before the year 2000 CE – from two relatively long marine-sediment records of the best quality available in the early 1970s. The authors demonstrated therewith that significant precessional and obliquity peaks near 20 and 40 kyr were present in this record’s spectral analysis; see Fig. 1. The power spectrum in the figure also made it quite clear that these peaks were superimposed on a continuous background – the stippled area in the figure – whose total variance much exceeded the sum of the variances present in the peaks.

The work of Hays et al. (1976) and of the subsequent CLIMAP and SPECMAP projects resulted in a much more detailed spatiotemporal mapping of the Quaternary and extended the belief in the pacemaking role of orbital variations into the more remote past. The spectral peaks near 20 and 40 kyr have been widely interpreted within the geological community as evidence for a linear response of the climate system to the orbital forcing (Imbrie and Imbrie, 1986). A third spectral peak at 100 kyr was, however, the most pronounced but much more difficult to reconcile with the orbital theory of Quaternary glaciations. Since no sufficiently pronounced counterpart can be found in the spectra of the seasonal insolation forcing, Hays et al. (1976) hypothesized a nonlinear response of the climate system in order to explain this dominant periodicity of the Late Pleistocene glacial-interglacial cycles. At the same time, the advent of higher-resolution marine cores and especially ice cores from both Greenland and Antarctica, led to the discovery of Heinrich events (Heinrich, 1988), Dansgaard-Oeschger (D-O) events (e.g., Dansgaard et al., 1993), and Bond cycles (Bond et al., 1997), which were hard to explain by orbital forcing given their shorter timescales.

In fact, interest in past climates was heightened not only by
these striking observational discoveries but also by the grow-
ing concerns about humanity's impact on the climate (SMIC,

1971; National Research Council, 1975). Given the declining
temperatures between the 1940s and 1970s, on the one hand,
as shown in Fig. 2 (see also Ghil and Vautard, 1991, Fig. 3),
and the substantial advances in the description of the Qua-
ternary glaciations, on the other hand, it is clear that interest
was mainly on the planet's falling into another ice age (e.g.,
National Research Council, 1975).

As a result of the twofold stimulation provided by data
about past glaciations and concern about future ones, a num-
ber of researchers in the early-to-mid 1970s worked on en-
ergy balance models (EBMs) of climate with multiple sta-
ble steady states (Held and Suarez, 1974; North, 1975; Ghil,
1976). Two such stable “equilibria” corresponded to the
present climate and to a “deep-freeze”, as it was called at
the time, i.e., to a totally ice-covered Earth. At the time there
was some disbelief about this second climate, as its calcu-
lated temperatures were much lower than those associated
with the Quaternary glaciations and incompatible with pale-
oclimatic evidence available in the 1970s.

New geochemical evidence, however, led in the early
1990s to the discovery of a “snowball” or, at least, “slush-
ball” Earth prior to the emergence of multicellular life, some-
time before 650 Myr b2k (Hoffman et al., 1998). It thus
turned out that this climate state – predicted by several
EBMs and confirmed by a general circulation model (GCM)
with much higher spatial resolution (Wetherald and Manabe,
1975) – had actually occurred and is now being modeled in

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

**Figure 1.** Power spectrum of a composite $\delta^{18}$O record using deep-sea cores RC11-120 and E49-18. This figure is based on the work of Hays et al. (1976), as presented by Imbrie and Imbrie (1986). Reprinted by permission from Springer Nature Customer Service Centre GmbH: Springer Nature. Topics in Geophysical Fluid Dynamics: Atmospheric Dynamics, Dynamo Theory, and Climate Dynamics by Ghil and Childress (1987), © 1987 by Springer Science+Business Media New York. All rights reserved. 1987.
<figure><img src="image_4.png" /></figure>

**Figure 2.** Comparison of six analyses of the annually and globally averaged surface temperature anomalies through 2018. The abscissa is time in years, and the ordinate is temperature anomalies in °C with respect to a 30-year climatological average for 1951–1980. NASA stands for the National Aeronautics and Space Administration; NOAA stands for the National Oceanic and Atmospheric Administration. Reprinted by permission from John Wiley & Sons Inc.: American Geophysical Union. Journal of Geophysical Research: Atmospheres. Improvements in the GISTEMP Uncertainty Model. Lenssen et al. (2019), © 2019. American Geophysical Union. All Rights Reserved. 2019.

much greater detail (Pierrehumbert, 2004; Ghil and Lucarini, 2020).

On the other hand, it also became clear that these early models, whose only stable solutions were stationary, could not reproduce the wealth of variability that the proxy records were describing very well, not even in the presence of stochastic forcing (e.g., Ghil, 1994). Certain theoretical paleoclimatologists therefore turned to coupling a “climate” equation, with temperature as its only dependent variable, with an ice sheet equation (Källén et al., 1979; Ghil and Le Treut, 1981) or a carbon dioxide equation (Saltzman et al., 1981; Saltzman and Maasch, 1988). These coupled climate models, albeit highly idealized, did produce oscillatory solutions that captured some of the features of the Quaternary glaciation cycles as known at that time. For instance, the models of Ghil and associates (Källén et al., 1979; Ghil and Le Treut, 1981) captured the phase differences between peak ice sheet extent and minimum temperatures in the North Atlantic suggested by Ruddiman and McIntyre (1981), while the work of Saltzman and associates (e.g., Saltzman and Maasch, 1988) captured the asymmetry of the glaciation cycles with their more rapid “terminations” (Broecker and Van Donk, 1970).

The stable self-sustained oscillations of these coupled models, however, were totally independent of any orbital or other time-dependent forcing, i.e., the solar input to their radiative budget was constant in time. Hence, they could not capture the wealth of spectral features, with their orbital and other peaks, of the paleo-records available by the 1980s. The basic quandary of the Quaternary glaciation cycles – at least

from the point of view of theoretical climate dynamics (Ghil and Childress, 1987, Part IV) – is formulated in Fig. 3 below; see also Ghil (1994): How does the quasi-periodic orbital forcing, with its relatively narrow spectral peaks, act on the climate’s internal variability to produce a response that is characterized by significant spectral peaks superimposed on a broad continuous background? In addition, how does the climate’s spectral peak at 100 kyr arise given its absence in the power spectrum of the forcing?

In this paper, we try to show a path toward resolving the four fundamental questions listed in Box 1. In the next section, we summarize existing results on how the climate system’s intrinsic variability arises on Quaternary timescales and on how this variability is modified by the time-dependent orbital forcing, which was added to the previously autonomous climate models as the next step in paleoclimate modeling evolution; see, for instance, Le Treut and Ghil (1983) and Le Treut et al. (1988) vs. Ghil and Le Treut (1981). In Sect. 3, we outline a more general framework for the study of such mechanisms, as given by the theory of non-autonomous and random dynamical systems (NDSs and RDSs), and sketch an application of this theory to other climate problems. An application to the problem at hand is proposed in Sect. 4 and conclusions follow in Sect. 5.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** The theoretical quandary of modeling the Quaternary glaciation cycles, illustrated here by schematic diagrams of the composite power spectral densities (PSD) of (a) the paleo-records and (b) the orbital forcing. In panel (a), the dominant peak for the Late Pleistocene is near 100 kyr, while in panel (b) eccentricity forcing is distributed over several spectral lines. The peaks at 6–8 and 1–2 kyr in panel (a) correspond to Bond cycles (Bond et al., 1992, 1993) and to the mean recurrence of D-O events, and they lack a match in the forcing lines of panel (b).</figcaption></figure>

Based on the discussion so far and on the discrepancy between Fig. 3a and b, we propose herein the following questions that still need to be resolved for a thorough understanding of Quaternary climate variability.

1. How does the dominant peak of the observed variability near 100 kyr arise, given the rather diffuse orbital forcing close to this periodicity?

2. What causes the continuous part of the observed spectrum, which contains most of the variance?

3. What gives rise to the high-frequency peaks due to the sudden warmings and to the approximate periodicities associated with Heinrich and D-O events, among others?

4. What are the contributions of the orbital forcing and the climate system's intrinsic variability to items (1)–(3), and how does the former factor modify the latter?

**Box 1.** Fundamental questions regarding the Quaternary glacial-interglacial cycles.

## 2 Self-sustained climate oscillators

### 2.1 A simple mechanism for climate oscillations

We follow Ghil (1994) in sketching the simplest physical mechanism for a self-sustained climate oscillation at fixed insolation forcing. Consider the Källén et al. (1979, hereafter KCG) oscillator, to the best of our knowledge the first such self-sustained climate oscillator. The model itself was built from the ground up, coupling a scalar version of the Ghil (1976) energy balance model (EBM) with a simplified, scalar version of the Weertman (1964, 1976) ice sheet model

(ISM). The model's details and further analyses of its ingredients and variants can be found in several references (e.g., Crafoord and Källén, 1978; Ghil and Tavantzis, 1983; Ghil, 1984; Bódai et al., 2015).

The basic workings of this climate oscillator can be represented by two coupled ordinary differential equations (ODEs), written symbolically as follows:

$$ \dot{T} \approx -V, \qquad (1a) $$

$$ \dot{V} \approx T. \qquad (1b) $$

Here $T$ stands for global temperature and $V$ for global ice volume, while Eq. (1a) is an EBM and Eq. (1b) is an ice sheet model (ISM). The “$\approx$” symbol stands for a binary relation of rough proportionality and is intended to neglect the details of the equation's right-hand side (RHS), including its nonlinearities. The EBM represents the well-known ice–albedo feedback used by both Budyko (1969) and Sellers (1969), while the ISM relies on the precipitation–temperature feedback postulated by KCG and used also by Ghil and Le Treut (1981), who coined the term.

The latter feedback can be better understood by writing the following equations:

$$ \dot{V} \approx p, \qquad (2a) $$

$$ p \approx p_{ac} - p_{ab}, \qquad (2b) $$

$$ p_{ac} \approx T. \qquad (2c) $$

Here $p$ is net precipitation on the single ice sheet of the globally integrated model, given by the difference in Eq. (2b) between the accumulation $p_{ac}$ and the ablation $p_{ab}$ (KCG).

As first observed by George C. Simpson – the meteorologist of Robert F. Scott's *Terra Nova* expedition to Antarctica in 1910–1912 and later the longest-serving Director of the U.K. Meteorological Office – warmer winters have more snow, and hence, at least in central Antarctica, the increase of $p_{ac}$ with $T$ exceeds the more obvious increase of $p_{ab}$ with
T. Hence, $p \simeq T$, and we have derived therewith Eq. (1b), i.e., $V \simeq p \simeq T$. For more recent studies of the precipitation–temperature feedback, see Tziperman and Gildor (2002).

More generally, the presence of feedbacks of opposite sign in a system of two linear coupled ODEs,

$$\dot{x} = y, \quad \dot{y} = -x,$$

leads to an oscillation, with the solution given by two trigonometric functions in quadrature with each other, $x(t) = \sin(t)$, $y(t) = \cos(t) = \sin(t + \pi/2)$, and the trajectory describing a circle in the $(x, y)$ phase plane, $x^2(t) + y^2(t) = 1$. In a nonlinear system, however – like the full KCG model or any other climate oscillator mentioned so far – the possibility of an oscillation, as indicated by the system (1), is actually realized in the explicit, full set of equations only for certain parameter values and not for others.

This can be understood by considering the so-called normal form of a Hopf bifurcation, which leads from a stable steady state, called a fixed point in dynamical systems theory, to a stable oscillatory solution, called a limit cycle. The easiest way to see this transition is by writing the normal form in polar coordinates, as in Arnold (1988) and in Ghil and Childress (1987, Sect. 12.2), namely,

$$\dot{z} = (\mu + i\omega)z + c(z\bar{z})z. \quad (3)$$

Here $z = x + i y$ is complex, where $i = \sqrt{-1}$ is the imaginary unit, while $\mu$ is a real bifurcation parameter, and $c, \omega$ are real and nonzero. Note that the KCG model per se is not in the normal form above, and we will discuss its bifurcation parameter $\mu_*$ in the next subsection.

A very natural transformation of variables,

$$\rho = z\bar{z} > 0, \quad z = \rho^{1/2} \exp(i\theta), \quad (4)$$

leads from the complex ODE (3) to the system of two real and decoupled ODEs,

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

$$\dot{\theta} = \omega. \quad (5b)$$

Equation (5b) simply provides an angular rotation around the origin $\rho = 0 = x = y$, since the complex exponential in Eq. (4) is periodic with period $2\pi$. Equation (5a) is quadratic in $\rho$, and thus it can have two real roots, $\rho = \rho_0 = 0$ and $\rho = \rho_* = -\mu/c$. But $\rho$ has to be positive, and thus in the case in which $c < 0$, the only possible solution for $\mu < 0$ is the origin, and it is stable since $\rho(\mu + c\rho)$ is negative for $\rho > 0$ in this case; hence, $\rho$ has to be monotonically decreasing, i.e., all the solutions of Eq. (5) spiral into the origin. The Hopf bifurcation from this stable steady state to a periodic solution, i.e., a limit cycle with radius $\rho^{1/2} = \rho_*^{1/2}$, occurs as $\mu$ crosses 0. Since now $\rho(\mu + c\rho) > 0$ for $0 < \rho < \rho_*$ and $\rho(\mu + c\rho) < 0$ for $\rho > \rho_*$, the limit cycle is stable and trajectories spiral out from inside this cycle and into it from outside; see Fig. 4 for the so-called supercritical (or soft) Hopf bifurcation case with $c < 0$.

## 2.2 Intrinsic climate oscillations and the mid-Pleistocene transition (MPT)

In this subsection, we present an argument for the role of intrinsic oscillations in the mid-Pleistocene transition (MPT). The first point to be made is that while orbital forcing clearly plays a major role in the power spectrum of the Quaternary’s climatic variability, it cannot be, in and of itself, the cause of the MPT. Indeed, changes in the solar system’s orbital periodicities only occur on much longer timescales than the entire Quaternary’s duration (Varadi et al., 2003; Laskar et al., 2004a). Our argument continues with a further analysis of the Hopf bifurcation presented in the previous subsection. Such an analysis was carried out for the KCG model by Ghil and Tavantzis (1983).

Physically speaking, the presence or absence of the regular, purely periodic oscillations obtained by KCG and illustrated in Ghil and Childress (1987, Fig. 12.6) depends on whether $c \ge 0$ in Eq. (5a). The KCG model’s bifurcation parameter is $\mu_* = c_T/c_L$, where $c_T$ is the heat capacity in its EBM, while $c_L$ is the “mass capacity” in its ISM (Ghil and Tavantzis, 1983). Large $\mu_*$ corresponds physically to a very small, possibly pre-Pleistocene ice cap (Ghil, 1984; Saltzman and Sutera, 1987). At these values of $\mu_*$, the KCG model’s nullclines and fixed points – the latter being given by the intersection of the former – are very different from those that are obtained for Quaternary-sized ice sheets, for which $c_L$ is comparable in value to $c_T$; see Ghil and Tavantzis (1983, Figs. 3–5). As $\mu_*$ decreases to $O(1)$, i.e., as we proceed from very small to more substantial ice sheets, the fixed point transfers its stability to a branch of periodic solutions via a subcritical Hopf bifurcation (Källén et al., 1979; Ghil and Tavantzis, 1983); see also Ghil and Childress (1987, Figs. 12.8 and 12.9).

To clarify the simple physical concepts that underlie subcritical and supercritical Hopf bifurcations, let us consider a purely mechanical oscillator with mass $m$, a spring $kx$, and a dashpot $\alpha\dot{x}$ (Landau and Lifshitz, 1960; Jordan and Smith, 1987),

$$m\ddot{x} = -\alpha(x)\dot{x} - k(x)x. \quad (6)$$

If $k = \text{const.}$ and $\alpha = \text{const.}$, we have the simplest, linear setup, but we will be interested in the nonlinear cases here. Normalizing by the mass $m$ and not changing notation otherwise, by rearranging terms and adding a periodic forcing we get the following equation:

$$\ddot{x} + \alpha(x)\dot{x} + k(x)x = F \cos(\omega t). \quad (7)$$

Two classical nonlinear cases are those of the Duffing (1918) equation, in which $k(x) = x^2$ and $\alpha = \text{const.}$, and of the Van der Pol (1926) equation, in which $k = \text{const.}$ and $\alpha(x) = \nu(x^2 - 1)$. The fully nonlinear case in which both the spring and the damping are nonlinear, with $k(x) = x^2$ and $\alpha(x) = \nu(x^2 - 1)$, is known as the Van der Pol–Duffing oscillator
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Supercritical Hopf bifurcation. (a) Vector field $\dot{z}(z)$ of Eq. (3) for the parameter values $\mu = -1$ and $c = -1$; $z = x+ty$. In this case, the origin constitutes the only stable fixed point, and all trajectories will spiral into this point, as illustrated by the single brown trajectory. (b) Vector field $\dot{z}(z)$ and three trajectories for $\mu = 1$. In this case, the origin is an unstable fixed point, while the limit cycle with radius $\rho_*^{1/2} = 1$ constitutes the only stable solution. Trajectories that start inside this limit cycle, with $\rho(t_0) < \rho_*$, tend to it by spiraling out – as illustrated by the magenta trajectory – while trajectories that start outside the limit cycle, with $\rho(t_0) > \rho_*$, approach it by spiraling inward, as illustrated by the gray trajectory. (c) Dependence of the stable solution of Eq. (3) on the parameter $\mu$, for $c = -1$. For $\mu \le 0$, the single stable solution is the fixed point located at the origin, $\rho \equiv 0$. For $\mu > 0$, the stable solution is the limit cycle given by $\rho = \mu$.</figcaption></figure>

(e.g., Jackson, 1991; Pierini et al., 2018). Note that all three types of nonlinear oscillators can exhibit chaotic behavior even in the presence of simple periodic forcing (e.g., Guckenheimer and Holmes, 1983; Pierini et al., 2018, and references therein). The idea of using such simple, classical oscillators in modeling Quaternary glaciation cycles goes back to Saltzman et al. (1981).

Jordan and Smith (1987, Sect. 5.6) specifically discuss the case of a soft and a hard spring for a generalized Duffing equation, with $k(x) = k_0 + \epsilon h(x)$, where $0 < |\epsilon| \ll 1$, $h(-x) = h(x)$, and $h'' > 0$. A spring is soft if it is sublinear, $\epsilon < 0$, and hard if it is superlinear, $\epsilon > 0$; see their Eq. (5.37) and Fig. 5.4a, b, with $h(x) = x^2$.

The supercritical Hopf bifurcation in the absence of forcing is analogous to the nonlinear response of a soft, sublinear spring to periodic forcing in which the oscillations in the position $x$ of the mass $m$ increase gradually in amplitude as the spring constant $k_0$ increases past a critical value $k_*$, while the subcritical case is analogous to the response of a hard, superlinear spring, for which the oscillations in $x$ jump suddenly from zero amplitude to a finite amplitude as the spring constant $k_0$ crosses the value $k_*$. Please compare the behavior of supercritical and subcritical Hopf bifurcations in Ghil and Childress (1987, Sect. 12.2 and Figs. 12.7–12.9) and see Jordan and Smith (1987, Fig. 5.7) for the change in the nonlinear response of a Duffing oscillator as its spring changes from soft, with $\epsilon < 0$, to hard, with $\epsilon > 0$.

There is a clear-cut analogy with the mid-Pleistocene transition, occurring at roughly 0.8 Ma b2k, at which small-amplitude climate variability with a dominant periodicity near 40 kyr becomes larger, dominated by a periodicity that is close to 100 kyr, as well as being more irregular (e.g., Huybers, 2009; Quinn et al., 2018; Rousseau et al., 2020). A fair number of distinct dynamical theories for this transition have been formulated (e.g., Maasch and Saltzman, 1990; Ghil, 1994; Crucifix, 2012; Ashwin and Ditlevsen, 2015; Daruka

and Ditlevsen, 2016; Omta et al., 2016; Ditlevsen and Ashwin, 2018). A rather obvious one is that a Hopf bifurcation occurs at that point, which leads to a more vigorous response to the multi-periodic orbital forcing; thus, the latter does not need to change in the least in order to explain the observed phenomena. In Saltzman and Sutera (1987), there is only a comment on the likely role of a Hopf bifurcation in the transition, but their Fig. 3 suggests that in their model such a bifurcation would have to be of the supercritical type and lead to a fairly gradual transition.

In contrast, the subcritical Hopf bifurcation of the KCG and Ghil and Le Treut (1981) oscillators would have to lead to a more abrupt transition, as suggested by Ghil (1984). Later, Crucifix (2012) showed that the models by Saltzman and Maasch (1990, 1991) exhibit MPT-like behavior via supercritical or subcritical Hopf bifurcations, depending on the parameter values. The existing $\delta^{18}$O and $\delta^{13}$C records might or might not have sufficient resolution back in time up to 1.2 Ma to settle this question about the abruptness of the transition. In case some of them do, an objective test of suddenness – as proposed by Bagniewski et al. (2021) for the high-resolution NGRIP record (North Greenland Ice Core Project members, 2004) – will have to be applied to such records.

# 3 Time-dependent forcing, NDSs, and RDSs

## 3.1 Orbital forcing of a climate oscillator

We start this section by describing some fairly simple ways in which the orbital forcing might have modified intrinsic climate variability, thus helping to solve the mismatch between Fig. 3a and b in Sect. 1. To explore this possibility, Le Treut and Ghil (1983) used somewhat simplified insolation forcing based on the calculations of Berger (1978) and applied it to a slightly modified version of the Ghil and Le Treut (1981) oscillator. These authors found that, as expected for a nonlinear oscillator, its internal frequency $f_0$ is affected strongly by
the forcing ones, {$f_1, \dots, f_5$}, resulting in both nonlinear resonance and combination tones (Landau and Lifshitz, 1960).

Depending on parameter values, the periodicity $P_0 = 1/f_0$ of the Ghil and Le Treut (1981) oscillator is $P_0 \approx 6-7$ kyr. The lines in the simplified insolation spectrum used by Le Treut and Ghil (1983) had the following periodicities: $P_1 = 19$ kyr, $P_2 = 23$ kyr, $P_3 = 41$ kyr, $P_4 = 100$ kyr, and $P_5 = 400$ kyr. These periodicities correspond to the two precessional ones, the obliquity one, and the two eccentricity ones. The actual celestial-mechanics calculations that Berger (1978) based his insolation calculations on have been substantially updated since (e.g., Varadi et al., 2003; Fienga et al., 2015). However, these advances have not modified the spectrum of the planetary-orbit solutions over the 2.6 Myr of the Quaternary very much, which is a rather short interval in celestial-mechanics terms.

The results of the Le Treut and Ghil (1983) model on the evolution of the primary climate variables *T* and *V* were converted to $\delta^{18}$O values in simulated isotopic records of marine sediment and ice cores by Le Treut et al. (1988); the spectra of the latter are plotted in Fig. 5. The values on the abscissa of Fig. 5a are values of the logarithm of frequency, while those in Fig. 5b are values of the frequency itself; the values on the ordinate of both panels are powers of 10. One refers to such figures as being in (a) log-log coordinates vs. (b) log-linear coordinates for short.

Aside from the spectral features noted in the figure caption and discussed in greater detail by Ghil (1994), it is important to realize (i) that the large continuous background in Fig. 5a is purely of deterministically chaotic origin since there is no stochastic element whatsoever in the Le Treut and Ghil (1983) model or in its forcing and (ii) that the dominant peak at 109 kyr is not directly forced by the $f_4 = 1/100$ kyr$^{-1}$ eccentricity line but rather it is due to the difference tone between the two precessional frequencies, $f_1$ and $f_2$. Finally, it is the nonlinear, broad resonance of the model's $f_0$ frequency with the quasi-periodic forcing that produces the bump in the spectrum of Fig. 5a to the right of the orbital frequencies.

In returning to the “fundamental question no. 2” in Box 1, one must recall that on the paleoclimatic timescales of interest – apart from deterministic chaos à la Lorenz (1963), as obtained by Le Treut and Ghil (1983) and Le Treut et al. (1988) and shown here in Fig. 5a – stochastic contributions à la Hasselmann (1976) to the continuous part of the spectrum must also play an important role. In fact, RDS theory, as outlined in the next subsection, provides an excellent framework for a “grand unification” of these two complementary points of view (Ghil, 2014, 2019). In the paleoclimatic context, Ditlevsen et al. (2020) have suggested that, aside from red noise processes, dating uncertainties in the proxy records from which the spectra are derived may contribute, in all likelihood, to this background; see also Boers et al. (2017a, b). In this context, Verbitsky and Crucifix (2020) also provide a simple theory that addresses scaling properties in the glacial

cycles and their spectra, based on the so-called Buckingham $\pi$-theorem (e.g., Barenblatt, 1996).

## 3.2 Basic facts of NDS and RDS life

The highly preliminary results summarized in Sect. 3.1 encourage us to pursue the effects of the orbital forcing on intrinsic climatic variability in a more systematic way, effects that may have contributed to generate the rich paleoclimate spectrum on Quaternary and longer timescales (e.g., Westerhold et al., 2020). In fact, several research groups in the climate sciences have carried out during the past 2 decades an important extension of the dynamical systems and model hierarchy framework presented by Ghil and Childress (1987) and by Ghil (2001), from deterministically autonomous to non-autonomous and random dynamical systems (NDSs and RDSs, e.g., Ghil et al., 2008; Chekroun et al., 2011; Bódai and Tél, 2012).

On the road to including deterministically time-dependent and random effects, one needs to realize first that the climate system – as well as any of its subsystems on any timescale – is not closed: it exchanges energy, mass, and momentum with its surroundings, whether this pertains to other subsystems or the interplanetary space and the solid earth. The typical applications of dynamical systems theory to climate variability until not so long ago have only taken into account exchanges that are constant in time, thus keeping the model – whether governed by ordinary, partial, or other differential equations – autonomous; i.e., the models had coefficients and forcings that were constant in time.

Alternatively, the external forcing or the parameters were assumed to change either much more slowly than a model's internal variability, meaning that the changes could be assumed to be quasi-adiabatic, or much faster, meaning that they could be approximated by stochastic processes. Some of these issues are covered in much greater detail by Ghil and Lucarini (2020, Sect. III.G). The key concepts and tools of NDSs and RDSs go beyond such approaches that rely in an essential way on a scale separation between the characteristic times of the forcing and the internal variability of a given system; such a separation is rarely if ever actually present in the climate sciences.

The presentation of the key NDS and RDS concepts and tools in this subsection is aimed at as large a readership as possible and follows Ghil (2014). Slightly more in-depth but still fairly expository presentations can be found in Crauel and Kloeden (2015) and in Caraballo and Han (2017). Readers who are less interested in this mathematical framework – which allows a truly thorough understanding of the way that orbital forcing acts on intrinsic climate variability on Quaternary timescales – may skip at a first reading the remainder of this section and continue with Sect. 4.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Power spectra of simulated (a) deep-sea and (b) ice core records for the forced climate oscillator of Le Treut et al. (1988) with idealized orbital forcing. Panel (a) is in log–log coordinates and clearly shows a dominant peak at 109 kyr and a very large amount of variance in the continuous spectrum, which has a roughly −2 power law slope. The orbital peaks at $P_1 = 19$ kyr, $P_2 = 23$ kyr and $P_3 = 41$ kyr are also present, along with peaks at 14.7 and 10.4 kyr, which correspond to the sum tones $f_2 + f_3$ and $f_1 + f_2$, respectively. Panel (b) is in log–linear coordinates and shows additional harmonics and sum tones, as well as the difference tone $f_1 - f_2$, which corresponds to the dominant peak at 109 kyr. Inside the figure, the notation ka instead of kyr has been kept unchanged from the original publication. Adapted by permission from John Wiley & Sons Inc.: American Geophysical Union. Journal of Geophysical Research: Atmospheres. Isotopic Modeling of Climatic Oscillations: Implications for a Comparative Study of Marine and Ice Core Records by Le Treut et al. (1988), © 1988 by the American Geophysical Union. 1988.</figcaption></figure>

### 3.2.1 Autonomous and non-autonomous systems

Succinctly, one can write an autonomous system as

$$ \dot{X} = F(X; \mu), \qquad (8) $$

where $X$ stands for any state vector or climate field. While $F$ is a smooth function of $X$ and of the parameter $\mu$, it does not depend explicitly on time. This autonomous character of Eq. (8) greatly facilitates the analysis of its solutions' properties.

For instance, two distinct trajectories, $X_1(t)$ and $X_2(t)$, of a well-behaved, smooth autonomous system cannot intersect – i.e., they cannot pass through the same point in phase space – because of the uniqueness of solutions. This property helps one draw the phase portrait of an autonomous system, as does the fact that we only need to consider the behavior of solutions $X(t)$ as time $t$ tends to $+\infty$. The sets of points so obtained are (possibly multiple) equilibria, periodic and quasi-periodic solutions, and chaotic sets. In the language of dynamical systems theory, these are called, respectively: fixed points, limit cycles, tori, and strange attractors.

We know only too well, however, that the seasonal cycle plays a key role in climate variability on interannual timescales, while orbital forcing is crucial on the Quaternary timescales of many millennia. In addition, more recently it has become obvious that anthropogenic forcing is of utmost importance on the interdecadal timescales between these two extremes.

How can one take into account these types of time-dependent forcings and analyze the non-autonomous systems

that they lead us to formulate? One succinctly writes such a system as follows:

$$ \dot{X} = F(X, t; \mu). \qquad (9) $$

In Eq. (9), the dependence of $F$ on $t$ may be periodic, $F(X, t + P) = F(X, t)$, as in various El Niño–Southern Oscillation (ENSO) models, where the period $P = 12$ months, or monotonic, $|F|(X, t + \tau) \ge |F|(X, t)$, as it is when studying scenarios of anthropogenic climate forcing. An even more general situation includes time dependence in one or more parameters $\{\mu_1, \dots, \mu_p\}$.

To illustrate the fundamental distinction between an autonomous system like Eq. (8) and a non-autonomous one like Eq. (9), consider the simple scalar version of these two equations:

$$ \dot{X} = -\beta X, \qquad (10a) $$

and

$$ \dot{X} = -\beta X + \gamma t, \qquad (10b) $$

respectively. We assume that both systems are dissipative, i.e., $\beta > 0$, and that the forcing is monotonically increasing, $\gamma \ge 0$, as would be the case for anthropogenic forcing in the industrial era. Lorenz (1963) pointed out the key role of dissipativity in giving rise to strange but attracting solution behavior, while Ghil and Childress (1987) emphasized its importance and pervasive character in climate dynamics. Clearly the only attractor for the solutions of Eq. (10a), given
any initial point $X(0) = X_0$, is the stable fixed point $X_* = 0$, attained as $t \to +\infty$.

In the case of Eq. (10b), however, this forward-in-time approach yields blow-up from $t \to +\infty$ for any initial point. To make sense of what happens in the case of time-dependent forcing, one instead introduces the pullback approach, in which solutions are allowed to still depend on the time $t$ at which we observe them but also on a time $s$ from which the solution is started, with $X(s) = X_0$ and $s \ll t$. With this little change of approach, one can easily verify that

$$|X(s, t; X_0) - A_t| \to 0 \quad \text{as} \quad s \to -\infty, \qquad (11)$$

for all $t$ and $X_0$, where $A_t = (\gamma/\beta)(t-1/\beta)$. We obtain therewith, in this pullback sense, the intuitively obvious result that the solutions, if we start them far enough in the past, all approach the family of attracting sets $A_t$; this family follows the forcing $\gamma t$, and it thus has a linear growth in time $t$. Hence, the fixed point $X_*$ of Eq. (10b) is in fact a moving target and it is given by $X_* = \gamma t/\beta$. Due to the system's inertia, the set $A_t$ that is approached by the trajectories lags this time-dependent fixed point by a constant offset of $\gamma/\beta^2$.

### 3.2.2 Pullback attractor (PBA)

Formally, the indexed family $\mathcal{A}$ of all pullback attracting sets $A_t$,

$$\mathcal{A} = \{A_t\}_{t \in \mathbb{R}}, \qquad (12)$$

is termed the pullback attractor (PBA) of an NDS if the following two conditions are fulfilled:

i. each snapshot $A_t$ is compact and the family of snapshots $\{A(t)\}_{t \in \mathbb{R}}$ is invariant with respect to the dynamics,

$$X(t, s; X_0) \in A_t \quad \forall s \le t \quad \text{and} \quad X_0 \in A_s, \qquad (13)$$

and

ii. the pullback attraction occurs for all times,

$$\lim_{s \to -\infty} |X(t, s; X_0) - A_t| = 0 \quad \forall t. \qquad (14)$$

To further improve the reader's intuition for PBAs, we provide a second illustrative example here. A system defined in polar coordinates by

$$\dot{\rho} = \alpha(\mu - \rho), \quad \dot{\phi} = \omega, \qquad (15)$$

with $\rho, \mu \in \mathbb{R}^+$, and $\phi \in \mathbb{R}/2\pi$ can easily be seen to exhibit a limit cycle in the $(x, y)$ plane with $(x = \rho \cos\phi, y = \rho \sin\phi)$. An initial deviation of $\rho$ from $\mu$ will decay exponentially, and the system converges to an oscillation of radius $\mu$ with the angular velocity $\omega$. Here, we transform this

autonomous dynamical system into a non-autonomous one by modulating the target radius $\mu$ with a sinusoidal forcing,

$$\mu \to \mu(t) = \mu_0 + \beta \sin(vt), \qquad (16)$$

where the modulation is moderate to guarantee that $\mu_0 + \beta \sin(vt) > 0$ for all $t$.

Since the dynamics of the phase $\phi$ and of the radius $\rho$ are decoupled, the corresponding equations can be solved and analyzed separately. While the temporal development of the phase is trivial, the pullback-invariant attracting set of the radius for the initial condition $\rho(s) = \rho_0$ is given by

$$A^{(\rho)}(t; \rho_0) = \lim_{s \to -\infty} \rho(t, s; \rho_0) = \alpha\beta \sin(vt + \vartheta) + \mu_0, \qquad (17a)$$

with

$$\vartheta = \arctan(-v/\alpha), \qquad (17b)$$

as shown in Appendix B. Note that in the limit $s \to -\infty$, the dependence on the initial value $\rho_0$ vanishes and the attracting set $A_t^{(\rho)}$ performs an oscillation of the same frequency as the forcing. It lags the phase of the time-dependent fixed point $\mu(t)$ by the constant $\vartheta$, while its amplitude is amplified by the factor $\alpha$. Since $\rho$ is restricted to positive values, this solution requires $\alpha\beta < \mu_0$.

The PBA with respect to the coordinate $\rho$ is comprised of the family of all the sets $A_t^{(\rho)}$ as defined in Eq. (18) and thus reads

$$\mathcal{A}^{(\rho)} = \{\alpha\beta \sin(vt + \vartheta) + \mu_0\}_{t \in \mathbb{R}}. \qquad (18)$$

Since the pullback limit for the phase $\phi$ does not exist, no constraints on it other than $\phi \in [0, 2\pi)$ are imposed by the dynamics. Hence, for the system (15) comprised of radius and phase, we find that the distance of any trajectory at time $t$ - i.e. $(\rho(t; t_0, \rho_0), \phi(t; t_0, \phi_0))$ - to the set $A_t = \{(\alpha\beta \sin(vt + \vartheta) + \mu_0, \varphi) : \varphi \in [0, 2\pi)\}$ tends to zero as we pullback the initial time $t_0$ to $-\infty$. The pullback attracting sets $A_t$ at time $t$ are circles in the $(x, y)$ plane with a radius that oscillates in time, and the system's PBA is given by the family of these circles

$$\mathcal{A} = \{(\alpha\beta \sin(vt + \vartheta) + \mu_0, \varphi) : \varphi \in [0, 2\pi)\}_{t \in \mathbb{R}}. \qquad (19)$$

Figure 6 shows trajectories of the system starting from different points in the past. In Fig. 6a, the trajectories are depicted in the three-dimensional (3-D) space spanned by the two Cartesian coordinates $(x, y)$ and the time $t$, where the usual transformation from polar to Cartesian coordinates was applied. The shaded surface in this panel represents the PBA of the system. Corresponding trajectories of $\rho(t)$ and their convergence to the PBA $A_t^{(\rho)}$ are shown Fig. 6b. Figure 6f shows a heat map (Wilkinson and Friendly, 2009) that approximates a portion of the PBA's invariant measure projected onto the $(x, y)$ plane. For a clean definition of such
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Trajectories and PBA of the system defined by Eqs. (15) and (16); in all the panels the parameter values are $\mu_0 = 15$, $\alpha = 10$, $\nu = 0.2$, and $\omega = 0.3$, unless stated otherwise. (a) Trajectories $(\rho(t), \phi(t))$ of the system starting from different times in the past in the 3-D space spanned by $(x, y)$ and time $t$; the system's PBA lies on the red-shaded surface. (b) Temporal evolution of the radius (solid colors) of the three trajectories shown in (a) together with its PBA (dotted red). (c–e) Trajectories integrated from $t_0 = -100$ to $t_f = 200$ for $\nu = 0.3, 0.2, 0.5^{1/2}$ in panels (c), (d), and (e), respectively. The values are chosen such that the ratio between $\nu$ and $\omega$ is $\omega/\nu = 1$, $\omega/\nu \in \mathbb{Z}\setminus\{1\}$ and $\omega/\nu \in \mathbb{R}\setminus\mathbb{Z}$ in panels (c) to (e). Correspondingly, in (c) the trajectory quickly converges to a circle, whose center is slightly shifted from the origin. In panel (d) a quasi-closed three-fold loop can be observed, since $\omega/\nu = 3/2$. In panel (e) the trajectories will densely fill the annular disk defined in the video supplement statement. (f) Heat map of numerous trajectories projected onto the $(x, y)$ plane. The trajectories start at random points in state space and are integrated from $t_i = -200$ to $t_f = 200$. A video of the heat map filling up as more and more trajectories with different initial conditions are added is provided in the video supplement. The heat map shown here is a snapshot of the video taken at time $t = 0.20$; see details in the video supplement.</figcaption></figure>

a measure in NDSs and RDSs, please see Caraballo and Han (2017), Chekroun et al. (2011), Crauel and Kloeden (2015), or Ghil et al. (2008). Essentially, the heat map here counts the number of times that 100 trajectories integrated from $t = -200$ to $t = 200$, with randomized initial conditions like the ones shown in Fig. 6a, cross small pixels in the $(x, y)$ plane.

Figure 6c–e demonstrate a particularity of this system, which is characteristic of dynamics confined to a torus. Namely, the structure of the system’s trajectories depends on the frequency ratio $\omega/\nu$, and three different cases must be distinguished. If the radius is modulated with the same frequency as the oscillation itself, i.e. $\omega = \nu$ as in panel (c), the forcing and the system have a fixed phase relation. That is,

for a given phase of the system, its radius is always attracted by the same fixed point. Hence, the system practically repeats its orbit after a short time. More precisely, the radius of the oscillation does differ from one “round trip” around the torus to the next, but this difference tends to zero as $\rho(t)$ approaches the PBA $A_t^{(\rho)}$.

If $\omega$ and $\nu$ are rationally related, i.e., $m\omega = n\nu$ with $n, m \in \mathbb{N}$, as in Fig. 6d, then – after $n$ periods of the radial modulation and $m$ periods of the system’s oscillation – the phase relation between the system and its forcing will repeat itself, and hence we observe the same quasi-repetition of the orbit after the time $n 2\pi/\nu = m2\pi/\omega$. That is, such a trajectory will appear as an $n$-fold quasi-closed loop.
Finally, if $\omega/v \notin \mathbb{Z}$, as in Fig. 6e, then a given phase of the system will never coincide with the same phase of the radius modulation more than once. Hence, the trajectory does not repeat itself but instead densely covers the annular disk $D = \{(\rho, \phi) : \rho \in [\mu_0 - \alpha\beta, \mu_0 + \alpha\beta] \text{ and } \phi \in [0, 2\pi)\}$.

### 3.2.3 Random attractor

Let us return now to the more general, nonlinear case of Eq. (9) and add not only deterministic time dependence $F(X, t)$ but also random forcing,

$$dX = F(X, t)dt + G(X)d\eta, \quad (20)$$

where $\eta = \eta(t, \omega)$ represents a Wiener process – with $d\eta$ commonly referred to as “white noise” – and $\omega$ now labels the particular realization of this random process. When $G = \text{const.}$, the noise is additive, while for $\partial G/\partial X \neq 0$ we speak of multiplicative noise. The distinction between $d\eta$ and $d\eta$ in the stochastic differential Eq. (20) is necessary since, roughly speaking and following the Einstein (1905) paper on Brownian motion, it is the variance of a Wiener process that is proportional to time and thus $d\eta \propto (dt)^{1/2}$. In Eq. (20), we dropped the dependence on a parameter $\mu$ for the sake of simplicity.

The noise processes may include “weather” and volcanic eruptions when $X(t)$ is “climate”, thus generalizing the linear model of Hasselmann (1976), or cloud processes when dealing with the weather itself: one person’s signal is another person’s noise, as the saying goes. In the case of random forcing of Eq. (20), the concepts introduced by the simple deterministic examples of Eq. (10b) and Eqs. (15) and (16) above can be illustrated by the random attractor $\mathcal{A}(\omega)$ in Fig. 7.

A key feature of the pullback point of view on noise-perturbed dynamical systems that characterizes RDS theory is the use of a single noise realization, as opposed to the traditional, forward viewpoint of the Fokker-Planck equation and associated concepts, in which multiple noise realizations play a role. For a precise definition of a random attractor – as well as the commonalities and differences between the deterministic and random cases of time-dependent forcing – please see Caraballo and Han (2017).

Chekroun et al. (2011) studied a specific case of such a random attractor for the paradigmatic, climate-related Lorenz (1963) convection model. The authors introduced multiplicative noise into each of the ODEs of the original, deterministically chaotic system, as shown below:

$$dX = P_r(Y - X)dt + \sigma X d\eta, \quad (21a)$$

$$dY = (rX - Y - XZ)dt + \sigma Y d\eta, \quad (21b)$$

$$dZ = (-bZ + XY)dt + \sigma Z d\eta, \quad (21c)$$

where $r = 28$, $P_r = 10$, and $b = 8/3$ are the standard parameter values for chaotic behavior in the absence of noise and $\sigma$ is a constant variance of the Wiener process that is not

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Schematic diagram of a random attractor $\mathcal{A}(\omega)$ and of the pullback attraction to it; here $\omega$ labels the particular realization of the random process $\theta(t)\omega$ that drives the system. We illustrate the evolution in time $t$ of the random process $\theta(t)\omega$ (light solid black line at the bottom), the random attractor $\mathcal{A}(\omega)$ itself (yellow band in the middle), with the snapshots $A_0(\omega) = \mathcal{A}(\omega; t = 0)$ and $\mathcal{A}(\omega; t)$ (the two vertical sections, heavy solid lines), and the flow of an arbitrary compact set $B$ from “pullback times” $t = -\tau_2$ and $t = -\tau_1$ onto the attractor (heavy blue arrows). See Ghil et al. (2008, Appendix A) for the requisite properties of the random process $\theta(t)\omega$ that drives the RDS formulated by Eq. (20). Adapted from Physica D: Nonlinear Phenomena, 237, Ghil et al. (2008), Climate dynamics and fluid mechanics: natural variability and related uncertainties, 2111–2126, © 2008 Elsevier B.V. All rights reserved. 2008, with permission from Elsevier.</figcaption></figure>

necessarily small; both the noise realization $\eta(t)$ and $\sigma$ are the same in all three equations. The well-known strange attractor of the deterministic case is replaced by the Lorenz model’s random attractor, dubbed LORA by the authors. Four snapshots $A_t(\omega)$ of LORA are plotted in Fig. 8 and a video of its evolution in time $\mathcal{A}(\omega) = \{A_t(\omega)\}_{t \in \mathbb{R}}$ is available as supplementary material to Chekroun et al. (2011) at https://doi.org/10.1016/j.physd.2011.06.005.

Charó et al. (2021) have further analyzed the striking effects of the noise on the nonlinear dynamics that are visible in Fig. 8 and in the video of Chekroun et al. (2011) and gathered further insights into the abrupt changes of the snapshots’ topology at critical points in time. These remarkable changes suggest the possibility of random processes giving rise to qualitative jumps, such as the MPT, in paleoclimatic variability.

## 3.3 Application to Dansgaard-Oeschger (D-O) events

Before discussing conceptual glacial cycle models, we take a little detour and introduce a simpler – yet interesting and at the same time highly instructive – application of NDS theory to another important climate phenomenon. During past glacial periods, Greenland experienced a series of sudden decadal-scale warming events that left a clear trace in ice core records (Dansgaard et al., 1993). These so-called D-O events were followed by intervals of steady moderate cooling, be-
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Heat maps of the time-dependent invariant measure $v_t(\omega)$ supported on four snapshots $A_t(\omega)$ of LORA. The values of the parameters $r, s$, and $b$ are the classical ones, where the variance of the noise $\sigma = 0.5$. The color bar is on a log scale, and it quantifies the probability of landing in a particular region of phase space; a projection is shown of the 3-D phase space $(X, Y, Z)$ onto the $(X, Z)$ plane. Note the complex, interlaced filament structures between highly populated regions (in yellow) and moderately populated ones (in red). Reprinted from Physica D: Nonlinear Phenomena, 240, Chekroun et al. (2011), Stochastic climate dynamics: Random attractors and time-dependent invariant measures, 1685–1700, © 2011 Elsevier B.V. All rights reserved. 2011, with permission from Elsevier.</figcaption></figure>

fore a short phase of enhanced cooling brought the temperatures back to their pre-event levels (e.g., Rasmussen et al., 2014). This pattern is very clearly apparent in NGRIP $\delta^{18}$O records (North Greenland Ice Core Project members, 2004), and it can be qualitatively simulated by the fast component of a FitzHugh–Nagumo (FHN) model (FitzHugh, 1961; Nagumo et al., 1962), as pointed out by Mitsui and Crucifix (2017), among others; see also Kwasniok (2013), Rial and Yang (2007), Roberts and Saha (2017), and Vettoretti et al. (2022).

We discuss the example of the FHN model at some length in order to illustrate how external forcing can act on a system’s internal variability and thereby give rise to more complex dynamics. This model’s concise mathematical formulation and its widespread application in paleoclimate modeling and other fields make it ideally suited for this goal. We start with a description of the autonomous model with no time-dependent forcing. Subsequently, we introduce a simple sinusoidal forcing and numerically compute the corresponding PBA. We then extend these consideration into the realm

of random dynamical systems by adding stochastic forcing and discuss the resulting random attractor. Finally, we replace the synthetic forcings by one that corresponds to a paleoclimate proxy record of past CO$_2$ concentrations retrieved from Antarctic ice cores (Bereiter et al., 2015a) and show that this setup brings the model’s trajectories into good qualitative agreement with the D-O patterns observed in $\delta^{18}$O records from Greenland ice cores. In doing so, we pay less attention to the physical interpretation of the model’s variables, while focusing on the detailed explanation of model behavior and on the role of the forcing in the resulting dynamics.

### 3.3.1 The FitzHugh–Nagumo (FHN) model of fast–slow oscillations

The FHN model consists of two coupled ODEs that govern behavior alternating between slow evolutions and fast transitions. Typically, the timescales of the two variables are separated by introducing the parameters $\tau_x$ and $\tau_y$, with $x(t)$ being the slow component and $y(t)$ the fast one:

$$ \dot{x} = \frac{1}{\tau_x}(y - \gamma), \qquad (22a) $$

$$ \dot{y} = \frac{1}{\tau_y}[\alpha(y - y^3) - x]. \qquad (22b) $$

In order to develop an understanding for the way that such a model can simulate the rapid D-O warmings, followed by slow coolings, we start by discussing its autonomous behavior for time-independent $\gamma$.

First, consider the case of large timescale separation

$$ \tau_x \gg \tau_y. \qquad (23) $$

This choice guarantees that the fast $y$ component adjusts adiabatically to quasi-static changes of the slow $x$ component. The time derivative of $y(t)$, as shown in Fig. 10a, exhibits either three real roots or a single one, depending on the value of $x(t)$, which shifts the graph of the cubic polynomial $P_3(x, y) = \alpha(y - y^3) - x$ globally upwards or downwards. Of the three potential roots, the outer two are stable fixed points for $y(t)$ at a fixed value of $x$, while the inner one is unstable. Note that the two stable fixed points are always located either left or right of the local minimum or maximum of $P_3(x, y)$, respectively. Accordingly, we label them $y_l$ and $y_r$. The positions of the local extrema, namely $y_{min} = -\sqrt{1/3}$ and $y_{max} = \sqrt{1/3}$, provide an upper and a lower bound for the left and right stable fixed point, respectively.

Now, let us investigate the coupled dynamics of the slow and fast variables $x(t)$ and $y(t)$. Assume we are in a state where $x > 2\alpha/\sqrt{27}$ such that $y_l$ is the only root of $P_3(x, y)$. Provided that $|\gamma| < \sqrt{1/3}$, the time derivative $\dot{x}$ and $x$ itself will have opposite signs, with the consequence that a slow adjustment process of $x(t)$ sets in. This will shift the polynomial $P_3(x, y)$ upwards in Fig. 10a, which in turn revives the other two roots of $P_3(x, y)$. The fast $y(t)$ closely follows
$y_1$, which is shifted accordingly. Once $x(t) = -2\alpha/\sqrt{27}$, the root $y_1$ stops existing and a fast transition to the neighborhood of the opposite root $y_r$ ensues. This transition switches the sign of $\dot{x}(t)$ and the slow adjustment of $x(t)$ now happens in the opposite direction.

In the $(x, y)$ plane, this behavior manifests itself as a stable limit cycle. However, any value $|\gamma| > \sqrt{1/3}$ prevents $\dot{x}(t)$ from switching sign and therefore interrupts the cyclic destruction and revival of the respective opposite root of $P_3(x, y)$. Instead, $|\gamma| > \sqrt{1/3}$ gives rise to a single stable fixed point for the entire system in the $(x, y)$ plane and in this case both variables relax towards it.

In fact, in the autonomous setting, the system's qualitative behavior is controlled by the value of the parameter $\gamma$, which decides between internal oscillations along a limit cycle and relaxation towards a fixed point in the $(x, y)$ plane. Note that previously we referred to the roots $y_1$ and $y_r$ also as stable fixed points of $y$ for a given value of $x$. Here, the term stable fixed point refers to the entire system defined by the coupled ODEs (22a) and (22b). Both $\gamma = \pm\sqrt{1/3}$ are critical values of $\gamma$ that give rise to supercritical Hopf bifurcations of the coupled system's fixed points; recall Fig. 4 of Sect. 2.1.

This behavior can be better understood by considering the nullclines of Eqs. (22b) and (22a) in the $(x, y)$ plane, as shown in Fig. 9. If the branches of the $y$ nullcline that correspond to $y_1$ and $y_r$, and thus to stable fixed points of Eq. (22b) for a given value of $x$, intersect with the $x$ nullcline given by $y = \gamma$, then this intersection constitutes a stable fixed point for the entire system.

If they do not, the system first relaxes along the fast direction toward the $y$ nullcline. Only then does the adjustment of the slow component start to drag the system along the $y$ nullcline in the direction where the distance to the $x$ nullcline decreases. At the point where the $y$ nullcline reverses, the fast component is immediately attracted by the other branch of the fast nullcline and the same process starts all over again.

So far we have described the formation of the limit cycle in the FHN model under the assumption of clear timescale separation and the independence of the $x$ nullcline $\{y = \gamma\}$ from $x$. See Roccoreanu et al. (2000) for the emergence of the limit cycle in the more general FHN model.

The highly nonlinear, two-time behavior of the FHN model somewhat modifies the way that stable limit cycles arise in it. While we saw the oscillation's radius grow with the square root of the bifurcation parameter in the case of the normal form given by Eqs. (3)–(5), in the case of the FHN model, the radius of the oscillations actually grows exponentially over a small range of $\gamma$ values right after the bifurcation point and then stabilizes. These exponentially growing stable limit cycles have been termed “canard cycles” (Benoît, 1983). Roberts and Saha (2017) have pointed out a possible link between canard cycles and D-O cycles, and they play a role in other excitable climate models; see Pierini and Ghil (2022) and references therein.

### 3.3.2 A pullback attractor of a periodically forced FHN model

Introducing a sinusoidal time dependence

$$ \gamma \to \gamma(t) \propto \sin(t/\tau_f) \quad (24) $$

into the slow Eq. (22a) makes the system non-autonomous, as discussed in general terms in Sect. 3.2, and it periodically switches the self-oscillatory behavior on and off. We consider here the case in which the variations of the external forcing $\gamma(t; \tau_f)$ occur on a slower timescale than the entire internal dynamics, i.e.,

$$ \tau_f > \tau_x > \tau_y. \quad (25) $$

Note that we give up here on the restriction of strict timescale separation between $x(t)$ and $y(t)$, as expressed by Eq. (23).

The trajectories plotted in Fig. 10c–e and f–h represent the solutions of such a periodically forced FHN model starting at different points in the past. The time units are arbitrary and the two sets of solutions are for two different forcing timescales $\tau_f$. These trajectories illustrate the applicability of the pullback perspective suggested by Sect. 3.2 and Fig. 6 to the periodically forced FHN model. In contrast to the illustrative examples of Sect. 3.2, an analytical solution is not available in this case. However, this small sample of numerically computed trajectories, together with the phenomenological discussion herein, is quite sufficient for an intuitive understanding of the system's pullback behavior.

For Fig. 10c–e, the forcing's timescale is much slower than the internal timescales. The amplitude of the forcing $\gamma(t)$ was chosen such that it repeatedly crosses the two thresholds $\pm\sqrt{1/3}$ and thus induces a sequence of Hopf bifurcations by switching between intervals of self-sustained oscillation and attraction to a stable fixed point. Crossing such a bifurcation point due to slow changes in the forcing is referred to as a bifurcation-induced tipping or “B-tipping” (Ashwin et al., 2012; Ghil, 2019; Ghil and Lucarini, 2020).

Strikingly, all trajectories converge to one another during non-oscillatory time intervals, when they are simultaneously attracted by the single existing fixed point. During oscillatory intervals, phase differences between individual trajectories may, in principle, persist. Still, convergence during a single non-oscillatory interval is so strong that after it the trajectories can no longer be discriminated visually. Numerically, however, the distance between trajectories only tends to zero but never reaches it. At the end of non-oscillatory intervals, the trajectories always re-enter the oscillatory regime from the same location in the $(x, y)$ plane – to within negligible numerical differences – and hence very nearly repeat themselves. Qualitatively speaking, the PBA $\mathcal{A}$ – i.e., the family of invariant snapshots $\{\mathcal{A}(t)\}_{t \in \mathbb{R}}$ of Eq. (13) – is an infinite repetition of the common trajectory structure that can be observed in Fig. 10d, e between –5000 and 15 000 time units. In the case at hand, each snapshot $\mathcal{A}(t)$ consists of a single point.
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Nullclines of the autonomous FHN model governed by Eq. (22) with large timescale separation, for $\alpha = 2$, $\tau_x = 100$, and $\tau_y = 10$; the nullclines of the fast y component are in blue, and those of the slow x component are in orange. (a–c) One illustrative trajectory (red dotted) for $\gamma = 0, 0.5$, and $1.0$, respectively. The upper branch and the lower branch of the y nullcline (solid blue) correspond to the roots $y_r$ and $y_l$, respectively, as defined in the text. For fixed values of x, they constitute stable fixed points for Eq. (22b). The middle branch (dashed blue) corresponds to the unstable fixed point for Eq. (22b) for a given x. The trajectories for panels (a) and (b) follow a limit cycle, as the x nullcline does not intersect with either of the stable branches of the y nullcline, given that $|\gamma| < \sqrt{1/3}$, while the trajectory in panel (c) approaches a single stable fixed point for the coupled system formed by the intersection of the x nullcline with the $y_l$ branch of the y nullcline.</figcaption></figure>

For Fig. 10f–h, the timescale separation between the forcing and the internal dynamics is reduced, resulting in a qualitatively different behavior of the non-autonomous system. The frequency of occurrence of B-tipping points is much higher, and hence the trajectories do not even execute a full oscillation during a single time interval that permits oscillations. As a result, two stable patterns of trajectories are formed. These two patterns can be brought into agreement by switching the sign of one pattern and shifting it in time by $\tau_f/2$. This symmetry reflects the symmetry of the stable nullcline of the fast system component, as shown in Fig. 9.

Again, the PBA of this non-autonomous system can be thought of as an infinite repetition of the common trajectory structure that can be observed in Fig. 10g, h between –5000 and 15 000 time units. In contrast to the slow-forcing case, each snapshot $A(t)$ now is comprised of two points in the $(x, y)$ plane. This example illustrates how the action of an external force on an autonomous system can give rise to considerably richer dynamics, which crucially depends on both the system’s internal variability and the nature of the forcing.

### 3.3.3 A random attractor of a periodically and stochastically forced FHN model

Based on the brief introduction to RDSs in Sect. 3.2, we take our investigation of the periodically forced FHN model one step further and include a random component into the external forcing, acting on the model’s fast y component:

$$dx = \frac{1}{\tau_x}(y - \gamma(t))dt, \quad (26a)$$

$$dy = \frac{1}{\tau_y} \left[ \alpha (y - y^3) - x \right] dt + \sigma d\eta. \quad (26b)$$

Here, $\eta$ denotes a Wiener process, as in Eq. (20) – i.e., a continuous stochastic process whose increments $d\eta$ are inde-

pendently and normally distributed, with mean zero and unit variance – and $\gamma(t)$ remains a periodic forcing of the x component proportional to $\sin(\tau_f t)$.

In order to study the random attractor of this system, we compute trajectories with random initial conditions over a time span long enough to reveal the asymptotic behavior, as shown in Fig. 11. For both the slower and faster deterministic forcing – i.e., $\tau_f = 1000$ and $\tau_f = 350$, respectively, as studied in the previous paragraphs – several random attractors are approximated for increasing noise variance values $\sigma$ in Eq. (26b). For each attractor approximation, we use 20 trajectories with the same noise realization. Each random attractor (red) is shown together with the corresponding PBA of the FHN system subject to purely periodic forcing (blue and green).

For the long periodic forcing with $\tau_f = 1000$, the trajectories in Fig. 11a converge fairly rapidly to a single one, as in the case with no noise in Fig. 10d, e. Furthermore, there is a clear similarity of pattern and proximity in phase between the PBA of the deterministic system and the random attractor. However, the deviations of the random attractor from the PBA increase in both pattern and phase as the noise variance increases. These deviations are most striking during the oscillatory intervals because the noise can induce phase shifts in the oscillations, which then persist for the duration of the oscillatory interval. During non-oscillatory intervals, the random attractor is less susceptible to the noise because the resulting perturbations decay in the presence of a stable fixed point of the underlying deterministic system.

For the shorter period forcing with $\tau_f = 350$, the PBA of the deterministic system features two stable branches, which could in principle persist in the random attractor. As for the case of $\tau_f = 1000$, a rapid convergence of the trajectories can be observed for all noise variances. However, two separate
<figure><img src="image_12.png" /><figcaption>**Figure 10.** FitzHugh-Nagumo (FHN) model with parameters $τ_x = 100$, $τ_y = 60$, and $α = 2$. (a) The cubic polynomial $P_3(x, y)$ of the fast derivative as a function of $y$ for $x = 0$ (solid blue line); the red lines point to the local maximal and minimal values of $P_3(x, y)$, namely $±2a/\sqrt{27}$, respectively – these are the maximal values by which $P_3$ can be shifted up or down, while maintaining all of its three roots; the dotted gray lines indicate the shifted function with $x = ±2a/\sqrt{27}$. The purple lines labeled $y_{min}$ and $y_{max}$ mark the right and left boundaries for the roots $y_1$ and $y_r$, respectively: $y_1$ and $y_r$ can never be located in between the two purple lines. (b) Trajectories of the non-autonomous model with $γ(t) = 0.8 \sin(t/τ_f)$ and $τ_f = 1000$, plotted in the $(x, y)$ phase plane; the trajectories are colored by their starting times $\{t_0 = -20000, t_1 = -16000, t_2 = -13000, t_3 = -7000\}$, and the initial positions were drawn from a standard Gaussian bivariate distribution. (c) The slow time-dependent forcing $γ(t) = 0.8 \sin(t/1000)$. (d, e) The same trajectories as in (b) but plotted in time as $y = y(t)$ and $x = x(t)$, respectively. Panels (f–h) are same as panels (c)–(e) but for the fast time-dependent forcing $γ(t) = 0.8 \sin(t/350)$. The gray shading in panels (c)–(h) indicates intervals during which $|γ| > \sqrt{1/3}$, and the internal oscillation is hence suppressed.</figcaption></figure>

bundles of trajectories persist, each associated with one of the two separate branches of the deterministic PBA. Only after some time do the two bundles merge and subsequently follow either one of the two branches of the PBA or the other. This phenomenon is called noise-induced synchronization (e.g., Arnold, 1998; Chekroun et al., 2011). Which of the two PBA branches the random attractor follows depends on the exact noise realization. The random attractor may also switch irregularly from one branch of the PBA to the other (not shown). Already in the simple setup explored here, one notices that the higher the noise variance, the faster the synchronization of the trajectories. This statement must be understood probabilistically: it may certainly happen for two given noise realizations that the random attractor with the

higher noise level takes longer to synchronize, as is the case for $σ = 0.7$ and $σ = 0.8$ in Fig. 11b.

The investigation carried out herein assesses a very special case of an FHN model’s random attractor. Random attractors of FHN-type models have been studied intensively (e.g., Wang, 2009; Yamakou et al., 2019). Most studies so far, however, have concentrated on the excitable regime of the stochastically forced FHN model, i.e., the deterministic model’s parameters were chosen so as to exhibit a single stable fixed point in the absence of the noise. The noise variance was then chosen sufficiently high to ensure that random fluctuations can push the system out of equilibrium and allow it to take a round trip on what would be the deterministic model’s limit cycle, given different parameter values.
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Random attractor of the periodically and stochastically forced FHN model governed by Eq. (26). (a) Results for slow periodic forcing with period $\tau_f = 1000$. The top graph shows the periodic forcing, with non-oscillatory intervals marked by gray shading. The next three graphs show the fast y component of the approximate random attractors, as per Eq. (26b); the noise variance $\sigma$ increases from top to bottom. Each random attractor is approximated by integrating 20 trajectories (solid red) with different initial conditions over time and using the same Wiener process as their common stochastic forcing; the corresponding deterministic PBA is shown in blue. The random attractors and the PBA are very similar for small noise variances, but they differ more and more as the noise variance increases. (b) Results for faster periodic forcing with period $\tau_f = 350$. The forcing is only shown for the first 8000 time units, with non-oscillatory intervals again shaded gray. The left part of panel (b) shows five approximate random attractors, computed as in panel (a), on a common time axis. The panel's right-hand side shows their continuations on individual time axes in order to display the moment when full noise-induced synchronization of the trajectories takes place. Prior to this point, the random attractor is split into two branches, which closely follow the PBA of the deterministic system (blue and green). For $\sigma = 0.7$, the synchronization already takes place during the first 8000 time units.</figcaption></figure>

The situation studied here is rather different, as shown across Figs. 9–11.

### 3.3.4 An FHN model of the NGRIP record

Readers who are familiar with the NGRIP $\delta^{18}$O record (North Greenland Ice Core Project members, 2004) might have already realized the qualitative resemblance between the proxy data and the fast component's trajectory of the periodically forced FHN model in Fig. 10d. In particular, the prominent sawtooth pattern of the data is satisfactorily captured by the fast and slow dynamics of the model.

Figure 12 shows a trajectory of the FHN model for which the sinusoidal forcing used in Fig. 10 was replaced by a rescaled time series of atmospheric $CO_2$ concentrations retrieved from Antarctic ice cores (Bereiter et al., 2015a):

$$ \gamma(t) \propto CO_2(t). \qquad (27) $$

Is it remarkable how well this simple forcing brings the oscillatory intervals of the FHN model into agreement with the time intervals of the record that are dominated by D-O cycles without any systematic tuning of the model parameters. Clearly, the $CO_2$-forced FHN model fails to reproduce the exact waiting times between D-O events. However, with these waiting times being at least in part stochastically determined (Ditlevsen et al., 2007), the purely deterministic FHN model is not meant to reproduce the exact pattern of D-O events. Vettoretti et al. (2022) have carried out a detailed
<figure><img src="image_14.png" /></figure>

**Figure 12.** FHN model fit to the NGRIP $\delta^{18}$O dataset of 20-year means, as published as a supplement to Seierstad et al. (2014) (see also Rasmussen et al., 2014). Originally, the North Greenland Ice Core Project members (2004) reported a 50-year-mean dataset. (a) The fast component $y_{FHN}(t)$ of an FHN model forced with historical CO$_2$ concentrations (orange), together with the observed $\delta^{18}$O record (blue) from the NGRIP ice core (Seierstad et al., 2014). (b) Rescaled atmospheric CO$_2$ concentrations from Antarctic ice cores in arbitrary units (au) (Bereiter et al., 2015a); the horizontal red lines indicate the lower and upper bounds, $y_{min} = -\sqrt{1/3}$ and $y_{max} = \sqrt{1/3}$, respectively, of the free-oscillation regime.

study on the use of a CO$_2$-forced FHN-type model to simulate D-O variability.

In fact, Rousseau et al. (2022, Fig. 6) describe in detail a somewhat more complex, proxy-record-based picture of the interaction between episodes full of D-O events, (Heinrich, 1988) events, and longer-term cooling trends. It is quite possible that a simple model like the one in this section but including explicitly continental ice sheets could capture such a detailed picture.

In the present framework, the FHN model's fast variable $y(t)$ may be interpreted as the intensity of the Atlantic meridional overturning circulation (AMOC), which switches between on and off states during self-oscillatory behavior; see, for instance, Henry et al. (2016), Ghil (1994, Table 5), and Ghil and Lucarini (2020, Table I). The slower $x$ variable that drives the transition between the on and off states of the AMOC may then be taken, for instance, as the waxing and waning of Northern Hemisphere ice sheets (e.g., Ghil et al., 1987), linked in turn to varying ice shelf extent (e.g., Boers et al., 2018), or as the weakening and strengthening of Antarctic bottom water production (Vettoretti et al., 2022).

The interaction between the fast variable and the slow one happens here in the presence of a climate forcing represented by CO$_2$ concentration. On the much slower timescales of Quaternary glaciations, an interplay between the CO$_2$ concentration and mean global temperature might also occur, as we shall see in the next section.

# 4 An NDS for the Quaternary glaciations

Apparently, it was Crucifix (2013) who first applied pullback ideas to the problems of Quaternary glaciations, independently of earlier work on the topic in the climate literature (Ghil et al., 2008; Chekroun et al., 2011; Bódai and Tél, 2012). His work concentrated mainly on the connection between the pacemaking role of the orbital forcing and the observed irregularity of the glacial terminations during the Late Pleistocene; cf. Broecker and Van Donk (1970) and Ghil and Childress (1987, Fig. 11.2).

Based on the considerable success of NDS and RDS applications to other climate problems – such as ENSO (Ghil et al., 2008; Ghil and Zaliapin, 2015; Chekroun et al., 2018; Marangio et al., 2019), the wind-driven ocean circulation (Pierini et al., 2016, 2018), or the evaluation of the ensemble simulations routinely performed in support of the Assessment Reports of the Intergovernmental Panel on Climate Change (Drótos et al., 2015; Vissio et al., 2020) – it would appear worthwhile to proceed further along these lines.

Box 2 summarizes open questions with respect to Quaternary glaciations whose investigation should be aided by NDS and RDS theory. Of course, each of these problems requires one or more distinct climate models, as well as very careful modeling of the kinds of time-dependent changes in forcing and parameters that are most enlightening, as well as most relevant and plausible. A good way would be to start testing ideas with relatively simple models and pursue the investigation systematically across a hierarchy of models – through intermediate ones and on to the most detailed ones – in order to further increase understanding of the climate system and
The concepts and tools summarized in Sect. 3.2 could shed further light on several problems that not yet quite solved concerning Quaternary glaciations and paleoclimate in general. A short and highly incomplete list of the interesting problems that might gain by such a treatment is given below.

1. What are the causes of the MPT, and what are the contributions of intrinsic climate variability and orbital forcing therein? See discussion in Sec. 2.2 for more details.

2. What causes the large continuous background of the observed spectrum: is it deterministic chaos, stochastic forcing, or both? See discussion in the caption of Fig. 4 for more details.

3. Particularly complex problems are associated with Heinrich and D-O events and their modulation over longer time intervals by orbital forcing and other slow-acting changes in atmospheric composition, geology, ocean circulation, and so on. See, for instance, Rousseau et al. (2021, and references therein) for more details.

4. It is well established by now that a snow-ball mode occurred several times in Earth's geologic history, but we still know fairly little about how the planet got into, and especially out of, such a mode (e.g., Pierrehumbert, 2004).

**Box 2.** Some open questions concerning Quaternary glaciations.

of its predictability on the various paleoclimatic time scales mentioned in the list above.

Such an approach can usefully complement the more common one of merely pushing onwards to higher and higher model resolution in order to achieve ever more detailed simulations of the system's behavior for a limited set of semi-empirical parameter values. Ghil (2001) and Held (2005), among others, have emphasized the need to pursue such a model hierarchy, as originally proposed by Schneider and Dickinson (1974), in order to balance the need to broaden the number of plausible hypotheses vs. the need for confronting them with spatiotemporal details derived from observations.

In this section, we illustrate how the PBA concept can help shed more light upon the dynamics of ice age models. As pointed out in Sect. 3.1 and elsewhere in this paper, there is a long history of modeling the climate of the Quaternary by means of conceptual models, and many non-autonomous models have been proposed to simulate glacial-interglacial cycles of the last 400 kyr to 2.6 Myr based on the orbital forcing. In Appendix A, we provide a long but still not exhaustive list of glacial-cycle models and specify some of their key characteristics, including the degree of their success at simulating the MPT; see also the discussion in Sect. 2.2.

Among these glacial-cycle models, the model of Daruka and Ditlevsen (2016, DD16 hereafter) belongs to the more abstract ones, as it is not derived from detailed physical considerations. Still, its concise form, interesting nonlinear dynamics, and ability to simulate glacial cycles, as well as the MPT, make the DD16 model well suited for our illustrative purposes. We first modify this model slightly from its original formulation. We do so mainly in order for the model to better approximate the benthic $\delta^{18}$O proxy reconstruction of glacial-interglacial cycles due to Lisiecki and Raymo (2005a), especially the timing of glacial terminations; compare our Fig. 13 with Fig. 1 in DD16. Thereafter, we compute the PBAs of the modified DD16 model (hereafter M-DD16) to investigate the dynamical stability of its glacial cycles over the past 2.6 Myr.

Our model's variables, following DD16, are a global temperature anomaly $y$ that is proportional to minus the global ice volume and an effective climatic memory term $x$ that represents the internal degrees of freedom. In the deterministic case, the governing equations of the M-DD16 model are given by

$$ \tau \dot{x} = \lambda y, \qquad (28a) $$

$$ \tau \dot{y} = -\alpha(t) + x - x^3 - \beta(t)F(t)x - \kappa y; \qquad (28b) $$

here $t$ is the time (in kyr) and $F(t)$ is the normalized 21 June insolation at 65° N, based on the calculations of Laskar et al. (2004a), as shown in Fig. 13a. The constant parameter values are chosen as $\kappa = 1$, $\tau = 100$, and $\lambda = 10$. Note that with this choice of $\lambda$ and unlike in the FHN model of D-O oscillations in Sect. 3.3, $x$ is the fast variable and $y$ is the slow one.

In the original DD16 model, MPT-like behavior was produced by a slow sigmoid variation of the parameter $\kappa$ in Eq. (28b),

$$ \kappa(t) = \kappa_1 + 0.5(\kappa_0 - \kappa_1) \left( 1.0 - \tanh \left( \frac{t - t_0}{t_s} \right) \right). \qquad (29) $$

In our M-DD16 model, we introduce instead a slow change in the parameters $\alpha(t)$ and $\beta(t)$ of Eq. (28b) as follows:

$$ \alpha(t) = 2.1 - 1.4 \tanh((t + 1100)/500), \qquad (30a) $$

$$ \beta(t) = 2.5 + 1.4 \tanh((t + 1100)/500). \qquad (30b) $$

The functions $\alpha(t)$ and $\beta(t)$ so defined are plotted in Fig. 13b, and they induce, as we shall see forthwith, a change in model behavior that not only resembles the MPT but also shows correct timings for most of the terminations. Moreover, to simulate $\delta^{18}$O$_{model}$, we add a linear trend to the slow variable $y$ to mimic the overall cooling at time scales of millions of years, and thus $\delta^{18}$O$_{model} = 4.3 - 1.4y + 0.0003t$. Equation (28) generalizes a first-order ODE system that is equivalent to the Duffing form of the nonlinear spring Eq (7) discussed in Sect. 2.2. The use of such classical first-order systems – like the Duffing and Van der Pol ones – in paleoclimate models was initiated by Saltzman et al. (1981); see also Sect. 2.2 herein.
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Glacial-interglacial cycles simulated by the modified Daruka and Ditlevsen model of Eqs. (28) and (30): (a) 21 June insolation $F(t)$ at 65° N, normalized to have mean zero and unit standard deviation over the last 1000 kyr, according to Laskar et al. (2004a). (b) Slowly changing parameters $\alpha(t)$ and $\beta(t)$ introduced to give rise to the MPT. (c) Simulated glacial-interglacial cycles $\delta^{18}\text{O}_{\text{model}}$ (red) in comparison with the benthic $\delta^{18}\text{O}$ data (blue); the latter dataset, provided by Lisiecki and Raymo (2005a), is also known as the LR04 stack.</figcaption></figure>

Figure 13c shows a time series of simulated glacial-interglacial changes $\delta^{18}\text{O}_{\text{model}}$ (red) in comparison with the benthic $\delta^{18}\text{O}$ (blue) of Lisiecki and Raymo (2005a). The model's initial condition is taken to be $x = -1$ and $y = 0$ at $t = 10000$ kyr b2k and, since the insolation forcing in Fig. 13a is prescribed as a time series with 100-year sampling, we solved Eq. (28) using Heun's predictor-corrector method (Isaacson and Keller, 2012, chap. 8) with a step size of 100 years. A large spin-up time is chosen to guarantee that transient effects caused by the initial conditions have abated by the year 2.6 Ma b2k, which is the starting point for the time interval under study. The correlation between the model simulation and the proxy record is 0.75 for the time interval from 2600 to 0 kyr b2k and 0.72 over the interval from 1000 to 0 kyr b2k. Varying the parameters slowly across the time interval of interest, as shown in Fig. 13b, leads to a change in the frequency – from a dominant 41 kyr periodicity prior to the MPT, at roughly 1.2 Ma kyr b2k, to a dominant 100 kyr periodicity after the MPT, at roughly 800 kyr b2k – and a substantial increase in the amplitude.

We next approximate the PBA by taking 40 random initial conditions at 10 Ma b2k and integrating the model of Eqs. (28) and (30) up to the present time. The PBA in this case is simply a moving fixed point, as plotted in Fig. 14a, since the model dynamics is predominantly stable in the long

time interval prior to the MPT. It would thus appear that the orbital forcing simply moves this fixed point around and fully determines Earth's climate. This agrees with the clear statement in DD16 that “first and foremost, our model does not have any internal periods of oscillation”.

However, when keeping the parameters $\alpha$ and $\beta$ fixed at their post-MPT values $\alpha = 0.7$ and $\beta = 3.9$ throughout the simulation interval and repeating the computation of the PBA, a more complex picture arises. In the latter case, Fig. 14b shows a bunching of trajectories into separate fuzzy clusters, subject to the quasi-periodic orbital forcing of Fig. 13a.

There are two interesting inferences to be drawn. First, post-MPT dynamics is much more irregular and unstable than the more stable dynamics prior to the MPT. The robustness of the 40 kyr glacial cycles and instability of 100 kyr glacial cycles against perturbations is in line with the conclusions of previous studies (Mitsui et al., 2015; Quinn et al., 2018). It also appears to be consistent with the Willeit et al. (2019) simulation of the last 3 Myr of Earth's history that used an Earth system model of intermediate complexity (CLIMBER-2) that included ice sheets and a carbon cycle, along with atmosphere-ocean dynamics. In their work, trajectories starting from different initial states tended to converge to a single attracting trajectory in the Early Pleistocene,
<figure><img src="image_16.png" /></figure>

**Figure 14.** PBA of the M-DD16 model governed by Eq. (28), approximated by 40 trajectories starting from random initial conditions in $(x, y) \in [-2, 2] \times [-2, 2]$ at $t = -10000$ kyr. (a) The PBA corresponding to Fig. 13, with the slowly changing parameters $\alpha(t)$ and $\beta(t)$ given by Eq. (30). (b) The PBA for $\alpha(t)$ and $\beta(t)$ kept constant at the post-MPT values of $\alpha = 0.7$ and $\beta = 3.9$. In panel (b), the PBA is shown over a shorter time interval of the last 1000 kyr so that the detailed structure is more clearly visible. Without changes in $\alpha$ and $\beta$ in time, the overall structure of the PBA is similar before and after 1000 kyr b2k.

while several distinct trajectories survived in the Late Pleistocene after the MPT.

Second, the separate bundles or “ropes” of trajectories in Fig. 14b seem to point to the type of generalized synchronization discussed in the paleoclimate context by De Saedeleer et al. (2013) and in the context of interannual and interdecadal climate variability by Pierini and Ghil (2022) and Vannitsem et al. (2021). Generalized synchronization in the strict sense of the existence of a map between a time-dependent control and the system’s asymptotic behavior has only been shown to hold for non-chaotic systems. Work is under way, however, to further generalize this concept to chaotic systems as well (e.g., Rulkov et al., 1995; Zhang et al., 2007).

# 5 Conclusions

In this review and research paper, we have covered the contributions of the 1970s to the rebirth of the Milankovitch (1920) theory of the ice ages in Sect. 1 and the 1980s advances in modeling the Quaternary climate's intrinsic variability in Sect. 2. In Sect. 3, we presented first results on the effects of

the orbital insolation forcing of the data discussed in Sect. 1 on the intrinsic variability of the data discussed in Sect. 2, and proceeded to introduce the novel concepts and tools of the theory of non-autonomous and random dynamical systems (NDSs and RDSs) that can help to better model and understand these effects. Section 3 concluded with the formulation and study of a FitzHugh-Nagumo (FHN)-type model of recurrent Dansgaard-Oeschger (D-O) events, in which historical CO₂ concentrations induced episodes of D-O events alternating with episodes of their absence in excellent qualitative agreement with NGRIP δ¹⁸O data; see Fig. 12.

Finally, in Sect. 4, we listed a number of open issues on Quaternary and longer paleoclimate timescales and proposed to address them by using the tools of Sect. 3.2. This approach was illustrated by a Duffing-type model of Daruka and Ditlevsen (2016), modified to include slow changes in the parameters that mimic such changes in the Earth system over the duration of the Quaternary period.

When the parameters are gradually changed in time so as to exhibit the mid-Pleistocene transition (MPT), the PBA is simply a moving fixed point. However, when the parameters are fixed at their post-MPT values, the PBA so obtained is chaotic and exhibits clusters of trajectories that we termed ropes. This suggests (a) that the stability of the system is gradually lost while crossing the MPT and (b) that the Late Pleistocene climate, albeit chaotic, may well be subject to a kind of generalized synchronization (cf. De Saedeleer et al., 2013; Pierini and Ghil, 2022; Vannitsem et al., 2021) with the orbital forcing that is illustrated in Fig. 5 of Sect. 3.1 herein. In the specific situation at hand, separate ropes may be associated with various combination tones of the forcing frequencies.

In a broader perspective – and leaving aside various finer points of the MPT conversation outlined in Sect. 2.2 – one can see the work that was reviewed and extended in this paper as a confirmation of the fine intuition of Emiliani and Geiss (1959), 6 decades ago, as summarized in and further expanded by Ghil and Childress (1987, Sect. 12, pp. 446–447).

Hence the following scenario (compare Emiliani and Geiss, 1959) suggests itself for the successive climatic transitions from Pliocene to Pleistocene and from Early to Late Pleistocene: As land masses moved towards more northerly positions, small ice caps formed on mountain chains and at high latitudes. These ice caps, due to their feedback on albedo, made climate more sensitive to insolation variations than it was in the total absence of ice. The response of the climatic system to such variations during the Early Pleistocene (2000 [kyr]–1000 [kyr] ago) was still relatively weak, of a fraction of a degree centigrade in global temperature perhaps, in agreement with the quasi-equilibrium results of Sect. 10.2.
As ice caps passed, about 1000 [kyr] ago, a certain critical size, the unforced system jumped from its stable equilibrium to its stable limit-cycle state (Figures 12.5 and 12.9), increasing dramatically the climate's total variability, to a few degrees centigrade in global temperature. Furthermore, resonant response became possible (see also Oerlemans (1984) [in Berger et al. (1984)] and Sergin (1979)), enhancing abruptly the amplitude of the peak at 100 kyr, among others.

The take-home message is that slow and fast processes, both intrinsic and extrinsic, interact on all paleoclimatic timescales and that we are mastering the art of modeling such interactions.

## Appendix A: Low-order dynamical-system models of glacial cycles

The dynamical modeling of glacial cycles dates back to the 1970s. Calder (1974) proposed a model of global ice volume changes that had different sensitivities to the insolation when ice sheets were waxing and waning, respectively. His model can be written as an NDS, according to Paillard (2001). Subsequently, conceptual dynamical models were further developed (e.g., Imbrie and Imbrie, 1980; Berger, 1999). Some of the more recent models simulate the proxy records of glacial cycles remarkably well (Paillard, 1998; Imbrie et al., 2011; Parrenin and Paillard, 2012).

Shortly after Calder (1974) presented his work, Weertman (1976) proposed a simple ice sheet model based on the flow law of a perfectly plastic solid. Next, researchers extended this simple ice sheet model by coupling it with an energy balance model (Källén et al., 1979) and further with the isostatic response of the underlying bedrock (Oerlemans, 1980; Ghil and Le Treut, 1981; Le Treut and Ghil, 1983; Pollard, 1983). Källén et al. (1979) found self-sustained oscillations in their simple coupled ice sheet-energy balance model. Le Treut and Ghil (1983) showed that the dominant 100 kyr periodicity of glacial cycles is generated – in their simple oscillator coupling ice sheet volume with the bedrock's isostatic rebound, on the one hand, and with the atmosphere and ocean's energy balance, on the other hand – via nonlinear resonance with the multi-periodic orbital forcing. More recently, Verbitsky et al. (2018) developed a simple physical model through a scaling argument that respects the underlying physics. Another branch of simple models explicitly includes the carbon cycle as an essential ingredient (Saltzman and Maasch, 1990; Paillard and Parrenin, 2004; Hogg, 2008; Toggweiler, 2008; Omta et al., 2016; Talento and Ganopolski, 2021).

A deeper understanding of glaciation cycles cannot be obtained without process-based models that focus on the detailed physics and biogeochemical phenomena involved (Berger et al., 1999; Ganopolski and Calov, 2011; Abe-Ouchi et al., 2013; Ganopolski and Brovkin, 2017; Willeit et al.,

2019). Still, simple dynamical models like the ones mentioned above, as well as models based on more mathematical considerations, are also useful for understanding the climate system's behavior and changes therein, since complex systems can sometimes exhibit familiar dynamics, regardless of the details (Nicolis and Nicolis, 2012; Crucifix, 2011). As Henri Poincaré pointed out, “mathematics is the art of giving the same name to different things”.

For example, coupled nonlinear oscillators frequently exhibit synchronization with simple frequency ratios, either with each other (Pikovsky et al., 2001) or with the forcing (Ghil and Childress, 1987). Thus, mathematical models that ignore many physical details may also help us elucidate emergent properties in paleoclimatic dynamics (Tziperman et al., 2006; Crucifix, 2011, 2012; De Saedeleer et al., 2013; Mitsui et al., 2015; Ashwin and Ditlevsen, 2015; Daruka and Ditlevsen, 2016) or test different orbital hypotheses (Huybers, 2011). Many low-order models of glacial dynamics are listed in Table A1, although it is by no means exhaustive.

## Appendix B: PBA for a limit cycle with sinusoidally modulated radius

Here we study the system of two formally decoupled ODEs

$$ 
\begin{aligned}
\frac{d\rho}{dt} &= \alpha(\mu + \beta \sin(vt) - \rho), \\
\frac{d\phi}{dt} &= \omega \quad \text{with} \quad \rho > 0, \mu > 0,
\end{aligned}
\qquad (B1) $$

that were introduced in Sect. 3.2 and analytically derive their invariant sets

$$ 
\begin{aligned}
\mathcal{A}(t) &= \{ (\alpha\beta \sin(vt + \vartheta) + \mu, \varphi) : \varphi \in [0, 2\pi) \} \\
\forall t \in \mathbb{R},
\end{aligned}
\qquad (B2) $$

as well as the corresponding pullback attractor (PBA). Following Crauel and Kloeden (2015), the PBA is given by the family

$$ \mathcal{A} = \{\mathcal{A}\}_{t:t \in \mathbb{R}}. \qquad (B3) $$

First, we define $\Delta\rho(t) = \rho(t) - \mu$, which gives rise to

$$ \frac{d\Delta\rho(t)}{dt} = -\alpha\Delta\rho(t) + \alpha\beta\sin(vt). \qquad (B4) $$

This is an inhomogeneous ODE and can thus be solved by the variation of parameters method (e.g., Boyce and DiPrima, 2005). The ansatz

$$ \Delta\rho(t) = c(t) e^{-\alpha(t-t_0)}, \qquad (B5) $$

yields

$$ \frac{d\Delta\rho(t)}{dt} = -\alpha\Delta\rho(t) + \frac{dc(t)}{dt} e^{-\alpha(t-t_0)}. \qquad (B6) $$
<table><caption>**Table A1.** List of simple conceptual glacial-cycle models with only 1–3 variables. Note that, while extensive, this list is not exhaustive.</caption><thead><tr><th>Author(s) (year)</th><th>Succinct description</th><th>Dynamical properties</th><th>Is the mid-Pleistocene transition (MPT)-like behavior generated?</th><th>Additional comments</th></tr></thead><tbody><tr><td colspan="5">Low-order dynamical systems models of glacial cycles</td></tr><tr><td>Calder (1974)</td><td>Model of the global ice volume with different sensitivities to the insolation during growth and retreat, respectively.</td><td>Non-chaotic. No self-sustained oscillation in the absence of forcing</td><td>Not specifically addressed.</td><td>His model can be written as a non-autonomous dynamical system (NDS) (Paillard, 2001).</td></tr><tr><td>Weertman (1976)</td><td>Simplified ice sheet model based on the flow law of a perfectly plastic solid.</td><td>Non-chaotic. No self-sustained oscillation in the absence of forcing.</td><td>Not specifically addressed.</td><td></td></tr><tr><td>Källén et al. (1979), Ghil and Le Treut (1981), Le Treut and Ghil (1983), Ghil (1994)</td><td>Ice sheet-bedrock-temperature coupled oscillator; see the text for details.</td><td>Chaotic in the presence of orbital forcing (Ghil, 1994). Exhibits self-sustained oscillations with a period of ≈ 10 kyr; the basic mechanism for getting ~ 100 kyr cycle is nonlinear resonance at a difference tone between the 19 and 23 kyr precession cycles (Le Treut and Ghil, 1983); see Sect. 2.1 herein.</td><td>Not explicitly. It can exhibit, however, the MPT in principle, since 41 kyr oscillations and dominant 100 kyr oscillations can be generated depending on the model's parameter values.</td><td>Källén et al. (1979) is an ice sheet-temperature coupled oscillator. Ghil and Le Treut (1981) introduced the bedrock and Le Treut and Ghil (1983) added the orbital forcing.</td></tr><tr><td>Imbrie and Imbrie (1980)</td><td>Piecewise linear model with two different timescales for ice sheet waxing and waning.</td><td>Non-chaotic. No self-sustained oscillation in the absence of forcing.</td><td>Not specifically addressed.</td><td>400 kyr periodicity is more dominant than 100 kyr periodicity.</td></tr><tr><td>Saltzman and Maasch (1988, 1990), Maasch and Saltzman (1990)</td><td>Coupled oscillator based on ice-CO2-ocean temperature (or NADW) coupling. Nonlinearity is only in the CO2 component.</td><td>Non-chaotic. After the MPT, it exhibits self-sustained oscillations with a period of ~ 100 kyr in the absence of forcing.</td><td>Yes, the MPT arises via a Hopf bifurcation in the underlying system; see also Cruci-fix (2012).</td><td>It has a strange non-chaotic attractor, as well as a chaotic one, depending on the parameter setting (Mitsui and Aihara, 2014).</td></tr><tr><td>Paillard (1998), Parrenin and Paillard (2003, 2012)</td><td>Hybrid dynamical system with discrete states that switch when conditions are satisfied.</td><td>Non-chaotic but the dynamics is sensitive to parameter changes near the switching boundaries (Paillard, 2001). No self-sustained oscillations in the absence of forcing.</td><td>Yes, in Paillard (1998).</td><td>The authors call these models relaxation oscillators, in spite of their discrete-state formulation.</td></tr><tr><td>Berger (1999)</td><td>Delayed differential equation for ice volume,</td><td>Robustness of trajectories against random perturbations is mentioned. It exhibits damped oscillations in the absence of forcing.</td><td>Not specifically addressed.</td><td></td></tr><tr><td>Paillard and Parrenin (2004)</td><td>Coupled oscillator based on Northern Hemisphere ice volume-Antarctic ice extent-CO2 coupling.</td><td>Non-chaotic. It exhibits self-sustained oscillations in the absence of forcing.</td><td>Yes, the MPT is induced by a slow drift in the bottom water formation efficiency around the Antarctic.</td><td>Ashwin et al. (2018) showed that this model can exhibit chaotic dynamics when it is slightly modified.</td></tr><tr><td>Ashkenazy and Tziperman (2004), Tziperman et al. (2006)</td><td>An ice mass model that simplifies the sea ice switch model of Gildor and Tziperman (2000).</td><td>It exhibits self-sustained oscillations in the absence of forcing. Non-chaotic under the orbital forcing.</td><td>Yes, in Ashkenazy and Tziperman (2004), where the maximum ice volume threshold is increased in time according to the regolith hypothesis.</td><td>Uses the temperature-precipitation feedback introduced by Källén et al. (1979).</td></tr><tr><td>Huybers (2007, 2011)</td><td>The ice mass grows monotonically and collapses to zero when it exceeds a threshold modulated by the obliquity (Huybers, 2007) or by a hybrid measure consisting of the obliquity and climatic precession (Huybers, 2011).</td><td>Non-chaotic. It exhibits self-sustained oscillations in the absence of forcing.</td><td>Yes, in Huybers (2007).</td><td></td></tr></tbody></table>
<table><caption>**Table A1. Continued.**</caption><thead><tr><td>Author(s) (year)</td><td>Succinct description</td><td>Dynamical properties</td><td>Is the mid-Pleistocene transition (MPT)-like behavior generated?</td><td>Additional comments</td></tr></thead><tbody><tr><td>Imbrie et al. (2011)</td><td>Data-based, phase space model for Pleistocene ice volume with thresholds.</td><td>Due to the presence of the threshold, the model is sensitive to changes in parameters or in its position in phase space.</td><td>Yes, MPT-like behavior is produced solely by changes in orbital parameters.</td><td></td></tr><tr><td>Crucifix (2012), De Saedeleer et al. (2013)</td><td>Forced Van der Pol (VdP) oscillator model; closely related to the FHN model used in Sect. 3.3 herein.</td><td>Non-chaotic (Crucifix, 2012) but sensitive to the noise. It exhibits self-sustained oscillations in the absence of forcing.</td><td>Yes, MPT-like behavior is generated via a Hopf bifurcation with an explosive character (Crucifix, 2012).</td><td>Ashwin et al. (2018) showed that this model has only a small parameter region corresponding to chaotic dynamics but may have a wider chaotic region when it is generalized to the VdP-Duffing system.</td></tr><tr><td>Mitsui et al. (2015)</td><td>One-dimensional phase oscillator model.</td><td>Non-chaotic.</td><td>Yes, MPT-like frequency change accompanies a smooth or non-smooth saddle-node bifurcation of tori.</td><td>It has a strange non-chaotic attractor or a quasi-periodic attractor in a classical sense, i.e., not in the pullback sense.</td></tr><tr><td>Ashwin and Ditlevsen (2015)</td><td>Two-dimensional forced limit cycle oscillator.</td><td>Non-chaotic; overall stability of simulated glacial cycles against dynamical noise is reported.</td><td>It exhibits MPT-like behavior via a trans-critical bifurcation of the slow manifold in the fast dynamics.</td><td></td></tr><tr><td>Omta et al. (2016)</td><td>Forced two-dimensional oscillator consisting of ocean alkalinity and calcifier population.</td><td>Non-chaotic or chaotic given periodic forcing, depending on the parameters. It exhibits self-sustained oscillations in the absence of forcing.</td><td>Yes.</td><td></td></tr><tr><td>Daruka and Ditlevsen (2016)</td><td>Forced Duffing oscillator-type model.</td><td>No self-sustained oscillations in the absence of forcing. Can be chaotic depending on parameter values; the authors call this the "climatic butterfly effect".</td><td>Yes, MPT-like behavior is induced by a slow change in the damping coefficient parameter κ.</td><td>See also Sect. 4 herein.</td></tr><tr><td>Huybers and Langmuir (2017)</td><td>Coupled system of ice volume, temperature, and atmospheric CO₂ incorporating a delayed CO₂ contribution from ocean ridge volcanism. It is reduced to a forced 1-D delay differential equation.</td><td>A phase-locking property is reported. The underlying system has two stable and one unstable stationary states. However, it is close to a Hopf bifurcation point and is excitable by the forcing.</td><td>Yes, the MPT is modeled as a switch from small-amplitude oscillations to large-amplitude ones, which is triggered by the amplitude modulation of obliquity cycles.</td><td></td></tr><tr><td>Quinn et al. (2018)</td><td>Scalar delay differential equation for ice volume, derived from the Saltzman and Maasch (1988) model via the linear chain approximation.</td><td>Non-chaotic before the MPT and temporarily chaotic after the MPT around 800 kyr b2k. In the absence of astronomical forcing, the delayed feedback leads to bistable behavior, in which stable large-amplitude oscillations and an equilibrium coexist.</td><td>Yes, MPT-like behavior is induced by summer insolation forcing as a transition from small-amplitude ~41 kyr cycles to large-amplitude ~100 kyr cycles.</td><td></td></tr><tr><td>Ashwin et al. (2018)</td><td>Study of several low-order dynamical systems models, including Paillard and Parrenin (2004), Crucifix (2012), and a generalized VdP-Duffing oscillator.</td><td>Chaotic or non-chaotic depending on the model and its parameter values.</td><td>Not explicitly discussed. The models can, however, exhibit MPT-like behavior in principle because 41 kyr oscillations and dominant 100 kyr oscillations are generated depending on parameter values.</td><td>The authors stress the possibility of chaotic dynamics occurring.</td></tr></tbody></table>
<table><caption>**Table A1. Continued.**</caption><thead><tr><th>Author(s) (year)</th><th>Succinct description</th><th>Dynamical properties</th><th>Is the mid-Pleistocene transition (MPT)-like behavior generated?</th><th>Additional comments</th></tr></thead><tbody><tr><td>Verbitsky et al. (2018)</td><td>Coupled model relying on ice sheet basal temperature-ocean temperature coupling, deduced from physical laws via a scaling analysis.</td><td>Non-chaotic. No self-sustained oscillations in the absence of forcing. The 100 kyr cycles are attributed to a period-doubling response to 41 kyr obliquity cycles.</td><td>Yes, MPT-like behavior is induced by an enhancement of positive feedbacks against negative feedbacks in the model.</td><td></td></tr><tr><td>Talento and Ganopolski (2021)</td><td>Ice volume-CO<sub>2</sub> coupled conceptual model.</td><td></td><td>Not specifically addressed.</td><td>This model was developed to simulate the last 800 kyr of glacial cycles and Earth's future climate on the million-year timescale.</td></tr><tr><td colspan="5">Low-order dynamical models in which stochastic processes are of the essence</td></tr><tr><td>Benzi et al. (1981, 1982), Nicolis (1981)</td><td>Two independent formulations of a stochastic-resonance model. Additive noise favors the bimodal response of a periodically modified double-well potential to global annual insolation variations.</td><td>Additive stochastic perturbations play a key role.</td><td>Not specifically addressed. An MPT-like transition could be obtained by subjecting both the depth of the deterministic double well and the pure periodicity of its modulation to a more or less gradual modification.</td><td>Benzi et al. (1981, 1982) and Nicolis (1981) proposed the stochastic-resonance idea for ice age cyclicity independently of each other.</td></tr><tr><td>Matteucci (1989)</td><td>Generalized stochastic resonance model with obliquity and precessional forcing.</td><td>Additive stochastic perturbations play a key role.</td><td>Not specifically addressed.</td><td></td></tr><tr><td>Pelletier (2003)</td><td>Model based on the coherence resonance mechanism (Pikovsky and Kurths, 1997) uses a single temperature ODE with a discrete delay.</td><td>Additive stochastic perturbations play a key role.</td><td>Yes.</td><td>The delay feedback in the temperature is based on ice sheet extent reconstructions and affects the model's albedo.</td></tr><tr><td>Ditlevsen (2010)</td><td>Generalized stochastic resonance model with obliquity and precessional forcing.</td><td>Additive stochastic perturbations play a key role.</td><td>Yes.</td><td>The model has an underlying bifurcation structure following Paillard (1998).</td></tr></tbody></table>

A comparison with Eq. (B4) requires

$$ \frac{d}{dt}c(t) = \alpha\beta \sin(vt)e^{+\alpha(t-t_0)}, \quad (B7) $$

and hence

$$ \begin{aligned} c(t) &= \int_{t_0}^{t} \alpha\beta \sin(vt')e^{+\alpha(t'-t_0)}dt' + \gamma \\ &= \left[ \alpha\beta \sin(vt') \frac{1}{\alpha} e^{+\alpha(t'-t_0)} \right]_{t_0}^{t} \\ &\quad - \int_{t_0}^{t} \frac{\nu}{\alpha} \alpha\beta \cos(vt')e^{+\alpha(t'-t_0)}dt' + \gamma. \end{aligned} $$

Repeated partial integration yields

$$ \begin{aligned} \int_{t_0}^{t} \frac{\nu}{\alpha} \alpha\beta \cos(vt')e^{+\alpha(t'-t_0)}dt' = & \\ & \left[ \frac{\nu}{\alpha} \alpha\beta \cos(vt') \frac{1}{\alpha} e^{+\alpha(t'-t_0)} \right]_{t_0}^{t} \\ & + \int_{t_0}^{t} \frac{\nu^2}{\alpha^2} \alpha\beta \sin(vt')e^{+\alpha(t'-t_0)}dt'. \end{aligned} \quad (B9) $$

Therefore, we find

$$ \begin{aligned} \left(1 + \frac{\nu^2}{\alpha^2}\right) \int_{t_0}^{t} \alpha\beta \sin(vt')e^{+\alpha(t'-t_0)}dt' = & \\ & \left[ \alpha\beta \sin(vt') \frac{1}{\alpha} e^{+\alpha(t'-t_0)} \right]_{t_0}^{t} \\ & - \left[ \frac{\nu}{\alpha} \alpha\beta \cos(vt') \frac{1}{\alpha} e^{+\alpha(t'-t_0)} \right]_{t_0}^{t}, \end{aligned} \quad (B10) $$
and finally

$$c(t) = \frac{1}{\left(1 + \frac{v^2}{\alpha^2}\right)} \left( \underbrace{\left[ \frac{\beta \sin(vt')e^{+\alpha(t'-t_0)}}{\beta[\sin(vt)e^{\alpha(t-t_0)} - \sin(vt_0)]} \right]_{t_0}}_{=\beta[\sin(vt)e^{\alpha(t-t_0)} - \sin(vt_0)]} \right. \\ - \underbrace{\left[ \frac{v\beta}{\alpha} \cos(vt')e^{+\alpha(t'-t_0)} \right]_{t_0}}_{=\frac{v\beta}{\alpha}[\cos(vt)e^{\alpha(t-t_0)} - \cos(vt_0)]} + \gamma. \quad (B11) \\ = \frac{v\beta}{\alpha} \left[ \cos(vt)e^{\alpha(t-t_0)} - \cos(vt_0) \right]$$

Plugging this result into the ansatz (Eq. B5) yields

$$\Delta\rho(t, t_0) = \frac{1}{\left(1 + \frac{v^2}{\alpha^2}\right)} \left( \beta \left[ \sin(vt) - \sin(vt_0)e^{-\alpha(t-t_0)} \right] \right. \\ \left. - \frac{v\beta}{\alpha} \left[ \cos(vt) - \cos(vt_0)e^{-\alpha(t-t_0)} \right] \right) \\ + \gamma e^{-\alpha(t-t_0)}, \quad (B12)$$

with the initial conditions

$$\Delta\rho(t_0, t_0) = \gamma. \quad (B13)$$

In the pullback limit, all the terms that carry a factor $e^{-\alpha(t-t_0)}$ vanish, and thus

$$\lim_{t_0 \to -\infty} \Delta\rho(t, t_0) = \frac{1}{\left(1 + \frac{v^2}{\alpha^2}\right)} \left( \beta \sin(vt) - \frac{v\beta}{\alpha} \cos(vt) \right) \\ = \alpha\beta \sin(vt + \vartheta), \quad (B14)$$

with $\vartheta = \arctan(-v/\alpha)$. For comparison, the modulation of the target radius itself was given by $\beta \sin(vt)$, and hence it is amplified by the factor of $\alpha$. Since $\rho$ is restricted to positive values, this solution requires $\alpha\beta < \mu$.

Since the evolution in time of the phase $\phi(t)$ is trivial, different initial conditions for the phase do not converge. Hence, the time-dependent sets that are invariant with respect to the dynamics of the system are

$$\mathcal{A}(t) = \{(\rho(t), \phi) : \phi \in [0, 2\pi)\} \\ = \{(\alpha\beta \sin(vt + \vartheta), \phi) : \phi \in [0, 2\pi)\}. \quad (B15)$$

Defined as the indexed family of all $\mathcal{A}(t)$, the system's PBA is comprised of the family of circles

$$\mathcal{A} = \{A(t)\}_{t \in \mathbb{R}} \\ = \{(\alpha\beta \sin(vt + \vartheta), \phi) : \phi \in [0, 2\pi)\}_{t \in \mathbb{R}}. \quad (B16)$$

**Code and data availability.** All code used to generate the figures presented in this article is available from the authors upon

request. The NGRIP $\delta^{18}$O – originally published in North Greenland Ice Core Project members (2004) – and the historical CO$_2$ data (Seierstad et al., 2014) shown in Fig. 12 are available as a supplement to Seierstad et al. (2014) or https://www.iceandclimate.nbi.ku.dk/data/ (last access: 13 April 2022) and from https://www.ncei.noaa.gov/access/paleo-search/study/17975 (last access: 13 April 2022, Bereiter et al., 2015b), respectively. Theenthic $\delta^{18}$O data (Lisiecki and Raymo, 2005a) shown in Fig. 13 can be obtained from https://doi.org/10.1594/PANGAEA.704257 (Lisiecki and Raymo, 2005b), and the 21 June insolation at 65° N (Laskar et al., 2004a) can be downloaded from http://vo.imcce.fr/insola/earth/online/earth/online/index.php (last access: 13 April 2022, Laskar et al., 2004b).

**Video supplement.** The video supplement to this article (https://doi.org/10.5281/zenodo.6346211, Riechers, 2022) illustrates the pullback attractor (PBA) associated with the simple system governed by Eqs. (15) and (16). It shows a heat map of the phase plane, derived from an increasing number of trajectories with common initial time $t_1 = -200$ and final time $t_f = 200$. The initial radius and phase of each trajectory are randomly sampled from Gaussian distributions centered at $\mu_\rho = 20$ and $\mu_\phi = 0$ with standard deviations of $\sigma_\rho = 5$ and $\sigma_\phi = 10$, respectively. Over the course of the video, 100 trajectories are continuously added to the heat map and the annular disk $D = \{(\rho, \phi) : \rho \in [\mu - \alpha\beta, \mu + \alpha\beta] \text{ and } \phi \in [0, 2\pi)\}$ fills up. The heat map in Fig. 6b is a snapshot from this video at time $t = 0.2$.

**Author contributions.** MG conceived and designed the study. KR and TM carried out the major part of the article's new research in close interaction with MG and NB. All authors interpreted and discussed the results and wrote the manuscript.

**Competing interests.** The contact author has declared that neither they nor their co-authors have any competing interests.

**Disclaimer.** Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Special issue statement.** This article is part of the special issue "A century of Milankovic's theory of climate changes: achievements and challenges (NPG/CP inter-journal SI)". It is a result of the conference "One Hundred Years of Milankovic's Theory of Climate Changes: synergy of the achievements and challenges of the next century", 17–18 November 2020.

**Acknowledgements.** We thank Andreas Groth for helpful comments on an earlier version of this manuscript. It is a pleasure to thank Tamás Bódai and two anonymous referees for their thorough and constructive comments. István Daruka's public comments in part motivated the addition of Appendix A and its Table A1, in order
to give a broader perspective of relevant work on the Quaternary's glacial cycles and the mid-Pleistocene transition (MPT). Takahito Mitsui and Niklas Boers acknowledge funding by the Volkswagen Foundation. The present work is TiPES contribution no. 52; the TiPES (Tipping Points in the Earth System) project has received funding from the European Union's Horizon 2020 research and innovation program under grant agreement no. 820970. Michael Ghil acknowledges support by the EIT Climate-KIC; EIT Climate-KIC is supported by the European Institute of Innovation & Technology (EIT), a body of the European Union.

**Financial support.** This research has been supported by the Horizon 2020 research and innovation program under grant agreement no. 820970 (TiPES), the Volkswagen Foundation, and the European Institute of Innovation & Technology via the EIT Climate-KIC.

The publication of this article was funded by the Open Access Fund of the Leibniz Association.

**Review statement.** This paper was edited by Marie-France Loutre and reviewed by Tamas Bodai and two anonymous referees.

## References

Abe-Ouchi, A., Saito, F., Kawamura, K., Raymo, M. E., Okuno, J., Takahashi, K., and Blatter, H.: Insolation-driven 100 000-year glacial cycles and hysteresis of ice-sheet volume, Nature, 500, 190–193, https://doi.org/10.1038/nature12374, 2013.

Arnold, L.: Random Dynamical Systems, 1st edn., Springer, Berlin, Heidelberg, https://doi.org/10.1007/978-3-662-12878-7, 1998.

Arnold, V.: Geometrical Methods in the Theory of Ordinary Differential Equations, Springer, New York, NY, https://doi.org/10.1007/978-1-4612-1037-5, 1988.

Ashkenazy, Y. and Tziperman, E.: Are the 41 kyr glacial oscillations a linear response to Milankovitch forcing?, Quaternary Sci. Rev., 23, 1879–1890, https://doi.org/10.1016/j.quascirev.2004.04.008, 2004.

Ashwin, P. and Ditlevsen, P.: The middle Pleistocene transition as a generic bifurcation on a slow manifold, Clim. Dynam., 45, 2683–2695, 2015.

Ashwin, P., Wieczorek, S., Vitolo, R., and Cox, P.: Tipping points in open systems: bifurcation, noise-induced and rate-dependent examples in the climate system, Philos. T. Roy. Soc. A, 370, 1166–1184, 2012.

Ashwin, P., David Camp, C., and von der Heydt, A. S.: Chaotic and non-chaotic response to quasiperiodic forcing: limits to predictability of ice ages paced by Milankovitch forcing, Dynamics and Statistics of the Climate System, 3, 1–20, https://doi.org/10.1093/climsys/dzy002, 2018.

Bagniewski, W., Ghil, M., and Rousseau, D. D.: Automatic detection of abrupt transitions in paleoclimate records, Chaos, 31, 113129, https://doi.org/10.1063/5.0062543, 2021.

Barenblatt, G. I.: Scaling, Self-similarity, and Intermediate Asymptotics: Dimensional Analysis and Inter-

mediate Asymptotics, Cambridge University Press, https://doi.org/10.1017/CBO9781107050242, 1996.

Benoit, E.: Systèmes lents-rapides dans R 3 et leurs canards, in: IIle rencontre de géométrie du Schnepfenried Volume 2 – 10–15 mai 1982, no. 109-110 in Astérisque, Société mathématique de France, http://www.numdam.org/item/AST_1983__109-110__159_0/ (last access: 13 April 2022), 1983.

Benzi, R., Sutera, A., and Vulpiani, A.: The mechanism of stochastic resonance, J. Phys. A-Math. Den., 14, L453, https://doi.org/10.1088/0305-4470/14/11/006, 1981.

Benzi, R., Parisi, G., Sutera, A., and Vulpiani, A.: Stochastic resonance in climatic change, Tellus, 34, 10–15, https://doi.org/10.3402/tellusa.v34i1.10782, 1982.

Bereiter, B., Eggleston, S., Schmitt, J., Nehrbass-Ahles, C., Stocker, T. F., Fischer, H., Kipfstuhl, S., and Chappellaz, J.: Revision of the EPICA Dome C CO₂ record from 800 to 600 kyr before present, Geophys. Res. Lett., 42, 542–549, 2015a.

Bereiter, B., Eggleston, S., Schmitt, J., Nehrbass-Ahles, C., Stocker, T. F., Fischer, H., Kipfstuhl, S., and Chappellaz, J. A.: Antarctic Ice Cores Revised 800KYr CO₂ Data, National Centers for Environmental Information, NESDIS, NOAA, U.S. Department of Commerce [data set], https://www.ncei.noaa.gov/access/paleo-search/study/17975 (last access: 13 April 2022), 2015b.

Berger, A.: Long-term variations of daily insolation and Quaternary climatic changes, J. Atmos. Sci., 35, 2362–2367, 1978.

Berger, A., Imbrie, J., Hays, J., Kukla, G., and Saltzman, B. (Eds.): Milankovitch and Climate: Understanding the Response to Astronomical Forcing, NATO ASI series. Series C, Mathematical and physical sciences; vol. 126, pts. 1–2, Springer, Dordrecht, https://doi.org/10.1007/978-94-017-4841-4, 1984.

Berger, A., Li, X., and Loutre, M.-F.: Modelling northern hemisphere ice volume over the last 3 Ma, Quaternary Sci. Rev., 18, 1–11, 1999.

Berger, W. H.: The 100 kyr ice-age cycle: Internal oscillation or inclinational forcing?, Int. J. Earth Sci., 88, 305–316, https://doi.org/10.1007/s005310050266, 1999.

Bódai, T. and Tél, T.: Annual variability in a conceptual climate model: Snapshot attractors, hysteresis in extreme events, and climate sensitivity, Chaos, 22, 023110, https://doi.org/10.1063/1.3697984, 2012.

Bódai, T., Lucarini, V., Lunkeit, F., and Boschi, R.: Global instability in the Ghil-Sellers model, Clim. Dynam., 44, 3361–3381, 2015.

Boers, N., Chekroun, M. D., Liu, H., Kondrashov, D., Rousseau, D.-D., Svensson, A., Bigler, M., and Ghil, M.: Inverse stochastic-dynamic models for high-resolution Greenland ice core records, Earth Syst. Dynam., 8, 1171–1190, https://doi.org/10.5194/esd-8-1171-2017, 2017a.

Boers, N., Goswami, B., and Ghil, M.: A complete representation of uncertainties in layer-counted paleoclimatic archives, Clim. Past, 13, 1169–1180, https://doi.org/10.5194/cp-13-1169-2017, 2017b.

Boers, N., Ghil, M., and Rousseau, D.-D.: Ocean circulation, ice shelf, and sea ice interactions explain Dansgaard-Oeschger cycles, P. Natl. Acad. Sci. USA, 115, E11005–E11014, https://doi.org/10.1073/pnas.1802573115, 2018.

Bond, G., Heinricht, H., Broecker, W., Labeyrie, L., Mcmanus, J., Andrews, J., Huonll, S., Jantschik, R., Clasen, S., Simet, C., Tedesco, K., Klas, M., Bonanitt, G., and Ivy, S.: Evi-
dence for massive discharges of icebergs into the North Atlantic ocean during the last glacial period, Nature, 360, 1668–1672, https://doi.org/10.1038/360245a0, 1992.

Bond, G., Broecker, W., Johnsen, S., McManus, J., Labeyrie, L., Jouzel, J., and Bonani, G.: Correlations between climate records from North Atlantic sediments and Greenland ice, Nature, 365, 143–147, https://doi.org/10.1038/365143a0, 1993.

Bond, G., Showers, W., Cheseby, M., Peter Almasi, R. L., deMenocal, P., Priore, P., Irka Hajdas, H. C., and Bonani, G.: A pervasive millennial-scale cycle in North Atlantic Holocene and glacial climates, Science, 278, 1257–1266, https://doi.org/10.1126/science.278.5341.1257, 1997.

Boyce, W. E. and DiPrima, R. C.: Elementary Differential Equations and Boundary Value Problems, 8th edn., John Wiley & Sons, ISBN 978-0-470-38334-6, 2005.

Broecker, W. S. and Van Donk, J.: Insolation changes, ice volumes, and the O¹⁸ record in deep-sea cores, Rev. Geophys., 8, 169–198, 1970.

Budyko, M. I.: The effect of solar radiation variations on the climate of the Earth, Tellus, 21, 611–619, 1969.

Calder, N.: Arithmetic of ice ages, Nature, 252, 216–218, https://doi.org/10.1038/252216a0, 1974.

Caraballo, T. and Han, X.: Applied Nonautonomous and Random Dynamical Systems: Applied Dynamical Systems, Springer, Cham, https://doi.org/10.1007/978-3-319-49247-6, 2017.

Charó, G. D., Chekroun, M. D., Sciamarella, D., and Ghil, M.: Noise-driven topological changes in chaotic dynamics, arXiv [preprint], arXiv:2010.09611v7, 2 August 2021.

Chekroun, M. D., Simonnet, E., and Ghil, M.: Stochastic climate dynamics: random attractors and time-dependent invariant measures, Physica D, 240, 1685–1700, https://doi.org/10.1016/j.physd.2011.06.005, 2011.

Chekroun, M. D., Simonnet, E., and Ghil, M.: Stochastic climate dynamics: random attractors and time-dependent invariant measures, Physica D, 240, 1685–1700, https://doi.org/10.1016/j.physd.2011.06.005, 2011.

Chekroun, M. D., Ghil, M., and Neelin, J. D.: Pullback attractor crisis in a delay differential ENSO model, in: Advances in Nonlinear Geosciences, edited by: Tsonis, A. A., Springer Science & Business Media, 1–33, https://doi.org/10.1007/978-3-319-58895-7, 2018.

Crafoord, C. and Källén, E.: A note on the condition for existence of more than one steady state solution in Budyko-Sellers type models, J. Atmos. Sci., 35, 1123–1125, 1978.

Crauel, H. and Kloeden, P. E.: Nonautonomous and random attractors, Jahresbericht der Deutschen Mathematiker-Vereinigung, 117, 173–206, 2015.

Crucifix, M.: How can a glacial inception be predicted?, Holocene, 21, 831–842, 2011.

Crucifix, M.: Oscillators and relaxation phenomena in Pleistocene climate theory, Philos. T. Roy. Soc. A, 370, 1140–1165, 2012.

Crucifix, M.: Why could ice ages be unpredictable?, Clim. Past, 9, 2253–2267, https://doi.org/10.5194/cp-9-2253-2013, 2013.

Dansgaard, W., Johnsen, S. J., Clausen, H. B., Dahl-Jensen, D., Gundestrup, N. S., Hammer, C. U., Hvidberg, C. S., Steffensen, J. P., Sveinbjörnsdottir, A. E., Jouzel, J., and Bond, G.: Evidence for general instability of past climate from a 250 kyr ice-core record, Nature, 364, 218–220, https://doi.org/10.1038/364218a0, 1993.

Daruka, I. and Ditlevsen, P. D.: A conceptual model for glacial cycles and the middle Pleistocene transition, Clim. Dynam., 46, 29–40, 2016.

De Saedeleer, B., Crucifix, M., and Wieczorek, S.: Is the astronomical forcing a reliable and unique pacemaker for climate? a conceptual model study, Clim. Dynam., 40, 273–294, 2013.

Ditlevsen, P., Mitsui, T., and Crucifix, M.: Crossover and peaks in the Pleistocene climate spectrum; understanding from simple ice age models, Clim. Dynam., 54, 1801–1818, 2020.

Ditlevsen, P. D.: Extension of stochastic resonance in the dynamics of ice ages, Chem. Phys., 375, 403–409, https://doi.org/10.1016/j.chemphys.2010.05.022, 2010.

Ditlevsen, P. D. and Ashwin, P. B.: Complex climate response to astronomical forcing: The middle-Pleistocene transition in glacial cycles and changes in frequency locking, AIP Conf. Proc., 6, 62, https://doi.org/10.3389/fphy.2018.00062, 2018.

Ditlevsen, P. D., Andersen, K. K., and Svensson, A.: The DO-climate events are probably noise induced: statistical investigation of the claimed 1470 years cycle, Clim. Past, 3, 129–134, https://doi.org/10.5194/cp-3-129-2007, 2007.

Drótos, G., Bódai, T., and Tél, T.: Probabilistic concepts in a changing climate: A snapshot attractor picture, J. Climate, 28, 3275–3288, 2015.

Duffing, G.: Erzwungene Schwingungen bei veränderlicher Eigenfrequenz und ihre technische Bedeutung, vol. 41/42 of Sammlung Vieweg, R. Vieweg & Sohn, Braunschweig, https://doi.org/10.1002/zamm.19210010109, 1918.

Einstein, A.: Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen, Annalen der Physik, 322, 549–560; reprinted in Investigations on the Theory of the Brownian Movement, five articles by A. Einstein, edited by: Furth, R., translated by: Cowper, A. D., 1956, Dover Publ., New York, 122 pp., https://doi.org/10.1002/andp.200590005, 1905.

Emiliani, C. and Geiss, J.: On glaciations and their causes, Geol. Rundsch., 46, 576–601, 1959.

Fienga, A., Laskar, J., Exertier, P., Manche, H., and Gastineau, M.: Numerical estimation of the sensitivity of INPOP planetary ephemerides to general relativity parameters, Celestial Mechanics and Dynamical Astronomy, 123, 325–349, 2015.

FitzHugh, R.: Impulses and physiological states in theoretical models of nerve membrane, Biophys. J., 1, 445–466, 1961.

Flint, R. F.: Glacial and Quaternary Geology, Wiley New York, ISBN 978-0471264354, 1971.

Ganopolski, A. and Brovkin, V.: Simulation of climate, ice sheets and CO₂ evolution during the last four glacial cycles with an Earth system model of intermediate complexity, Clim. Past, 13, 1695–1716, https://doi.org/10.5194/cp-13-1695-2017, 2017.

Ganopolski, A. and Calov, R.: The role of orbital forcing, carbon dioxide and regolith in 100 kyr glacial cycles, Clim. Past, 7, 1415–1425, https://doi.org/10.5194/cp-7-1415-2011, 2011.

Ghil, M.: Climate stability for a Sellers-type model, J. Atmos. Sci., 33, 3–20, 1976.

Ghil, M.: Climate sensitivity, energy balance models, and oscillatory climate models, J. Geophys. Res.-Atmos., 89, 1280–1284, 1984.

Ghil, M.: Cryothermodynamics: the chaotic dynamics of paleoclimate, Physica D, 77, 130–159, 1994.
Ghil, M.: Hilbert problems for the geosciences in the 21st century, Nonlin. Processes Geophys., 8, 211–211, https://doi.org/10.5194/npg-8-211-2001, 2001.

Ghil, M.: Climate variability: Nonlinear and random aspects, in: Encyclopedia of Atmospheric Sciences, 2nd edn., edited by G. R. North, J. P. and Zhang, F., vol. 2, 38–46, Elsevier, ISBN 9780123822253, 2014.

Ghil, M.: A century of nonlinearity in the geosciences, Earth Space Sci., 6, 1007–1042, https://doi.org/10.1029/2019EA000599, 2019.

Ghil, M. and Childress, S.: Topics in Geophysical Fluid Dynamics: Atmospheric Dynamics, Dynamo Theory, and Climate Dynamics, Springer Science+Business Media, Berlin/Heidelberg, https://doi.org/10.1007/978-1-4612-1052-8, 1987.

Ghil, M. and Le Treut, H.: A climate model with cryodynamics and geodynamics, J. Geophys. Res.-Oceans, 86, 5262–5270, 1981.

Ghil, M. and Lucarini, V.: The physics of climate variability and climate change, Rev. Mod. Phys., 92, 035002, https://doi.org/10.1103/RevModPhys.92.035002, 2020.

Ghil, M. and Tavantzis, J.: Global Hopf bifurcation in a simple climate model, SIAM J. Appl. Math., 43, 1019–1041, https://doi.org/10.1137/0143067, 1983.

Ghil, M. and Vautard, R.: Interdecadal oscillations and the warming trend in global temperature time series, Nature, 350, 324–327, 1991.

Ghil, M. and Zaliapin, I.: Understanding ENSO variability and its extrema: A delay differential equation approach, in: Extreme Events: Observations, Modeling and Economics, Geophysical Monograph 214, edited by: Chavez, M., Ghil, M., and Urrutia-Fucugauchi, J., Wiley Online Library, 63–78, https://doi.org/10.1002/9781119157052.ch6, 2015.

Ghil, M., Mullhaupt, A., and Pestiaux, P.: Deep water formation and Quaternary glaciations, Clim. Dynam., 2, 1–10, 1987.

Ghil, M., Chekroun, M. D., and Simonnet, E.: Climate dynamics and fluid mechanics: natural variability and related uncertainties, Physica D, 237, 2111–2126, https://doi.org/10.1016/j.physd.2008.03.036, 2008.

Gildor, H. and Tziperman, E.: Sea ice as the glacial cycles’ climate switch: Role of seasonal and orbital forcing, Paleooceanography, 15, 605–615, 2000.

Guckenheimer, J. and Holmes, P. J.: Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields, Applied Mathematical Sciences, Springer Science & Business Media, https://doi.org/10.1007/978-1-4612-1140-2, 1983.

Hasselmann, K.: Stochastic climate models. I: Theory, Tellus, 28, 473–485, 1976.

Hays, J. D., Imbrie, J., and Shackleton, N. J.: Variations in the Earth’s orbit: pacemaker of the ice ages, Science, 194, 1121–1132, 1976.

Heinrich, H.: Origin and consequences of cyclic ice rafting in the Northeast Atlantic Ocean during the past 130 000 years, Quaternary Res., 29, 142–152, https://doi.org/10.1016/0033-5894(88)90057-9, 1988.

Held, I. M.: The gap between simulation and understanding in climate modeling, B. Am. Meteorol. Soc., 86, 1609–1614, https://doi.org/10.1175/bams-86-11-1609, 2005.

Held, I. M. and Suarez, M. J.: Simple albedo feedback models of the ice caps, Tellus, 26, 613–629, 1974.

Henry, L. G., McManus, J. F., Curry, W. B., Roberts, N. L., Piotrowski, A. M., and Keigwin, L. D.: North Atlantic ocean circulation and abrupt climate change during the last glaciation, Science, 353, 470–474, https://doi.org/10.1126/science.aaf5529, 2016.

Hoffman, P. F., Kaufman, A. J., Halverson, G. P., and Schrag, D. P.: A Neoproterozoic snowball earth, Science, 281, 1342–1346, 1998.

Hogg, A. M. C.: Glacial cycles and carbon dioxide: A conceptual model, Geophys. Res. Lett., 35, 1–5, https://doi.org/10.1029/2007GL032071, 2008.

Huybers, P.: Glacial variability over the last two million years: an extended depth-derived agemodel, continuous obliquity pacing, and the Pleistocene progression, Quaternary Sci. Rev., 26, 37–55, https://doi.org/10.1016/j.quascirev.2006.07.013, 2007.

Huybers, P.: Pleistocene glacial variability as a chaotic response to obliquity forcing, Clim. Past, 5, 481–488, https://doi.org/10.5194/cp-5-481-2009, 2009.

Huybers, P.: Combined obliquity and precession pacing of late Pleistocene deglaciations, Nature, 480, 229–232, https://doi.org/10.1038/nature10626, 2011.

Huybers, P. and Langmuir, C. H.: Delayed CO₂ emissions from mid-ocean ridge volcanism as a possible cause of late-Pleistocene glacial cycles, Earth Planet. Sc. Lett., 457, 238–249, 2017.

Imbrie, J. and Imbrie, J. Z.: Modeling the Climatic Response to Orbital Variations, Science, 207, 943–953, 1980.

Imbrie, J. and Imbrie, K. P.: Ice Ages: Solving the Mystery, 2nd edn., Harvard University Press, ISBN 9780674440753, 1986.

Imbrie, J. Z., Imbrie-Moore, A., and Lisiecki, L. E.: A phase-space model for Pleistocene ice volume, Earth Planet. Sci. Lett., 307, 94–102, https://doi.org/10.1016/j.epsl.2011.04.018, 2011.

Isaacson, E. and Keller, H. B.: Analysis of numerical methods, Dover Publications, Inc., New York, NY, ISBN 9780486137988, 2012.

Jackson, E. A.: Perspectives of Nonlinear Dynamics, Cambridge University Press, New York, ISBN 9780198596219, 1991.

Jordan, D. W. and Smith, P.: Nonlinear Ordinary Differential Equations – An Introduction for Scientists and Engineers, 2nd edn., Oxford University Press, Oxford/New York, ISBN 0-19-859657-X, 1987.

Källén, E., Crafoord, C., and Ghil, M.: Free oscillations in a climate model with ice-sheet dynamics, J. Atmos. Sci., 36, 2292–2303, 1979.

Kwasniok, F.: Analysis and modelling of glacial climate transitions using simple dynamical systems, Philos. T. Roy. Soc. A, 371, https://doi.org/10.1098/rsta.2011.0472, 2013.

Landau, L. D. and Lifshitz, E. M.: Mechanics, vol. I of Course on Theoretical Physics, Pergamon Press, Oxford, 1960.

Laskar, J., Robutel, P., Joutel, F., Gastineau, M., Correia, A., and Levrard, B.: A long-term numerical solution for the insolation quantities of the Earth, Astron. Astrophys., 428, 261–285, 2004a.

Laskar, J., Robutel, P., Joutel, F., Gastineau, M., Correia, A. C. M., and Levrard, B.: Computation of various insolation quantities for Earth, IMCCE, Observatoire de Paris/CNRS [data set], http://vo.imcce.fr/insola/earth/online/earth/online/index.php (last access: 13 April 2022), 2004b.

Lenssen, N. J. L., Schmidt, G. A., Hansen, J. E., Menne, M. J., Persin, A., Ruedy, R., and Zyss, D.: Improvements in the GIS-
TEMP uncertainty model, J. Geophys. Res.-Atmos., 124, 6307–6326, 2019.

Le Treut, H. and Ghil, M.: Orbital forcing, climatic interactions, and glaciation cycles, J. Geophys. Res.-Oceans, 88, 5167–5190, 1983.

Le Treut, H., Portes, J., Jouzel, J., and Ghil, M.: Isotopic modeling of climatic oscillations: Implications for a comparative study of marine and ice core records, J. Geophys. Res.-Atmos., 93, 9365–9383, 1988.

Lisiecki, L. E. and Raymo, M. E.: A Pliocene-Pleistocene stack of 57 globally distributed benthic $δ^{18}$O records, Paleooceanography, 20, PA1003, https://doi.org/10.1029/2004PA001071, 2005a.

Lisiecki, L. E. and Raymo, M. E.: Pliocene-Pleistocene stack of globally distributed benthic stable oxygen isotope records, PANGAEA [data set], https://doi.org/10.1594/PANGAEA.704257, 2005b.

Lorenz, E. N.: Deterministic nonperiodic flow, J. Atmos. Sci., 20, 130–141, 1963.

Maasch, K. A. and Saltzman, B.: A Low-Order Dynamical Model of Global Climatic Variability Over the Full Pleistocene, J. Geophys. Res., 95, 1955–1963, 1990.

Marangio, L., Sedro, J., Galatolo, S., Di Garbo, A., and Ghil, M.: Arnold maps with noise: Differentiability and non-monotonicity of the rotation number, J. Stat. Phys., 179, 1–31, https://doi.org/10.1007/s10955-019-02421-1, 2019.

Matteucci, G.: Orbital forcing in a stochastic resonance model of the Late-Pleistocene climatic variations, Clim. Dynam., 3, 179–190, https://doi.org/10.1007/BF01058234, 1989.

Milankovitch, M.: Théorie mathématique des phénomènes thermiques produits par la radiation solaire, Gauthier-Villars, Paris, 1920.

Mitsui, T. and Aihara, K.: Dynamics between order and chaos in conceptual models of glacial cycles, Clim. Dynam., 42, 3087–3099, https://doi.org/10.1007/s00382-013-1793-x, 2014.

Mitsui, T. and Crucifix, M.: Influence of external forcings on abrupt millennial-scale climate changes: a statistical modelling study, Clim. Dynam., 48, 2729–2749, 2017.

Mitsui, T., Crucifix, M., and Aihara, K.: Bifurcations and strange nonchaotic attractors in a phase oscillator model of glacial-interglacial cycles, Physica D, 306, 25–33, 2015.

Nagumo, J., Arimoto, S., and Yoshizawa, S.: An active pulse transmission line simulating nerve axon, Proceedings of the IRE, 50, 2061–2070, 1962.

National Research Council: Understanding Climatic Change, a Program for Action, National Academy of Sciences, Washington, DC, 239 pp., ISBN 978-0309023238, 1975.

Nicolis, C.: Solar variability and stochastic effects on climate, Sol. Phys., 74, 473–478, https://doi.org/10.1007/BF00154530, 1981.

Nicolis, G. and Nicolis, C.: Foundations of Complex Systems, World Scientific, 2nd edn., https://doi.org/10.1142/8260, 2012.

North, G. R.: Analytical solution to a simple climate model with diffusive heat transport, J. Atmos. Sci., 32, 1301–1307, 1975.

North Greenland Ice Core Project members: High-resolution record of the Northern Hemisphere climate extending into the last interglacial period, Nature, 431, 147–151, 2004.

Oerlemans, J.: Model experiments on the 100 000-year glacial cycle, Nature, 287, 430–432, 1980.

Oerlemans, J.: On the origin of the ice ages, in: Milankovitch and Climate: Understanding the Response to Astronomical Force-

ing, vols. I & II, edited by: Berger, A., Imbrie, J., Hays, J., Kukla, G., and Saltzman, B., D. Reidel Publ. Co., 607–611, https://doi.org/10.1007/978-94-017-4841-4, 1984.

Omta, A. W., Kooi, B. W., van Voorn, G. A., Rickaby, R. E., and Follows, M. J.: Inherent characteristics of sawtooth cycles can explain different glacial periodicities, Clim. Dynam., 46, 557–569, 2016.

Paillard, D.: The timing of Pleistocene glaciations from a simple multiple-state climate model, Nature, 391, 378–381, 1998.

Paillard, D.: Glacial cycles: Toward a new paradigm, Rev. Geophys., 39, 325–346, https://doi.org/10.1029/2000RG000091, 2001.

Paillard, D. and Parrenin, F.: The Antarctic ice sheet and the triggering of deglaciations, Earth Planet. Sci. Lett., 227, 263–271, https://doi.org/10.1016/j.epsl.2004.08.023, 2004.

Parrenin, F. and Paillard, D.: Amplitude and phase of glacial cycles from a conceptual model, Earth Planet. Sci. Lett., 214, 243–250, https://doi.org/10.1016/S0012-821X(03)00363-7, 2003.

Parrenin, F. and Paillard, D.: Terminations VI and VIII (~530 and ~720 kyr BP) tell us the importance of obliquity and precession in the triggering of deglaciations, Clim. Past, 8, 2031–2037, https://doi.org/10.5194/cp-8-2031-2012, 2012.

Pelletier, J. D.: Coherence resonance and ice ages, J. Geophys. Res., 108, 1–14, https://doi.org/10.1029/2002jd003120, 2003.

Pierini, S. and Ghil, M.: Climate tipping points induced by parameter drift: an excitable system study, Sci. Rep., in press, 2022.

Pierini, S., Ghil, M., and Chekroun, M. D.: Exploring the pullback attractors of a low-order quasigeostrophic ocean model: The deterministic case, J. Climate, 29, 4185–4202, 2016.

Pierini, S., Chekroun, M. D., and Ghil, M.: The onset of chaos in nonautonomous dissipative dynamical systems: a low-order ocean-model case study, Nonlin. Processes Geophys., 25, 671–692, https://doi.org/10.5194/npg-25-671-2018, 2018.

Pierrehumbert, R. T.: High levels of atmospheric carbon dioxide necessary for the termination of global glaciation, Nature, 429, 646–649, https://doi.org/10.1038/nature02640, 2004.

Pikovsky, A., Rosenblum, M. G., and Kurths, J.: Synchronization, A Universal Concept in Nonlinear Sciences, Cambridge University Press, Cambridge, https://doi.org/10.1017/CBO9780511755743, 2001.

Pikovsky, A. S. and Kurths, J.: Coherence resonance in a noise-driven excitable system, Phys. Rev. Lett., 78, 775–778, https://doi.org/10.1103/physrevlett.78.775, 1997.

Poincaré, H.: Méthodes nouvelles de la Mécanique céleste, vols. I–III, Gauthier-Villars, https://doi.org/10.3931/e-rara-421, 1892–1899.

Pollard, D.: A coupled climate-ice sheet model applied to the Quaternary ice ages, J. Geophys. Res.-Oceans, 88, 7705–7718, 1983.

Quinn, C., Sieber, J., Von Der Heydt, A. S., and Lenton, T. M.: The Mid-Pleistocene Transition induced by delayed feedback and bistability, Dynamics and Statistics of the Climate System, 3, dzy005, https://doi.org/10.1093/climsys/dzy005, 2018.

Rasmussen, S. O., Bigler, M., Blockley, S. P., Blunier, T., Buchardt, S. L., Clausen, H. B., Cvijanovic, I., Dahl-Jensen, D., Johnsen, S. J., Fischer, H., Gkinis, V., Guillevic, M., Hoek, W. Z., Lowe, J. J., Pedro, J. B., Popp, T., Seierstad, I. K., Steffensen, J. P., Svensson, A. M., Vallelonga, P., Vinther, B. M., Walker, M. J., Wheatley, J. J., and Winstrup, M.: A stratigraphic framework for abrupt climatic changes during the Last Glacial period based on
three synchronized Greenland ice-core records: refining and ext-
tending the INTIMATE event stratigraphy, Quaternary Sci. Rev.,
106, 14–28, 2014.

Rial, J. A. and Yang, M.: Is the frequency of abrupt climate change modulated by the orbital insolation?, Geophys. Monog. Ser., 173, 167–174, https://doi.org/10.1029/173GM12, 2007.

Riechers, K.: kriechers/Orbital-Insolation-Variations-Intrinsic-
Climate-Variability-and-Quaternary-Glaciations: Video Su-
plement to: Orbital Insolation Variations, Intrinsic Climate
Variability, and Quaternary Glaciations (1.0), Zenodo [video],
https://doi.org/10.5281/zenodo.6346211, 2022.

Roberts, A. and Saha, R.: Relaxation oscillations in an ideal-
ized ocean circulation model, Clim. Dynam., 48, 2123–2134,
https://doi.org/10.1007/s00382-016-3195-3, 2017.

Rocsoreanu, C., Georgescu, A., and Giurgiteanu, N.: The FitzHugh-
Nagumo model: bifurcation and dynamics, 1st edn., Springer,
Dordrecht, https://doi.org/10.1007/978-94-015-9548-3, 2000.

Rousseau, D.-D., Antoine, P., Boers, N., Lagroix, F., Ghil, M., Lo-
max, J., Fuchs, M., Debret, M., Hatté, C., Moine, O., Gauthier,
C., Jordanova, D., and Jordanova, N.: Dansgaard–Oeschger-
like events of the penultimate climate cycle: the loess point
of view, Clim. Past, 16, 713–727, https://doi.org/10.5194/cp-16-
713-2020, 2020.

Rousseau, D.-D., Bagniewski, W., and Ghil, M.: Abrupt cli-
mate changes and the astronomical theory: are they related?,
Clim. Past, 18, 249–271, https://doi.org/10.5194/cp-18-249-
2022, 2022.

Ruddiman, W. F. and McIntyre, A.: The North Atlantic Ocean dur-
ing the last deglaciation, Palaeogeogr. Palaeocl., 35, 145–214,
1981.

Rulkov, N. F., Sushchik, M. M., Tsimring, L. S., and Abar-
banel, H. D. I.: Generalized synchronization of chaos in direc-
tionally coupled chaotic systems, Phys. Rev. E, 51, 980–994,
https://doi.org/10.1103/physreve.51.980, 1995.

Saltzman, B. and Maasch, K. A.: Carbon cycle instability as a cause of the late Pleistocene ice age oscillations: modeling the asym- metric response, Global Biogeochem. Cy., 2, 177–185, 1988.

Saltzman, B. and Maasch, K. A.: A first-order global model of late Cenozoic climatic change, Transactions of the Royal Society of Edinburgh: Earth Sciences, 81, 315–325, https://doi.org/10.1017/S0263593300020824, 1990.

Saltzman, B. and Maasch, K. A.: A first-order global model of late Cenozoic climatic change II. Further analysis based on a simplification of CO₂ dynamics, Clim. Dynam., 5, 201–210, https://doi.org/10.1007/BF00210005, 1991.

Saltzman, B. and Sutera, A.: The mid-Quaternary climatic transi-
tion as the free response of a three-variable dynamical model, J.
Atmos. Sci., 44, 236–241, 1987.

Saltzman, B., Sutera, A., and Evenson, A.: Structural stochastic sta-
bility of a simple auto-oscillatory climatic feedback system, J.
Atmos. Sci., 38, 494–503, 1981.

Schneider, S. H. and Dickinson, R. E.: Climate modelling, Rev.
Geophys. Space Ge., 25, 447–493, 1974.

Seierstad, I. K., Abbott, P. M., Bigler, M., Blunier, T., Bourne,
A. J., Brook, E., Buchardt, S. L., Buizert, C., Clausen,
H. B., Cook, E., Dahl-Jensen, D., Davies, S. M., Guille-
vic, M., Johnsen, S. J., Pedersen, D. S., Popp, T. J.,
Rasmussen, S. O., Severinghaus, J. P., Svensson, A., and
Vinther, B. M.: Consistently dated records from the Green-

land GRIP, GISP2 and NGRIP ice cores for the past 104 ka
reveal regional millennial-scale δ¹⁸O gradients with possi-
ble Heinrich event imprint, Quaternary Sci. Rev., 106, 29–
46, https://doi.org/10.1016/j.quascirev.2014.10.032, 2014 (data
available at: https://www.iceandclimate.nbi.ku.dk/data/, last ac-
cess: 3 April 2022).

Sellers, W. D.: A global climatic model based on the energy balance of the Earth atmosphere, J. Appl. Meteorol., 8, 392–400, 1969.

Sergin, V. Y.: Numerical modeling of the glaciers-ocean-
atmosphere global system, J. Geophys. Res.-Oceans, 84, 3191–
3204, 1979.

SMIC: Inadvertent Climate Modification: Report of the Study of Man's Impact on Climate, The MIT Press, Cambridge, Mass., 308 pp., ISBN 9780262191012, 1971.

Talento, S. and Ganopolski, A.: Reduced-complexity model for the impact of anthropogenic CO₂ emissions on future glacial cycles, Earth Syst. Dynam., 12, 1275–1293, https://doi.org/10.5194/esd-12-1275-2021, 2021.

Toggweiler, J.: Origin of the 100,000-year timescale in Antarc-
tic temperatures and atmospheric CO₂, Paleooceanography, 23,
PA2211, https://doi.org/10.1029/2006PA001405, 2008.

Tziperman, E. and Gildor, H.: The stabilization of the thermohaline circulation by the temperature-precipitation feedback, J. Phys. Oceanogr., 32, 2707–2714, 2002.

Tziperman, E., Raymo, M. E., Huybers, P., and Wunsch, C.: Consequences of pacing the Pleistocene 100 kyr ice ages by nonlinear phase locking to Milankovitch forcing, Paleooceanography, 21, 1–11, https://doi.org/10.1029/2005PA001241, 2006.

Van der Pol, B.: On relaxation-oscillations, The London, Edinburgh and Dublin Phil. Mag. and J. Sci., 2, 978–992, 1926.

Vannitsem, S., Demaeyer, J., and Ghil, M.: Extratropical low-frequency variability with ENSO forcing: A reduced-order coupled model study, J. Adv. Model. Earth Syst., 13, e2021MS002530, https://doi.org/10.1029/2021MS002530, 2021.

Varadi, F., Runnegar, B., and Ghil, M.: Successive refinements in long-term integrations of planetary orbits, Astrophys. J., 592, 620–630, 2003.

Verbitsky, M. Y. and Crucifix, M.: π-theorem generalization of the ice-age theory, Earth Syst. Dynam., 11, 281–289, https://doi.org/10.5194/esd-11-281-2020, 2020.

Verbitsky, M. Y., Crucifix, M., and Volobuev, D. M.: A theory of Pleistocene glacial rhythmicity, Earth Syst. Dynam., 9, 1025–1043, https://doi.org/10.5194/esd-9-1025-2018, 2018.

Vettoretti, G., Ditlevsen, P., Jochum, M., and Rasmussen, S.
O.: Atmospheric CO₂ control of Spontaneous Millennial-Scale
Ice Age Climate Oscillations, Nat. Geosci., 15, 300–306,
https://doi.org/10.1038/s41561-022-00920-7, 2022.

Vissio, G., Lembo, V., Lucarini, V., and Ghil, M.: Evaluat-
ing the performance of climate models based on Wasser-
stein distance, Geophys. Res. Lett., 47, e2020GL089385,
https://doi.org/10.1029/2020GL089385, 2020.

Wang, B.: Random attractors for the stochastic FitzHugh-Nagumo system on unbounded domains, Nonlinear Anal.-Theor., 71, 2811–2828, https://doi.org/10.1016/j.na.2009.01.131, 2009.

Weertman, J.: Rate of growth or shrinkage of non-equilibrium ice-sheets, J. Glaciol., 6, 145–158, 1964.

Weertman, J.: Milankovitch solar radiation variations and ice-age ice-sheet sizes, Nature, 261, 17–20, 1976.
Westerhold, T., Marwan, N., Drury, A. J., Liebrand, D., Agnini, C., Anagnostou, E., Barnet, J. S. K., Bohaty, S. M., Vleeschouwer, D. D., Florindo, F., Frederichs, T., Hodell, D. A., Holbourn, A. E., Kroon, D., Lauretano, V., Littler, K., Lourens, L. J., Lyle, M., Pálike, H., Röhl, U., Tian, J., Wilkens, R. H., Wilson, P. A., and Zachos, J. C.: An astronomically dated record of Earth’s climate and its predictability over the last 66 million years, Science, 369, 1383–1387, 2020.

Wetherald, R. T. and Manabe, S.: The effects of changing the solar constant on the climate of a general circulation model, J. Atmos. Sci., 32, 2044–2059, 1975.

Wilkinson, L. and Friendly, M.: The History of the Cluster Heat Map, Am. Stat., 63, 179–184, https://doi.org/10.1198/tas.2009.0033, 2009.

Willeit, M., Ganopolski, A., Calov, R., and Brovkin, V.: Mid-Pleistocene transition in glacial cycles explained by declining CO₂ and regolith removal, Sci. Adv., 5, eaav7337, https://doi.org/10.1126/sciadv.aav7337, 2019.

Yamakou, M. E., Tran, T. D., Duc, L. H., and Jost, J.: The stochastic Fitzhugh–Nagumo neuron model in the excitable regime embeds a leaky integrate-and-fire model, J. Math. Biol., 79, 509–532, https://doi.org/10.1007/s00285-019-01366-z, 2019.

Zhang, G., Liu, Z., and Ma, Z.: Generalized synchronization of different dimensional chaotic dynamical systems, Chaos Soliton. Fract., 32, 773–779, https://doi.org/10.1016/j.chaos.2005.11.099, 2007.