## Abstract

Ice sheets appeared in the northern hemisphere around 3 Ma (million years) ago and glacial–interglacial cycles have paced Earth's climate since then. Superimposed on these long glacial cycles comes an intricate pattern of millennial and sub-millennial variability, including Dansgaard–Oeschger and Heinrich events. There are numerous theories about these oscillations. Here, we review a number of them in order to draw a parallel between climatic concepts and dynamical system concepts, including, in particular, the relaxation oscillator, excitability, slow–fast dynamics and homoclinic orbits. Namely, almost all theories of ice ages reviewed here feature a phenomenon of synchronization between internal climate dynamics and astronomical forcing. However, these theories differ in their bifurcation structure and this has an effect on the way the ice age phenomenon could grow 3 Ma ago. All theories on rapid events reviewed here rely on the concept of a limit cycle excited by changes in the surface freshwater balance of the ocean. The article also reviews basic effects of stochastic fluctuations on these models, including the phenomenon of phase dispersion, shortening of the limit cycle and stochastic resonance. It concludes with a more personal statement about the potential for inference with simple stochastic dynamical systems in palaeoclimate science.

## 1. Introduction

The Pliocene and the Pleistocene cover approximately the past 5 Myr. The climatic fluctuations that characterized this period may be reconstructed from numerous natural archives, including marine, continental and ice core records. These archives show a complex climate history. Ice sheets appeared in the northern hemisphere around 3–3.5 Ma (million years) ago [1,2]. The volume of these ice sheets fluctuated with the variations of the seasonal and spatial distributions of incoming solar radiation (insolation), which are induced by changes in the geometry of the Earth's orbit and the angle (obliquity) between Earth's equator and the ecliptic [3–5]. This is called the astronomical forcing.^{1} Glacial cycles had an average duration of about 40 000 years [6] until about 800 000 years ago. The dominant period of glacial cycles increased around 800 000 years ago and this is referred to as the Middle Pleistocene Transition. Data and models about the Middle Pleistocene Transition are reviewed by Clark *et al*. [7]. Time-series analyses based on band-pass filtering provide further evidence of the nonlinear nature of the climate response to the astronomical forcing, from about 1.4 Ma ago [8]. The latest four glacial cycles, in particular, are distinguished by a pronounced saw-tooth time structure: ice accumulates over the continents during about 80 000 years and then melts in about 10 000 years (figure 1).

Superimposed on these long glacial cycles comes a complex pattern of millennial and sub-millennial variability [11]. For example, the Greenland record features at least 20 events of abrupt rise and slower decline in oxygen isotopic ratio (a proxy for temperature) [12,13] and methane [14] during the latest glacial epoch. These events are known as Dansgaard–Oeschger events. They were found to occur from at least the last glacial inception [15] and the Antarctic ice core records provide evidence that they are characteristic of Pleistocene glacial climates [16]. Some of these events follow pulses of iceberg discharges into the North Atlantic Ocean, called Heinrich events [17–19]. Heinrich events and Dansgaard–Oeschger events have left climatic footprints all over the globe [20], including in Antarctica [16]. The current interglacial period is referred to as the Holocene. It is also characterized by millennial and centennial variability, mainly observed in the North Atlantic [21–23], but of a much weaker amplitude than during the preceding glacial period.

This paper reviews attempts to explain these fluctuations with concepts that originate in dynamical system theory. These are the concepts of limit cycle, synchronization and excitability. The central message of the paper is that current theories of ice ages and rapid events may often be interpreted in terms of generic deterministic models, which are also used in other areas of science such as biology and ecology. However, stochastic parametrizations are an essential part of any complex system model, and their effects on climatic oscillations have to be taken into account.

Dynamical system theory entered palaeoclimate science with idealized models representing the response of ice sheets to the astronomical forcing. These models were directly derived from the physics of the ice sheet–atmosphere system [24–27]. Ghil & Childress [28], in particular, insisted on the interest of analysing such models in terms of bifurcation theory. For modelling the complex carbon cycle response, authors sometimes adopted a more heuristic approach by considering simple models and confronting the results of palaeoclimate evidence [29].

Nowadays, climate research is largely oriented towards large climate simulators (typically general circulation models), which are developed to include as many climate processes as possible. However, thinking in terms of dynamical system theory remains insightful. Indeed, the behaviour of a complex system at a certain spatio-temporal scale is in practice often dominated by a few leading modes, of which the dynamics may be captured fairly convincingly with a low-order dynamical system. Climate scientists are increasingly using this property. For example, they formulate simpler models to explain the seemingly complex behaviours observed in ocean–atmosphere simulators. Examples have been provided in recent years by focusing on interannual [30], centennial [31] and millennial [32,33] variability. In parallel, so-called hysteresis experiments, which aim at identifying the number of stable states in individual components of the climate system such as the ocean circulation [34] or ice sheets [35], contribute to a dynamical system-founded understanding of the climate system. This approach may also help us to predict and communicate about the proximity of bifurcations, which may result in catastrophic climatic changes. Timmermann & Jin [36] termed our ability to anticipate bifurcation phenomena as *predictability of the third kind*, by reference to the predictabilities of the *first* and *second* kind originally introduced by Lorenz [37].

The article is structured as follows. Section 2 reviews some of the basic concepts of oscillator theory. This is no substitute for proper textbook reading, but the reader will find essential notions and definitions needed to understand the remainder. Section 3 reviews how these concepts enter theories of ice ages and rapid events. Section 4 discusses effects of stochastic fluctuations and, finally, §5 is a more personal statement about the potential for inference with simple stochastic dynamical systems in palaeoclimate science.

## 2. Vocabulary and elementary notions

The reader will find an accessible introduction to dynamical system theory and concepts in Strogatz [38]. More formal background on oscillator theory, albeit a bit dated, is available in Guckenheimer & Holmes [39]. Bifurcation and oscillator theory is explicitly connected to climate theory in Ghil & Childress [28], in particular, ch. 12 and Saltzman [40], ch. 7. Background on synchronization and an introduction to the phenomenon of excitability is available in Pikovski *et al*. [41]. Finally, the Scholarpaedia peer-reviewed website is an increasingly rich and authoritative source of information on dynamical systems. Only the notions essential for the present article are summarized here.

*Oscillator*. The oscillator is a dynamical system that has a globally attracting limit cycle. In more simple terms, it oscillates even in the absence of an external drive. Here, we are interested in oscillators to describe climate phenomena, which involve dissipation of energy. The minimal model for a *dissipative* oscillator includes two ordinary differential equations, of which at least one is nonlinear.

*Relaxation oscillator*. The relaxation oscillator is a particular kind of oscillator featuring an interplay between relaxation dynamics (generally fast) and a destabilization process (generally slow). The relaxation is the process by which the system is attracted to a region of the phase space. This evokes the relaxation of a spring. In a relaxation oscillator, the system continues to evolve slowly after the relaxation phase. During this slow evolution phase, the system's stability diminishes gradually until the system is ejected out of its relaxation state, either towards another relaxation state or to the same relaxation state via a dissipative loop. In this review, we will encounter three kinds of relaxation oscillators (figure 2): relaxation founded on slow–fast dynamics (involving a slow manifold); relaxations structured by a homoclinic orbit (involving only one relaxation state); and relaxations structured around a focus. More details are provided in the caption of figure 2.

*Excitability*. An excitable system has a globally attracting fixed point (it does not oscillate spontaneously). However, an external perturbation may have the effect of *exciting* it. During this excitation, the system is being ejected far from its fixed point and then returns to it.

*Link between relaxation dynamics and excitability*. In practice, it is often found that a relaxation oscillator may be transformed into an excitable system by a mere change in parameter, and vice versa. The reason is the following. A relaxation oscillation is often structured globally in the phase space; for example, by a slow manifold (figure 2*a*) or by one or several saddle points (figure 2*b*). Suppose now that the oscillation displayed by such a system ceases because a parameter has been changed. The system is then no longer an oscillator, but the ‘backbone’ of the oscillation dynamics are still latent in the phase space because the elements that structured the limit cycle (the slow manifold or the saddle points) have not disappeared. Consequently, the system may be run on a trajectory close to the defunct limit cycle if it is being pushed by some external force (the excitation) into the region of the phase space previously occupied by this limit cycle. This point is illustrated on the basis of slow–fast dynamics in figure 2*d*, but similar excitation dynamics generally occur near any kind of ‘explosive bifurcation’, that is, bifurcations that give rise rapidly to a fully developed limit cycle. This includes homoclinic, heteroclinic and certain Hopf bifurcations (two examples follow and are illustrated in figure 6).

## 3. Oscillators, relaxation and excitability in palaeoclimates

### (a) Models of ice ages

#### (i) The Saltzman et al. models

Saltzman established a theory in which ice ages are interpreted as a limit cycle synchronized on the astronomical forcing. Saltzman *et al.* [42] wrote a series of articles on the subject, starting with the introduction of the limit cycle idea [43] and synchronization hypothesis, the interpretation of the Middle Pleistocene Transition as a bifurcation [44], and the more complete models in the mid-1990s [45]. The full theory is developed in a book [40]. Here, we concentrate on two intermediate models [46,47]. They are called SM90 and SM91, by reference to the authors (Saltzman & Maasch) and the year of publication. The variables *I*, *μ* and *θ* are the continental ice mass, CO_{2} concentration and deep-ocean temperature, respectively. The reader is referred to the original publications for the meaning and value of the different parameters. They are not crucial here; it suffices to know that they are all positive.
SM90
and
SM91

In both models, the first equation describes the ice mass response to changes in CO_{2} (*μ*) and the astronomical forcing (*F*_{I}(*t*)); Saltzman adopts the so-called Milankovitch view^{2} that an increase in insolation causes a decrease in ice mass. Increases in CO_{2} or in ocean temperature have the same effects.

The other two equations describe the dynamics of CO_{2} and the response of deep-ocean temperature to changes in ice volume. It is further assumed that the mean state of the climate varied slowly throughout the Pliocene–Pleistocene; in particular, in response to a ‘tectonically driven’ decline in the average concentration in CO_{2}, consistently with an earlier proposal [49]. This tectonically driven decline is modelled here as a slow decrease in the forcing term *F*_{μ}(*t*) throughout the Pleistocene.

Consider the bifurcation diagrams of the SM90 and SM91 models with respect to *F*_{μ}, assuming no astronomical forcing (*F*_{I}=0) (figure 3). The systems are then said to be free or autonomous. Depending on *F*_{μ}, both models show regimes with a stable fixed point, and regimes for which the fixed point is unstable, so that the system orbits along a limit cycle.

These considerations led Saltzman to interpret the Middle Pleistocene Transition as a bifurcation between a ‘quasi-linear’ response regime to the astronomical forcing (in the fixed-point regime) to a regime of nonlinear synchronization (resonance) on the astronomical forcing. He concluded that ice ages would occur today even in the absence of astronomical forcing. The main effect of the astronomical forcing is to control the timing of glaciations.

Saltzman's theory is seductive because it explains in a consistent framework several aspects of the Pleistocene climate history, including the change from linear to nonlinear regime [8], the presence of 100 000 year periodicity in climate records [51], the lack of a 400 000 year spectral peak in the ice volume record (such a peak appears in the simple piece-wise linear model devised by Imbrie & Imbrie [52], owing to rectification of the precession signal), the synchronization of deglaciations on the astronomical forcing [53,54], and the occurrence of large climatic transitions even when eccentricity, which modulates the effect of precession on insolation, is at its lowest.

The difficulty for accepting Saltzman's models as a definitive theory lies in the physical interpretation of the CO_{2} equation. This equation encapsulates all the interesting dynamics of the system and it is thus crucial to the theory. Some semi-empirical justification for the CO_{2} equation is given in Saltzman & Sutera [44], but the form of this equation has undergone somewhat ad hoc adjustments in SM90. The form present in SM91 is again different, with important effects on the bifurcation structure, while the authors did not justify this latter change based on physical or biogeochemical considerations.

To better appreciate the structural differences between the two models, let us return to the bifurcation diagrams. Consider SM90. As the forcing is decreased, the fixed point gives rise to a locally unstable limit cycle. The system must, therefore, find a stable limit cycle further away from the fixed point, but in this case not much further. This stable limit cycle is under the influence of the unstable fixed point and, in particular, the system slows down when it passes near it (figure 4). This scenario is depicted in figure 2*c*. The limit cycle evolves as the tectonic forcing is further decreased, until it shrinks smoothly around a perpetually glaciated state.

In SM91, the bifurcation induced by the decrease in tectonic forcing is much more explosive. The system lands on a stable limit cycle that turns out to be little affected by the position of the unstable point. The cycle dynamics do not show clear phases of acceleration and the system cannot be regarded as a relaxation system. The limit cycle disappears abruptly as the tectonic forcing is further decreased, through a phenomenon called a saddle bifurcation of cycles. The consequences of the difference between the bifurcation structures of SM90 and SM91 may be further appreciated in the transient experiments shown in figure 5.

#### (ii) Paillard's (1998) ice age model (P98)

Paillard & Labeyrie [55] have been advocating the concept of relaxation for understanding palaeoclimate dynamics, both ice ages and the more abrupt events, since the publication of a seminal paper in 1994. We return to this article later on, and concentrate on another article published in 1998, in which Paillard [56] introduces a conceptual model of ice ages. Ice volume dynamics respond to an ordinary differential equation,
P98
In this equation, the ice volume *x* is linearly *relaxed* to *x*_{R} with characteristic relaxation time *τ*_{R}. This relaxation process is further perturbed by the astronomical forcing *F*(*t*) with a characteristic time *τ*_{f}. Such a system is said to be hybrid [57] because the relaxation equation involves a discrete state variable, here denoted by *y*. Its state may be ‘deep glacial’ (*G*), ‘mild glacial’ (*g*) or ‘interglacial’ (*i*). The numerical values of *x*_{R} and *τ*_{R} depend on this climate state. Climate states *y* follow a sequence according to a set of conditions formulated on the level of glaciation *x* and insolation. Namely, the transition is triggered when the forcing *F*(*t*) exceeds a certain threshold. Occurrence of *G* drives climate quickly into an interglacial state *i* because *x*_{R}(*G*) and *τ*_{R}(*G*) are specified in the model to be low.

Paillard is not very specific about the physical meaning of the discrete variable, but it accommodates the paradigm that the Atlantic Ocean circulation has gone through three different states during the latest glacial period: intermediate circulation, shut down of the circulation and modern, deep-sinking circulation. The system (P98) features the concept of slow–fast relaxation dynamics. However, this is *not* an oscillator because the shift from *g* to *G* is determined by the course of the external forcing. The Middle Pleistocene Transition is induced in (P98) in a fashion similar to that in Saltzman, and on the basis of similar physical assumptions (tectonically driven decline in CO_{2}). The drift in climatic conditions induced by tectonics is accounted for by a term added to the astronomical forcing. In a later review, Paillard [58] further emphasizes empirical evidence for the relevance of the relaxation concept in the phenomenon of deglaciation.

#### (iii) The Gildor–Tziperman model

Gildor & Tziperman [59] take a moderate step towards higher model complexity by considering a slightly more explicit representation of atmosphere, ocean, sea ice and land ice dynamics. Namely, the ocean is divided into eight boxes, and the atmosphere into four boxes. The sea ice fraction responds to standard energy balance equations. More crucially, land ice growth is influenced by a somewhat controversial feedback between sea ice and precipitation. The feedback is controversial because it is assumed that cold climate results in a *reduction* in ice volume: sea ice growth causes a reduction in precipitation in ice-covered areas and, by this mechanism, almost suppresses accumulation of snow on ice sheets. The latter then no longer compensates for ice ablation and ice volume shrinks.

A free oscillation arises from the fact that the ice volume thresholds for switching sea ice cover ‘on’ and ‘off’ differ. In other words, sea ice displays a hysteresis response to variations in ice volume. This is exactly the principle of the slow–fast relaxation oscillator depicted in figure 2*a*: the curve of equilibrium of sea ice with respect to ice volume is the slow manifold, and ice volume integrates the state of sea ice in time. In turn, this oscillation can be synchronized on the astronomical forcing.

The Gildor–Tziperman model is coupled to a biogeochemical cycle in a companion paper [60], but the essential dynamics of the glacial oscillation are unchanged. Tziperman *et al.* [61] further comment on the model and its property of synchronization on the astronomical forcing, and find that its behaviour is essentially reducible to a hybrid dynamical system.

#### (iv) The Paillard–Parrenin model

Paillard & Parrenin [62] proposed yet another relaxation model in 2004 (PP04). The prognostic variables are ice volume *I*, the area of the Antarctic continental ice sheet *A* and the atmospheric concentration in CO_{2} (*μ*) (*a*…*j* are parameters):
PP04

As in the other ice age models, ice volume is a slow variable driven by the astronomical forcing. It is here coupled to a variable with a similar time scale (*τ*_{A}∼*τ*_{I}) and a faster one (*τ*_{μ}=*τ*_{I}/3). The term *H*(−*D*), where *H* is the Heaviside function, represents the ventilation of the Southern Ocean. CO_{2} is released into the atmosphere when the Southern Ocean is ventilated (*D*<0), which drives deglaciation. Ice then grows slowly, until a Southern Ocean ventilation flush sends the system back to interglacial conditions. Ocean ventilation is thus the fast process in this model and it is the only nonlinear process accounted for. However, contrary to the Gildor–Tziperman model, it *does not* present a hysteresis behaviour. Consequently, the glacial cycles featured by this model cannot be interpreted in terms of shifts between the branches of a slow manifold.

To better understand the dynamics of glacial cycles in this model, we consider the bifurcation diagram along typical solutions in the phase space for the free (i.e. unforced) system (figure 6). The parameter *g* is taken in this example as the control parameter, in order to preserve Saltzman's idea that ice age cycles appear as the consequence of a slow perturbation of the carbon cycle. As in SM90, PP04 exhibits a limit cycle arising from a subcritical Hopf bifurcation. The dynamics along the limit cycle close to the bifurcation point are strongly influenced by the presence of the unstable focus. This is the configuration shown in figure 2*c*. Depending on *g*, the focus is either on the low-ice-volume side of the limit cycle (i.e. the system spends most of its time with high CO_{2}) or on the high-volume side of the limit cycle (i.e. the system spends most of its time in low CO_{2}). Parrenin and Paillard estimate that we are currently in the second configuration.

#### (v) A minimal model of ice ages

It has been claimed [61,63] that any model that has some form of 100 000 year internal periodicity could be used to reproduce the course of ice volume over the last 800 000 years. Taking the argument at face value, Crucifix [64] used one of the simplest possible slow–fast oscillators—the van der Pol (VDP) oscillator—with minimal modifications to account for the astronomical forcing and the asymmetry between the phase of ice build-up and melt during the Late Pleistocene (*α*, *β*, *γ* and *τ* are parameters; *F*(*t*) is the astronomical forcing),
VDP

The system dynamics are determined by the structure of the slow manifold *x*=*y*^{3}/3−*y*. The parameter *β* controls the position of the fixed point on the slow manifold and, consequently, the ratio of times spent by the system in the two branches (‘glacial’ and ‘interglacial’) of the slow manifold. The ice age curve can be captured with some tuning (figure 7), although it is fair to add that a small change in parameters may shift the timing of one or several ice age cycles. This minimal model was used to challenge intuitive arguments about the predictability of ice ages [64].

### (b) Models for millennial climate variability

#### (i) Dansgaard–Oeschger events as relaxation oscillations

Welander [65] introduced the concept of relaxation oscillations in the context of ocean dynamics. He described a *heat–salt oscillator* involving exchanges of heat and salt within a single oceanic column, coupled to a phenomenon of surface temperature relaxation. The destabilization process needed for the relaxation oscillation to appear is here related to diffusion between the deep ocean and the mixed layer. The system dynamics are further controlled by the mean freshwater flux at the top of the ocean column. It determines the transitions from a regime characterized by perpetual convection in the oceanic column, to a regime with intermittent convection (oscillation) and finally to a regime with no convection [66]. The bifurcations between the different regimes bear the character of *global* bifurcations, with the oscillation period approaching infinity near the bifurcation points (in particular, the second bifurcation bears the character of a homoclinic bifurcation). The heat–salt oscillator belongs thus to the class shown in figure 2*b*.

Welander [67] and Winton & Sarachik [68] later introduce the concept of another kind of relaxation oscillator in the ocean. It involves the meridional structure of the ocean thermohaline circulation, and the key nonlinear process is the meridional advection of heat and salt. The oscillations featured by this model are termed ‘deep-decoupling’ oscillations [68]. Given that the slow process now relates to heat accumulating in the global ocean, the characteristic time of deep-decoupling oscillations is of the order of 1000 years. The net flux of freshwater delivered to the North Atlantic acts as a bifurcation parameter controlling the transition between non-oscillating and oscillating regimes in the Winton–Sarachik model [69].

Millennial oscillations have since been observed across a hierarchy of ocean models, including three-dimensional ocean models with prescribed freshwater flux and restoring conditions to surface temperature [70], and three-dimensional models coupled to a simple atmosphere [71,72]. Sakai & Peltier [73] proposed that millennial deep-decoupling oscillations could explain Dansgaard–Oeschger events. Colin de Verdière *et al.* [33,74,75] complement this early proposal with a fairly complete theory based on ocean circulation model experiments. The oscillations described by Colin de Verdière *et al*. involve the processes of turbulent vertical mixing (neglected in Winton & Sarachik [68]), advection and convection, which unify the salt oscillator with the deep-decoupling oscillation model. Incidentally, Colin de Verdière *et al*. [74] dismiss the nonlinearity of the equation of state as the cause of the oscillations.

There is, across the model hierarchy, consistency about the fact that the transition between the oscillating circulation regime and the so-called diffusive, haline regime (without deep convection) is associated with a homoclinic bifurcation [33,76]. The nature of the bifurcation between the convective regime and the oscillation is more model-dependent. Timmermann *et al*. [76], based on experiments with the eight-ocean box Gildor–Tziperman model, find a Hopf bifurcation; salt-conserving experiments with a two-dimensional ocean model show a transition towards a finite-period cycle, but of increasing period as the bifurcation is approached; experiments with a more idealized model, formulated as a two-equation dynamical system, reveal the signature of an infinite-period bifurcation [33].^{3} The latter implies that Dansgaard–Oeschger events, at the time when they appear soon after the glacial inception process, should be very long but of a similar amplitude to the Dansgaard–Oeschger events coming later in the glacial cycle. This feature is consistent with the Greenland ice core record (figure 1). More specifically, the first Dansgaard–Oeschger cycles that appeared at the beginning of the glacial era were characterized by a long ‘plateau’ phase (also called interstadial) during which the thermohaline circulation was certainly very active [15]. In the Colin de Verdière *et al*. [74] theory, the plateau phase is the phase of the trajectory influenced by the ‘ghost of the saddle point’.

#### (ii) Dansgaard–Oeschger cycles as the manifestation of an excitable system

Given the explosive nature of the bifurcations involved in ocean dynamics it is no surprise to find excitability properties in ocean models. Weaver & Hughes [70] discuss this effect in salt-conserving experiments with an idealized-geometry, ocean model. The ocean–atmosphere model of intermediate complexity, CLIMBER (CLIMate BiosphERe model), was shown to exhibit excitability properties when boundary conditions are set to be typical of the latest glacial era [77]. The ocean circulation has then one stable state, with a moderate Atlantic overturning and a ‘quasi-stable state’ with more intense overturning. The conceptual sketch of the excitation cycles shown by Ganopolski & Rahmstorf [78], fig. 1 can be interpreted in terms of slow–fast dynamics, in which the different states of the ocean circulation constitute the different branches of a slow manifold. The intense overturning state, which is the ‘plateau’ phase of the Dansgaard–Oeschger event, may thus be viewed as the repelling branch of the slow manifold in the excitable regime (figure 2*d*). The excitable Dansgaard–Oeschger hypothesis was used as a possible basis to explain how a weak forcing, exogenous to the system, could explain the observed 1500 year periodicity of Dansgaard–Oeschger cycles (on this periodicity refer to earlier studies [79,80], but see the other view in the study of Ditlevsen & Ditlevsen [81]). Two such theories were developed on the basis of experiments with CLIMBER. One suggests that Dansgaard–Oeschger events are excited by stochastic fluctuations, modulated by a weak, hypothetical solar periodic forcing [82] (more on the effects of stochastic fluctuations in §4). The alternative theory suggests that the excitation is induced by the interference between two solar forcings with periods close to 1470/7(=210) and 1470/17(≈87) years [83], possibly combined with noise [84].

#### (iii) Heinrich cycles as a relaxation oscillation

MacAyeal [85] proposed an ice binge–purge theory to explain Heinrich events. The theory rests on experiments with a one-spatial direction model of ice flow dynamics. Suppose, as a starting point, that ice volume grows in response to net accumulation of snow. The growth continues until the accumulated effect of geothermal heat flux causes basal sliding. A volume of ice is then released into the ocean (this is the ‘purge’), causing the release of icebergs characteristic of Heinrich events. Ice volume thus decreases, until ice accumulation wins over so that ice volume can grow again. The ice binge–purge model is thus a relaxation oscillator combining a slow integrating process (ice mass accumulation) with a fast lateral discharge process.

#### (iv) Coupling between Heinrich and Dansgaard–Oeschger events

To what extent may Heinrich events interfere with Dansgaard–Oeschger dynamics? Paillard [55,86] investigated this question by coupling the MacAyeal ice model—but reduced to ordinary differential equations by Galerkin truncation—with a three-ocean box model. The coupling simply assumes that ice released into the ocean causes a net freshening of the surface of the North Atlantic that alters the deep-ocean circulation. Paillard realized that this coupling could lead to fairly non-intuitive and complex effects, such as the succession of Dansgaard–Oeschger events of decreasing amplitude between Heinrich events. This succession is known in the literature on palaeoclimate records as *Bond cycles* [18]. Paillard also found that the oscillations are aperiodic in this model under certain parameter configurations.

The issue is further explored in Schulz *et al.* [69], based on the Winton–Sarachik ocean model, and in Timmermann *et al*. [76], based on the slightly more sophisticated Gildor–Tziperman [59] ocean model. The objective was to study the response of deep-decoupling ocean oscillations to prescribed Heinrich cycles. Schulz *et al.* [69] noted that deep-decoupling oscillations could be synchronized on the Heinrich cycles. Timmermann *et al.* [76] then proposed, on the basis of numerical experiments with a fairly idealized model, that Heinrich events excite Dansgaard–Oeschger cycles because the variation in ice volume caused by a Heinrich event modifies slowly the amount of net freshwater released in the ocean. In turn, they suggested, Dansgaard–Oeschger may have a control on ice volume growth. This yields a two-way coupling between Dansgaard–Oeschger and Heinrich events.

Experiments with more comprehensive models of the ocean–ice sheet–atmosphere system [87,88] generally support the idea that the different water and heat fluxes involved in the different phases of ice build-up and iceberg release are quantitatively sufficient to support a coupling between ice sheets and ocean circulation during the latest glacial era. However, it was also noted that ‘three-dimensional thermomechanical ice-sheet models are unable to satisfactorily reproduce the binge–purge mechanism without an ad hoc basal parameterisation’ [89].

To address this difficulty, a theory in which Dansgaard–Oeschger events trigger Heinrich events was recently proposed [89]. The ice shelf plays a key role, in blocking the ice stream flow from the ice sheet to the oceans. Heinrich events occur when this ice shelf is broken; for example, under the influence of ocean sub-surface warming associated with a Dansgaard–Oeschger event. The resulting model is a system displaying a slow ice build-up—the Heinrich release cycle *excited* by fluctuations in ocean sub-surface temperature.

#### (v) Holocene oscillations and relationship with Dansgaard–Oeschger events

The much smaller ocean oscillations that characterized the Holocene period may also be a relaxation phenomenon. Schulz *et al.* [31] observed oscillations in the atmosphere–ocean model of intermediate complexity ECBILT-CLIO. These oscillations are related to the convective activity in the Labrador Sea.

Schulz *et al.* [90] considered the existence of such an oscillator in an earlier reference and speculated on the possible interactions between the centennial oscillations, millennial oscillations and Heinrich cycles. They considered a model in which each of these three kinds of oscillations is modelled as a Morris–Lecar relaxation oscillator. Their working hypothesis is that glacial conditions induce a coupling between these oscillators. They then observed that a very stable 1500 year oscillation appears, which they interpreted as a model equivalent of Dansgaard–Oeschger events.

## 4. Stochastic effects

The myriad of chaotic motions that characterize the dynamics of the ocean and the atmosphere may be taken into account in the form of parametrizations involving stochastic time processes. The method was introduced in climatology in the 1970s [91] and the theoretical justifications, which allow one to model chaotic motions as a (linear) stochastic process, are reviewed by Penland [92,93]. In a statistical inferential framework, the stochastic parametrizations may also be viewed as a way to account for the distance necessarily existing between the concepts and dynamics featured by the model, and the complex system being observed.

The effects of stochastic processes on relaxation oscillators and excitable systems are generally well documented in the literature because this is a topic of general interest [94]. Here, we review some of them in the specific context of palaeoclimate dynamics.

### (a) Stochastic effects on ice age dynamics

#### (i) Phase dispersion

One of the basic effects of noise on oscillators is the phenomenon of *phase dispersion*: a weak stochastic forcing on an oscillator causes a fading out of the memory of the exact initial conditions, even though the gross structure of the oscillation visualized in the phase space is conserved. The phenomenon is well known and it is an immediate consequence of the neutral stability of the phase of a free oscillator with respect to fluctuations. It was previously suggested that this phenomenon of phase dispersion may concern ice ages [95], but it is more commonly believed that ice ages are phase locked on the astronomical forcing. This phase locking should act against dispersion and permit a very long predictability horizon of ice ages. However, a phenomenon of phase dispersion may happen in oscillators that are locked on a periodic forcing. A stochastic fluctuation may momentarily cause a burst of desynchronization, called phase slip, during which the system is unhooked from its corresponding deterministic trajectory and attracted to another trajectory, which leads or lags the original one by one forcing period [41], §3.1.3. The difference between phase diffusion in a free system and in a periodic forcing-driven oscillator is that the diffusion effect has, in the latter, a quantum nature. More formally, it is said that the stochastic forcing disperses the system states around the different attractors that are compatible with the forcing. In a work in preparation, we suggest that the astronomically forced climate system may satisfy the conditions for a similar phenomenon of phase dispersion to occur (B. De Saedeleer, M. Crucifix & S. Wieczorek 2011, unpublished data). Given that the astronomical forcing is aperiodic the description of the phenomenon requires a suitable theoretical framework, which relies on the notion of a ‘local *pullback* attractor’. The equivalent of a phase slip is, in the aperiodic forcing context, a stochastic shift from one of the deterministic pullback attractors to another one. The phenomenon is illustrated based on experiments with the VDP model in figure 8.

#### (ii) Reduction of period

Additive fluctuations generally reduce the period of relaxation oscillators. In an oscillator presenting a homoclinic orbit such as the Duffing oscillator, additive fluctuations reduce the time spent near the unstable focus [96]. This implies that, even at the corresponding bifurcation point in the deterministic system, the return time of oscillations in the stochastically perturbed system remains finite. In a slow–fast oscillator such as the VDP oscillator, additive fluctuations generally result in early escapes of the branch of the slow manifold on which the system lies (figure 9). The period of the oscillator is thus affected by a correction that increases approximately linearly with the noise variance in the slow–fast VDP oscillator [97].^{4} This property was used at least once in Pleistocene theory, in the silicate weathering hypothesis advanced by Toggweiler [99]. Additive fluctuations reduce the limit cycle period from 800 000 to about 100 000 years. The reasons for the period reduction being so dramatic are left for another study.

### (b) Stochastic effects on Dansgaard–Oeschger dynamics

#### (i) Stochastic excitation and resonance

Noise may naturally act as an excitation agent in an excitable system (figure 9). The topic is extensively reviewed by Lindner *et al*. [94]. Excitation loops are sporadic if the noise amplitude is weak, in which case the *recurrence time* of the events is set by the noise amplitude, whereas the *amplitude of the events* is set by the structure of the deterministic vector field. The frequency of excitation loops increases with the noise amplitude, until the system behaviour is qualitatively similar to a limit cycle regime. This is the coherent resonance regime. For yet higher noise amplitudes, the limit cycle structure is destroyed.

The concept of stochastic excitation has been considered several times in Dansgaard–Oeschger theories. The idea is introduced based on experiments with an ocean general circulation model with idealized geometry and forcing [70]. The effect of stochastic fluctuations is not only to excite oscillations in a system normally at rest, but also to reduce the period of these oscillations when the system is in oscillatory regime.

A phenomenon of non-autonomous stochastic resonance may occur if the noise is superimposed on a weak external drive. For this to happen, the autonomous system needs to be stable but excitable. The external forcing must be too weak to cause excitation by itself. The role of the noise is to provide the additional power to induce excitation. The timing of the excitation is then related to the phase of the external forcing. The mechanism was proposed several times [76,82] to explain the 1500 year recurrence time of Dansgaard–Oeschger events. The idea remains questioned, either on the grounds that the 1500 year recurrence time observed in palaeoclimate records is coincidental [81] or on the grounds that the 1500 year external forcing is unidentified [33,83]. A more subtle case of stochastic resonance involves the combination of noise with two solar cycles of 210 and 87 years, which yields the concept of ‘ghost resonance’ [100] for which some support, albeit not conclusive, is found in the observations [80].

#### (ii) Decreased sensitivity to noise in resonant oscillators

Coupled oscillators may exhibit, collectively, a resonance period that is more robust to external fluctuations than the uncoupled oscillators. Schulz *et al*. [90] used this property to explain the stability of the Dansgaard–Oeschger recurrence period of 1500 years in the presence of random fluctuations, without having to invoke an external forcing.

#### (iii) Pseudo-oscillations in two-well systems

Finally, a behaviour reminiscent of oscillations may occur in a system that is neither oscillating nor excitable, but which presents several stable states. Noise then simply induces jumps between these different states. The simplest mathematical model is the Langevin equation and this is on this basis that Schulz *et al.* [31] interpret the Holocene oscillations observed in the ECBILT-CLIO climate model.

## 5. Concluding discussion: can dynamical systems be used for inference?

The review has shown that relaxation oscillations are a popular and powerful model to explain oscillations observed in the Pleistocene record. The concept of *relaxation* implies some form of slow–fast separation, in the sense that at least one component of the system spends most of its time in ‘quasi-equilibrium’ states (this may be a ‘slow manifold branch’ or a region influenced by a saddle point, depending on the system structure), with acceleration phases.

Some of these models were constructed following a fairly careful procedure of truncation of a system of partial differential equations, which describes some of the fluid dynamics of the climate system. Others were proposed on a more conceptual basis, the idea being precisely to test a hypothesis based on palaeoclimate observations. The latter approach is sometimes criticized, on the grounds that box models, for example, cannot reasonably be taken as an adequate representation of the complex dynamics of the oceans [101].

This leads us to the last question of this review: can dynamical systems be used for inference on palaeoclimates? Inference implies that something is being learned by confronting a model with observations. This inference process may take the form of a calibration procedure (update our knowledge on parameters on the basis of observations) or a model selection procedure (which model, among different alternatives, explains the observations best).

The position taken here is that there is not such a thing as an ‘attractor’ of the climate system that is to be ‘discovered’. The hope is that some of its modes of behaviour are sufficiently decoupled from the rest of the variability to justify the fact that simple dynamical systems may capture the fundamental dynamical properties of these modes, and we want to learn about these modes from palaeoclimate observations.

The programme is challenging. Indeed, it was emphasized that different physical assumptions may lead to dynamical systems with dynamical properties that are similar enough to produce a convincing visual fit on palaeoclimate data [61]. The message is largely echoed in this review. The modeller's challenge is, therefore, to operate a model selection on more stringent criteria than just fitting some standard time series. For example, palaeoclimate observations may yield constraints on the bifurcation structure of the system. The Middle Pleistocene Transition is an attractive test case in this respect.

In a statistical inference process, the observations should be a *plausible* outcome or realization of the model. This makes sense only if the model has a stochastic component, which describes its uncertainties, limitations, and the noise that emerges from the chaotic motions of the atmosphere and oceans.

Stochastic dynamical systems are beginning to be used for inference on palaeoclimate time series. In a method called ‘potential analysis’, the climate system is modelled as a Langevin equation, that is, the combination of a down-gradient drift with a Wiener additive process and inference is made on the number of wells of the potential function [102]. The method was applied on Pleistocene climate records, yielding the conclusion that the number of wells increased from two to three over the course of the Pleistocene [103].

However, our position so far has been to favour a Bayesian methodology, because it allows one to encode physical constraints in the form of prior distributions on model parameters. The Bayesian formalism is also naturally designed for model calibration, selection and probabilistic predictions.

The fact is that Bayesian methods for selection and calibration of dynamical systems on noisy observations are only emerging. In a recent attempt, we considered a particle filter for parameter and state estimation [104]. To be honest, there is ample room for progress. Whether the process of inference with simple dynamical systems on palaeoclimate data will lead new insight in this context still needs to be demonstrated.

## Acknowledgements

This review benefited from stimulating discussions during the Isaac Newton Programme *Mathematical and Statistical Approaches to Climate Modelling and Prediction* held in Cambridge during the summer–autumn 2010. Guillaume Lenoir, Bernard De Saedeleer and Jonathan Rougier provided helpful comments on a first version of the article. Thanks are also due to Didier Paillard, Olivier Arzel, Alain Colin de Verdière, André Paul and Sebastian Wieczorek for e-mail correspondence during this review. The author is Research Associate with the Belgian National Fund of Scientific Research. This research is supported by the FP7-ERC starting grant ITOP ERG-StG-2009-239604. The editor (Jan Sieber) and two reviewers are acknowledged.

## Footnotes

One contribution of 13 to a Theme Issue ‘Climate predictions: the influence of nonlinearity and randomness’.

↵1 The astronomical forcing will generally be taken into account here in the form of a normalized measure of insolation during the month or on the day of summer solstice at a northerly latitude, typically 60 or 65

^{°}N. This is a fairly complex, aperiodic signal, with dominant harmonics corresponding to the phenomena of precession (23 716, 22 428 and 18 976 years), and obliquity (41 000 years) [5].↵2 In fact, the view is introduced by Murphy [48] but it is developed mathematically in the ‘canon of insolation’ authored by Milankovitch [4].

↵3 This particular case is not illustrated in figure 2. There is no saddle point along or near the orbit, but there is a combination of parameters for which a fixed point appears on the limit cycle. Beyond this particular parameter value, this fixed point splits into a saddle and a node. This particular parameter value corresponds to an ‘infinite-period’ bifurcation. In practice, as long as the limit cycle exists, the trajectory slows down near the point where the saddle–node will appear. Some authors then refer to the influence of the ‘ghost’ of the saddle point [38].

↵4 This result is established assuming fluctuations added to the slow variable. A much more general theory, suitable for different slow manifolds and additive fluctuations to slow and fast variables, is now available [98].

- This journal is © 2012 The Royal Society

This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.