Orbital insolation variations, intrinsic climate variability, and Quaternary glaciations

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 . Based on his own calculations and on insightful suggestions from Wladimir Köppen and Alfred Wegener , 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 . 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 .

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 and composition . This state of affairs led, after the initial success of the theory of the ice ages, to severe criticism of the temporal mismatch between insolation minima and glaciation maxima e.g.,.

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

Figure 1

Power spectrum of a composite ${\mathit{\delta }}^{\mathrm{18}}$O record using deep-sea cores RC11-120 and E49-18. This figure is based on the work of , as presented by . 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 , © 1987 by Springer Science+Business Media New York. All rights reserved. 1987.

[Figure omitted. See PDF]

The work of 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 . 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, 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 , Dansgaard–Oeschger (D-O) events e.g.,, and Bond cycles , which were hard to explain by orbital forcing given their shorter timescales.

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 ${}^{\circ }$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. , © 2019. American Geophysical Union. All Rights Reserved. 2019.

[Figure omitted. See PDF]

In fact, interest in past climates was heightened not only by these striking observational discoveries but also by the growing concerns about humanity's impact on the climate . Given the declining temperatures between the 1940s and 1970s, on the one hand, as shown in Fig.  (see also , Fig. 3), and the substantial advances in the description of the Quaternary glaciations, on the other hand, it is clear that interest was mainly on the planet's falling into another ice age e.g.,.

As a result of the twofold stimulation provided by data about past glaciations and concern about future ones, a number of researchers in the early-to-mid 1970s worked on energy balance models (EBMs) of climate with multiple stable steady states . 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 calculated temperatures were much lower than those associated with the Quaternary glaciations and incompatible with paleoclimatic evidence available in the 1970s.

New geochemical evidence, however, led in the early 1990s to the discovery of a “snowball” or, at least, “slushball” Earth prior to the emergence of multicellular life, sometime before 650 Myr b2k . 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 – had actually occurred and is now being modeled in much greater detail .

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.,. Certain theoretical paleoclimatologists therefore turned to coupling a “climate” equation, with temperature as its only dependent variable, with an ice sheet equation or a carbon dioxide equation . 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 captured the phase differences between peak ice sheet extent and minimum temperatures in the North Atlantic suggested by , while the work of Saltzman and associates e.g., captured the asymmetry of the glaciation cycles with their more rapid “terminations” .

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 and to the mean recurrence of D-O events, and they lack a match in the forcing lines of panel (b).

[Figure omitted. See PDF]

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 Part IV – is formulated in Fig.  below; see also : 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 . 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, and vs. . In Sect. , 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.  and conclusions follow in Sect. .

Box 1

Fundamental questions regarding the Quaternary glacial–interglacial cycles.

[Figure omitted. See PDF]

2 Self-sustained climate oscillators

2.1 A simple mechanism for climate oscillations

We follow in sketching the simplest physical mechanism for a self-sustained climate oscillation at fixed insolation forcing. Consider the 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 energy balance model (EBM) with a simplified, scalar version of the ice sheet model (ISM). The model's details and further analyses of its ingredients and variants can be found in several references e.g.,.

The basic workings of this climate oscillator can be represented by two coupled ordinary differential equations (ODEs), written symbolically as follows: $$\begin{array}{}\text{1a}& {\displaystyle }\dot{T}\simeq -V,\text{1b}& {\displaystyle }\dot{V}\simeq T.\end{array}$$ Here $T$ stands for global temperature and $V$ for global ice volume, while Eq. () is an EBM and Eq. () is an ice sheet model (ISM). The “$\simeq$” 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 and , while the ISM relies on the precipitation–temperature feedback postulated by KCG and used also by , who coined the term.

The latter feedback can be better understood by writing the following equations: $$\begin{array}{}\text{2a}& {\displaystyle }\dot{V}\simeq p,\text{2b}& {\displaystyle }p\simeq p_{\text{ac}}-p_{\text{ab}},\text{2c}& {\displaystyle }{p}_{\text{ac}}\simeq T.\end{array}$$ Here $p$ is net precipitation on the single ice sheet of the globally integrated model, given by the difference in Eq. () between the accumulation $p_{\text{ac}}$ and the ablation $p_{\text{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_{\text{ac}}$ with $T$ exceeds the more obvious increase of $p_{\text{ab}}$ with $T$. Hence, $p\simeq T$, and we have derived therewith Eq. (), i.e., $V\simeq p\simeq T$. For more recent studies of the precipitation–temperature feedback, see .

More generally, the presence of feedbacks of opposite sign in a system of two linear coupled ODEs, $$\dot{x}=y,\dot{y}=-x,$$ leads to an oscillation, with the solution given by two trigonometric functions in quadrature with each other, $x\left(t\right)=\sin \left(t\right)$, $y\left(t\right)=\cos \left(t\right)=\sin (t+\mathit{\pi }/\mathrm{2})$, and the trajectory describing a circle in the $(x,y)$ phase plane, $x^{\mathrm{2}}\left(t\right)+y^{\mathrm{2}}\left(t\right)=\mathrm{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 (), 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 and in Sect. 12.2, namely,

3$$\dot{z}=(\mathit{\mu }+\mathit{ı}\mathit{\omega })z+c\left(z\stackrel{\mathrm{‾}}{z}\right)z.$$ Here $z=x+\mathit{ı}y$ is complex, where $\mathit{ı}=\sqrt{-\mathrm{1}}$ is the imaginary unit, while $\mathit{\mu }$ is a real bifurcation parameter, and $c,\mathit{\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 ${\mathit{\mu }}_{\ast }$ in the next subsection.

Figure 4

Supercritical Hopf bifurcation. (a) Vector field $\dot{z}\left(z\right)$ of Eq. () for the parameter values $\mathit{\mu }=-\mathrm{1}$ and $c=-\mathrm{1}$; $z=x+\mathit{ı}y$. 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}\left(z\right)$ and three trajectories for $\mathit{\mu }=\mathrm{1}$. In this case, the origin is an unstable fixed point, while the limit cycle with radius ${\mathit{\rho }}_{*}^{\mathrm{1}/\mathrm{2}}=\mathrm{1}$ constitutes the only stable solution. Trajectories that start inside this limit cycle, with $\mathit{\rho }\left(t_{\mathrm{0}}\right)<{\mathit{\rho }}_{*}$, tend to it by spiraling out – as illustrated by the magenta trajectory – while trajectories that start outside the limit cycle, with $\mathit{\rho }\left(t_{\mathrm{0}}\right)>{\mathit{\rho }}_{*}$, approach it by spiraling inward, as illustrated by the gray trajectory. (c) Dependence of the stable solution of Eq. () on the parameter $\mathit{\mu }$, for $c=-\mathrm{1}$. For $\mathit{\mu }\le \mathrm{0}$, the single stable solution is the fixed point located at the origin, $\mathit{\rho }\equiv \mathrm{0}$. For $\mathit{\mu }>\mathrm{0}$, the stable solution is the limit cycle given by $\mathit{\rho }=\mathit{\mu }$.

[Figure omitted. See PDF]

A very natural transformation of variables, 4$$\mathit{\rho }=z\stackrel{\mathrm{‾}}{z}>\mathrm{0},z={\mathit{\rho }}^{\mathrm{1}/\mathrm{2}}\exp \left(\mathit{ı}\mathit{\theta }\right),$$ leads from the complex ODE () to the system of two real and decoupled ODEs, $$\begin{array}{}\text{5a}& {\displaystyle }\dot{\mathit{\rho }}=\mathrm{2}\mathit{\rho }(\mathit{\mu }+c\mathit{\rho }),\text{5b}& {\displaystyle }\dot{\mathit{\theta }}=\mathit{\omega }.\end{array}$$ Equation () simply provides an angular rotation around the origin $\mathit{\rho }=\mathrm{0}=x=y$, since the complex exponential in Eq. () is periodic with period $\mathrm{2}\mathit{\pi }$. Equation () is quadratic in $\mathit{\rho }$, and thus it can have two real roots, $\mathit{\rho }={\mathit{\rho }}_{\mathrm{0}}=\mathrm{0}$ and $\mathit{\rho }={\mathit{\rho }}_{\ast }=-\mathit{\mu }/c$. But $\mathit{\rho }$ has to be positive, and thus in the case in which $c<\mathrm{0}$, the only possible solution for $\mathit{\mu }<\mathrm{0}$ is the origin, and it is stable since $\mathit{\rho }(\mathit{\mu }+c\mathit{\rho })$ is negative for $\mathit{\rho }>\mathrm{0}$ in this case; hence, $\mathit{\rho }$ has to be monotonically decreasing, i.e., all the solutions of Eq. () spiral into the origin. The Hopf bifurcation from this stable steady state to a periodic solution, i.e., a limit cycle with radius ${\mathit{\rho }}^{\mathrm{1}/\mathrm{2}}={\mathit{\rho }}_{\ast }^{\mathrm{1}/\mathrm{2}}$, occurs as $\mathit{\mu }$ crosses $\mathrm{0}$. Since now $\mathit{\rho }(\mathit{\mu }+c\mathit{\rho })>\mathrm{0}$ for $\mathrm{0}<\mathit{\rho }<{\mathit{\rho }}_{\ast }$ and $\mathit{\rho }(\mathit{\mu }+c\mathit{\rho })<\mathrm{0}$ for $\mathit{\rho }>{\mathit{\rho }}_{\ast }$, the limit cycle is stable and trajectories spiral out from inside this cycle and into it from outside; see Fig.  for the so-called supercritical (or soft) Hopf bifurcation case with $c<\mathrm{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 . 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 .

Physically speaking, the presence or absence of the regular, purely periodic oscillations obtained by KCG and illustrated in Fig. 12.6 depends on whether $c\mathit{\gtrless }\mathrm{0}$ in Eq. (). The KCG model's bifurcation parameter is ${\mathit{\mu }}_{\ast }=c_T/c_L$, where $c_T$ is the heat capacity in its EBM, while $c_L$ is the “mass capacity” in its ISM . Large ${\mathit{\mu }}_{\ast }$ corresponds physically to a very small, possibly pre-Pleistocene ice cap . At these values of ${\mathit{\mu }}_{\ast }$, 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 Figs. 3–5. As ${\mathit{\mu }}_{\ast }$ decreases to $O\left(\mathrm{1}\right)$, 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 ; see also Figs. 12.8 and 12.9.

To clarify the simple physical concepts that underlie subcritical and supercritical Hop bifurcations, let us consider a purely mechanical oscillator with mass $m$, a spring $kx$, and a dashpot $\mathit{\alpha }\dot{x}$ ,

6$$m\ddot{x}=-\mathit{\alpha }\left(x\right)\dot{x}-k\left(x\right)x.$$ If $k$ $=$ const. and $\mathit{\alpha }$ $=$ 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: 7$$\ddot{x}+\mathit{\alpha }\left(x\right)\dot{x}+k\left(x\right)x=F\cos \left(\mathit{\omega }t\right).$$ Two classical nonlinear cases are those of the equation, in which $k\left(x\right)=x^{\mathrm{2}}$ and $\mathit{\alpha }$ = const., and of the equation, in which $k$ $=$ const. and $\mathit{\alpha }\left(x\right)=\mathit{\nu }(x^{\mathrm{2}}-\mathrm{1})$. The fully nonlinear case in which both the spring and the damping are nonlinear, with $k\left(x\right)=x^{\mathrm{2}}$ and $\mathit{\alpha }\left(x\right)=\mathit{\nu }(x^{\mathrm{2}}-\mathrm{1})$, is known as the Van der Pol–Duffing oscillator e.g.,. Note that all three types of nonlinear oscillators can exhibit chaotic behavior even in the presence of simple periodic forcing e.g.,and references therein. The idea of using such simple, classical oscillators in modeling Quaternary glaciation cycles goes back to .

Sect. 5.6 specifically discuss the case of a soft and a hard spring for a generalized Duffing equation, with $k\left(x\right)=k_{\mathrm{0}}+\mathit{ϵ}h\left(x\right)$, where $\mathrm{0}<\left|\mathit{ϵ}\right|\ll \mathrm{1}$, $h(-x)=h\left(x\right)$, and $h^{\prime \prime }>\mathrm{0}$. A spring is soft if it is sublinear, $\mathit{ϵ}<\mathrm{0}$, and hard if it is superlinear, $\mathit{ϵ}>\mathrm{0}$; see their Eq. (5.37) and Fig. 5.4a, b, with $h\left(x\right)=x^{\mathrm{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_{\mathrm{0}}$ increases past a critical value $k_{\ast }$, 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_{\mathrm{0}}$ crosses the value $k_{\ast }$. Please compare the behavior of supercritical and subcritical Hopf bifurcations in Sect. 12.2 and Figs. 12.7–12.9 and see Fig. 5.7 for the change in the nonlinear response of a Duffing oscillator as its spring changes from soft, with $\mathit{ϵ}<\mathrm{0}$, to hard, with $\mathit{ϵ}>\mathrm{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.,. A fair number of distinct dynamical theories for this transition have been formulated e.g.,. 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 , 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 oscillators would have to lead to a more abrupt transition, as suggested by . Later, showed that the models by exhibit MPT-like behavior via supercritical or subcritical Hopf bifurcations, depending on the parameter values. The existing ${\mathit{\delta }}^{\mathrm{18}}$O and ${\mathit{\delta }}^{\mathrm{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 for the high-resolution NGRIP record – 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. a and b in Sect. . To explore this possibility, used somewhat simplified insolation forcing based on the calculations of and applied it to a slightly modified version of the oscillator. These authors found that, as expected for a nonlinear oscillator, its internal frequency $f_{\mathrm{0}}$ is affected strongly by the forcing ones, $\mathit{\{}{f}_{\mathrm{1}},\mathrm{\dots },f_{\mathrm{5}}\mathit{\}}$, resulting in both nonlinear resonance and combination tones .

Depending on parameter values, the periodicity $P_{\mathrm{0}}=\mathrm{1}/f_{\mathrm{0}}$ of the oscillator is $P_{\mathrm{0}}\simeq \mathrm{6}$–7 kyr. The lines in the simplified insolation spectrum used by had the following periodicities: $P_{\mathrm{1}}=\mathrm{19}\text{kyr}$, $P_{\mathrm{2}}=\mathrm{23}\text{kyr}$, $P_{\mathrm{3}}=\mathrm{41}\text{kyr}$, $P_{\mathrm{4}}=\mathrm{100}\text{kyr}$, and $P_{\mathrm{5}}=\mathrm{400}\text{kyr}$. These periodicities correspond to the two precessional ones, the obliquity one, and the two eccentricity ones. The actual celestial-mechanics calculations that based his insolation calculations on have been substantially updated since e.g.,. 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.

Figure 5

Power spectra of simulated (a) deep-sea and (b) ice core records for the forced climate oscillator of 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 $-\mathrm{2}$ power law slope. The orbital peaks at $P_{\mathrm{1}}=\mathrm{19}$ kyr, $P_{\mathrm{2}}=\mathrm{23}$ kyr and $P_{\mathrm{3}}=\mathrm{41}$ kyr are also present, along with peaks at 14.7 and 10.4 kyr, which correspond to the sum tones $f_{\mathrm{2}}+f_{\mathrm{3}}$ and $f_{\mathrm{1}}+f_{\mathrm{2}}$, respectively. Panel (b) is in log–linear coordinates and shows additional harmonics and sum tones, as well as the difference tone $f_{\mathrm{1}}-f_{\mathrm{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 , © 1988 by the American Geophysical Union. 1988.

[Figure omitted. See PDF]

The results of the model on the evolution of the primary climate variables $T$ and $V$ were converted to ${\mathit{\delta }}^{\mathrm{18}}$O values in simulated isotopic records of marine sediment and ice cores by ; the spectra of the latter are plotted in Fig. . The values on the abscissa of Fig. a are values of the logarithm of frequency, while those in Fig. b 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 , it is important to realize (i) that the large continuous background in Fig. a is purely of deterministically chaotic origin since there is no stochastic element whatsoever in the model or in its forcing and (ii) that the dominant peak at 109 kyr is not directly forced by the $f_{\mathrm{4}}=\mathrm{1}/\mathrm{100}$ kyr${}^{-\mathrm{1}}$ eccentricity line but rather it is due to the difference tone between the two precessional frequencies, $f_{\mathrm{1}}$ and $f_{\mathrm{2}}$. Finally, it is the nonlinear, broad resonance of the model's $f_{\mathrm{0}}$ frequency with the quasi-periodic forcing that produces the bump in the spectrum of Fig. a 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 , as obtained by and and shown here in Fig. a – stochastic contributions à la 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 . In the paleoclimatic context, 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 . In this context, also provide a simple theory that addresses scaling properties in the glacial cycles and their spectra, based on the so-called Buckingham $\mathit{\pi }$-theorem e.g.,.

3.2 Basic facts of NDS and RDS life

The highly preliminary results summarized in Sect.  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.,. 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 and by , from deterministically autonomous to non-autonomous and random dynamical systems NDSs and RDSs, e.g.,.

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 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 . Slightly more in-depth but still fairly expository presentations can be found in and in . 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. .

3.2.1 Autonomous and non-autonomous systems

Succinctly, one can write an autonomous system as

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

For instance, two distinct trajectories, ${\mathit{X}}_{\mathrm{1}}\left(t\right)$ and ${\mathit{X}}_{\mathrm{2}}\left(t\right)$, 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 $\mathit{X}\left(t\right)$ as time $t$ tends to $+\mathrm{\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: 9$$\dot{\mathit{X}}=\mathit{F}(\mathit{X},t;\mathit{\mu }).$$ In Eq. (), the dependence of $\mathit{F}$ on $t$ may be periodic, $\mathit{F}(\mathit{X},t+P)=\mathit{F}(\mathit{X},t)$, as in various El Niño–Southern Oscillation (ENSO) models, where the period $P=\mathrm{12}$ months, or monotonic, $\left|\mathit{F}\right|(\mathit{X},t+\mathit{\tau })\ge \left|\mathit{F}\right|(\mathit{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 $\mathit{\{}{\mathit{\mu }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\mu }}_p\mathit{\}}$.

To illustrate the fundamental distinction between an autonomous system like Eq. () and a non-autonomous one like Eq. (), consider the simple scalar version of these two equations: 10a$$\dot{X}=-\mathit{\beta }X,$$ and 10b$$\dot{X}=-\mathit{\beta }X+\mathit{\gamma }t,$$ respectively. We assume that both systems are dissipative, i.e., $\mathit{\beta }>\mathrm{0}$, and that the forcing is monotonically increasing, $\mathit{\gamma }\ge \mathrm{0}$, as would be the case for anthropogenic forcing in the industrial era. pointed out the key role of dissipativity in giving rise to strange but attracting solution behavior, while emphasized its importance and pervasive character in climate dynamics. Clearly the only attractor for the solutions of Eq. (), given any initial point $X\left(\mathrm{0}\right)=X_{\mathrm{0}}$, is the stable fixed point $X_{\ast }=\mathrm{0}$, attained as $t\to +\mathrm{\infty }$.

In the case of Eq. (), however, this forward-in-time approach yields blow-up from $t\to +\mathrm{\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\left(s\right)=X_{\mathrm{0}}$ and $s\ll t$. With this little change of approach, one can easily verify that 11$$\left|X\right(s,t;X_{\mathrm{0}})-{\mathcal{A}}_t|\to \mathrm{0}\text{as}s\to -\mathrm{\infty },$$ for all $t$ and $X_{\mathrm{0}}$, where ${\mathcal{A}}_t=(\mathit{\gamma }/\mathit{\beta })(t-\mathrm{1}/\mathit{\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 ${\mathcal{A}}_t$; this family follows the forcing $\mathit{\gamma }t$, and it thus has a linear growth in time $t$. Hence, the fixed point $X_{\ast }$ of Eq. () is in fact a moving target and it is given by $X_{\ast }=\mathit{\gamma }t/\mathit{\beta }$. Due to the system's inertia, the set ${\mathcal{A}}_t$ that is approached by the trajectories lags this time-dependent fixed point by a constant offset of $\mathit{\gamma }/{\mathit{\beta }}^{\mathrm{2}}$.

3.2.2 Pullback attractor (PBA)

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

12$$\colorbox[rgb]{}{$$}=\mathit{\{}{\mathcal{A}}_t{\mathit{\}}}_{t\in \mathbb{R}},$$ is termed the pullback attractor (PBA) of an NDS if the following two conditions are fulfilled:

• i.

  each snapshot ${\mathcal{A}}_t$ is compact and the family of snapshots $\mathit{\left\{}\mathcal{A}\right(t){\mathit{\}}}_{t\in \mathbb{R}}$ is invariant with respect to the dynamics,13$$\mathit{X}(t,s;{\mathit{X}}_{\mathrm{0}})\in {\mathcal{A}}_t\forall s\le t\mathrm{and}{\mathit{X}}_{\mathrm{0}}\in A_s,$$and

• ii.

  the pullback attraction occurs for all times,14$$\underset{s\to -\mathrm{\infty }}{lim}\left|\mathit{X}\right(t,s;{\mathit{X}}_{\mathrm{0}})-{\mathcal{A}}_t|=\mathrm{0}\forall t.$$

To further improve the reader's intuition for PBAs, we provide a second illustrative example here. A system defined in polar coordinates by 15$$\dot{\mathit{\rho }}=\mathit{\alpha }(\mathit{\mu }-\mathit{\rho }),\dot{\mathit{\varphi }}=\mathit{\omega },$$ with $\mathit{\rho },\mathit{\mu }\in {\mathbb{R}}^{+}$, and $\mathit{\varphi }\in \mathbb{R}/\mathrm{2}\mathit{\pi }$ can easily be seen to exhibit a limit cycle in the $(x,y)$ plane with $(x=\mathit{\rho }\cos \mathit{\varphi },y=\mathit{\rho }\sin \mathit{\varphi })$. An initial deviation of $\mathit{\rho }$ from $\mathit{\mu }$ will decay exponentially, and the system converges to an oscillation of radius $\mathit{\mu }$ with the angular velocity $\mathit{\omega }$. Here, we transform this autonomous dynamical system into a non-autonomous one by modulating the target radius $\mathit{\mu }$ with a sinusoidal forcing, 16$$\mathit{\mu }\to \mathit{\mu }\left(t\right)={\mathit{\mu }}_{\mathrm{0}}+\mathit{\beta }\sin \left(\mathit{\nu }t\right),$$ where the modulation is moderate to guarantee that ${\mathit{\mu }}_{\mathrm{0}}+\mathit{\beta }\sin \left(\mathit{\nu }t\right)>\mathrm{0}$ for all $t$.

Since the dynamics of the phase $\mathit{\varphi }$ and of the radius $\mathit{\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 $\mathit{\rho }\left(s\right)={\mathit{\rho }}_{\mathrm{0}}$ is given by 17a$${\mathcal{A}}^{\left(\mathit{\rho }\right)}(t;{\mathit{\rho }}_{\mathrm{0}})=\underset{s\to -\mathrm{\infty }}{lim}\mathit{\rho }(t,s;{\mathit{\rho }}_{\mathrm{0}})=\mathit{\alpha }\mathit{\beta }\sin (\mathit{\nu }t+\mathit{\vartheta })+{\mathit{\mu }}_{\mathrm{0}},$$ with 17b$$\mathit{\vartheta }=\arctan (-\mathit{\nu }/\mathit{\alpha }),$$ as shown in Appendix . Note that in the limit $s\to -\mathrm{\infty }$, the dependence on the initial value ${\mathit{\rho }}_{\mathrm{0}}$ vanishes and the attracting set ${\mathcal{A}}_t^{\left(\mathit{\rho }\right)}$ performs an oscillation of the same frequency as the forcing. It lags the phase of the time-dependent fixed point $\mathit{\mu }\left(t\right)$ by the constant $\mathit{\vartheta }$, while its amplitude is amplified by the factor $\mathit{\alpha }$. Since $\mathit{\rho }$ is restricted to positive values, this solution requires $\mathit{\alpha }\mathit{\beta }<{\mathit{\mu }}_{\mathrm{0}}$.

The PBA with respect to the coordinate $\mathit{\rho }$ is comprised of the family of all the sets ${\mathcal{A}}_t^{\left(\mathit{\rho }\right)}$ as defined in Eq. (18) and thus reads 18$$\colorbox[rgb]{}{$$}=\mathit{\left\{}\mathit{\alpha }\mathit{\beta }\sin \right(\mathit{\nu }t+\mathit{\vartheta })+{\mathit{\mu }}_{\mathrm{0}}{\mathit{\}}}_{t\in \mathbb{R}}.$$ Since the pullback limit for the phase $\mathit{\varphi }$ does not exist, no constraints on it other than $\mathit{\varphi }\in [\mathrm{0},\mathrm{2}\mathit{\pi })$ are imposed by the dynamics. Hence, for the system () comprised of radius and phase, we find that the distance of any trajectory at time $t$ – i.e. $\left(\mathit{\rho }(t;t_{\mathrm{0}},{\mathit{\rho }}_{\mathrm{0}}),\mathit{\varphi }(t;t_{\mathrm{0}},{\mathit{\varphi }}_{\mathrm{0}})\right)$ – to the set ${\mathcal{A}}_t=\mathit{\{}\left(\mathit{\alpha }\mathit{\beta }\sin (\mathit{\nu }t+\mathit{\vartheta })+{\mathit{\mu }}_{\mathrm{0}},\mathit{\phi }\right):\mathit{\phi }\in [\mathrm{0},\mathrm{2}\mathit{\pi }\left)\mathit{\right\}}$ tends to zero as we pullback the initial time $t_{\mathrm{0}}$ to $-\mathrm{\infty }$. The pullback attracting sets ${\mathcal{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 19$$\colorbox[rgb]{}{$$}=\mathit{\{}\left(\mathit{\alpha }\mathit{\beta }\sin (\mathit{\nu }t+\mathit{\vartheta })+{\mathit{\mu }}_{\mathrm{0}},\mathit{\phi }\right):\mathit{\phi }\in [\mathrm{0},\mathrm{2}\mathit{\pi }){\mathit{\}}}_{t\in \mathbb{R}}.$$

Figure 6

Trajectories and PBA of the system defined by Eqs. () and (); in all the panels the parameter values are ${\mathit{\mu }}_{\mathrm{0}}=\mathrm{15}$, $\mathit{\alpha }=\mathrm{10}$, $\mathit{\nu }=\mathrm{0.2}$, and $\mathit{\omega }=\mathrm{0.3}$, unless stated otherwise. (a) Trajectories $\left(\mathit{\rho }\right(t),\mathit{\varphi }(t\left)\right)$ 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_{\mathrm{0}}=-\mathrm{100}$ to $t_{\text{f}}=\mathrm{200}$ for $\mathit{\nu }=\mathrm{0.3},\mathrm{0.2},{\mathrm{0.5}}^{\mathrm{1}/\mathrm{2}}$ in panels (c), (d), and (e), respectively. The values are chosen such that the ratio between $\mathit{\nu }$ and $\mathit{\omega }$ is $\mathit{\omega }/\mathit{\nu }=\mathrm{1}$, $\mathit{\omega }/\mathit{\nu }\in \mathbb{Z}\backslash \mathit{\left\{}\mathrm{1}\mathit{\right\}}$ and $\mathit{\omega }/\mathit{\nu }\in \mathbb{R}\backslash \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 $\mathit{\omega }/\mathit{\nu }=\mathrm{3}/\mathrm{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_{\mathrm{i}}=-\mathrm{200}$ to $t_{\mathrm{f}}=\mathrm{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=\mathrm{0.20}$; see details in the video supplement.

[Figure omitted. See PDF]

Figure  shows trajectories of the system starting from different points in the past. In Fig. a, 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 $\mathit{\rho }\left(t\right)$ and their convergence to the PBA ${\mathcal{A}}_t^{\left(\mathit{\rho }\right)}$ are shown Fig. b. Figure f shows a heat map that approximates a portion of the PBA's invariant measure projected onto the $(x,y)$ plane. For a clean definition of such a measure in NDSs and RDSs, please see , , , or . Essentially, the heat map here counts the number of times that 100 trajectories integrated from $t=-\mathrm{200}$ to $t=\mathrm{200}$, with randomized initial conditions like the ones shown in Fig. a, cross small pixels in the $(x,y)$ plane.

Figure c–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 $\mathit{\omega }/\mathit{\nu }$, and three different cases must be distinguished. If the radius is modulated with the same frequency as the oscillation itself, i.e. $\mathit{\omega }=\mathit{\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 $\mathit{\rho }\left(t\right)$ approaches the PBA ${\mathcal{A}}_t^{\left(\mathit{\rho }\right)}$.

If $\mathit{\omega }$ and $\mathit{\nu }$ are rationally related, i.e., $m\mathit{\omega }=n\mathit{\nu }$ with $n,m\in \mathbb{N}$, as in Fig. d, 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\mathrm{2}\mathit{\pi }/\mathit{\nu }=m\mathrm{2}\mathit{\pi }/\mathit{\omega }$. That is, such a trajectory will appear as an $n$-fold quasi-closed loop.

Finally, if $\mathit{\omega }/\mathit{\nu }\notin \mathbb{Z}$, as in Fig. e, 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 $\mathcal{D}=\mathit{\left\{}\right(\mathit{\rho },\mathit{\varphi }):\mathit{\rho }\in [{\mathit{\mu }}_{\mathrm{0}}-\mathit{\alpha }\mathit{\beta },{\mathit{\mu }}_{\mathrm{0}}+\mathit{\alpha }\mathit{\beta }]\text{and}\mathit{\varphi }\in [\mathrm{0},\mathrm{2}\mathit{\pi }\left)\mathit{\right\}}$.

3.2.3 Random attractor

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

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

The noise processes may include “weather” and volcanic eruptions when $\mathit{X}\left(t\right)$ is “climate”, thus generalizing the linear model of , 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. (), the concepts introduced by the simple deterministic examples of Eq. () and Eqs. () and () above can be illustrated by the random attractor $\left(\mathit{\omega }\right)$ in Fig. .

Figure 7

Schematic diagram of a random attractor $\left(\mathit{\omega }\right)$ and of the pullback attraction to it; here $\mathit{\omega }$ labels the particular realization of the random process $\mathit{\theta }\left(t\right)\mathit{\omega }$ that drives the system. We illustrate the evolution in time $t$ of the random process $\mathit{\theta }\left(t\right)\mathit{\omega }$ (light solid black line at the bottom), the random attractor $\left(\mathit{\omega }\right)$ itself (yellow band in the middle), with the snapshots ${\mathcal{A}}_{\mathrm{0}}\left(\mathit{\omega }\right)=\mathcal{A}(\mathit{\omega };t=\mathrm{0}))$ and $\mathcal{A}(\mathit{\omega };t)$ (the two vertical sections, heavy solid lines), and the flow of an arbitrary compact set $\mathcal{B}$ from “pullback times” $t=-{\mathit{\tau }}_{\mathrm{2}}$ and $t=-{\mathit{\tau }}_{\mathrm{1}}$ onto the attractor (heavy blue arrows). See Appendix A for the requisite properties of the random process $\mathit{\theta }\left(t\right)\mathit{\omega }$ that drives the RDS formulated by Eq. (). Adapted from Physica D: Nonlinear Phenomena, 237, , Climate dynamics and fluid mechanics: natural variability and related uncertainties, 2111–2126, © 2008 Elsevier B.V. All rights reserved. 2008, with permission from Elsevier.

[Figure omitted. See PDF]

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 .

studied a specific case of such a random attractor for the paradigmatic, climate-related convection model. The authors introduced multiplicative noise into each of the ODEs of the original, deterministically chaotic system, as shown below: $$\begin{array}{}\text{21a}& {\displaystyle }\text{d}X=P_r(Y-X)\text{d}t+\mathit{\sigma }X\text{d}\mathit{\eta },\text{21b}& {\displaystyle }\text{d}Y=(rX-Y-XZ)\text{d}t+\mathit{\sigma }Y\text{d}\mathit{\eta },\text{21c}& {\displaystyle }\text{d}Z=(-bZ+XY)\text{d}t+\mathit{\sigma }Z\text{d}\mathit{\eta },\end{array}$$ where $r=\mathrm{28}$, $P_r=\mathrm{10}$, and $b=\mathrm{8}/\mathrm{3}$ are the standard parameter values for chaotic behavior in the absence of noise and $\mathit{\sigma }$ is a constant variance of the Wiener process that is not necessarily small; both the noise realization $\mathit{\eta }\left(t\right)$ and $\mathit{\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 ${\mathcal{A}}_t\left(\mathit{\omega }\right)$ of LORA are plotted in Fig.  and a video of its evolution in time $\left(\mathit{\omega }\right)=\mathit{\left\{}{\mathcal{A}}_t\right(\mathit{\omega }){\mathit{\}}}_{t\in \mathbb{R}}$ is available as supplementary material to at 10.1016/j.physd.2011.06.005.

Figure 8

Heat maps of the time-dependent invariant measure ${\mathit{\nu }}_t\left(\mathit{\omega }\right)$ supported on four snapshots $\left(\mathit{\omega }\right)$ of LORA. The values of the parameters $r$, $s$, and $b$ are the classical ones, where the variance of the noise $\mathit{\sigma }=\mathrm{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, , Stochastic climate dynamics: Random attractors and time-dependent invariant measures, 1685–1700, © 2011 Elsevier B.V. All rights reserved. 2011, with permission from Elsevier.

[Figure omitted. See PDF]

have further analyzed the striking effects of the noise on the nonlinear dynamics that are visible in Fig.  and in the video of 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 . These so-called D-O events were followed by intervals of steady moderate cooling, before a short phase of enhanced cooling brought the temperatures back to their pre-event levels e.g.,. This pattern is very clearly apparent in NGRIP ${\mathit{\delta }}^{\mathrm{18}}$O records , and it can be qualitatively simulated by the fast component of a FitzHugh–Nagumo (FHN) model , as pointed out by , among others; see also , , , and .

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${}_{\mathrm{2}}$ concentrations retrieved from Antarctic ice cores and show that this setup brings the model's trajectories into good qualitative agreement with the D-O patterns observed in ${\mathit{\delta }}^{\mathrm{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 ${\mathit{\tau }}_x$ and ${\mathit{\tau }}_y$, with $x\left(t\right)$ being the slow component and $y\left(t\right)$ the fast one: $$\begin{array}{}\text{22a}& {\displaystyle }\dot{x}={\displaystyle \frac{\mathrm{1}}{{\mathit{\tau }}_x}}(y-\mathit{\gamma }),\text{22b}& {\displaystyle }\dot{y}={\displaystyle \frac{\mathrm{1}}{{\mathit{\tau }}_y}}\left[\mathit{\alpha }\right(y-y^{\mathrm{3}})-x].\end{array}$$ 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 $\mathit{\gamma }$.

First, consider the case of large timescale separation

23$${\mathit{\tau }}_x\gg {\mathit{\tau }}_y.$$ This choice guarantees that the fast $y$ component adjusts adiabatically to quasi-static changes of the slow $x$ component. The time derivative of $y\left(t\right)$, as shown in Fig. a, exhibits either three real roots or a single one, depending on the value of $x\left(t\right)$, which shifts the graph of the cubic polynomial $P_{\mathrm{3}}(x,y)=\mathit{\alpha }(y-y^{\mathrm{3}})-x$ globally upwards or downwards. Of the three potential roots, the outer two are stable fixed points for $y\left(t\right)$ 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_{\mathrm{3}}(x,y)$, respectively. Accordingly, we label them $y_{\mathrm{l}}$ and $y_{\mathrm{r}}$. The positions of the local extrema, namely $y_{\text{min}}=-\sqrt{\mathrm{1}/\mathrm{3}}$ and $y_{\text{max}}=\sqrt{\mathrm{1}/\mathrm{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\left(t\right)$ and $y\left(t\right)$. Assume we are in a state where $x>\mathrm{2}\mathit{\alpha }/\sqrt{\mathrm{27}}$ such that $y_{\mathrm{l}}$ is the only root of $P_{\mathrm{3}}(x,y)$. Provided that $\left|\mathit{\gamma }\right|<\sqrt{\mathrm{1}/\mathrm{3}}$, the time derivative $\dot{x}$ and $x$ itself will have opposite signs, with the consequence that a slow adjustment process of $x\left(t\right)$ sets in. This will shift the polynomial $P_{\mathrm{3}}(x,y)$ upwards in Fig. a, which in turn revives the other two roots of $P_{\mathrm{3}}(x,y)$. The fast $y\left(t\right)$ closely follows $y_{\mathrm{l}}$, which is shifted accordingly. Once $x\left(t\right)=-\mathrm{2}\mathit{\alpha }/\sqrt{\mathrm{27}}$, the root $y_{\mathrm{l}}$ stops existing and a fast transition to the neighborhood of the opposite root $y_{\mathrm{r}}$ ensues. This transition switches the sign of $\dot{x}\left(t\right)$ and the slow adjustment of $x\left(t\right)$ now happens in the opposite direction.

In the $(x,y)$ plane, this behavior manifests itself as a stable limit cycle. However, any value $\left|\mathit{\gamma }\right|>\sqrt{\mathrm{1}/\mathrm{3}}$ prevents $\dot{x}\left(t\right)$ from switching sign and therefore interrupts the cyclic destruction and revival of the respective opposite root of $P_{\mathrm{3}}(x,y)$. Instead, $\left|\mathit{\gamma }\right|>\sqrt{\mathrm{1}/\mathrm{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 $\mathit{\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_{\mathrm{l}}$ and $y_{\mathrm{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 () and (). Both $\mathit{\gamma }=\pm \sqrt{\mathrm{1}/\mathrm{3}}$ are critical values of $\mathit{\gamma }$ that give rise to supercritical Hopf bifurcations of the coupled system's fixed points; recall Fig.  of Sect. .

Figure 9

Nullclines of the autonomous FHN model governed by Eq. () with large timescale separation, for $\mathit{\alpha }=\mathrm{2}$, ${\mathit{\tau }}_x=\mathrm{100}$, and ${\mathit{\tau }}_y=\mathrm{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 $\mathit{\gamma }=\mathrm{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_{\mathrm{r}}$ and $y_{\mathrm{l}}$, respectively, as defined in the text. For fixed values of $x$, they constitute stable fixed points for Eq. (). The middle branch (dashed blue) corresponds to the unstable fixed point for Eq. () 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 $\left|\mathit{\gamma }\right|<\sqrt{\mathrm{1}/\mathrm{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_{\mathrm{l}}$ branch of the $y$ nullcline.

[Figure omitted. See PDF]

Figure 10

FitzHugh–Nagumo (FHN) model with parameters ${\mathit{\tau }}_x=\mathrm{100}$, ${\mathit{\tau }}_y=\mathrm{60}$, and $\mathit{\alpha }=\mathrm{2}$. (a) The cubic polynomial $P_{\mathrm{3}}(x,y)$ of the fast derivative as a function of $y$ for $x=\mathrm{0}$ (solid blue line); the red lines point to the local maximal and minimal values of $P_{\mathrm{3}}(x,y)$, namely $\pm \mathrm{2}a/\sqrt{\mathrm{27}}$, respectively – these are the maximal values by which $P_{\mathrm{3}}$ can be shifted up or down, while maintaining all of its three roots; the dotted gray lines indicate the shifted function with $x=\pm \mathrm{2}\mathit{\alpha }/\sqrt{\mathrm{27}}$. The purple lines labeled $y_{min}$ and $y_{max}$ mark the right and left boundaries for the roots $y_{\mathrm{l}}$ and $y_{\mathrm{r}}$, respectively: $y_{\mathrm{l}}$ and $y_{\mathrm{r}}$ can never be located in between the two purple lines. (b) Trajectories of the non-autonomous model with $\mathit{\gamma }\left(t\right)=\mathrm{0.8}\sin (t/{\mathit{\tau }}_{\text{f}})$ and ${\mathit{\tau }}_{\text{f}}=\mathrm{1000}$, plotted in the $(x,y)$ phase plane; the trajectories are colored by their starting times $\mathit{\{}{t}_{\mathrm{0}}=-\mathrm{20}\mathrm{000},t_{\mathrm{1}}=-\mathrm{16}\mathrm{000},t_{\mathrm{2}}=-\mathrm{13}\mathrm{000},t_{\mathrm{3}}=-\mathrm{7000}\mathit{\}}$, and the initial positions were drawn from a standard Gaussian bivariate distribution. (c) The slow time-dependent forcing $\mathit{\gamma }\left(t\right)=\mathrm{0.8}\sin (t/\mathrm{1000})$. (d, e) The same trajectories as in (b) but plotted in time as $y=y\left(t\right)$ and $x=x\left(t\right)$, respectively. Panels (f–h) are same as panels (c)–(e) but for the fast time-dependent forcing $\mathit{\gamma }\left(t\right)=\mathrm{0.8}\sin (t/\mathrm{350})$​​​​​​​. The gray shading in panels (c)–(h) indicates intervals during which $\left|\mathit{\gamma }\right|>\sqrt{\mathrm{1}/\mathrm{3}}$, and the internal oscillation is hence suppressed.

[Figure omitted. See PDF]

This behavior can be better understood by considering the nullclines of Eqs. () and () in the $(x,y)$ plane, as shown in Fig. . If the branches of the $y$ nullcline that correspond to $y_{\mathrm{l}}$ and $y_{\mathrm{r}}$, and thus to stable fixed points of Eq. () for a given value of $x$, intersect with the $x$ nullcline given by $y=\mathit{\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 $\mathit{\{}y=\mathit{\gamma }\mathit{\}}$ from $x$. See 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. ()–(), in the case of the FHN model, the radius of the oscillations actually grows exponentially over a small range of $\mathit{\gamma }$ values right after the bifurcation point and then stabilizes. These exponentially growing stable limit cycles have been termed “canard cycles” . have pointed out a possible link between canard cycles and D-O cycles, and they play a role in other excitable climate models; see and references therein.

3.3.2 A pullback attractor of a periodically forced FHN model

Introducing a sinusoidal time dependence

24$$\mathit{\gamma }\to \mathit{\gamma }\left(t\right)\propto \sin (t/{\mathit{\tau }}_{\text{f}})$$ into the slow Eq. () makes the system non-autonomous, as discussed in general terms in Sect. , and it periodically switches the self-oscillatory behavior on and off. We consider here the case in which the variations of the external forcing $\mathit{\gamma }(t;{\mathit{\tau }}_{\text{f}})$ occur on a slower timescale than the entire internal dynamics, i.e., 25$${\mathit{\tau }}_{\text{f}}>{\mathit{\tau }}_x>{\mathit{\tau }}_y.$$ Note that we give up here on the restriction of strict timescale separation between $x\left(t\right)$ and $y\left(t\right)$, as expressed by Eq. ().

The trajectories plotted in Fig. c–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 ${\mathit{\tau }}_{\text{f}}$. These trajectories illustrate the applicability of the pullback perspective suggested by Sect.  and Fig.  to the periodically forced FHN model. In contrast to the illustrative examples of Sect. , 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. c–e, the forcing's timescale is much slower then the internal timescales. The amplitude of the forcing $\mathit{\gamma }\left(t\right)$ was chosen such that it repeatedly crosses the two thresholds $\pm \sqrt{\mathrm{1}/\mathrm{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” .

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 – i.e., the family of invariant snapshots $\mathit{\left\{}\mathcal{A}\right(t){\mathit{\}}}_{t\in \mathbb{R}}$ of Eq. () – is an infinite repetition of the common trajectory structure that can be observed in Fig. d, e between $-\mathrm{5000}$ and 15 000 time units. In the case at hand, each snapshot $\mathcal{A}\left(t\right)$ consists of a single point.

For Fig. f–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 ${\mathit{\tau }}_{\text{f}}/\mathrm{2}$. This symmetry reflects the symmetry of the stable nullcline of the fast system component, as shown in Fig. .

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. g, h between $-\mathrm{5000}$ and 15 000 time units. In contrast to the slow-forcing case, each snapshot $\mathcal{A}\left(t\right)$ 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. , 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: $$\begin{array}{}\text{26a}& {\displaystyle }\text{d}x={\displaystyle \frac{\mathrm{1}}{{\mathit{\tau }}_x}}(y-\mathit{\gamma }(t\left)\right)\text{d}t,\text{26b}& {\displaystyle }\text{d}y={\displaystyle \frac{\mathrm{1}}{{\mathit{\tau }}_y}}\left[\mathit{\alpha }\left(y-y^{\mathrm{3}}\right)-x\right]\text{d}t+\mathit{\sigma }\text{d}\mathit{\eta }.\end{array}$$ Here, $\mathit{\eta }$ denotes a Wiener process, as in Eq. () – i.e., a continuous stochastic process whose increments $\text{d}\mathit{\eta }$ are independently and normally distributed, with mean zero and unit variance – and $\mathit{\gamma }\left(t\right)$ remains a periodic forcing of the $x$ component proportional to $\sin \left({\mathit{\tau }}_{\text{f}}t\right)$.

Figure 11

Random attractor of the periodically and stochastically forced FHN model governed by Eq. (). (a) Results for slow periodic forcing with period ${\mathit{\tau }}_{\text{f}}=\mathrm{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. (); the noise variance $\mathit{\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 ${\mathit{\tau }}_{\text{f}}=\mathrm{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 $\mathit{\sigma }=\mathrm{0.7}$, the synchronization already takes place during the first 8000 time units.

[Figure omitted. See PDF]

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. . For both the slower and faster deterministic forcing – i.e., ${\mathit{\tau }}_{\text{f}}=\mathrm{1000}$ and ${\mathit{\tau }}_{\text{f}}=\mathrm{350}$, respectively, as studied in the previous paragraphs – several random attractors are approximated for increasing noise variance values $\mathit{\sigma }$ in Eq. (). 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 ${\mathit{\tau }}_{\text{f}}=\mathrm{1000}$, the trajectories in Fig. a converge fairly rapidly to a single one, as in the case with no noise in Fig. d, 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 ${\mathit{\tau }}_{\text{f}}=\mathrm{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 ${\mathit{\tau }}_{\text{f}}=\mathrm{1000}$, a rapid convergence of the trajectories can be observed for all noise variances. However, two separate 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.,. 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 $\mathit{\sigma }=\mathrm{0.7}$ and $\mathit{\sigma }=\mathrm{0.8}$ in Fig. b.

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.,. 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. The situation studied here is rather different, as shown across Figs. –.

3.3.4 An FHN model of the NGRIP record

Readers who are familiar with the NGRIP ${\mathit{\delta }}^{\mathrm{18}}$O record 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. d. In particular, the prominent sawtooth pattern of the data is satisfactorily captured by the fast and slow dynamics of the model.

Figure 12

FHN model fit to the NGRIP ${\mathit{\delta }}^{\mathrm{18}}$O dataset of 20-year means, as published as a supplement to (see also ). Originally, the reported a 50-year-mean dataset. (a) The fast component $y_{\text{FHN}}\left(t\right)$ of an FHN model forced with historical CO${}_{\mathrm{2}}$ concentrations (orange), together with the observed ${\mathit{\delta }}^{\mathrm{18}}$O record (blue) from the NGRIP ice core . (b) Rescaled atmospheric CO${}_{\mathrm{2}}$ concentrations from Antarctic ice cores in arbitrary units (au) ; the horizontal red lines indicate the lower and upper bounds, $y_{\text{min}}=-\sqrt{\mathrm{1}/\mathrm{3}}$ and $y_{\text{max}}=\sqrt{\mathrm{1}/\mathrm{3}}$, respectively, of the free-oscillation regime.

[Figure omitted. See PDF]

Figure  shows a trajectory of the FHN model for which the sinusoidal forcing used in Fig.  was replaced by a rescaled time series of atmospheric CO${}_{\mathrm{2}}$ concentrations retrieved from Antarctic ice cores : 27$$\mathit{\gamma }\left(t\right)\propto {\mathrm{CO}}_{\mathrm{2}}\left(t\right).$$ 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${}_{\mathrm{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 , the purely deterministic FHN model is not meant to reproduce the exact pattern of D-O events. have carried out a detailed study on the use of a CO${}_{\mathrm{2}}$-forced FHN-type model to simulate D-O variability.

In fact, Fig. 6 describe in detail a somewhat more complex, proxy-record-based picture of the interaction between episodes full of D-O events, 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\left(t\right)$ 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, , Table 5, and 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.,, linked in turn to varying ice shelf extent e.g.,, or as the weakening and strengthening of Antarctic bottom water production .

The interaction between the fast variable and the slow one happens here in the presence of a climate forcing represented by CO${}_{\mathrm{2}}$ concentration. On the much slower timescales of Quaternary glaciations, an interplay between the CO${}_{\mathrm{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 who first applied pullback ideas to the problems of Quaternary glaciations, independently of earlier work on the topic in the climate literature . 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.  and Fig. 11.2.

Based on the considerable success of NDS and RDS applications to other climate problems – such as ENSO , the wind-driven ocean circulation , or the evaluation of the ensemble simulations routinely performed in support of the Assessment Reports of the Intergovernmental Panel on Climate Change – it would appear worthwhile to proceed further along these lines.

Box  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 of its predictability on the various paleoclimatic time scales mentioned in the list above.

Box 2

Some open questions concerning Quaternary glaciations.

[Figure omitted. See PDF]

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. and , among others, have emphasized the need to pursue such a model hierarchy, as originally proposed by , 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.  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 , 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. .

Figure 13

Glacial–interglacial cycles simulated by the modified Daruka and Ditlevsen model of Eqs. () and (): (a) 21 June insolation $F\left(t\right)$ at 65${}^{\circ }$ N, normalized to have mean zero and unit standard deviation over the last 1000 kyr, according to . (b) Slowly changing parameters $\mathit{\alpha }\left(t\right)$ and $\mathit{\beta }\left(t\right)$ introduced to give rise to the MPT. (c) Simulated glacial–interglacial cycles ${\mathit{\delta }}^{\mathrm{18}}$O${}_{\text{model}}$ (red) in comparison with the benthic ${\mathit{\delta }}^{\mathrm{18}}$O data (blue); the latter dataset, provided by , is also known as the LR04 stack.

[Figure omitted. See PDF]

Among these glacial-cycle models, the model of 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 ${\mathit{\delta }}^{\mathrm{18}}$O proxy reconstruction of glacial–interglacial cycles due to , especially the timing of glacial terminations; compare our Fig.  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 $$\begin{array}{}\text{28a}& {\displaystyle }\mathit{\tau }\dot{x}=\mathit{\lambda }y,\text{28b}& {\displaystyle }\mathit{\tau }\dot{y}=-\mathit{\alpha }\left(t\right)+x-x^{\mathrm{3}}-\mathit{\beta }\left(t\right)F\left(t\right)x-\mathit{\kappa }y;\end{array}$$ here $t$ is the time (in kyr) and $F\left(t\right)$ is the normalized 21 June insolation at 65${}^{\circ }$ N, based on the calculations of , as shown in Fig. a. The constant parameter values are chosen as $\mathit{\kappa }=\mathrm{1}$, $\mathit{\tau }=\mathrm{100}$, and $\mathit{\lambda }=\mathrm{10}$. Note that with this choice of $\mathit{\lambda }$ and unlike in the FHN model of D–O oscillations in Sect. , $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 $\mathit{\kappa }$ in Eq. (), 29$$\mathit{\kappa }\left(t\right)={\mathit{\kappa }}_{\mathrm{1}}+\mathrm{0.5}\left({\mathit{\kappa }}_{\mathrm{0}}\mathrm{-}{\mathit{\kappa }}_{\mathrm{1}}\right)\left(\mathrm{1.0}\mathrm{-}\tanh \left({\displaystyle \frac{t\mathrm{-}{t}_{\mathrm{0}}}{t_s}}\right)\right).$$ In our M-DD16 model, we introduce instead a slow change in the parameters $\mathit{\alpha }\left(t\right)$ and $\mathit{\beta }\left(t\right)$ of Eq. () as follows: $$\begin{array}{}\text{30a}& {\displaystyle }\mathit{\alpha }\left(t\right)=\mathrm{2.1}-\mathrm{1.4}\tanh \left((t+\mathrm{1100})/\mathrm{500}\right),\text{30b}& {\displaystyle }\mathit{\beta }\left(t\right)=\mathrm{2.5}+\mathrm{1.4}\tanh \left((t+\mathrm{1100})/\mathrm{500}\right).\end{array}$$ The functions $\mathit{\alpha }\left(t\right)$ and $\mathit{\beta }\left(t\right)$ so defined are plotted in Fig. b, 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 ${\mathit{\delta }}^{\mathrm{18}}$O${}_{\text{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 ${\mathit{\delta }}^{\mathrm{18}}$O${}_{\text{model}}=\mathrm{4.3}-\mathrm{1.4}y+\mathrm{0.0003}t$. Equation () generalizes a first-order ODE system that is equivalent to the Duffing form of the nonlinear spring Eq () discussed in Sect. . The use of such classical first-order systems – like the Duffing and Van der Pol ones – in paleoclimate models was initiated by ; see also Sect.  herein.

Figure c shows a time series of simulated glacial–interglacial changes ${\mathit{\delta }}^{\mathrm{18}}$O${}_{\text{model}}$ (red) in comparison with the benthic ${\mathit{\delta }}^{\mathrm{18}}$O (blue) of . The model's initial condition is taken to be $x=-\mathrm{1}$ and $y=\mathrm{0}$ at $t=\mathrm{10}\mathrm{000}$ kyr b2k and, since the insolation forcing in Fig. a is prescribed as a time series with 100-year sampling, we solved Eq. () using Heun's predictor–corrector method chap. 8 with a step size of $\mathrm{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. b, 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.

Figure 14

PBA of the M-DD16 model governed by Eq. (), approximated by 40 trajectories starting from random initial conditions in $(x,y)\in [-\mathrm{2},\mathrm{2}]\times [-\mathrm{2},\mathrm{2}]$ at $t=-\mathrm{10}\mathrm{000}$ kyr. (a) The PBA corresponding to Fig. , with the slowly changing parameters $\mathit{\alpha }\left(t\right)$ and $\mathit{\beta }\left(t\right)$ given by Eq. (). (b) The PBA for $\mathit{\alpha }\left(t\right)$ and $\mathit{\beta }\left(t\right)$ kept constant at the post-MPT values of $\mathit{\alpha }=\mathrm{0.7}$ and $\mathit{\beta }=\mathrm{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 $\mathit{\alpha }$ and $\mathit{\beta }$ in time, the overall structure of the PBA is similar before and after 1000 kyr b2k.

[Figure omitted. See PDF]

We next approximate the PBA by taking 40 random initial conditions at 10 Ma b2k and integrating the model of Eqs. () and () up to the present time. The PBA in this case is simply a moving fixed point, as plotted in Fig. a, 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 $\mathit{\alpha }$ and $\mathit{\beta }$ fixed at their post-MPT values $\mathit{\alpha }=\mathrm{0.7}$ and $\mathit{\beta }=\mathrm{3.9}$ throughout the simulation interval and repeating the computation of the PBA, a more complex picture arises. In the latter case, Fig. b shows a bunching of trajectories into separate fuzzy clusters, subject to the quasi-periodic orbital forcing of Fig. a.

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 . It also appears to be consistent with the 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, while several distinct trajectories survived in the Late Pleistocene after the MPT.

Second, the separate bundles or “ropes” of trajectories in Fig. b seem to point to the type of generalized synchronization discussed in the paleoclimate context by and in the context of interannual and interdecadal climate variability by and . 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.,.

5 Conclusions

In this review and research paper, we have covered the contributions of the 1970s to the rebirth of the theory of the ice ages in Sect.  and the 1980s advances in modeling the Quaternary climate's intrinsic variability in Sect. . In Sect. , we presented first results on the effects of the orbital insolation forcing of the data discussed in Sect.  on the intrinsic variability of the data discussed in Sect. , 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  concluded with the formulation and study of a FitzHugh–Nagumo (FHN)-type model of recurrent Dansgaard–Oeschger (D-O) events, in which historical CO${}_{\mathrm{2}}$ concentrations induced episodes of D-O events alternating with episodes of their absence in excellent qualitative agreement with NGRIP ${\mathit{\delta }}^{\mathrm{18}}$O data; see Fig. .

Finally, in Sect. , we listed a number of open issues on Quaternary and longer paleoclimate timescales and proposed to address them by using the tools of Sect. . This approach was illustrated by a Duffing-type model of , 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. with the orbital forcing that is illustrated in Fig.  of Sect.  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.  – one can see the work that was reviewed and extended in this paper as a confirmation of the fine intuition of , 6 decades ago, as summarized in and further expanded by 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 [in ] and ), 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. 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 . Subsequently, conceptual dynamical models were further developed e.g.,. Some of the more recent models simulate the proxy records of glacial cycles remarkably well .

Shortly after presented his work, 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 and further with the isostatic response of the underlying bedrock . found self-sustained oscillations in their simple coupled ice sheet–energy balance model. 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, 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 .

A deeper understanding of glaciation cycles cannot be obtained without process-based models that focus on the detailed physics and biogeochemical phenomena involved . 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 . As Henri Poincaré pointed out, “mathematics is the art of giving the same name to different things”.

Table A1

List of simple conceptual glacial-cycle models with only 1–3 variables. Note that, while extensive, this list is not exhaustive.

Author(s) (year)	Succinct description	Dynamical properties	Is the mid-Pleistocene transition (MPT)-like behavior generated?	Additional comments
Low-order dynamical systems models of glacial cycles
	Model of the global ice volume with different sensitivities to the insolation during growth and retreat, respectively.	Non-chaotic. No self-sustained oscillation in the absence of forcing	Not specifically addressed.	His model can be written as a non-autonomous dynamical system (NDS) .
	Simplified ice sheet model based on the flow law of a perfectly plastic solid.	Non-chaotic. No self-sustained oscillation in the absence of forcing.	Not specifically addressed.
, , ,	Ice sheet–bedrock–temperature coupled oscillator; see the text for details.	Chaotic in the presence of orbital forcing . Exhibits self-sustained oscillations with a period of $\approx$ 10 kyr; the basic mechanism for getting $\sim$ 100 kyr cycle is nonlinear resonance at a difference tone between the 19 and 23 kyr precession cycles ; see Sect.  herein.	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.	is an ice sheet–temperature coupled oscillator. introduced the bedrock and added the orbital forcing.
	Piecewise linear model with two different timescales for ice sheet waxing and waning.	Non-chaotic. No self-sustained oscillation in the absence of forcing.	Not specifically addressed.	400 kyr periodicity is more dominant than 100 kyr periodicity.
(, ), ()	Coupled oscillator based on ice–CO${}_{\mathrm{2}}$–ocean temperature (or NADW) coupling. Nonlinearity is only in the CO${}_{\mathrm{2}}$ component.	Non-chaotic. After the MPT, it exhibits self-sustained oscillations with a period of $\sim$ 100 kyr in the absence of forcing.	Yes, the MPT arises via a Hopf bifurcation in the underlying system; see also .	It has a strange non-chaotic attractor, as well as a chaotic one, depending on the parameter setting .
, Parrenin and Paillard (, )	Hybrid dynamical system with discrete states that switch when conditions are satisfied.	Non-chaotic but the dynamics is sensitive to parameter changes near the switching boundaries . No self-sustained oscillations in the absence of forcing.	Yes, in .	The authors call these models relaxation oscillators, in spite of their discrete-state formulation.
	Delayed differential equation for ice volume,	Robustness of trajectories against random perturbations is mentioned. It exhibits damped oscillations in the absence of forcing.	Not specifically addressed.
()	Coupled oscillator based on Northern Hemisphere ice volume–Antarctic ice extent–CO${}_{\mathrm{2}}$ coupling.	Non-chaotic. It exhibits self-sustained oscillations in the absence of forcing.	Yes, the MPT is induced by a slow drift in the bottom water formation efficiency around the Antarctic.	showed that this model can exhibit chaotic dynamics when it is slightly modified.
,	An ice mass model that simplifies the sea ice switch model of .	It exhibits self-sustained oscillations in the absence of forcing. Non-chaotic under the orbital forcing.	Yes, in , where the maximum ice volume threshold is increased in time according to the regolith hypothesis.	Uses the temperature–precipitation feedback introduced by .
	The ice mass grows monotonically and collapses to zero when it exceeds a threshold modulated by the obliquity or by a hybrid measure consisting of the obliquity and climatic precession .	Non-chaotic. It exhibits self-sustained oscillations in the absence of forcing.	Yes, in .

Table A1

Continued.

Author(s) (year)	Succinct description	Dynamical properties	Is the mid-Pleistocene transition (MPT)-like behavior generated?	Additional comments
	Data-based, phase space model for Pleistocene ice volume with thresholds.	Due to the presence of the threshold, the model is sensitive to changes in parameters or in its position in phase space.	Yes, MPT-like behavior is produced solely by changes in orbital parameters.
,	Forced Van der Pol (VdP) oscillator model; closely related to the FHN model used in Sect.  herein.	Non-chaotic but sensitive to the noise. It exhibits self-sustained oscillations in the absence of forcing.	Yes, MPT-like behavior is generated via a Hopf bifurcation with an explosive character .	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.
	One-dimensional phase oscillator model.	Non-chaotic.	Yes, MPT-like frequency change accompanies a smooth or non-smooth saddle-node bifurcation of tori.	It has a strange non-chaotic attractor or a quasi-periodic attractor in a classical sense, i.e., not in the pullback sense.
()	Two-dimensional forced limit cycle oscillator.	Non-chaotic; overall stability of simulated glacial cycles against dynamical noise is reported.	It exhibits MPT-like behavior via a trans-critical bifurcation of the slow manifold in the fast dynamics.
	Forced two-dimensional oscillator consisting of ocean alkalinity and calcifier population.	Non-chaotic or chaotic given periodic forcing, depending on the parameters. It exhibits self-sustained oscillations in the absence of forcing.	Yes.
()	Forced Duffing oscillator–type model.	No self-sustained oscillations in the absence of forcing. Can be chaotic depending on parameter values; the authors call this the “climatic butterfly effect”.	Yes, MPT-like behavior is induced by a slow change in the damping coefficient parameter $\mathit{\kappa }$.	See also Sect.  herein.
()	Coupled system of ice volume, temperature, and atmospheric CO${}_{\mathrm{2}}$ incorporating a delayed CO${}_{\mathrm{2}}$ contribution from ocean ridge volcanism. It is reduced to a forced 1-D delay differential equation.	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.	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.
	Scalar delay differential equation for ice volume, derived from the model via the linear chain approximation.	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.	Yes, MPT-like behavior is induced by summer insolation forcing as a transition from small-amplitude $\sim$ 41 kyr cycles to large-amplitude $\sim$ 100 kyr cycles.
	Study of several low-order dynamical systems models, including , , and a generalized VdP–Duffing oscillator.	Chaotic or non-chaotic depending on the model and its parameter values.	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.	The authors stress the possibility of chaotic dynamics occurring.

Table A1

Continued.

Author(s) (year)	Succinct description	Dynamical properties	Is the mid-Pleistocene transition (MPT)-like behavior generated?	Additional comments
	Coupled model relying on ice sheet basal temperature–ocean temperature coupling, deduced from physical laws via a scaling analysis.	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.	Yes, MPT-like behavior is induced by an enhancement of positive feedbacks against negative feedbacks in the model.
()	Ice volume–CO${}_{\mathrm{2}}$ coupled conceptual model.		Not specifically addressed.	This model was developed to simulate the last 800 kyr of glacial cycles and Earth's future climate on the million-year timescale.
Low-order dynamical models in which stochastic processes are of the essence
,	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.	Additive stochastic perturbations play a key role.	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.	and proposed the stochastic-resonance idea for ice age cyclicity independently of each other.
	Generalized stochastic resonance model with obliquity and precessional forcing.	Additive stochastic perturbations play a key role.	Not specifically addressed.
	Model based on the coherence resonance mechanism uses a single temperature ODE with a discrete delay.	Additive stochastic perturbations play a key role.	Yes.	The delay feedback in the temperature is based on ice sheet extent reconstructions and affects the model's albedo.
	Generalized stochastic resonance model with obliquity and precessional forcing.	Additive stochastic perturbations play a key role.	Yes.	The model has an underlying bifurcation structure following .

For example, coupled nonlinear oscillators frequently exhibit synchronization with simple frequency ratios, either with each other or with the forcing . Thus, mathematical models that ignore many physical details may also help us elucidate emergent properties in paleoclimatic dynamics or test different orbital hypotheses . Many low-order models of glacial dynamics are listed in Table , 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 B1$$\begin{array}{rl}& {\displaystyle \frac{\mathrm{d}\mathit{\rho }}{\mathrm{d}t}}=\mathit{\alpha }(\mathit{\mu }+\mathit{\beta }\sin (\mathit{\nu }t)-\mathit{\rho }),\\ & {\displaystyle \frac{\mathrm{d}\mathit{\varphi }}{\mathrm{d}t}}=\mathit{\omega }\text{with}\mathit{\rho }>\mathrm{0},\mathit{\mu }>\mathrm{0},\end{array}$$ that were introduced in Sect.  and analytically derive their invariant sets B2$$\begin{array}{rl}& \mathcal{A}\left(t\right)=\mathit{\{}\left(\mathit{\alpha }\mathit{\beta }\sin (\mathit{\nu }t+\mathit{\vartheta })+\mathit{\mu },\mathit{\phi }\right):\mathit{\phi }\in [\mathrm{0},\mathrm{2}\mathit{\pi }\left)\mathit{\right\}}\\ & \forall t\in \mathbb{R},\end{array}$$ as well as the corresponding pullback attractor (PBA). Following , the PBA is given by the family B3$$\colorbox[rgb]{}{$$}=\mathit{\{}\mathcal{A}{\mathit{\}}}_{t:t\in \mathbb{R}}.$$

First, we define $\mathrm{\Delta }\mathit{\rho }\left(t\right)=\mathit{\rho }\left(t\right)-\mathit{\mu }$, which gives rise to B4$${\displaystyle \frac{\text{d}\mathrm{\Delta }\mathit{\rho }\left(t\right)}{\text{d}t}}=-\mathit{\alpha }\mathrm{\Delta }\mathit{\rho }\left(t\right)+\mathit{\alpha }\mathit{\beta }\sin \left(\mathit{\nu }t\right).$$ This is an inhomogeneous ODE and can thus be solved by the variation of parameters method e.g.,. The ansatz B5$$\mathrm{\Delta }\mathit{\rho }\left(t\right)=c\left(t\right)e^{-\mathit{\alpha }(t-t_{\mathrm{0}})},$$ yields B6$${\displaystyle \frac{\text{d}\mathrm{\Delta }\mathit{\rho }\left(t\right)}{\text{d}t}}=-\mathit{\alpha }\mathrm{\Delta }\mathit{\rho }\left(t\right)+{\displaystyle \frac{\text{d}c\left(t\right)}{\text{d}t}}{e}^{-\mathit{\alpha }(t-t_{\mathrm{0}})}.$$ A comparison with Eq. () requires B7$${\displaystyle \frac{\text{d}}{\text{d}t}}c\left(t\right)=\mathit{\alpha }\mathit{\beta }\sin \left(\mathit{\nu }t\right)e^{+\mathit{\alpha }(t-t_{\mathrm{0}})},$$ and hence B8$$\begin{array}{rl}c\left(t\right)& =\underset{t_{\mathrm{0}}}{\overset{t}{\int }}\mathit{\alpha }\mathit{\beta }\sin \left(\mathit{\nu }{t}^{\prime }\right)e^{+\mathit{\alpha }(t^{\prime }-t_{\mathrm{0}})}\text{d}{t}^{\prime }+\mathit{\gamma }\\ & ={\left[\mathit{\alpha }\mathit{\beta }\sin \left(\mathit{\nu }{t}^{\prime }\right){\displaystyle \frac{\mathrm{1}}{\mathit{\alpha }}}{e}^{+\mathit{\alpha }(t^{\prime }-t_{\mathrm{0}})}\right]}_{t_{\mathrm{0}}}^t\\ & -\underset{t_{\mathrm{0}}}{\overset{t}{\int }}{\displaystyle \frac{\mathit{\nu }}{\mathit{\alpha }}}\mathit{\alpha }\mathit{\beta }\cos \left(\mathit{\nu }{t}^{\prime }\right)e^{+\mathit{\alpha }(t^{\prime }-t_{\mathrm{0}})}\text{d}{t}^{\prime }+\mathit{\gamma }.\end{array}$$

Repeated partial integration yields B9$$\begin{array}{rl}\underset{t_{\mathrm{0}}}{\overset{t}{\int }}& {\displaystyle \frac{\mathit{\nu }}{\mathit{\alpha }}}\mathit{\alpha }\mathit{\beta }\cos \left(\mathit{\nu }{t}^{\prime }\right)e^{+\mathit{\alpha }(t^{\prime }-t_{\mathrm{0}})}\text{d}{t}^{\prime }=\\ & {\left[{\displaystyle \frac{\mathit{\nu }}{\mathit{\alpha }}}\mathit{\alpha }\mathit{\beta }\cos \left(\mathit{\nu }{t}^{\prime }\right){\displaystyle \frac{\mathrm{1}}{\mathit{\alpha }}}{e}^{+\mathit{\alpha }(t^{\prime }-t_{\mathrm{0}})}\right]}_{t_{\mathrm{0}}}^t\\ & +\underset{t_{\mathrm{0}}}{\overset{t}{\int }}{\displaystyle \frac{{\mathit{\nu }}^{\mathrm{2}}}{{\mathit{\alpha }}^{\mathrm{2}}}}\mathit{\alpha }\mathit{\beta }\sin \left(\mathit{\nu }{t}^{\prime }\right)e^{+\mathit{\alpha }(t^{\prime }-t_{\mathrm{0}})}\text{d}{t}^{\prime }.\end{array}$$ Therefore, we find B10$$\begin{array}{rl}\left(\mathrm{1}+{\displaystyle \frac{{\mathit{\nu }}^{\mathrm{2}}}{{\mathit{\alpha }}^{\mathrm{2}}}}\right)& \underset{t_{\mathrm{0}}}{\overset{t}{\int }}\mathit{\alpha }\mathit{\beta }\sin \left(\mathit{\nu }{t}^{\prime }\right)e^{+\mathit{\alpha }(t^{\prime }-t_{\mathrm{0}})}\text{d}{t}^{\prime }=\\ & {\left[\mathit{\alpha }\mathit{\beta }\sin \left(\mathit{\nu }{t}^{\prime }\right){\displaystyle \frac{\mathrm{1}}{\mathit{\alpha }}}{e}^{+\mathit{\alpha }(t^{\prime }-t_{\mathrm{0}})}\right]}_{t_{\mathrm{0}}}^t\\ & -{\left[{\displaystyle \frac{\mathit{\nu }}{\mathit{\alpha }}}\mathit{\alpha }\mathit{\beta }\cos \left(\mathit{\nu }{t}^{\prime }\right){\displaystyle \frac{\mathrm{1}}{\mathit{\alpha }}}{e}^{+\mathit{\alpha }(t^{\prime }-t_{\mathrm{0}})}\right]}_{t_{\mathrm{0}}}^t,\end{array}$$

and finally B11$$\begin{array}{rl}c\left(t\right)& ={\displaystyle \frac{\mathrm{1}}{\left(\mathrm{1}+\frac{{\mathit{\nu }}^{\mathrm{2}}}{{\mathit{\alpha }}^{\mathrm{2}}}\right)}}\left(\underset{=\mathit{\beta }\left[\sin \left(\mathit{\nu }t\right)e^{\mathit{\alpha }(t-t_{\mathrm{0}})}-\sin \left(\mathit{\nu }{t}_{\mathrm{0}}\right)\right]}{\underbrace{{\left[\mathit{\beta }\sin \left(\mathit{\nu }{t}^{\prime }\right)e^{+\mathit{\alpha }(t^{\prime }-t_{\mathrm{0}})}\right]}_{t_{\mathrm{0}}}^t}}\right.\\ & \left.-\underset{=\frac{\mathit{\nu }\mathit{\beta }}{\mathit{\alpha }}\left[\cos \left(\mathit{\nu }t\right)e^{\mathit{\alpha }(t-t_{\mathrm{0}})}-\cos \left(\mathit{\nu }{t}_{\mathrm{0}}\right)\right]}{\underbrace{{\left[{\displaystyle \frac{\mathit{\nu }\mathit{\beta }}{\mathit{\alpha }}}\cos \left(\mathit{\nu }{t}^{\prime }\right)e^{+\mathit{\alpha }(t^{\prime }-t_{\mathrm{0}})}\right]}_{t_{\mathrm{0}}}^t}}+\mathit{\gamma }\right).\end{array}$$

Plugging this result into the ansatz (Eq. ) yields B12$$\begin{array}{rl}\mathrm{\Delta }\mathit{\rho }(t,t_{\mathrm{0}})& ={\displaystyle \frac{\mathrm{1}}{\left(\mathrm{1}+\frac{{\mathit{\nu }}^{\mathrm{2}}}{{\mathit{\alpha }}^{\mathrm{2}}}\right)}}\left(\mathit{\beta }\left[\sin \left(\mathit{\nu }t\right)-\sin \left(\mathit{\nu }{t}_{\mathrm{0}}\right)e^{-\mathit{\alpha }(t-t_{\mathrm{0}})}\right]\right.\\ & -\left.{\displaystyle \frac{\mathit{\nu }\mathit{\beta }}{\mathit{\alpha }}}\left[\cos \left(\mathit{\nu }t\right)-\cos \left(\mathit{\nu }{t}_{\mathrm{0}}\right)e^{-\mathit{\alpha }(t-t_{\mathrm{0}})}\right]\right)\\ & +\mathit{\gamma }{e}^{-\mathit{\alpha }(t-t_{\mathrm{0}})},\end{array}$$ with the initial conditions B13$$\mathrm{\Delta }\mathit{\rho }(t_{\mathrm{0}},t_{\mathrm{0}})=\mathit{\gamma }.$$ In the pullback limit, all the terms that carry a factor $e^{-\mathit{\alpha }(t-t_{\mathrm{0}})}$ vanish, and thus B14$$\begin{array}{rl}\underset{t_{\mathrm{0}}\to -\mathrm{\infty }}{lim}\mathrm{\Delta }\mathit{\rho }(t,t_{\mathrm{0}})& ={\displaystyle \frac{\mathrm{1}}{\left(\mathrm{1}+\frac{{\mathit{\nu }}^{\mathrm{2}}}{{\mathit{\alpha }}^{\mathrm{2}}}\right)}}\left(\mathit{\beta }\sin \left(\mathit{\nu }t\right)-{\displaystyle \frac{\mathit{\nu }\mathit{\beta }}{\mathit{\alpha }}}\cos \left(\mathit{\nu }t\right)\right)\\ & =\mathit{\alpha }\mathit{\beta }\sin (\mathit{\nu }t+\mathit{\vartheta }),\end{array}$$ with $\mathit{\vartheta }=\arctan (-\mathit{\nu }/\mathit{\alpha })$. For comparison, the modulation of the target radius itself was given by $\mathit{\beta }\sin \left(\mathit{\nu }t\right)$, and hence it is amplified by the factor of $\mathit{\alpha }$. Since $\mathit{\rho }$ is restricted to positive values, this solution requires $\mathit{\alpha }\mathit{\beta }<\mathit{\mu }$.

Since the evolution in time of the phase $\mathit{\varphi }\left(t\right)$ 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 B15$$\begin{array}{rl}\mathcal{A}\left(t\right)& =\mathit{\left\{}\right(\mathit{\rho }\left(t\right),\mathit{\varphi }):\mathit{\varphi }\in [\mathrm{0},\mathrm{2}\mathit{\pi }\left)\mathit{\right\}}\\ & =\mathit{\left\{}\right(\mathit{\alpha }\mathit{\beta }\sin (\mathit{\nu }t+\mathit{\vartheta }),\mathit{\varphi }):\mathit{\varphi }\in [\mathrm{0},\mathrm{2}\mathit{\pi }\left)\mathit{\right\}}.\end{array}$$ Defined as the indexed family of all $\mathcal{A}\left(t\right)$, the system's PBA is comprised of the family of circles B16$$\begin{array}{rl}\colorbox[rgb]{}{$$}& =\mathit{\left\{}A\right(t){\mathit{\}}}_{t\in \mathbb{R}}\\ & =\mathit{\left\{}\right(\mathit{\alpha }\mathit{\beta }\sin (\mathit{\nu }t+\mathit{\vartheta }),\mathit{\varphi }):\mathit{\varphi }\in [\mathrm{0},\mathrm{2}\mathit{\pi }){\mathit{\}}}_{t\in \mathbb{R}}.\end{array}$$

Code and data availability

All code used to generate the figures presented in this article is available from the authors upon request. The NGRIP ${\mathit{\delta }}^{\mathrm{18}}$O – originally published in – and the historical CO${}_{\mathrm{2}}$ data shown in Fig.  are available as a supplement to 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, ), respectively. The benthic ${\mathit{\delta }}^{\mathrm{18}}$O data shown in Fig.  can be obtained from 10.1594/PANGAEA.704257 , and the 21 June insolation at 65${}^{\circ }$ N can be downloaded from http://vo.imcce.fr/insola/earth/online/earth/online/index.php (last access: 13 April 2022, ).

Video supplement

The video supplement to this article (10.5281/zenodo.6346211, ) illustrates the pullback attractor (PBA) associated with the simple system governed by Eqs. () and (). It shows a heat map of the phase plane, derived from an increasing number of trajectories with common initial time $t_{\text{i}}=-\mathrm{200}$ and final time $t_{\text{f}}=\mathrm{200}$. The initial radius and phase of each trajectory are randomly sampled from Gaussian distributions centered at ${\mathit{\mu }}_{\mathit{\rho }}=\mathrm{20}$ and ${\mathit{\mu }}_{\mathit{\varphi }}=\mathrm{0}$ with standard deviations of ${\mathit{\sigma }}_{\mathit{\rho }}=\mathrm{5}$ and ${\mathit{\sigma }}_{\mathit{\varphi }}=\mathrm{10}$, respectively. Over the course of the video, 100 trajectories are continuously added to the heat map and the annular disk $\mathcal{D}=\mathit{\left\{}\right(\mathit{\rho },\mathit{\varphi }):\mathit{\rho }\in [\mathit{\mu }-\mathit{\alpha }\mathit{\beta },\mathit{\mu }+\mathit{\alpha }\mathit{\beta }]\text{ and }\mathit{\varphi }\in [\mathrm{0},\mathrm{2}\mathit{\pi }\left)\mathit{\right\}}$ fills up. The heat map in Fig. b is a snapshot from this video at time $t=\mathrm{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.