Terminations VI and VIII ( (cid:24) 530 and (cid:24) 720 kyr BP) tell us the importance of obliquity and precession in the triggering of deglaciations

. The main variations of ice volume of the last million years can be explained from orbital parameters by assuming climate oscillates between two states: glaciations and deglaciations (Parrenin and Paillard, 2003; Imbrie et al., 2011) (or terminations). An additional combination of ice volume and orbital parameters seems to form the trigger of a deglaciation, while only orbital parameters seem to play a role in the triggering of glaciations. Here we present an optimized conceptual model which realistically reproduce ice volume variations during the past million years and in particular the timing of the 11 canonical terminations. We show that our model looses sensitivity to initial conditions only after ∼ 200 kyr at maximum: the ice volume observations form a strong attractor. Both obliquity and precession seem necessary to reproduce all 11 terminations and both seem to play approximately the same role. More precisely, obliquity plays a fundamental role in the triggering of termination VI ( ∼ 530 kyr BP),

Abstract. The main variations of ice volume of the last million years can be explained from orbital parameters by assuming climate oscillates between two states: glaciations and deglaciations (Parrenin and Paillard, 2003;Imbrie et al., 2011) (or terminations). An additional combination of ice volume and orbital parameters seems to form the trigger of a deglaciation, while only orbital parameters seem to play a role in the triggering of glaciations. Here we present an optimized conceptual model which realistically reproduce ice volume variations during the past million years and in particular the timing of the 11 canonical terminations. We show that our model looses sensitivity to initial conditions only after ∼ 200 kyr at maximum: the ice volume observations form a strong attractor. Both obliquity and precession seem necessary to reproduce all 11 terminations and both seem to play approximately the same role. More precisely, obliquity plays a fundamental role in the triggering of termination VI (∼530 kyr BP), while precession plays a fundamental role in the triggering of termination VIII (∼720 kyr ago).

Introduction
Understanding past climates could help us to improve our predictions of future climatic variations. The reconstructions of Earth's climate over the past million years from either ice cores  or marine cores (Lisiecki and Raymo, 2005) show a succession of long glaciations and short deglaciations (or terminations), with a period of ∼ 100 kyr and known as glacial-interglacial cycles. Changes in in-coming solar insolation (Loutre, 1993) due to changes in Earth orbital parameters (Berger, 1978;Laskar et al., 2004) is the only known major external forcing of the climate system at these time scales. Investigation in the frequency domain suggests that variations in Earth's orbit indeed paced the observed changes (Hays et al., 1976), but the amplitude and saw-tooth shape of climatic variations implies that amplifications (e.g. through ice/snow albedo and greenhouse gases changes) and non-linearities exist. To trigger a complete deglaciation, it seems that a large ice volume and appropriate orbital parameters are necessary (Raymo, 1997), while only orbital forcing (and amplifying feedbacks) seems important to trigger a glaciation (Paillard, 1998;Khodri et al., 2001). Paillard (1998) used such ideas to build a 3-states climate system which correctly simulates (within a few kyr, 1000 yr) the timing of terminations (except termination VI) during the last million years. This idea has been further developed by assuming an additional combination of ice volume and orbital parameters forms the trigger of a termination (Parrenin and Paillard, 2003;Imbrie et al., 2011). It has also been suggested (Parrenin and Paillard, 2003) that ice volume influences the exact timing of deglaciation, and the prediction of this model has been later checked with accurate chronologies of the Dome Fuji and Vostok ice cores . This is a element of proof that such conceptual models are not useless and that they can have some predictive capabilities. More recently, based on slightly different evolution equations but on a similar deglaciation threshold than Parrenin and Paillard (2003), a 2-states conceptual climate model was proposed (Imbrie et al., 2011) (Laskar et al., 2004), ice volume data (Bintanja et al., 2005)  that the increase in 100 kyr power during the mid-Pleistocene transition could be caused by changes in Earth orbital parameters rather than by changes within the climate system. Two debates still take place in the circle of Quaternary climate scientists. The first question is whether the climate system is deterministic or stochastic, in other words is climate state predictable or not (see Crucifix and Rougier (2009) for a review of this question)? The second question concerns the timing of terminations, which have been proposed to be linked either to obliquity (Huybers and Wunsch, 2005) or to precession (Tzedakis et al., 2012) or to both (Huybers, 2011).
In this study we develop a 2-states (glaciation, deglaciation) conceptual model of Quaternary ice volume similar to Parrenin and Paillard (2003) and tune it to observed ice volume. We first determine the sensitivity of our model to initial conditions to determine whether the climate system is deterministic or stochastic. Second, we examine the relative role of obliquity and precession in the triggering of deglaciations.

Forward model description
Our study is based on a conceptual model of Quaternary climate which simulates the ice volume from the orbital parameters (Paillard et al., 1996;Laskar et al., 2004). As for the model of Imbrie et al. (2011), this model takes as inputs 3 functions of the orbital parameters which are normalized to zero mean and unit variance on the last million years: Esi ∼ e sin ω (precession, with ω the precession angle taken from the vernal equinox), Eco ∼ e cos ω (phase-shifted precession) and O ∼ ε (obliquity). Insolation at most latitudes and seasons can be represented quite accurately by a linear combination of these three orbital functions (Loutre, 1993). The model of Parrenin and Paillard (2003) was taking as input obliquity and summer solstice insolation at 65 • N and this modification allows to disentangle the role of obliquity and precession. For the rest, the model is similar to the one in Parrenin and Paillard (2003) except that the predicted ice volume is now dimensional which makes the comparison with the observed ice volume easier. It has two different states of evolution: the glaciation state g and the deglaciation state d and the evolution of ice volume v (expressed in m sea level) in these states is simply described by two linear equations: where "O" is obliquity normalized to zero mean and unit variance and Esi tr and Eco tr are respectively calculated from Esi and Eco the precession parameters using a truncation function: if x>0 : (where a is a constant) and then normalized to zero mean and unit variance. This truncation is similar to the one used by Paillard (1998) and appears necessary to simulate the lower ice volume sensitivity to precession during cold periods than during warm periods. We now need to define when the model jumps from one state to the other. For this we define thresholds on a linear combination of orbital parameters (and ice volume for the g-to-d transition): We start with v = v init in the S init state at t = 1000 kyr BP and we solve the evolution of v with a Runge-Kutta 4th order method and with a time step of 100 yr.

Monte Carlo fitting of parameters
To infer the value of the parameters of this model, we fit it to an ice volume reconstruction (Bintanja et al., 2005) based on the LR04 marine isotopic stack (Lisiecki and Raymo, 2005).
In the time space, the good agreement of the LR04 with the independent glaciological EDC3 age scale  for the last 400 kyr suggests that the age scale errors do not exceed ∼ 3 kyr. In the ice-volume space, the errors of the Bintanja et al. (2005) reconstruction are taken into account as follows. We define the density of probability of an ice volume simulation v(t): where k is a multiplicative constant and where: is the transposed residual vector with t 0 = 1000, t 1 = 999, . . . t N = 0 kyr BP, v d (t) is the observed sea level over the last million years (Bintanja et al., 2005;Lisiecki and Raymo, 2005) and C, the covariance matrix, takes into account both the modeling and data errors (Tarantola, 1987). We assume that the errors are independent (C is diagonal) and with a confidence interval σ d = 20 m. We then explore the parameters space using a 1 000 000 experiments random walk based on the Metropolis-Hastings algorithm (Metropolis et al., 1953;Hastings, 1970)  one. This calibration procedure is similar to the one employed by Hargreaves and Annan (2002).

Simulation of the last million years ice volume
Our optimized model (experiment called Best, see Table 1 and Fig. 1) is in good agreement with the data (standard deviation of the residuals of 12.5 m) and in particular reproduces correctly all the 11 terminations of the last million years. The agreement with the data is slightly better than for the model described by Parrenin and Paillard (2003) because of the use of phase-shifted precession and of the Monte Carlo optimization algorithm. The time periods when the model deviates most from Bintanja's sea level reconstruction are marine isotope stages (MIS) 13 and 7, corresponding to terminations III and VI. There, the model also simulates larger terminations with deglaciated MIS 13.3 and 7.5. The model then glaciates too much for the following stages (MIS 13.1 and 7.3-7.1). This model-data mismatch is not solved when using other sea level proxies for MIS 7 (Dutton et al., 2009).
It thus seem that terminations III and VI happen during two precessional cycles and it is where our simple modeling formulation shows its limits (this is also true for terminations VII and VIII but with a lower extend). Note that the model by Imbrie et al. (2011) does reproduce high sea levels during MIS 13.1 and 7.3-7.1 in better agreement with the data, by simulating terminations at the onset of these time periods. However, to do this, this model simulates for previous MIS 13.2 and 7.4 time periods very glaciated states, which is not in agreement with the data. Table 2 gives the timing of the onsets and ends of terminations (as defined in the model by the threshold crossing) in the Best experiment as well as their durations, which are also plotted in Fig. 2. TV is by far the longest termination with 21.6 kyr. TVIII is the shortest with 6.1 kyr. All the other terminations have a duration ranging between 8.8 and 14.6 kyr. We do not find two groups of terminations based on their duration as was recently suggested for interglacials (Tzedakis et al., 2012). Figure 3 also illustrates the timing of the closest precession and obliquity maxima with respect to our termination onsets. All maxima happen after the termination onset, except for the obliquity maxima at TIII and the precession maxima at TV. Not surprisingly, the precession and obliquity maxima timings are anti-correlated which is a consequence of our Eq. (4). One can also see in Fig. 3 that these timings are strongly varying: precession timings have a standard deviation of ∼ 4 kyr while obliquity timings have a standard deviation of ∼ 6 kyr. This shows that, even in a orbitally-forced model, the phase relationship between climate and orbital parameters may not be constant. This is also in contradiction with a recent study suggesting that onset of interglacials are always near a precession maximum (Tzedakis et al., 2012), also we are not speaking of exactly the same thing here (onset of terminations vs onset of interglacials).

Complexity reduction
The data have ∼ 90 extremas (two per precessional cycle), i.e. ∼ 90 degrees of freedom (the degree of a fitting polynomial). The model has 14 parameters. Among those, two are boundary conditions (S init and v init ) that could be fixed according to the data. Two others are optional: a similar agreement with the data is found with α Eco = 0 and κ Eco = 0 (see Supplement). So the complexity of the system is significantly reduced from 90 to 10 degrees of freedom. Considering the timing of the 11 terminations, the 3 parameters that matter are κ Esi , κ O and v 0 , which represent a reduction of the complexity from 11 to 3 degrees of freedom. This complexity reduction suggests that our model does capture the main structure of the Quaternary climate variations, although we cannot exclude that there is no other satisfying modeling based on different concepts and although our equations need to be linked to precise physical mechanisms.

Sensitivity to initial conditions
To answer the question whether Quaternary climate is stochastic or deterministic in our conceptual model, we perform a sensitivity experiment with respect to the initial condition at time 1000 kyr BP. All model parameters are kept as in the Best model experiment except v init which is given the values 15 m, 30 m, 45 m, 60 m and 75 m. The results are displayed in Fig. 4. One can see that 4 initial values (15 m, 45 m, 60 m and 75 m) give scenarios close to the observed ice volume. Only the 30 m initial value gives a different scenario for the first 200 kyr, but then converges again to the data. From these observations, we conclude that the ice volume data form a strong attractor in our conceptual model.

Importance of obliquity and precession in the triggering of deglaciation
The question now is how necessary are precession and obliquity to simulate the last million year ice volume termina- tions? In other words, using our model formulation, what are the acceptable values for the κ P = √ (α 2 Esi + α 2 Eco ) and κ O parameters?
In our Control experiment, obliquity and precession play approximately the same role in the deglaciation trigger with a marginally more important role for obliquity: κ P ≈15 m and κ O ≈18 m. Figure 5 displays κ O vs. κ P and the κ P /κ O ratio for all experiments selected by the Monte Carlo algorithm. One can see that κ P is slightly better constrained than κ O and that the κ P /κ O ratio is well constrained to 0.87 ± 0.07. This means that in all selected experiments, obliquity and precession approximately play the same role in the deglaciation trigger. Therefore, both obliquity and precession seem necessary to explain the last million years terminations.
To reinforce this affirmation, we apply our Monte Carlo optimization algorithm assuming κ O = 0 and find a most probable experiment called Best-wo (see Table 1 and Fig. 6).
In this experiment, all the terminations are correctly placed except termination VI. It is not surprising that the model fails to reproduce termination VI, since it has devaforable precession (low Esi maximum) and ice volume (low ice volume at MIS14.2) configurations and only a favorable obliquity configuration.
If we try to force more terminations by decreasing the deglaciation threshold v 0 to 98 m in Eq. (4) (experiment Test-wo, see supplement), deglaciations appear at ∼ 270 and 810 kyr BP, periods with a precession parameter larger than that of termination VI, but a lower obliquity. The fact that some terminations appear both before and after termination VI exclude a long-term trend in one of the model's parameters (Paillard, 1998) to explain the observations.
We now test our model without precession influence on deglaciation threshold and apply our Monte Carlo optimization algorithm with κ Esi = 0 and κ Eco = 0. We find a bestguess experiment called Best-wp (see Table 1 and Fig. 6). Interestingly, this experiment simulates a lot more terminations than the Best experiment, but some do not last very long. All the 11 canonical terminations are reproduced, except TVIII which is shifted ∼ 20 kyr toward younger ages. It is not surprising that the model fails to reproduce termination VIII, since it has a unfavorable obliquity configuration but a favorable precession configuration.
If, here also, we try to force more terminations by decreasing the deglaciation threshold v 0 to 80 m (experiment Testwp, see Supplement), some terminations appear at ∼ 760 kyr or 180 kyr BP before termination VIII is reached.

Conclusions
A conceptual model of continental ice volume, taking as input only the Earth orbital parameters, has been built and successfully simulates the observations within their confidence interval. This model suggests strongly varying duration and timing of terminations with respect to orbital parameters. Our model is not very sensitive to initial conditions, which suggests that climate evolution was deterministic during the Quaternary. A study of the relative role of obliquity and precession in the triggering of deglaciation suggests that both orbital parameters are necessary to explain all the 11 observed terminations of the last million years and that both played approximately the same quantitative role. More precisely, termination VI cannot be explained without any influence of obliquity while termination VIII cannot be explained without any influence of precession. TVI and the precursor of TVIII have already been emphasized by their lag of CO 2 to global ice volume (Lisiecki, 2010). By this article, we hope to stimulate further studies focused on these terminations.