How can a glacial inception be predicted?
Contribution to the special issue of the Holocene on the Early Anthropogenic Hypothesis
note : minor modifications may have been introduced at the proofediting stage
Abstract
The Early Anthropogenic Hypothesis considers that greenhouse gas concentrations should have declined during the Holocene in absence of humankind activity, leading to glacial inception around the present. It partly relies on the fact that present levels of northern summer incoming solar radiation are close to those that, in the past, preceded a glacial inception phenomenon, associated to declines in greenhouse gas concentrations. However, experiments with various numerical models of glacial cycles show that next glacial inception may still be delayed by several ten thousands of years, even with the assumption of greenhouse gas concentration declines during the Holocene. Furthermore, as we show here, conceptual models designed to capture the gross dynamics of the climate system as a whole suggest also that small disturbances may sometimes cause substantial delays in glacial events, causing a fair level of unpredictability on ice age dynamics. This suggests the need of a validated mathematical description of the climate system dynamics that allows us to quantify uncertainties on predictions. Here, it is proposed to organise our knowledge about the physics and dynamics of glacial cycles through a Bayesian inference network. Constraints on the physics and dynamics of climate can be encapsulated into a stochastic dynamical system. These constraints include, in particular, estimates of the sensitivity of the components of climate to external forcings, inferred from plans of experiments with large simulators of the atmosphere, oceans and ice sheets. On the other hand, palaeoclimate observations are accounted for through a process of parameter calibration. We discuss promises and challenges raised by this programme.
1 Introduction
The Early Anthropogenic Hypothesis considers that the natural evolution of climate consisted in a decline in greenhouse gas concentrations throughout the Holocene, leading today to conditions favourable to accumulation of ice in the Northern Hemisphere (Ruddiman et al. (2011) and references therein). This hypothesis supposes an important premise: it is possible to predict the slow evolution of climate several millennia ahead. Indeed, suppose that climatologist lived 8,000 years ago. What the Early Anthropogenic Hypothesis says is that a forecast for the 8,000 years to come made by this early climatologists would have been a decline in greenhouse gas concentrations and ultimately, glacial inception.
Throughout this paper we will consider that this premise should not be taken for granted. The general problem of predicting the evolution of the climate system several millennia ahead may be tackled from different perspectives. One method of prediction consists in seeking episodes in the past during which the climatic conditions and forcings were analogous to those of the Early Holocene, and observe the subsequent climate evolution. The other method consists in using models, generally numerical models, which account for known physical constraints on the dynamics of the climate system. The two methods are briefly reviewed and discussed in the next section of this article.
Most investigators generally appreciate that the two methods need, in practice, to be combined. The question, addressed in section 3, is how observations and models may be combined in a way that is transparent, optimal, and avoids the dangers of a circular reasoning. The methodology proposed here consists in (a) formulating stochastic dynamical systems capturing palaeoclimate dynamics. These may be viewed as generators of palaeoclimatic histories, designed to encapsulate information inferred from experiments with large numerical models and/or other hypotheses about climate system dynamics; (b) calibrate the parameters and initial conditions of these dynamical systems on palaeoclimate observations. We explain why the Bayesian statistical framework is adapted to this programme.
2 Traditional approaches to predict the slow evolution of climate
2.1 The analogues method and its limitations
Among the solutions envisioned by Edwar Lorenz (1969) to forecast weather, was an empirical approach also known the ‘analogues method’. Considering that similar causes should lead to similar effects at least within a finite time horizon, weather predictions may be attempted by seeking analogous synoptic situations to today’s in the archives of meteorological services. Likewise, predicting the natural climate evolution during the Holocene may be addressed by identifying in palaeoclimate archives situations for which the state of climate and the forcing conditions were similar to those experienced during the Holocene.
At these times scales climate is influenced by the slow changes in incoming solar radiation, and the latter is a function of the Earth orbital elements (characterised by the eccentricity and the true solar longitude of the perihelion ) and the tilt—or obliquity—of the equator on the ecliptic, denoted . These three elements have all specific and wellknown effects on the seasonal and spatial distributions of incoming solar radiation (insolation). They also have their own periodicities: You never encounter twice a same orbital configuration (see, on this topic, the early works of Berger (1977) and Berger (1979)).
Palaeoclimatologists therefore tend to consider specific measures of insolation, supposed to summarise in one quantity the effect of insolation on climate, and look at past times where that specific measure of insolation was the same as during the Holocene. However, the choice of an ‘analog’ depends on how these measures are constructed (Loutre et al., 2008).
To illustrate this point consider three measures of insolation classically referred to in the tradition of Milankovitch’s theory (Figure 1): daily mean insolation at the summer solstice (21th June), calculated according to Berger (1978a); (b) insolation integrated over the caloric summer (the 180 days of highest insolation, assuming a 360day year) (Berger, 1978b)—this measure is favoured by Ruddiman (2007)— (c) insolation integrated over all days for which insolation exceeds 300 W/m (after a proposal by Huybers and Tziperman (2008)). All are calculated here for a latitude of 65 N. It is easily shown that these measures of insolation are very well approached by a linear combination of climatic precession () and obliquity (); they are here ordered by increasing influence of obliquity vs climatic precession.
Circles on Figure 1 indicate times at which the different insolation measures equal their presentday values. Call these ‘insolation matches’. Those occurring during a climate environment analogous to the present Holocene (arbitrarily defined as a CO concentration larger than 250 ppm and benthic less than 3.8 per mil) are highlighted by vertical bars.
The timing of insolation matches differ by about 1 ka (thousands of years) depending on which insolation is considered: for example the time of insolation match at marine isotopic stage (MIS) 11 ranges between ka (daily mean at solstice) and ka (Huybers’ thresholded sum). In general, though, insolation levels similar to the present day slightly precede or coincide with times during which ice volume markedly increased. Bearing in mind the caveats introduced by uncertainties on the chronology of the palaeoclimate records, this observation suggests that today’s insolation levels are close to those that prevailed at the end of previous interglacial periods. However, observe also that presentday insolations are at a minimum on their secular course (daily mean solstice) or close to it (seasonintegrated insolations). This remark invites us to be prudent about the prediction of a natural end to the Holocene around the present, again assuming no anthropogenic perturbation.
More crucial differences among insolation measures lie in their dynamics. Climate does not respond instantaneously to the astronomical forcing. Ice buildup and ocean alkalinity adjustments involve time scales of the order of 10,000 years at least. Predicting the evolution of climate therefore requires one to consider the history of astronomical forcing over long enough a timewindow, corresponding to the time needed by the different climate system components to adjust dynamically to that specific astronomical history. Ruddiman (2007) and Loutre et al. (2008) noted that the astronomical history preceding ka (marine isotopic stage 9) admittedly looks as a reasonable candidate to the last 20 ka by reference to seasonintegrated insolations (consider the ’shapes’ of the caloric summer insolation before ka); though, it is not straightforward to demonstrate that this really is the ‘best possible analog’ quantitatively, in a way that is convincingly robust and objective (Loutre, pers. comm. and further tests made for this article).
There are further caveats to a straightforward use of insolation analogues to predict climate change. First, which insolation (daily mean solstice or seasonal integration) is the best predictor of climate change may depend on the level of glaciation (Rohling et al., 2010). Second, the analogues method supposes as a premise that one single possible climate history may correspond to a given history of astronomical forcing. In the formal language of dynamical system theory this property is called general synchronisation of climate on astronomical forcing (see Rulkov et al. (1995) for a formal definition). It is shown in the next section that this might not be the case. Third, climate cycles evidently changed in shape and variance over the last million years due to slow changes in the environment. For example, the midPleistocene revolution, when glacial cycles shifted from a period of 40 ka to a period of the order of 100 ka is explained by Saltzman and Maasch (1990) as a bifurcation related to a slow, tectonicallydriven decline in background levels of CO. Clark et al. (2006) review a number of alternative explanations. The midBrünhes transition about 400 ka ago, after which interglacial periods became characterised by higher greenhouse gas concentrations than before this transition, is still unexplained (see the review of Tzedakis et al. (2009)). Consequently, one cannot be assured that the levels of insolation leading to glacial inception determined on the basis of an event 400 ka ago still hold valid today.
2.2 Using climate simulators to predict the glacial inception
The other approach to the problem of predicting an ice age proposed by Edwar Lorenz relies on a mathematical model of the Earth’s climate. In a visionary paper entitled “climate modelling as a mathematical problem” (Lorenz, 1970), he expresses the hope that “in the 21st century”, large mathematical models accounting for “every feature of the atmosphere and its environment” as well as ocean circulation and vegetation could be run and predict ice ages.
The ‘mathematical’ models envisioned by Lorenz are known today as ‘general circulation models’ or ‘global climate models’. They reproduce many features of the climate system visible at the global scale (the ocean circulation, modes of variability such as ElNiño) based on equations that are expressed at the level of a grid cell. We prefer to call them ’simulators’ in order to avoid confusion with other meanings of the word ‘model’.
Accounting for physical constraints is naturally considered to strengthen the reliability of a prediction, especially if this prediction involves a situation that is unprecedented. Lorenz’s proposal is therefore a generally well accepted way of dealing with the problems of climate predictability at all time scales. However, it has a number of drawbacks that need to be emphasised in the present context. An oftenheard objection to Lorenz’ proposal is computing time. Stateoftheart simulators resolving the dynamics of the ocean and the atmosphere are designed to be run on time scales of the order of a century. One experiment of thousand years with such a simulator may take months of computing time when it is possible. Experiments over several millennia therefore require changes in the simulator design. Modellers have generally chosen to degrade, sometimes very drastically, the resolution of the atmosphere and ocean and simplify or suppress many of the parametrisations in order to gain computing time. Cloud cover, for example, may simply be considered as constant. In the meantime, they attempted to account for the dynamics of other components of the climate system, such as the ice sheets and the carbon cycle, which react on the longer time scales. The simulators satisfying these criteria are known as Earth Models of Intermediate Complexity, or simply EMICS (Claussen et al., 2002).
EMIC experiments designed to test the early anthropogenic hypothesis were presented soon after the original publication of Ruddiman (2003). An EMIC called CLIMBERSICOPOLIS ((Petoukhov et al., 2000; Calov and Ganopolski, 2005)) was used to estimate the effects of declining CO concentrations during the Holocene from about 265 ppm to about 235 ppm (a similar decline was imposed on methane concentrations). With that scenario, CLIMBER does not accumulate ice any time during the Holocene, while it correctly reproduces the increase in ice area during the Eemian / Weschelian transition (115 ka ago) Claussen et al. (2005). Similar experiments were carried out with the LLN2D simulator of Gallée et al. (1992). Consistent with the experiments with CLIMBER, reasonable noantropogenic scenarios with CO declining to about 240 cause almost no accumulation before several tens of thousand years in the future. Only very drastic declines in CO concentrations, such as to reach a CO concentration of 225 ppm at present and 206 ppm 4 ka in the future, yield some ice accumulation (about 27 m of equivalent eustatic sealevel in 22 ka).
Results of experiments with EMICS may however be questioned because they misrepresent (if at all) features of the climate system that may be important to understand and predict its slow evolution, such as monsoons, western boundary currents, modes of variability like ENSO and the details of the ocean circulation in the Southern Ocean. For these reasons a number of experiments were published, based on simulators with more stateoftheart representations of ocean and atmosphere dynamics. These simulators are not designed to compute the accumulation of ice but they are useful to indicate spots where snow is likely to accumulate from one year to the next, and/or to analyse the conditions propitious for accumulation. They also provide information about the sensitivity of the ocean surface and circulation to changes in the astronomical forcing and greenhouse gas concentrations.
In the present special issue, three articles discuss experiments with the Community Climate Model in different configurations: Vavrus et al. (2011) examine the influence of topography representation on the snow accumulation process and find that accounting for higher resolution topography increases the sensitivity of snow accumulation on the external forcing; Kutzbach (2011) examine the influence of decreasing greenhouse gas concentrations on seasurface temperatures and attempts to quantify the effects on the carbon cycle; they suggest that the feedback of the surface cooling on the carbon cycle is substantial enough to accommodate Ruddiman’s suggestion of a natural amplification of the natural perturbation Ruddiman (2007). Vettoretti and Peltier (2011) uses the oceanatmosphere version of the Community Climate Model and compares the effects of decreasing CO concentrations with those of the orbital forcing on snow accumulation and the abyssal circulation in the Atlantic. They come to a somewhat challenging conclusion for the Early Anthropogenic Hypothesis, that is, astronomical forcing is a more important driver of ice accumulation than CO. Earlier, Vettoretti and Peltier (2004) presented a series of eight experiments with the Canadian Climate Centre Atmospheric Model GCM2 to identify the conditions propitious to yeartoyear snow accumulation on ice nucleation sites (see also references to closely related work in that article). They found that obliquity greatly influences snow accumulation in their model.
Here, I would like to point to perhaps a more fundamental challenge to Lorenz’s vision. Even the most sophisticated simulators are an incomplete (truncated) representation of reality. Some of the truncated processes are represented with semiempirical formulae, called parameterisations (see Palmer (2005) for a very accessible overview). A fraction of this uncertainty—called the parametric uncertainty—can be quantified by considering the sensitivity of predictions on parameters involved in the representation of the truncated processes. Pioneering works of Murphy et al. (2004) on this subject lead to great research activity. However, it will never be possible to guarantee that all the physical and biogeochemical processes are correctly and accurately represented in a numerical model of the climate system. There is structural uncertainty. In particular, remember that CO is active not only as a greenhouse gas, but also as a fertiliser and an ocean acidifier; changes in CO concentrations at glacialinterglacial time scales have innumerable consequences on life, and hence on climate.
Saltzman already noted Saltzman (2001)(sect. 5.1) that expecting from a climate simulator to reproduce ice ages without any prior information about the timing of glacialinterglacial cycles is illusory, so huge would the accuracy needed on the snow accumulation imbalance and carbon dioxide fluxes to capture the rate and timing of glacial cycles. Given that glacial cycles cannot be simulated ab initio, past observations have to be somehow incorporated into the simulator. Very often, the procedure is quite implicit: The uncertain parameters of the climate simulator are varied until the simulated climate (or climate changes) is reasonably realistic. This is called tuning. One must be clear about the fact that published predictions of the next glacial inception are all obtained with climate simulators ’tuned’ on past climate observations. For example, the LLN2D simulator, used in Berger and Loutre (2002), Loutre and Berger (2003) and Loutre (2003) was tuned to capture the timing and rate of the latest glacial inception, about 115 ka ago (André Berger, pers. comm.). Confidence in the simulator was gained by the fact that once tuned, it reproduces the last glacialinterglacial cycle (Gallée et al., 1992) reasonably well, as well as the dynamics of previous interglacial periods (Loutre and Berger, 2003) and, admittedly with more difficulties, the succession of marine isotopic stages 121110 (400,000 years ago) (Loutre, 2003).
The problem is that tuning a simulator is a fairly nontransparent and nonoptimal way of using past climate information for the prediction. It is nontransparent in the sense that it is not always clear which information has been incorporated and which one was effectively predicted (or ’hindcasted’); it is nonoptimal in the sense that only a fraction of the palaeoclimate information has been used in the inference process, and that the simulator can never be perfectly calibrated. Furthermore, the ‘satisfaction’ criteria which preside the tuning process are not explicitly defined and the consequences of the deviations between simulator results and observations on predictions are not explicitly quantified.
3 Dynamical systems as models of glacial cycles
3.1 What do we learn from simple dynamical systems?
The complexity and multitude of mechanisms involved in glacialinterglacial cycles may leave the modeller hopeless. Fortunately, this is a fact of nature that very complex systems may exhibit fairly regular and structured dynamics. This phenomenon, called ’emergence’, is related to the fact there are constrains on the system taken as a whole. The dynamics of the system may therefore be inferred and understood by identifying the constraints that are relevant at the scale of interest without having to consider all the interactions between the individual components of the system.
This argument justifies the interest of scientists for ‘conceptual’ models of glacial cycles. These models do have educational interest and they also have predictive power if they are correctly formulated (see examples in Paillard (2001); for more general background about complex systems consider the monographs of Haken (2006) and Nicolis and Nicolis (2007)).
For illustration consider the following minimal conceptual setting. Suppose that the state of the climate system may be summarised by two variables. Call them and . Variable represents the total continental ice volume. It accumulates variations in time due to the astronomical forcing plus a drift term associated with the variable , which is defined later. Mathematically this reads, if is a variation in over a fixed time interval :
(1a)  
Parameter is introduced later. is a characteristic time that controls the rate at which ice volume grows or decays, and which is tuned on observations. Variable represents the state of a climatic component that responds with a hysteresis behaviour to changes in ice volume. The variations in surface pCO associated with ocean circulation changes may present such a behaviour (e.g. (Gildor and Tziperman, 2001; Paillard and Parrenin, 2004)). This is mathematically modelled as follows:  
(1b) 
where is another parameter that controls the rapidity of the transition between the two branches of the hysteresis (to understand the principle of this second equation, consider the situation . There are then two possible values of that would be in equilibrium with (i.e.: ): or . Whichever is approached depends on the system history, i.e.: there is hysteresis). Equations (1a,b) constitute the dynamical core to several conceptual models of ice ages available in the literature, including models by Saltzman and Maasch (1991); Gildor and Tziperman (2001) and Paillard and Parrenin (2004). Mathematicians will recognize a particular formulation of the celebrated FitzHughNagumo model (e.g. Izhikevich and FitzHugh (2006)) originally developed to study neuronal responses to electrical stimuli (the term ’+1’ added at the end of equation (1b) plays no other role than constraining to be almost always positive, which is more physical for an ice volume). In the absence of astronomical forcing the behaviour of is mainly dictated by . It may converge to a fixed point close to zero (if ), much above zero () or oscillate between the two (if ). Within the oscillation regime controls the asymmetry of the spontaneous oscillations and , the sharpness of the shifts between the glaciation and deglaciation regimes. Here we chose to capture the notion of slow glacial buildup and fast deglaciation, and for rapid shifts in compared to variations in . The oscillations are also known as relaxation oscillations, first documented by Van der Pol in 1921 (Kanamaru, 2007).
For forcing the system we use a linear combination of climatic precession and obliquity in order to capture the notion that summer insolation controls ice ages (], where is the departure to the mean). With this scaling the forcing amplitude is approximately of the order of and thus of the same order as . There are already two lessons to be drawn from this simple experiment. Simulations with this model are shown on Figure 2. First, observe that the simulated future climate history is an ice volume decline up to 20 ka into the future, in spite of the fact that the model correctly reproduces the increases in benthic trends at the end of marine isotopic stages 11 (400 ka ago), 9 (325 ka ago, although the model spuriously produces a subsequent drop in ice volume), 7 (200 ka ago) and 5e (115 ka ago, but as for stage 9 the subsequent drop towards stage 5c is too strong). Hence, the argument that present ice volume should be growing (assuming no anthropogenic perturbation) because it did so during analogous situations in past interglacials is not sufficient. It is necessary to demonstrate, in addition, that this simple model is inappropriate. Second, observe that shifts to negative values at the end of the deglaciations and, in this model, this shift preconditions glacial inception. This may justify why, rather than considering insolation, it has sometimes been decided to seek an analogue to the present by considering the timing of terminations, because they are a proxy for the state of system variables that cannot be immediately observed. Namely, (EPICA community members, 2004) chose to align the previous deglaciation on the deglaciation that lead to stage 11. Based on the observed trends in the D, they estimated that glacial inception should not occur before more than ten thousand years from now. This alignment was debated on the ground that it violates the insolation alignment (Crucifix and Berger, 2006; Ruddiman, 2007) (further informal debate and correspondence followed). However, they may be a more fundamental reason why any ‘alignment’, i.e., any blind application of the analogues method, may yield a wrong or at least overconfident prediction.
To see this slightly modify (1b) to add a random number of standard deviation 0.05 every 1 ka (assume normally distributed) :
(1b’) 
The random (stochastic) process parameterises unpredictable oceanatmosphere disturbances such as DansgaardOeschger events (cf. Ditlevsen and Ditlevsen (2009)), or volcanic eruptions, which impact on climate dynamics. The theory that allows us to represent nonresolved dynamical processes with stochastic processes was introduced in climatology by Hasselmann (1976). Nicolis (2005) discuss to what extent such parameterisations may also be used to estimate the effect of accumulating model errors, and a full volume entirely devoted to stochastic parameterisations was recently published (Palmer and Williams, 2010).
Together, equations (1a) and (1b’) form a stochastic dynamical system. Two trajectories simulated with two different realisations of the random numbers are shown on (Figure 3). In principle the standard deviation of the noise (0.05) is an uncertain quantity that should be inferred from realistic computer experiments and / or calibrated on observations. Given that the present purpose is merely illustrative we are satisfied with the fact that the noise added to the evolution of and looks reasonable compared to typical benthic and ice core records. Indeed, it is seen that the stochastic disturbances are here too small to compromise the shape of glacial cycles and the action of astronomical forcing as a pacemaker, but they are large enough to alter the timing of glacial events and shift the succession of glacial cycles. Observe that the climatic history shown in black is similar to the deterministic simulation shown in Figure 2. The red one undergoes a shift around ka. The length of the simulated interglacials around ka are differ in the two experiments, in spite of the fact that the preceding deglaciations are in phase and of similar amplitude.
A similar phenomenon was is fact visible in the experiments with the conceptual model of Paillard (1998), who showed that two future climate scenarios, one with immediate glacial inception and one with a longdelayed glacial inception, with only tiny changes in model parameters. It is possible to explain these shifts with arguments of dynamical system theory (article in preparation). They are related to the fact that the nonperturbed dynamical system may accommodate different plausible climate histories, characterised by different series of deglaciation timings. At certain times, random perturbations (or small parameter variations) may cause a shift from one possible history to another one. The (unsolved) mathematical problem is to identify and characterise the moments at which such phase shifts are the most likely to occur. The possibility of such shifts contrast with the earlier conclusions of Tziperman et al. (2006). Whether phase shifts effectively occurred in the true climatic history is still to be elucidated. However, if this phenomenon is confirmed, it implies that the predictive horizon of glacial cycles is intrinsically limited and that the analogues method may lead us to be overconfident about what can actually be predicted.
3.2 Developing dynamical systems of glacial cycles.
The above example suggests that the gross dynamics of glacialinterglacial cycles may be characterised by a simple dynamical system. However, Tziperman et al. (2006) previously warned us about a form of ambiguity related with the interpretation of such simple dynamical systems of ice ages. Indeed, we have not been specific at all about the meaning of the variable . In the model of Paillard and Parrenin (2004), the role of is played by a variable quantifying the overturning of the southern ocean circulation. Saltzman and Maasch (1990) choose to introduce a third variable in order to distinguish the ocean circulation from the carbon cycle. In their model, the nonlinear terms appear in the carbondioxide equations. More interpretations of the meaning of and are possible. Ambiguity per se is not necessarily a bad thing. It is useful to explore and understand the system dynamics without having to enquire about the details of the mechanisms involved. However, it is desirable to attach a more concrete meaning to climatic variables in order to explore mechanisms and provide credible predictions outside the range of observations. There are two complementary ways of doing this. The first one is to verify that the relationships between the different variables are compatible with the information inferred from experiments with numerical simulators. For example, eq. (1a) encodes a linear relationship between the ice mass balance and insolation. How does it compare with climate simulators? The second one is to connect the variables with actual palaeoclimate observations. If and represent ice volume and, say, the North Atlantic Overturning cell, their simulated variations should be consistent with the current interpretations of planctonic and benthic foraminifera records.
In this section we focus on how we include information inferred from numerical simulators in simple dynamical systems (Figure 4 supports the exposé). Our ambition is to reach a level of understanding of the simulator response that is more general than mere predictions of the future ice volume evolution obtained with this simulator. To this end we propose to examine the bifurcation structure of the simulator in response to varying input parameters.
To this end, consider again LLN2D. This model is an implementation of equations of quasigeostrophic motion in the atmosphere and rheological equations for icesheet dynamics assembled with the purpose of studying glacial cycles. As we noted above the level of sophistication of LLN2D is largely superseded by modern simulators. Nevertheless it is a valid example because LLN2D was used in the debate over the early anthropogenic (Crucifix et al., 2005; Loutre et al., 2007) .
In this example, we attempt to determine the number and stability of steadystates obtained with LLN2D, as a function of the precession parameter. Mathematicians devised very ingenious and sophisticated ways of exploring the dynamical structure of large mathematical systems but here we follow a very simple and intuitive method, previously used to analyse simulators of the oceanatmosphere Rahmstorf et al. (2005) and oceanatmosphereicesheet systems Calov and Ganopolski (2005). It is known as the ‘hysteresis experiment’. This is achieved as follows: consider eccentricity of 0.025, obliquity of , a CO concentration of 210 ppm, and run the simulator until equilibrium with (perihelion on the 21st June). Then, slowly vary along the trigonometrical circle (the experiment takes the model equivalent of 400,000 years) and record the total continental ice volume as varies. The resulting curve (Figure 5) depicts an ensemble of quasiequilibrium points. The first observation is that there is at least one stable equilibrium response for all the forcings in the range explored. Secondly, there are two stable states for within about and . A similar curve was calculated for a CO concentration of 260 ppm. It is similar in shape but the twostablestate regime now ranges from to and the maximum ice volume is ). Similar results were previously obtained with the CLIMBERSICOPOLIS simulator (Calov and Ganopolski, 2005).
The bifurcation structure that emerges from LLN2D was not apparent until we did the simulations. We have therefore learned a new bit of information that we are ready to trust, or at least to test against the observations. This information may be encoded as a simple equation, which is known to have a similar bifurcation structure as the one found in LLN2D:
(2) 
where is a variation in ice volume over a timestep , is a function of ice volume, astronomical forcing and CO and other parameters (summarised by ) (the full expression is given in Appendix 1).
Equation 2 is a surrogate of the LLN2D model: it is a simple, fast process that behaves dynamically in a way that is similar to LLN2D. Here, we calibrated the parameter such that the surrogate reproduces approximately the behaviour of LLN2D over the last 800 ka. The calibration did not yield a simple best value of , but rather a distribution of plausible values of the parameter that reproduce reasonably well the model output. The resulting fit is shown on Figure 6.
By finding a surrogate for LLN2D we have progressed on two fronts. First, we have distilled the vast network of partial differential equations implemented in LLN2D into a much simpler structure. This structure can now be compared with canonical dynamical systems in order to comprehend its important dynamical characteristics (Ditlevsen (2009) precisely studied an equation similar to (2) in the context of glacial cycles). Secondly, it is possible to efficiently express our beliefs and uncertainties about this simulator. Namely, uncertainty about the position of tipping points may be expressed as uncertainty on the coefficients (cf. Appendix 1).
Unfortunately the LLN2D is a very rough simulator by today’s standards. One should therefore attempt to correct or augment eq. (2) based on the knowledge gained from simulators more specifically designed to study the response of the atmosphere and oceans to changes in forcing and boundary conditions. This is where welldesigned experiments with general circulation simulators of the ocean and the atmosphere, such as those reviewed in the previous section, may be useful. Such plans of experiments may help us to locate points in the parameter space from which snow may accumulate from one year to the next. More generally, they should allow us to explore the (possibly nonlinear) response structure of the atmosphereocean system to changes in forcing and boundary conditions. To this end, statisticians have developed two concepts that should be useful:
 The experiment design theory.

This theory provides practical guidance to efficiently sample the space of possible inputs (forcing and parameters) of a numerical model in order to learn as much as possible about the simulator with as few experiments as possible (Santner et al., 2003).
 The Emulator
Just as the LLN2D surrogate, the emulator provides us with a fast, efficient, highlevel description of the general circulation simulator response. It may constitute the basis for nonlinear parameterisations of, the northern hemisphere snow accumulation balance response to obliquity, precession, CO and ice volume. In turn, these parameterisations can be included into a dynamical model of ice ages.
At last, it happens that certain processes supposed to relevant at the glacial time scales are difficult to simulate reliably with numerical simulators. For example, calculating the response of the ocean circulation near the Antarctic ice shelves to changes in sealevel is particularly challenging because it involves connexions between processes at very different spacial scales. Yet, Paillard and Parrenin (2004) speculated that these circulation changes may play a role on carbon cycle dynamics in glacial cycles, and they support their claim with a number qualitative arguments. Likewise, Ruddiman (2006) suggests a direct effect of precession on the carbon cycle through tropical dynamics and, possible, southern ocean temperatures. In such cases where numerical simulation are difficult or considered to be too unreliable, the framework of the climate stochastic dynamical model gives us the freedom to formulate hypotheses about the role of such processes in the form of simple equations, and then to test them against observations.
3.3 A statistical framework to accommodate observations
Many of the ideas expressed above were previously formulated by Barry Saltzman (2001). The present programme for investigation and prediction of ice ages proposed here goes further. It follows a Bayesian methodology because uncertainties are expressed by means of probability distribution functions and they are updated using observations.
Haslett et al. (2006) depicts the problem quite efficiently: the system defined by eq. (1a), (1b’) may be interpreted as a generator of a priori plausible climate histories. They are a priori plausible in the sense that they are generated from a system of equations and probability distributions of parameters that express, transparently, our beliefs about the dynamics of the climate system.
In practice, the calibration of this ‘generator’ consists in rating these climate histories against palaeoclimate observations. This operation consists in attributing more or less weight to the different possible model parameters, depending on the likelihood of the climate histories that they generate. This weighting process constitutes the mechanism by which the a priori uncertainty on model parameters is updated by accounting for observations. The more general principle of model selection consists, in the presence of two different models, in deciding which one produces the climatic histories that are the most likely given the observations. The concept of Bayesian model calibration was first applied to the problem of glacial cycle dynamics by Hargreaves and Annan (2002).
The likelihood is, in that framework, an important quantity. It is formally defined as the probability of observing what was observed, if the simulated palaeoclimate history was correct. Estimating the likelihood therefore requires to express carefully the connections between modelled climate histories (the ’model’), the real world, and observations. These connections may be visualised by means of a Bayesian graphical network, an illustrative example of which is given in Figure 7. The graphic is greatly inspired from that published by Haslett et al. (2006) but it is here adapted to the glacial cycles problem. Arrows indicate causal dependencies. Namely, our knowledge of benthic sampled at a certain depth depends on the that was actually preserved at that depth (after, e.g., bioturbation), which itself depends on the age corresponding to that depth, as well as on the trapped by the foraminifera at the time corresponding to that age. Further up, the generated climate histories depend on the parameters of the climate stochastic model and on the external forcing.
The climate stochastic dynamical system used to generate ‘palaeoclimate histories’ is thus part of a large framework generating potential ’palaeoclimate records’, chronologies and sedimentation histories. Haslett and Parnell (2008) provide an algorithm to generate sediment chronologies compatible with time markers (radiocarbon dates or, in our case, tephra layers and/or magnetic inversions).
Figure 7 is a simplified diagram; the full network should also include parameters controlling the response of proxies on climate. Nonetheless, it is good enough to clarify some of our ideas about palaeoclimate how inferences are drawn. For example, palaeoclimate scientists often take advantage of astronomical forcing to constrain the chronology of climate records Imbrie et al. (1984). The Bayesian graphical network shown here makes it clear that this operation requires a climate model (the white ellipse representing climate histories generated by the climate model lies between the forcing and observations).
Bayesian methodologies have become increasingly popular in climate science because they provide transparent ways of expressing our uncertainties, modelling choices (Annan and Hargreaves, 2006; Rougier, 2007; Rougier and Sexton, 2007; Sansó et al., 2008) and the distance estimated between the model and reality (on this specific point see Goldstein and Rougier (2009)).
Calibrating dynamical systems is however a very difficult problem. In general, the number of possible climate histories compatible with a model is so big that fairly complex algorithms are needed to reduce as much as possible the number of experiments required to efficiently update prior distributions. Palaeoclimates present additional specific challenges related to dating uncertainties and the complexity of the environmental factors affecting palaeoclimate archives. Hargreaves and Annan (2002) proposed a straightforward application of an algorithm well known by statisticians called ‘MarkovChain Metropolis Hastings’. The algorithm works well as long as the calibration period is relatively short and that climate trajectories all cluster around a same common history. Unfortunately this may not be the case if local instabilities occur, as in the example shown on Figure 3. Furthermore, the algorithm demands large computing resources. Crucifix and Rougier (2009) then proposed a solution based on a statistical algorithm called ‘particle filter for parameter and state estimation’ (Liu and West, 2001). This is a sequential filter: it ‘sweeps’ observations forward in time from the past until the present in order to generate a posteriori distributions of parameters and model states at the end of the run. The filter seemed to be a powerful solution to the problem posed by the existence of local instabilities. Unfortunately (again), further sensitivity studies performed after publication lead us to tone down our enthusiasm. The filter performs poorly on the model selection problem because it fails to discriminate models on the basis of their longterm dynamics. For example, some of the models selected by the filter no longer exhibit glacialinterglacial cycles .
Hope resides in more advanced Bayesian methods, which combine the MetropolisHastings strategy with the particle filter (Andrieu et al., 2010). An alternative solution is to calibrate the parameters on the basis of invariant summary statistics Wood (2010) using a method known an ’Approximate Bayesian Computation’ Sisson et al. (2007); Wilkinson (2008). Such statistics allow one to characterise a climatic trajectory in a way that is not sensitive to its initial conditions, nor to the exact timing of climatic events. For example, one may attempt to calibrate the dynamical system based on the average duration of a glacial cycle, or in the ‘skewness’ or ‘asymmetry’ of these cycles (King, 1996). More complex metrics may be envisaged to capture complex features associated with the nonlinear dynamics of the system (Marwan et al., 2009; Kantz and Schreiber, 2003, and. ref. therein). The use of summary statistics, however, should still respect the inference process outlined on Figure 7. An easy mistake is to apply summary statistics on climate model outputs and compare them straightforwardly with palaeoclimate data, without regard to the data preservation, sampling and possibly interpolation processes. Good summary statistics should allow us to efficiently discriminate between climate models and, at the same time, they should be robust to hypotheses about the data preservation and sampling processes. The works by Mudelsee and Stattegger (1994); Mudelsee (2001); Mudelsee et al. (2009) and ref. therein constitute important references in this respect.
4 Conclusion
The Pleistocene record suggests that modern insolation is not much above that needed for glacial inception. However, the complexity of the climate response to the complex astronomical forcing implies that further theoretical elaborations are needed to transform this statement into a reliable prediction about the slow evolution of climate like: ‘greenhouse gas should have declined during the Holocene, leading to glacial inception’. In particular, there is no perfect insolation analogue to the Holocene in the past and, would there be one, it cannot be guaranteed that climate would behave exactly the same way as during that hypothetical analogue. We mentioned several reasons for this. One is the possibility that small disturbances may, at strategic times, delay glacial events by several thousands of years. Locating such ’highsensitivity’ times is an attractive challenge. If the Holocene was such a time, any small disturbance, including the preindustrial anthropogenic activity, may have durably inflected the course of climate.
The dynamical system approach proposed here offers opportunities to tackle this problem in the tradition of hypotheticdeductive scientific investigation, subject to the usual criteria of model parsimony, prediction skill and accommodation of available observations. We have seen that very simple dynamical systems may already roughly reproduce the benthic curve; on the other hand, the useful output of more sophisticated numerical simulators may also be captured by fairly simply structured surrogates or emulators. These considerations encourage our quest for fairly compact and efficient predictive models of glacialinterglacial cycles. They are our main hope to credibly accommodate the fairly low amount of information present in palaeoclimate records with the flow of information generated by the complex simulators of the climate system. Indeed, reduced dynamical models designed to emulate complex numerical simulators contain a much smaller number of adjustable parameters and they are way more computationally efficient. They are therefore easier to calibrate on palaeoclimate observations. Furthermore, they make it possible to infer very general statements about climate predictability based on the large body of literature available on dynamical system theory.
Nevertheless, one should not be deceived by the apparent success our simple model and its close parents. Finding a dynamical system that convincingly accommodates more than one palaeoclimate data series is not trivial. A convincing model of glacial cycles should arguably accommodate observations about the oxygen isotopic composition of foraminifera, greenhouse gases, and carbon isotopic data, both in the oceans and in the atmosphere. It should explain spectral peaks, stochastic aspects, the asymmetry of cycles and slow changes in the system dynamics like those that seemingly occurred over the last million years. So far, this target “is still elusive” (Raymo and Huybers, 2008).
Bayesian inference constitutes a powerful framework to handle knowledge and uncertainties associated with palaeoclimate dynamics. This framework will lead us to formulate probabilistic statements about the next glacial inception and, more generally, assist our investigations about the climate system mechanics and dynamics.
Such a programme for statisticaldynamical modelling of (palaeo)climates is in its beginnings and serious challenges are faced. The more general problem of statistical calibration of dynamical systems is already quite difficult. Palaeoclimate data bear further difficulties, both at the climate system level (intertwining of dynamical structures at different time scales like DansgaardOeschger events and glacial cycles) and the observation level (time uncertainties, sediment disturbances, uncertain relationships between climate state and proxy…) It is therefore my opinion that we can’t yet express our predictions about natural trends in ice volume and greenhouse gases by probabilistic statements, and hence, that natural archives of previous interglacials do not falsify natural explanations for increasing greenhouse gas trends during the Holocene.
Appendix : the LLN2D surrogate
The form adopted for the surrogate is:
(3) 
where is ice volume, the ice volume increment over ka.; is the CO concentration, is a random number with Gaussian distribution, are regression coefficients. is further forced to be positive by imposing a lower bound of on .
Priors on regression coefficients are normal as follows (volumes are in , angles in radians, carbon dioxide concentrations in ppm and times in ka): , , , , , (conventionally, is a normal distribution of centre and standard deviation ). Posteriors are estimated by running a particle filter for state and parameter estimation (Liu and West, 2001) on an experiment spanning the interval forced by astronomical forcing (Berger, 1978a) and CO data by EPICA (Luethi et al., 2008) and Vostok (Petit et al., 1999). The likelihood function used to run the filter assumes that the errors of the surrogate every 1kyr time step are independent and normally distributed with a standard deviation equal to .
The red curve on Figure 6 is a climate simulation with the LLN2D model, while the blue curve is computed using the surrogate with the central estimates for . Grey curves are trajectories obtained with parameters sampled from the posterior distribution. Note that the timing of the deglaciations is here constrained by the history of CO imposed in these experiments.
Acknowledgements
Material present in this paper was greatly influenced by discussions at the SUPRAnet workshops organised by Caitlin Buck (University of Sheffield), in particular with Jonty Rougier, John Haslett and Andrew Parnell, as well as by further discussions with Jonty Rougier made possible with the exchange programme between the French Community of Belgium and the British council. The editor (Bill Ruddiman) and two anonymous reviewers provided particularly constructive comments. The author is funded by the Belgian Fund for Scientific Research. This is a contribution to the European Research Council starting grant ‘Integrating Theory and Observations over the Pleistocene’ nr ERC2009Stg 239604ITOP.
References
 Andrieu et al. (2010) Andrieu, C., Doucet, A., and Holenstein, R. 2010: Particle Markov chain Monte Carlo methods. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 72, 269–342.
 Annan and Hargreaves (2006) Annan, J. and Hargreaves, J.C. 2006: Using multiple observationallybased constraints to estimate climate sensitivity. Geophys. Res. Lett. 33, L06704.
 Berger (1977) Berger, A. 1977: Longterm variations of the Earth’s orbital elements. Celes. Mech. 15, 53–74.
 Berger (1979) — 1979: Insolation signatures of Quaternary climatic changes. Il Nuovo Cimento 2C, 63–87.
 Berger and Loutre (2002) Berger, A. and Loutre, M.F. 2002: An exceptionally long interglacial ahead? Science 297, 1287–1288.
 Berger (1978a) Berger, A.L. 1978a: Longterm variations of daily insolation and Quaternary climatic changes. J. Atmos. Sci. 35, 2362–2367.
 Berger (1978b) Berger, A. 1978b: Longterm variations of caloric insolation resulting from the earth’s orbital elements. Quaternary Research 9, 139 – 167.
 Calov and Ganopolski (2005) Calov, R. and Ganopolski, A. 2005: Multistability and hysteresis in the climatecryosphere system under orbital forcing. Geophysical Research Letters 32.
 Clark et al. (2006) Clark, P.U., Archer, D., Pollard, D., Blum, J.D., Rial, J.A., Brovkin, V., Mix, A.C., Pisias, N.G., and Roy, M. 2006: The middle Pleistocene transition: characteristics. mechanisms, and implications for longterm changes in atmospheric PCO2. Quat. Sci. Rev. 25, 3150–3184.
 Claussen et al. (2005) Claussen, M., Brovkin, V., Calov, R., Ganopolski, A., and Kubatzki, C. 2005: Did humankind prevent a Holocene glaciation ? Clim. Change 69, 409–417.
 Claussen et al. (2002) Claussen, M., Mysak, L., Weaver, A., Crucifix, M., Fichefet, T., Loutre, M.F., Weber, S., Alcamo, J., Alexeev, V., Berger, A., Calov, R., Ganopolski, A., Goosse, H., Lohmann, G., Lunkeit, F., Mokhov, I., Petoukhov, V., Stone, P., and Wang, Z. 2002: Earth system models of intermediate complexity: closing the gap in the spectrum of climate system models. Clim. Dyn. 18, 579–586.
 Crucifix and Berger (2006) Crucifix, M. and Berger, A. 2006: How long wil out interglacial be? Eos, Trans. Am. Geophys. Union 87, 352.
 Crucifix et al. (2005) Crucifix, M., Loutre, M.F., and Berger, A. 2005: Commentary on ”the anthropogenic greenhouse era began thousands of years ago”. Clim. Change 69, 419–426.
 Crucifix and Rougier (2009) Crucifix, M. and Rougier, J. 2009: On the use of simple dynamical systems for climate predictions: A Bayesian prediction of the next glacial inception. European Physics Journal  Special Topics 174, 11–31.
 Ditlevsen (2009) Ditlevsen, P.D. 2009: Bifurcation structure and noiseassisted transitions in the Pleistocene glacial cycles. Paleoceanography 24.
 Ditlevsen and Ditlevsen (2009) Ditlevsen, P.D. and Ditlevsen, O.D. 2009: On the Stochastic Nature of the Rapid Climate Shifts during the Last Ice Age. Journal of Climate 22, 446–457.
 EPICA community members (2004) EPICA community members 2004: Eight glacial cycles from an Antarctic ice core. Nature 429, 623–628.
 Gallée et al. (1992) Gallée, H., van Ypersele, J.P., Fichefet, T., Marsiat, I., Tricot, C., and Berger, A. 1992: Simulation of the last glacial cycle by a coupled, sectorially averaged climateice sheet model. Part II : Response to insolation and CO variation. J. Geophys. Res. 97, 15713–15740.
 Gildor and Tziperman (2001) Gildor, H. and Tziperman, E. 2001: Physical mechanisms behind biogeochemical glacialinterglacial CO variations. Geophysical Research Letters 28, 2421–2424.
 Goldstein and Rougier (2009) Goldstein, M. and Rougier, J. 2009: Reified Bayesian modelling and inference for physical systems. Journal of Statistical Planning and Inference 139, 1221 – 1239.
 Haken (2006) Haken, H. 2006: Information and SelfOrganization: A Macroscopic Approach to Complex Systems. Springer Series in Synergetics. Springer.
 Hargreaves and Annan (2002) Hargreaves, J.C. and Annan, J.D. 2002: Assimilation of paleodata in a simple Earth system model. Clim. Dyn. 19, 371–381.
 Haslett et al. (2006) Haslett, J., Whiley, M., Bhattacharya, S., SalterTownshend, M., Wilson, S.P., Allen, J.R.M., Huntley, B., and Mitchell, F.J.G. 2006: Bayesian palaeoclimate reconstruction. Journal of the Royal Statistical Society: Series A (Statistics in Society) 169, 395–438.
 Haslett and Parnell (2008) Haslett, J. and Parnell, A. 2008: A simple monotone process with application to radiocarbondated depth chronologies. Journal of the Royal Statistical Society: Series C (Applied Statistics) 57, 399–418.
 Hasselmann (1976) Hasselmann, K. 1976: Stochastic climate models. Part I: Theory. Tellus 28, 473–485.
 Huybers and Tziperman (2008) Huybers, P. and Tziperman, E. 2008: Integrated summer insolation forcing and 40,000year glacial cycles: The perspective from an icesheet/energybalance model. Paleoceanography 23.
 Imbrie et al. (1984) Imbrie, J.J., Hays, J.D., Martinson, D.G., McIntyre, A., Mix, A.C., Morley, J.J., Pisias, N.G., Prell, W.L., and Shackleton, N.J. 1984: The orbital theory of Pleistocene climate: Support from a revised chronology of the marine Orecord. In A. Berger, J. Imbrie, J. Hays, J. Kukla, and B. Saltzman, editors, Milankovitch and Climate, Part I, pages 269–305. D. Reidel, Norwell, Mass.
 Izhikevich and FitzHugh (2006) Izhikevich, E.M. and FitzHugh, R. 2006: FitzHughNagumo model. Scholarpedia 1, 1349.
 Kanamaru (2007) Kanamaru, T. 2007: Van der Pol oscillator. Scholarpedia 2, 2202.
 Kantz and Schreiber (2003) Kantz, H. and Schreiber, T. 2003: Nonlinear Time Series Analysis. Cambridge University Press.
 Kennedy and O’Hagan (2001) Kennedy, M. and O’Hagan, A. 2001: Bayesian calibration of computer models. J. Roy. Stat. Soc. B, 425–464.
 King (1996) King, T. 1996: Quantifying nonlinearity and geometry in time series of climate. Quaternary Science Reviews 15, 247 – 266.
 Kutzbach (2011) Kutzbach, J.E. 2011: Comparisons of coupled atmosphereocean simulations of greenhouse gasinduced climate change for preindustrial and hypothetical ‘noanthropogenic’ radiative forcing, relative to present day. The Holocene 21, in press.
 Lisiecki and Raymo (2005) Lisiecki, L.E. and Raymo, M.E. 2005: A PliocenePleistocene stack of 57 globally distributed benthic O records. Paleoceanogr. 20, PA1003.
 Liu and West (2001) Liu, J. and West, M. 2001: Combined parameter and state estimation in simulationbased filtering. In A. Doucet, N. de Freitas, and N. Gordon, editors, Sequential Monte Carlo Methods in Practice. Springer.
 Lorenz (1969) Lorenz, E.N. 1969: Three approaches to atmospheric predictability. Bull. Am. Meteorol. Soc. 50, 345–351.
 Lorenz (1970) — 1970: Climate change as a mathematical problem. J. Appl. Meteor. 9, 325–329.
 Loutre (2003) Loutre, M.F. 2003: Clues from MIS11 to predict the future climate. A modelling point of view. Earth Planet. Sci. Lett. 212, 213–234.
 Loutre and Berger (2003) Loutre, M.F. and Berger, A. 2003: Marine Isotope Stage 11 as an analogue for the present interglacial. Glob. Plan. Change 36, 209–217.
 Loutre et al. (2007) Loutre, M.F., Berger, A., Crucifix, M., Desprat, S., and SánchezGoñi, M.F. 2007: Interglacials as simulated by the LLN2D NH and MoBidiC climate models. In F. Sirocko, M. Claussen, M.F. Sánchez Goñi, and T. Litt, editors, The climate of past interglacials, volume 7 of Developments in Quaternary Science, pages 547–582. Elsevier.
 Loutre et al. (2008) Loutre, M.F., Crucifix, M., and Berger, A. 2008: Chasing an Analogue for the Holocene : The astronomical forcing. Eos, Trans. Am. Geophys. Union 89, Fall Meet. Suppl., Abstract U33B–01.
 Luethi et al. (2008) Luethi, D., Le Floch, M., Bereiter, B., Blunier, T., Barnola, J.M., Siegenthaler, U., Raynaud, D., Jouzel, J., Fischer, H., Kawamura, K., and Stocker, T.F. 2008: Highresolution carbon dioxide concentration record 650,000800,000 years before present. Nature 453, 379–382.
 Marwan et al. (2009) Marwan, N., Donges, J.F., Zou, Y., Donner, R.V., and Kurths, J. 2009: Complex network approach for recurrence analysis of time series. Physics Letters A 373, 4246 – 4254.
 Mudelsee (2001) Mudelsee, M. 2001: The phase relations among atmospheric CO content, temperature and global ice volume over the past 420 ka. Quaternary Science Reviews 20, 583–589.
 Mudelsee et al. (2009) Mudelsee, M., Scholz, D., Röthlisberger, R., Fleitmann, D., Mangini, A., and Wolff, E.W. 2009: Climate spectrum estimation in the presence of timescale errors. Nonlinear Processes in Geophysics 16, 43–56.
 Mudelsee and Stattegger (1994) Mudelsee, M. and Stattegger, K. 1994: Application of the Grassberger. Procaccia algorithm to the 6180 record from ODP site 659: selected methodical aspects. In J.H. Kruhl, editor, Fractals and dynamics systems in geosciences, pages 399–413. SpringerVerlag.
 Murphy et al. (2004) Murphy, J.M., Sexton, D.M.H., Barnett, D.N., Jones, G.S., Webb, M.J., Collins, M., and Stainforth, D.A. 2004: Quantification of modelling uncertainties in a large ensemble of climate change simulations. Nature 430, 768–772.
 Nicolis (2005) Nicolis, C. 2005: Can error source terms in forecasting models be represented as Gaussian Markov noises? Quarterly Journal of the Royal Meteorological Society 131, 2151–2170.
 Nicolis and Nicolis (2007) Nicolis, G. and Nicolis, C. 2007: Foundations of complex systems: nonlinear dynamics, statistical physics, information and prediction. World Scientific.
 Paillard (1998) Paillard, D. 1998: The timing of Pleistocene glaciations from a simple multiplestate climate model. Nature 391.
 Paillard (2001) — 2001: Glacial Cycles: Toward a New Paradigm. Rev. Geophys. 39.
 Paillard and Parrenin (2004) Paillard, D. and Parrenin, F. 2004: The Antarctic ice sheet and the triggering of deglaciations. Earth Planet. Sc. Lett. 227, 263–271.
 Palmer (2005) Palmer, T. 2005: Global warming in a nonlinear climate  Can we be sure? Europhysics News 36, 42–46.
 Palmer and Williams (2010) Palmer, T. and Williams, P. 2010: Stochastic Physics and Climate Modeling. Camb. Univ. Press.
 Petit et al. (1999) Petit, J.R., Jouzel, J., Raynaud, D., Barkov, N.I., Barnola, J.M., Basile, I., Bender, M., Chappellaz, J., Davis, M., Delaygue, G., Delmotte, M., Kotlyakov, V.M., Legrand, M., Lipenkov, V.Y., Lorius, C., Pepin, L., Ritz, C., Saltzman, E., and Stievenard, M. 1999: Climate and atmospheric history of the past 420, 000 years from the Vostok ice core, Antarctica. Nature 399, 429–436.
 Petoukhov et al. (2000) Petoukhov, V., Ganopolski, A., Brovkin, V., Claussen, M., Eliseev, A., Kubatzki, C., and Rahmstorf, S. 2000: CLIMBER2: a climate system model of intermediate complexity. Part I: model description and performance for present climate. Clim. Dyn. 16, 1–17.
 Rahmstorf et al. (2005) Rahmstorf, S., Crucifix, M., Ganopolski, A., Goosse, H., Kamenkovich, I., Knutti, R., Lohmann, G., Marsh, R., Mysak, L.A., Wang, Z., and Weaver, A.J. 2005: Thermohaline circulation hysteresis: A model intercomparison. Geophys. Res. Lett. 32, L23605.
 Raymo and Huybers (2008) Raymo, M.E. and Huybers, P. 2008: Unlocking the mysteries of the ice ages. Nature 451, 284–285.
 Rohling et al. (2010) Rohling, E., Braun, K., Grant, K., Kucera, M., Roberts, A., Siddall, M., and Trommer, G. 2010: Comparison between Holocene and Marine Isotope Stage11 sealevel histories. Earth and Planetary Science Letters 291, 97 – 105.
 Rougier (2007) Rougier, J. 2007: Probabilistic inference for future climate using an ensemble of climate model evaluations. Climatic Change 81, 247–264.
 Rougier and Sexton (2007) Rougier, J. and Sexton, D.M.H. 2007: Inference in ensemble experiments. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 365, 2133–2143.
 Rougier et al. (2009) Rougier, J., Sexton, D.M.H., Murphy, J.M., and Stainforth, D. 2009: Analyzing the Climate Sensitivity of the HadSM3 Climate Model Using Ensembles from Different but Related Experiments. Journal of Climate 22, 3540–3557.
 Ruddiman (2003) Ruddiman, W.F. 2003: The anthropogenic greenhouse era began thousands of years ago. Clim. Change 61, 261–293.
 Ruddiman (2006) — 2006: Orbital changes and climate. Quat. Sci. Rev. 25, 3092–3112.
 Ruddiman (2007) — 2007: The early anthropogenic hypothesis a year later: challenges and responses. Rev. Geophys. 45, RG4001.
 Ruddiman et al. (2011) Ruddiman, W.F., Kutzbach, J.E., and Vavrus, S. 2011: Can natural or anthropogenic explanations of late Holocene CO and CH increases be falsified? The Holocene 21.
 Rulkov et al. (1995) Rulkov, N.F., Sushchik, M.M., Tsimring, L.S., and Abarbanel, H.D.I. 1995: Generalized synchronization of chaos in directionally coupled chaotic systems. Phys. Rev. E 51, 980–994.
 Saltzman (2001) Saltzman, B. 2001: Dynamical paleoclimatology: Generalized Theory of Global Climate Change (International Geophysics), volume 80 of International Geophysics Series. Academic Press.
 Saltzman and Maasch (1990) Saltzman, B. and Maasch, K.A. 1990: A firstorder global model of late Cenozoic climate. Trans. R. Soc. Edinburgh Earth Sci 81, 315–325.
 Saltzman and Maasch (1991) — 1991: A firstorder global model of late Cenozoic climate. II Further analysis based on a simplification of the CO dynamics. Clim. Dyn. 5, 201–210.
 Sansó et al. (2008) Sansó, B., Forest, C.E., and Zantedeschi, D. 2008: Inferring Climate System Properties Using a Computer Model. Bayesian Analysis 3, 1–38.
 Santner et al. (2003) Santner, T., Williams, B., and Notz, W. 2003: The Design and Analysis of Computer Experiments. New York: Springer.
 Sisson et al. (2007) Sisson, S.A., Fan, Y., and Tanaka, M.M. 2007: Sequential Monte Carlo without likelihoods. Proceedings of the National Academy of Sciences 104, 1760–1765.
 Tzedakis et al. (2009) Tzedakis, P.C., Raynaud, D., McManus, J.F., Berger, A., Brovkin, V., and Kiefer, T. 2009: Interglacial diversity. Nature Geosci 2, 751–755.
 Tziperman et al. (2006) Tziperman, E., Raymo, M.E., Huybers, P., and Wunsch, C. 2006: Consequences of pacing the Pleistocene 100 kyr ice ages by nonlinear phase locking to Milankovitch forcing. Paleoceanography 21, PA4206.
 Vavrus et al. (2011) Vavrus, S., Kutzbach, J.E., and Ruddiman, W.F. 2011: The Role of GCM Resolution in Simulating Glacial Inception. The Holocene 21, in press.
 Vettoretti and Peltier (2004) Vettoretti, G. and Peltier, W.R. 2004: Sensitivity of glacial inception to orbital and greenhouse gas climate forcing. Quaternary Science Reviews 23, 499–519.
 Vettoretti and Peltier (2011) — 2011: The Impact of Insolation, Greenhouse Gas Forcing and Ocean Circulation Changes on Glacial Inception. The Holocene 21, in revision.
 Wilkinson (2008) Wilkinson, R.D. 2008: Approximate Bayesian computation (ABC) gives exact results under the assumption of model error. ArXiv eprints page http://arxiv.org/abs/0811.3355v1.
 Wood (2010) Wood, S.N. 2010: Statistical inference for noisy nonlinear ecological dynamic systems. Nature 466, 1102–1104.
List of Figures
 1 Three measures of insolation covering the period ka to ka (departure to the present) calculated using the Berger (1978a) solution, aligned with a stack of 57 benthic foraminifera (Lisiecki and Raymo, 2005) and a reconstruction of atmospheric CO concentration (Luethi et al., 2008; Petit et al., 1999). Circles mark insolations equal to the present. Times corresponding to CO concentrations above 250 ppm and benthic below 3.8 per mil are highlighted to subjectively mark interglacial conditions.
 2 Timeseries generated with a twovariable model, equations (1a) and (1b). The system is forced by a synthetic mix of obliquity and precession. The variable integrates the forcing with a correction . shifts between two states according to a hysteresis behaviour controlled by . In this conceptual model may be interpreted as the total volume of continental ice, and is often associated to the state of the ocean circulation. is here superimposed to the benthic stack by Lisiecki and Raymo (2005), which is appropriately scaled to highlight the similarity between the simulation and the palaeoclimate record.
 3 Same model as Figure 2 but with a small stochastic forcing added to (equation 1b’). The two timeseries shown correspond to two realisations of the stochastic forcing. Observe the phase shift aroud ka. It corresponds to temporary loss of synchronisation of the system on its (astronomical) forcing.
 4 Route map for the design of a stochastic generator of climate histories. Statistical approximations of large numerical simulators (surrogates or emulators) are developed on the basis of appropriate plans of experiments. These mathematical objects are much smaller and tractable than the original simulators, and may be coupled with each other and/or with more conceptual models. Expert elicitation is required for assessing the uncertainties associated with numerical simulators and formulate of prior probability distributions of parameters. The resulting stochastic dynamical system may then be used to generate ensembles of palaeoclimate histories compatible with this model.
 5 Simulated volume of ice over the northern hemisphere near steadystate as a function of precession using the LLN2D simulator by Gallée et al. (1992), and given fixed eccentricity, obliquity and CO concentration. The experiment consists in varying the precession parameter slowly all along the trigonometric cycle in 400 ka of model time.
 6 Simulation of ice volume with the LLN2D model (red) forced by astronomical forcing and variations in CO (Petit et al., 1999; Luethi et al., 2008), compared with 50 realisations of the emulator (eq. (2)). The blue trajectory is obtained with the emulator using the central estimates of regression coefficients. All experiments are deterministic.
 7 Simplified Bayesian inference network for selection and calibration of stochastic dynamical systems of palaeoclimates (cf. Haslett et al. (2006) for more theoretical background). Rectangles correspond to observations or certain quantities. Ellipses are uncertain quantities. Continuous functions of time are symbolically marked with ; functions of depth by and discrete series are marked with . The sequence from left to right corresponds to forcing, climate model, natural archiving of the climate information (blue), sampling and measure (yellow). Modelling consists in defining the relationships symbolised by the arrows and attributing prior probability distributions to uncertain quantities at the top of the chain (e.g.: model parameters). Model validation and selection are judged based on whether actual observations are compatible the model (measured by the likelihood); statistical calibration consists in propagating the information backwards, from the observations to parameters in order to make predictions (green). ’Summary statistics’ are variables such as spectral power or other integrated measures which summarise several observations.