Climate of the Past The role of orbital forcing , carbon dioxide and regolith in 100 kyr glacial cycles

The origin of the 100 kyr cyclicity, which dominates ice volume variations and other climate records over the past million years, remains debatable. Here, using a comprehensive Earth system model of intermediate complexity, we demonstrate that both strong 100 kyr periodicity in the ice volume variations and the timing of glacial terminations during past 800 kyr can be successfully simulated as direct, strongly nonlinear responses of the climate-cryosphere system to orbital forcing alone, if the atmospheric CO 2 concentration stays below its typical interglacial value. The existence of long glacial cycles is primarily attributed to the North American ice sheet and requires the presence of a large continental area with exposed rocks. We show that the sharp, 100 kyr peak in the power spectrum of ice volume results from the long glacial cycles being synchronized with the Earth’s orbital eccentricity. Although 100 kyr cyclicity can be simulated with a constant CO 2 concentration, temporal variability in the CO2 concentration plays an important role in the amplification of the 100 kyr cycles.


Introduction
Although it is generally accepted that, as postulated by the Milankovitch theory (Milankovitch, 1941), Earth's orbital variations play an important role in Quaternary climate dynamics, the nature of glacial cycles still remains poorly understood.One of the major challenges to the classical Milankovitch theory is the presence of 100 kyr cycles that dominate global ice volume and climate variability over the past million years (Hays et al., 1976;Imbrie et al., 1993;Paillard, 2001).This periodicity is practically absent in the principal "Milankovitch forcing" -variations of summer insolation Correspondence to: A. Ganopolski (andrey@pik-potsdam.de) at high latitudes of the Northern Hemisphere (NH).The eccentricity of Earth's orbit does contain periodicities close to 100 kyr and the robust phase relationship between glacial cycles and 100-kyr eccentricity cycles has been found in the paleoclimate records (Hays et al., 1976;Berger et al., 2005;Lisiecki, 2010).However, the direct effect of the eccentricity on Earth's global energy balance is very small.Moreover, eccentricity variations are dominated by a 400 kyr cycle which is also seen in some older geological records (e.g.Zachos et al., 1997), but is practically absent in the frequency spectrum of the ice volume variations for the last million years.In view of this long-standing problem, it was proposed that the 100 kyr cycles do not originate directly from the orbital forcing but rather represent internal oscillations in the climate-cryosphere (Gildor and Tziperman, 2001) or climate-cryosphere-carbonosphere system (e.g.Saltzman and Maasch, 1988;Paillard and Parrenin, 2004), which can be synchronized (phase locked) to the orbital forcing (Tziperman et al., 2006).Alternatively, it was proposed that the 100 kyr cycles result from the terminations of ice sheet buildup by each second or third obliquity cycle (Huybers and Wunsch, 2005) or each fourth or fifth precessional cycle (Ridgwell et al., 1999) or they originate directly from a strong, nonlinear, climate-cryosphere system response to a combination of precessional and obliquity components of the orbital forcing (Paillard, 1998).Since a number of conceptual models based on fundamentally different assumptions were able to reproduce reconstructed ice volume variations with similar skill, it became clear that a further advance in understanding of 100 kyr cyclicity requires physically-based models.
Simulations with simplified climate-cryosphere models (Pollard, 1983;Deblonde and Peltier, 1991;Berger et al., 1999;Crowley and Hyde, 2008) have shown that 100 kyr cyclicity does appear in ice volume variations driven by orbital variations alone.However, in most cases, the simulated 100 kyr cycles were weaker than in the paleoclimate records A. Ganopolski and R. Calov: The role of orbital forcing, carbon dioxide and regolith in 100 kyr glacial cycles and were additionally accompanied by pronounced variability at another eccentricity frequency -400 kyr -which is not seen in the spectra of reconstructed ice volume.It was only when realistic CO 2 forcing was applied in addition to orbital forcing that realistic simulations of the glacial cycles became possible (Berger et al., 1998).The notable exception is the work by Pollard (1983) wherein after adding several nonlinear process, the model forced by orbital variations alone simulates strong 100 kyr cycles in agreement with the ice volume reconstructions available at that time.It is interesting to note that the agreement is even more impressive when Pollard's modelling results are compared to the most recent reconstructions of ice volume.
Although simplified climate-cryosphere models demonstrate the possibility of the appearance of the 100 kyr cycle as a direct response of the climate-cryosphere system to the orbital forcing, due to their simplicity (usually these models were based on a one-dimensional ice sheet model and an energy balance atmosphere model), doubt remains whether these results are fully applicable to the real world.Moreover, the presence of 400 kyr cycles in many simulations remains an obvious problem.Therefore, it is crucial to corroborate earlier results and further advance the understanding of glacial cycles by using more physically based and geographically explicit climate-cryosphere models.While coupled GCMs still remain too expensive for simulating glacial cycles, models of intermediate complexity (EMICs, Claussen et al., 2002) can be coupled to 3-D ice sheet models and are sufficiently computationally efficient to perform simulations of glacial cycles.Using CLIMBER-2 coupled to different ice sheet models, Bonelli et al. (2009) and Ganopolski et al. (2010) performed simulations of the last glacial cycles; Calov and Ganopolski (2005) analysed the stability of the climate-cryosphere system in the phase space of Milankovitch forcing and Bauer and Ganopolski (2010) reported simulations of the last four glacial cycles.Here we will present a large suite of simulations for the last 800 kyr, a period of time dominated by 100 kyr cyclicity.

Model description and experimental setup
The model used in the study is the most recent version of the Earth system model of intermediate complexity CLIMBER-2 (Petoukhov et al., 2000;Ganopolski et al., 2001;Brovkin et al., 2002), which includes a 3-D polythermal ice sheet model (Greve, 1997).The ice sheet model is only applied to the Northern Hemisphere and is coupled to the climate component via a high-resolution, physically-based surface energy and mass balance interface (Calov et al., 2005), which explicitly accounts for the effect of aeolian dust deposition on snow albedo.Here we use the same model and experimental setup as Ganopolski et al. (2010) but extended it to simulate glacial cycles over the past 800 000 years.
In all experiments the equilibrium state of the climatecryosphere system obtained for present-day conditions was used as the initial condition and the model was run from 860 kyr BP until the present.The first 60 000 years, representing the model spin-up, were not used for further analysis.
In the Baseline Experiment (referred to hereafter as BE), we prescribed variations in orbital parameters following Berger (1978) and the equivalent CO 2 concentration, which accounts for the radiative forcing of three major greenhouse gases -carbon dioxide, methane and nitrous oxide.Their concentrations were derived from the Antarctic ice cores (Petit et al., 1999;EPICA, 2004).The method used to calculate the equivalent CO 2 concentration is described by Ganopolski et al. (2010).A continuous record of N 2 O is not available for the last 800 kyr, but existing data suggest that, to the first approximation, the N 2 O concentration has a temporal dynamic similar to CO 2 .Therefore, we assumed that the radiative forcing of N 2 O (relative to preindustrial) is 20 % of that of CO 2 during the whole simulated period, i.e. the ratio between radiative forcings of N 2 O and CO 2 is the same as at the LGM.
Although concentrations of GHGs from the ice cores are only available for the last 800 000 years, the time 800 kyr BP is not the best choice for the beginning of the simulations because it was close to a glacial maximum and therefore would require an initialization of the large continental ice sheets in the Northern Hemisphere.For this reason, we begin our simulations at 860 kyr BP, which corresponds to the MIS 21 interglacial for which we can use the equilibrium present-day climate state as initial conditions.However, this choice of the initial state requires prescription of the equivalent CO 2 concentration for the time interval when reliable data for GHGs concentration are not yet available.To extend the time series of equivalent CO 2 concentration beyond 800 kyr BP, we made use of a close correlation between the total radiative forcing of GHGs and benthic δ 18 O c stack (Lisiecki and Raymo, 2005) observed for the last 800 kyr.By using a simple linear regression, we calculated equivalent CO 2 concentration for these initial 60 kyr.Since this period was considered as the model spin-up and was not used for further analysis, the accuracy of this reconstruction of the equivalent CO 2 is not crucial for the results presented in the paper.
In addition to the BE, we performed a large suite of experiments with constant CO 2 , modified orbital forcing, and terrestrial sediment mask.These experiments are summarized in Table 1.concentration. Figure 1c and d show simulated Greenland, east Antarctic and tropical SST in comparison with corresponding ice core records (NGRIP, 2004;Jouzel et al., 2007) and tropical SST stack (Herbert et al., 2010).The model correctly simulates both the magnitude and temporal dynamics of temperature variations associated with glacial cycles.The Greenland temperature also contains strong millennial scale variability associated with simulated reorganizations of the Atlantic thermohaline circulation (Ganopolski et al., 2010).

Baseline experiment
The magnitude of glacial-interglacial Greenland temperature changes is about 20 • C, which is close to the estimates for the last glacial maximum.The variations in Antarctic temperature are roughly twice as small as those in Greenland.
It is noteworthy that, although the temporal dynamics of Antarctic temperature variations are simulated rather realistically, the magnitude of warming during recent interglacials (MIS 5a, 7, 9 and 11) is underestimated.While in the model they are only slightly warmer than the present interglacial, the reconstructions suggest that previous interglacials were warmer by about 4 • C (Jouzel et al., 2007).Much warmer conditions in Antarctica during previous interglacials can be attributed to the disappearance of the West Antarctic ice sheet (Holden et al., 2010) and the initial temperature overshoot due to the mechanism of the bipolar seesaw (Ganopolski and Roche, 2009).The first mechanism is not taken into account in our simulations since we did not include changes in the Antarctic ice sheet.It is also noteworthy that simulated Antarctic interglacial temperatures prior to MIS 11 are in much better agreement with reconstructions.The agreement between simulated tropical SST with the stack in general is rather good.This is why the much stronger simulated cooling compared to empirical data during MIS 6 and 16 is difficult to explain.
Figure 2 shows that the model successfully simulates the waxing and waning of the ice sheets with a dominant 100 kyr periodicity and a pronounced asymmetry of the glacial cycles.For the second half of the run (Fig. 2d), modeling results agree favorably with reconstructed variations of global sea level by Waelbroeck et al. (2002).For the earlier part of the modeled period, reliable reconstructions of global ice volume are lacking, and the benthic δ 18 O c stack by Lisiecki and Raymo (2005) was used for comparison.Since benthic δ 18 O c is not an accurate proxy for ice volume, we computed the model's equivalent of δ 18 O c from the simulated global ice volume and deep ocean temperature using a simple relationship between the three (Duplessy et al., 1991).In addition, based on the results of simulations of the Antarctic ice sheet evolution during the last glacial cycle (Huybrechts, 2002), we assume that the Southern Hemisphere contributed an additional 10 % to global ice volume variations.Computed in this way, the modelled δ 18 O c agrees reasonably well with the empirical benthic δ 18 O c stack (Fig. 2c).It is also not surprising that the agreement is better for the second part of the model run than for the first one (prior to MIS 11).Indeed, the model used in this study was tuned for the last glacial cycles and, since several previous cycles resemble the last one in many respects, it is not surprising that the model results compare well with the proxy data also for these cycles.However, prior to MIS 11, the magnitude of δ 18 O c and CO 2 variations was smaller, which is likely due to difference in the Earth's topography and/or sediment distribution, which gradually evolved over geological time scales.Since these changes were not taken into account, degrading of the model agreement with the data further into the past is natural.The largest discrepancy between simulated and empirical δ 18 O c stack occurs during MIS 14 (ca.550 ka).Whether this is an indication of model deficiencies or dating problems of the proxy records is not clear.
The frequency spectra of modeled and empirical δ 18 O c are also in good agreement (Fig. 3a).All three major peaks -100, 41 and 23 kyr are reproduced with the dominance of 100 kyr cycle and a weaker precessional cycle, even though the modeled δ 18 O c contains more spectral power in the precessional band than the empirical spectrum.
Figure 2e shows the simulated volume of the North American and Eurasian ice sheets.Similarly to the last glacial cycle (Ganopolski et al., 2010), the global ice volume is dominated by the North American ice sheet, which agrees with Bintanja and van de Wal (2008).It also has stronger 100 kyr cyclicity whilst the Eurasian ice sheet reveals more variability at the precessional time scale.Simulated maxima in ice volume and ice area (not shown) are rather similar for both ice sheets for the LGM and penultimate glacial maximum.This is explained by the fact that the history and the magnitude of the orbital and GHG forcings prior to both glacial maxima were rather similar and, as a result, the simulated temperature over the NH continents at the end of MIS 6 and MIS 2 are rather similar as well.At the same time, paleoclimate reconstructions suggest (Svendsen et al., 2004) that the maximum extent of the Eurasian ice sheet in Russia and Siberia during the penultimate glaciation was significantly larger than at the LGM.This implies that our model underestimates the sensitivity of the Eurasian ice sheet area to the orbital forcings.In particular, it is possible that the background dust deposition taken from GCM simulations for the LGM (Mahowald et al., 2006) is too high over this area, which restricts the eastward advance of the Eurasian ice sheet, thus making its area less sensitive to climate forcing than it is in reality.In the future, we are planning to avoid this problem by using a fully interactive dust model instead of the prescribed GCM output (Bauer and Ganopolski, 2010).
As discussed by Ganopolski et al. (2010), both the magnitude and temporal dynamics of the simulated glacial cycles are very sensitive to the choice of model parameters, in particular, related to surface mass balance, basal sliding, and glaciogenic dust.It was also shown that without enhanced dust deposition, complete deglaciation is not achieved in the model.Figure 4b shows the simulated dust deposition rate over the southern margin of the North American ice sheet.Each glacial termination is associated with a large increase in dust deposition, which reduces surface albedo and enhances ablation.There are two main reasons for the increase in the dust deposition rate during glacial termination: (i) a considerable portion of the North American ice sheet at the glacial maxima spreads over the area covered by thick terrestrial sediments and (ii) most of the ice sheet base over this area is at the pressure melting point.Both of these factors enable fast sliding of the ice sheets and a large sediment transport towards the ice margins which, in turn, lead to enhanced glaciogenic dust production and dust deposition over the ice sheets.This amplifies the direct effect of rising summer insolation and GHGs concentration.
The realistic simulation of climate-cryosphere dynamics over past 800 kyr in the BE performed with prescribed orbital and GHG forcings represents an important test for the model.However, such an experiment does not answer the question about the origin of the strong 100 kyr cycles, since it is possible that the dominant 100 kyr periodicity and the correct timing of glacial terminations are solely attributed to the prescribed GHG forcing, the temporal dynamics of which strongly resembles the ice volume.

Experiments with constant CO 2
To clarify whether the 100 kyr cycles directly originate from the orbital forcing, we performed a set of additional experiments (referred to hereafter as CCn, where n is the prescribed CO 2 concentration in ppm, also see Table 1) with the same orbital forcing as in the BE described above but maintaining a constant CO 2 concentration in time.We performed five experiments with the CO 2 concentration ranging from 200 to 280 ppm (for every 20 ppm). Figure 5a shows a representative subset of these simulations, while Fig. 6a shows the results of all CCn experiments for the whole range of CO 2 concentrations from 200 to 280 ppm.Quasi-regular glacial cycles are simulated for all CO 2 concentrations, with the magnitude of the ice volume variations increasing for decreasing CO 2 .For a CO 2 concentration of 280 ppm, simulated glacial cycles are dominated by obliquity and precession, but for lower CO 2 concentrations, the model simulates long and asymmetric glacial cycles with a strong peak in the 100 kyr band in the frequency spectra (Figs.6a and 3c).It is noteworthy that, together with 100 kyr peak another one, although weaker, appears in the 400 kyr band.This peak is also seen in the results of previous simulations of the glacial cycles and coincides with another eccentricity period.Therefore, the presence of both 100 and 400 kyr periodicities strongly www.clim-past.net/7/1415/2011/Clim.Past, 7, 1415-1425, 2011 indicates a direct relationship of the long glacial cycles with eccentricity variations.Unlike (Crowley and Hyde, 2008) who found the existence of 100 kyr cycles in a narrow range of CO 2 concentrations, in our simulations 100 kyr cycles are robust over a broad range of CO 2 concentrations (200-260 ppm).
It is important to note that in the simulations with a constant CO 2 concentration below 260 ppm, not only the ice volume changes are dominated by the 100 kyr cycles, but simulated glacial terminations also occur at the same time as in the experiment with prescribed time-dependent CO 2 and, within the dating accuracy, are in good agreement with paleoclimate reconstructions.The only exception is for MIS 11 (around 400 kyr BP), when complete deglaciation of the Northern Hemisphere does not occur in the experiments with low CO 2 concentrations.The latter fact is not surprising since the orbital forcing was weak during MIS 11 due to low eccentricity.Whether this problem implies that for this specific termination the role of CO 2 is more important than for the others or that the model is still not sufficiently non-linear to remove the ice under MIS 11 orbital forcing cannot be answered within the context of this study.
Although 860 kyr BP represents a convenient time to start simulations of the last glacial cycles, since it corresponds to an interglacial state and therefore does not require initialisation of the continental ice sheets, it is theoretically possible that this choice can be crucial for the timing of simulated glacial terminations and therefore a good agreement between simulated and real glacial terminations would be accidental.To show that this is not the case and that the timing of the glacial terminations is solely controlled by the orbital forcing, we performed an additional set of model simulations for constant CO 2 concentration equal to 200 ppm, where we began the model runs at the different astronomical times between 800 and 900 kyr BP using the same (interglacial) initial conditions.These experiments are referred to as CC200/m, where m denotes the timing of the start of the experiment.Figure 6b shows that in all CC200/m experiments the simulated ice volume converged within one glacial cycle to the same solution.Therefore, the temporal dynamics and timings of glacial terminations after the model spin-up are not sensitive to the choice of the beginning of the models runs.
As was shown by Ganopolski et al. (2010), simulation of a complete glacial termination, even with the prescribed GHG forcing, is only possible when the deposition of glaciogenic dust is taken into account.This is even more important in the case of a constant CO 2 concentration.Without glaciogenic dust, long glacial cycles are not simulated in our model.Similarly to the BE, in the experiments with constant CO 2 (Fig. 4c), the rate of dust deposition over the southern flank of the North American ice sheet increases significantly during glacial terminations, which leads to enhanced ablation and facilitates the ice sheet's retreat.Similarly to the last glacial cycle (see Fig. 11c in Ganopolski et al., 2010), simulated glacial cycles over the entire 800 kyr period are rather sensitive to the poorly constrained parameters of the dust deposition scheme.Nonetheless, Fig. 6c shows that both doubling and halving of the glaciogenic dust deposition rate does not change the dominant 100 kyr cyclicity or the timing of major glacial terminations.Therefore, although details of the simulated glacial cycles depend on the dust parameterization, the simulated 100 kyr cyclicity is robust with respect to the strength of the dust feedback.

Sensitivity of glacial cycles to different components of the orbital forcing
To find which component of the orbital forcing is responsible for the existence of 100 kyr cyclicity, we performed a suite of additional experiments in which, similar to the CCn set, the CO 2 concentration was held constant but the orbital forcing was modified.In the first set of experiments (referred as COBn), we removed the effect of obliquity variations by setting obliquity constant in time and equal to its average value over the simulated period (Fig. 7b).As shown in Fig. 5b, the removal of obliquity variations does not qualitatively affect the simulated glacial cycles.For sufficiently low CO 2 concentrations, long glacial cycles with a sharp maximum in the frequency spectra at 100 kyr periodicity are simulated (Fig. 3d).However, fixing obliquity results in a narrower range of CO 2 concentrations for which 100 kyr cyclicity dominates the frequency spectra.In addition, fixing obliquity makes glacial terminations less robust for the periods of low eccentricity, particularly during the most recent termination.At the same time, all terminations that occurred during high eccentricity are correctly simulated in the COBn experiments.
In the complimentary set of experiments (referred to as CECn/e), we modified the orbital forcing by setting eccentricity constant in time with a value in the range 0-0.05, which covers real variations of eccentricity (Fig. 7c).Variations of obliquity and precession were the same as in reality.For eccentricity lower than 0.02, the model failed to simulate pronounced glacial cycles and for eccentricity of 0.04 and higher, the ice volume variations are dominated by precession (not shown).However, for intermediate values of eccentricity (0.02 and 0.03) and sufficiently low CO 2 concentrations, the model simulates long glacial cycles (Fig. 5c).Moreover, many (but not all) simulated glacial terminations occur in the CEC220/0.02experiments at the right time.However, the frequency spectra of ice volume variations in CECn/e experiments lack a sharp peak in the 100 kyr band (Fig. 3e).In fact, the shape of the frequency spectrum in these experiments is very sensitive to the CO 2 level and the maximum of spectral power tends to occur at one or multiples of obliquity periods, rather than at 100 kyr.Therefore, in our experiments obliquity variations themselves are not responsible for the existence of the 100 kyr cycles, but they do contribute to the robustness of the long glacial cycles.It is noteworthy that the experiments described above essentially repeat the work by Pollard (1983).Although the principal mechanism of glacial terminations is different in our and in Pollard's models, the main results are very similar.

The role of regolith
To gain insight into the possible origin of the mid-Pleistocene transition (around 1 million years ago), when the dominant periodicity of the glacial cycles changed from 40 kyr to 100 kyr (Ruddiman et al., 1989), we performed a set of model experiments identical to the CCn set, except for the spatial distribution of terrestrial sediments.We prescribed namely the presence of a thick terrestrial sediment layer for all continental grid cells, while in the previous experiments we used a realistic distribution of terrestrial sediments.This set of experiments is referred to as REGn.With the continents completely covered by sediments, the 100 kyr cyclicity is absent in the simulated ice volume for the whole range of prescribed CO 2 concentrations and the ice volume variations are of a smaller magnitude than in the CCn experiments (Fig. 5d).This result lends support to the hypothesis that the gradual removal of the sediments from the large area of North America stabilized the North American ice sheet and led to the regime change in glacial cycles (Clark and Pollard, 1998).It is important to note that the presence of terrestrial sediments in our model has a dual effect on the ice sheets: (i) it enhances the velocity of ice sheet sliding in the areas where the ice base is at the pressure melting point, and (ii) it increases production of glaciogenic dust around the margins of the ice sheets (Fig. 4d), which affects surface albedo and facilitates surface melt (Ganopolski et al., 2010).Both factors affect the stability of the ice sheets by making them more sensitive to changes in orbital forcing.Only when a large area of North America is free of sediments, can the ice sheet survive several precessional cycles before it expands well into the area covered by the sediments, making the ice sheet more sensitive to summer insolation changes.
Glacial cycles simulated in the REG experiments have, for the same CO 2 concentrations, considerably smaller amplitude than in the CC experiments.Even for a CO 2 concentration of 200 ppm, which is likely near the lower limit of CO 2 concentrations during the 40-kyr world (Hönisch et al., 2009), the maximum simulated NH ice volume for most of glacial cycles is within the range 30-40 m of sea level.This is significantly smaller then simulated by Bintanja and van de Wal (2008) but more in line with Siddall et al. (2010).Such a small magnitude of ice volume variations in the REG experiments has an important implication for the interpretation of benthic δ 18 O records, since it implies a much smaller relative contribution of the NH ice volume to the glacial-interglacial δ 18 O variations for the 40-kyr world compared to the 100kyr world.At the same time, the maximum southward extent of the Laurentide ice sheet in the REG200 experiment (Fig. 8) is similar to that simulated in the BE and CC200 experiments.This is consistent with paleoclimate evidence of a similar southward extent of the Laurentide ice sheet in the 40-and 100-kyr worlds reviewed by Clark et al. (2006).The large difference in volume for ice sheets with a similar area is explained by the fact that in the REG experiments, due to strong sliding over a broader sediment-covered area, the maximum thickness of the central part of the Laurentide ice sheet is much smaller than at the end of the long glacial cycles simulated in BE or CC200 experiments.

Discussion
The results presented above demonstrate that our model simulates long asymmetric glacial cycles both with realistic (time-dependent) and constant CO 2 concentrations for a broad range of CO 2 concentrations and model parameters.The simulated long glacial cycles can be understood in the framework of the stability diagram of the climate-cryosphere system in the phase space of the orbital forcing (Calov and Ganopolski, 2005).According to this analysis, the interglacial state becomes unstable when summer insolation falls below a certain threshold value and the system enters the glacial state.Although the ice sheets experience both growth and decay during different phases of orbital forcing, the system remains in the domain of attraction of the glacial state most of time and, therefore, the ice volume has a general tendency to grow.Such a stability diagram, however, does not explain why glacial cycles are rapidly terminated ca. 100 kyr after their start.The explanation of glacial termination requires an additional strong nonlinear mechanism, which, in our case, is the dust feedback.This feedback is activated after the ice sheets spread well into the area covered by thick terrestrial sediments.High rates of dust deposition over the ice sheets reduce their albedo, which enhances ablation and thus amplifies the ice sheet response to rising insolation.Note that this mechanism was already proposed by Peltier and Marshall (1995), but they treat this in a rather simplistic manner, namely, they prescribed two to three times higher ablation rates during termination only.
This mechanism related to dust feedback also explains the difference between 100-kyr and 40-kyr worlds.If the Northern Hemisphere continents are entirely covered by thick sediments, the strong dust feedback already starts to operate with the first significant rise of insolation, and glacial cycles are terminated each 20-40 kyr, depending on magnitude of eccentricity and the relative phase of precession and obliquity.
If large portions of the northern continents are cleared from terrestrial sediments, as is the case at present, the ice sheets can survive several insolation maxima before they spread into the area covered by the sediments and the mechanism described above is activated.This usually occurs after the minimum of eccentricity is reached, but the precise timing of glacial termination depends on the phase between obliquity and precession.For this reason, the timing of glacial terminations relative to the eccentricity minima can vary in a broad range (Huybers and Wunsch, 2005).However, even this rather general tendency of the NH ice sheets to grow during periods of low eccentricity and to terminate during the rise of eccentricity is sufficient to produce a sharp peak in the frequency spectrum of ice volume in the 100 kyr band.
Although in our experiments the effects of fast ice sheet sliding and dust deposition are sufficient to terminate long glacial cycles, it is likely that other nonlinear processes also play an important role.Several such mechanisms have been proposed: accelerated flow and wastage at the southern tip of the ice sheet due to proglacial lakes and/or marine incursions (Pollard, 1983), large-scale instabilities of marine ice sheets (Peltier and Marshall, 1995).Another factor that contributes to the termination of glacial cycles is the rise of CO 2 during glacial terminations (Wolff et al., 2009).However, the latter is rather a positive feedback in the Earth system than a cause of glacial termination.
Although the mechanisms of 100 kyr cyclicity presented here invoke the concept of phase locking, it is fundamentally different from the mechanism of phase locking of internal oscillations in the climate-cryosphere or climate-cryospherecarbonosphere system to the orbital forcing (e.g.Saltzman and Maasch, 1988;Paillard and Parrenin, 2004;Tziperman et al., 2006).In the latter case, glacial cycles with a typical periodicity of 100 kyr exist in the system even without orbital forcing (i.e. with constant orbital parameters).The orbital forcing only sets the precise timing of glacial terminations.In our case, glacial cycles are the response of the climatecryosphere system to the orbital forcing and both 100 kyr cyclicity and the timing of terminations are set by the orbital forcing.If orbital parameters are fixed, our model does not simulate glacial cycles.More detailed comparison of these two concepts of glacial cycles will be given in forthcoming work.

Conclusions
Results of our experiments support the notion that 100 kyr cycles represent a direct, strongly nonlinear response of the climate-cryosphere system to orbital forcing and they are directly related to the corresponding eccentricity period.In terms of nonlinear dynamics, this link can be interpreted as the phase-locking of the long glacial cycles to the shortest (100 kyr) eccentricity cycles.Physically, this phase-locking is explained by the fact that the ice sheets tend to grow monotonously during periods of low eccentricity and reach their critical size (volume) around the minimum of eccentricity.When eccentricity starts to grow, the first sufficiently large positive anomaly in orbital forcing can lead to the rapid and irreversible meltback of the Northern Hemisphere ice sheets.This mechanism requires the existence of long glacial cycles that, in turn, require sufficiently low CO 2 concentrations and a large area of the continents to be free of sediment.The CO 2 concentration not only determines the dominant regime of glacial variability, but also strongly amplifies 100 kyr cycles.Therefore, realistic simulations of the glacial cycles require comprehensive Earth system models that include both physical and bio-geochemical components of the Earth system.

Figure
Figure 1a-b show two principal forcings -orbital (illustrated by the maximum summer insolation at 65 • N) and the global mean radiative forcing of prescribed equivalent CO 2

Fig. 4 .
Fig. 4. (a) Maximum summer insolation at 65 • N; (b) simulated volume of the North American ice sheet in msl (black line) and dust deposition rate (orange) at the southern margin of the ice sheet in the Baseline experiment, (c) the same as (b) but for the CC200 experiment, (d) the same as (b) but for REG200 experiment.The dust deposition rate is sampled at the location (45 • N, 100 • E).

Fig. 5 .
Fig. 5. Simulated ice volume variations in a subset of the experiments with constant CO 2 (colored lines) versus the Baseline experiment (grey shading).(a) CCn experiments, (b) COBn experiments (constant obliquity), (c) CECn/0.02experiments (constant eccentricity equal to 0.02), (d) REGn experiments (with continents completely covered by thick sediment layer).Black lines correspond to a CO 2 concentration of 200 ppm, blue -220 ppm and red -280 ppm.Dashed line at the top shows eccentricity.

Fig. 7 .
Fig. 7. Maximum summer insolation at 65 • N in (a) the BE, CCn and REGn experiments, (b) the experiment with constant obliquity (COBn) and (c) constant eccentricity (e = 0.02, CECn/0.02).The dashed line in panel (b) shows the amplitude modulation of the precessional cycle by eccentricity.

Fig. 8 .
Fig. 8. Thickness of the simulated ice sheets when the maximum area of the entire simulation was achieved in experiment CC200 (a) and REG200 (b), respectively.

Table 1 .
List of model experiments.