Mechanisms and time scales of glacial inception simulated with an Earth system model of intermediate complexity

Abstract. We investigate glacial inception and glacial thresholds in the climate-cryosphere system utilising the Earth system model of intermediate complexity CLIMBER-2, which includes modules for atmosphere, terrestrial vegetation, ocean and interactive ice sheets. The latter are described by the three-dimensional polythermal ice-sheet model SICOPOLIS. A bifurcation which represents glacial inception is analysed with two different model setups: one setup with dynamical ice-sheet model and another setup without it. The respective glacial thresholds differ in terms of maximum boreal summer insolation at 65° N (hereafter referred as Milankovitch forcing (MF)). The glacial threshold of the configuration without ice-sheet dynamics corresponds to a much lower value of MF compared to the full model. If MF attains values only slightly below the aforementioned threshold there is fast transient response. Depending on the value of MF relative to the glacial threshold, the transient response time of inland-ice volume in the model configuration with ice-sheet dynamics ranges from 10 000 to 100 000 years. Due to these long response times, a glacial threshold obtained in an equilibrium simulation is not directly applicable to the transient response of the climate-cryosphere system to time-dependent orbital forcing. It is demonstrated that in transient simulations just crossing of the glacial threshold does not imply large-scale glaciation of the Northern Hemisphere. We found that in transient simulations MF has to drop well below the glacial threshold determined in an equilibrium simulation to initiate glacial inception. Finally, we show that the asynchronous coupling between climate and inland-ice components allows one sufficient realistic simulation of glacial inception and, at the same time, a considerable reduction of computational costs.


Introduction
Simulation of glacial inception with numerical models different in complexity and detail of description of climate components has already a long tradition. An early approach has been undertaken by Weertman (1964) with an ice-sheet model based on the assumption of perfect plasticity in ice flow. Andrews and Mahaffy (1976) utilised an, at that time, relatively sophisticated two dimensional longitude-latitude vertically integrated ice-sheet model to simulate the Laurentide ice sheet dynamics. Common in these approaches is an extremely simplified representation of the surface mass balance by mainly prescribing the movement of a so-called equilibrium line (height of zero mass balance). There are many attempts to model glacial cycles, in particular glacial inception, using conceptional models (Paillard, 2001;Saltzman, 2002) which try to describe the entire climate system with a minimum set of ordinary differential equations. The advantage of these models is that they are computational very efficient. Due to large computational costs, simulations with GCMs (general circulation models of the atmosphere) so far were restricted to time slices experiments (e.g. Dong and Valdes, 1995;Yoshimori et al., 2002;Vettoretti and Peltier, 2004) or, in case of transient simulations, tied up to relative short integration time or asynchronous coupling (Pollard and Thompson, 1997;Vizcaíno et al., 2008). Earth system models of intermediate complexity (EMICS;Claussen et al., 2002) are a compromise, because of their sufficient grade of Published by Copernicus Publications on behalf of the European Geosciences Union. complexity paired with adequate detail in description of the climate subsystems. These models allow full synchronous transient simulations over the time of decreasing boreal summer insolation during glacial inceptions including all relevant components of the climate system. In contrast to time slice simulations, the transient approach captures adequately the inertia of the climate-cryosphere system, which will shown to be important in this paper.
One of the first EMICs including ice-sheet dynamics was presented by Gallée et al. (1991). Wang and Mysak (2002) published a simulation of glacial inception with a geographical explicit EMIC including interactive ice sheets (the McGill Paleoclimate Model). A number of feedbacks and mechanism in the context of glacial inception have been investigated with EMICs. One common results of different simulations is that the terrestrial vegetation amplifies glacial inception (Kageyama et al., 2004;Wang et al., 2005;Kubatzki et al. 2006), mainly due the southward retreat of boreal forest (Calov et al., 2005b). A further, even more important process is the reduction of SSTs and expansion of sea ice, which lead to an additional strengthening of glacial inception Calov et al., 2005b). Dust deposition on the snow lowers the albedo of snow considerably (Calov et al., 2005a). This effect plays an important role for the ice cover in Alaska and north-east Eurasia, regions which mainly remained ice-free during the Quaternary. The importance of dust deposition for the extent of the ice sheets during last glacial inception was demonstrated by Calov et al. (2005b). GCM simulations which treated dust in snow-albedo parameterisation resulted in a strong increase of snow free days per year over Alaska and Northern Asia (Krinner et al., 2006), what supports the results yielded with an EMIC by Calov et al. (2005b). The disregarding of dust deposition could be the reason why many other modellers have a too expanded ice cover over Alaska and north-east Eurasia. Calov et al. (2005a) showed that glacial inception is primarily triggered by a drop in boreal summer insolation and is linked to a bifurcation from interglacial to glacial state caused by snow-albedo feedback. They explicitly demonstrated the existence of two equilibria for the same forcing: one interglacial and one glacial equilibrium (bistability). Finally, Calov and Ganopolski (2005) found hysteresis in the climate-cryosphere system by tracing the space of Milankovitch forcing. In this paper, we closer inspect glacial thresholds for model configurations with and without icesheet dynamics.
By Ruddiman (2003) and in further publications on this topic (e.g. Ruddiman, 2007) it was postulated that early Holocene emission of greenhouse gases (GHGs) prevented late Holocene glacial inception. Claussen et al. (2005) commented on Ruddimans's hypothesis by comparing a transient simulation through last glacial inception with a Holocene simulation where Ruddiman's postulated natural drop in atmospheric CO 2 applied. Here, we are widening Claussen's approach by treating Ruddiman's hypothesised drop in atmo-spheric CH 4 . We mimic the CH 4 radiative effect with an additional drop in prescribed atmospheric CO 2 . But our results do still not support Ruddiman's hypothesis in the sense that even an appreciable lowering of CO 2 does not lead to noticeable glaciation of the Northern Hemisphere during the late Holocene. The current minimum MF is close to the glacial threshold and lowering of CO 2 does affect this threshold. But still MF does not drop far enough below the glacial threshold to cause a late Holocene growth of the ice sheets over the continents of the Northern Hemisphere.
Further on, our analysis will support efforts to model glacial inception with state-of-the-art Earth system models. Due to high computation cost, GCMs so far can afford only rather short integration times (some 1000 model years). Therefore, the use of asynchronous coupling is the only way to perform transient simulations of the whole glacial inception. Below, we investigate the influence of different asynchronous coupling strategies on the accuracy of transient simulation of glacial inception.
The paper starts with the model description (Sect. 2) and an analysis of glacial thresholds (Sect. 3). The implication for the transient cryosphere-climate system if being subject to a perturbation distant from the thresholds will be inspected in Sect. 4. The investigation of transient response during different glacial inceptions with a particular link to Ruddiman's hypothesis is included in Sect. 5. An analysis of the temporal coupling properties of the climate and ice-sheet modules will follow in Sect. 6. We close our paper with a discussion of our findings (Sect. 7) and conclusions (Sect. 8).

The model
The Earth system model of intermediate complexity CLIMBER-2 consists of an atmospheric part based on the statistical-dynamical approach (Petoukhov et al., 2000). The atmospheric component has a resolution of 10 • in latitude and about 51 • in longitude. The atmosphere is coupled to a three-basin ocean model (similar to Stocker et al., 1992), which is run on a latitudinal resolution of 2.5 • and resolves 20 unequally thick layers in the vertical. The three basins are linked in the circumpolar oceans. A thermodynamic sea-ice model is included. The dynamical global vegetation model VECODE (Brovkin et al., 2002) describes potential fractions of tree, grass, and bare soil coverage. Land-surface processes are simulated on the spatial resolution of the atmospheric model. The original CLIMBER-2 model has been applied to many types of climate studies, including the simulation of different palaeoclimates (e.g. Ganopolski et al., 1998a, b;Brovkin et al., 1999;Claussen et al., 1999;Rahmstorf and Ganopolski, 1999;Kubatzki et al., 2000;Ganopolski et al., 2001).
In this paper, we utilise the upgraded model version of CLIMBER-2 (Calov et al., 2005a), which includes the three-dimensional polythermal ice-sheet model SICOPOLIS (Greve, 1997). The domain of SICOPOLIS is restricted to the Northern Hemisphere extra-tropics in the given application. SICOPOLIS is run on a resolution of 0.75 • in latitude and 1.50 • in longitude. It treats inland-ice extent, thickness, temperature, velocity, and water content. An interface module has been developed for direct bi-directional coupling between SICOPOLIS and the coarse-resolution climate component of CLIMBER-2 described above. This coupler, SEMI, provides SICOPOLIS with annual values of surface temperature and mass balance. SEMI interpolates the characteristics it needs to determine the snow cover from the coarse atmospheric grid to the fine grid of the ice-sheet model taking into account the high-resolution orography. It calculates snow thickness and surface temperature from this high-resolution information using the same equations that are used in the coarse-grid land-surface model. SEMI determines the annual surface mass balance of the ice sheets and passes it together with the surface temperature to the ice-sheet model. The inland-ice model returns the updated values of inlandice fraction, surface elevation, and sea-level back to the climate model. For the calculation of snow albedo, variations of the deposition of mineral dust were taken into account. For details see Calov et al. (2005a). In this work, we also use a model configuration without the ice-sheet model, but with the area of inland ice diagnosed by SEMI on the spatial grid of the ice-sheet model.
CLIMBER-2 including the ice-sheet model has been utilised for the simulation of Heinrich events  and transient simulations of last glacial inception (Calov et al., 2005a, b;Kubatzki et al., 2006). Calov and Ganopolski (2005) used this model to perform a stability analysis of the climate-cryosphere system in the phase space of MF, whilst Claussen et al. (2005) used this model to evaluate a hypothesis by Ruddiman (2003) on the possible prevention of glacial inception during the late Holocene by early anthropogenic emission of CO 2 and CH 4 .

Stability analysis
In order to better understand the mechanisms of glacial inception, we performed a stability analysis of the climatecryosphere system in the phase space of MF similar to that carried out by Calov and Ganopolski (2005). We use two different model settings for our stability analysis. In the first setting, we use the full CLIMBER-2 model including the icesheet model (ice-sheet dynamics "switched on" referred to as ID-on configurations hereafter). In the second setting, the module which describes the ice-sheet dynamics is "switched off" (ID-off configuration). In the ID-off configuration, the surface energy and mass balance interface SEMI is still enabled and all land grid points where SEMI simulates positive surface mass balance are treated as permanent snow cover area. The area of permanent snow cover is then averaged over each grid cell of the coarse resolution climate model and fed back to the climate model as a time-dependent area of "inland ice", one of six surface types considered in the CLIMBER-2 model. This surface type is characterised by high surface albedo equal to the maximum albedo of snow. Therefore, in the ID-off configuration the positive ice-albedo feedback operates the same way as in the ID-on configuration, but additional feedbacks associated with ice sheet lateral spreading and elevation changes are ignored. The ID-off configuration is intended to emulate simulations of glacial inception with an AGCM (atmosphere only GCM) or an AOGCM (GCM of atmosphere and ocean) in absence of an ice-sheet model: a situation which often appeared for time-slice simulation of glacial inception performed in the past (Syktus et al., 1994;Dong and Valdes, 1995;Gallimore and Kutzbach, 1996;de Noblet et al., 1996;Kageyama et al., 1999;Khodri et al., 2001;Yoshimori et al., 2002;Meissner et al., 2003;Vettoretti and Peltier, 2004;Kaspar et al., 2007). For both configurations, we use equilibrium interglacial climate as initial condition. While in the ID-on configuration the Greenland ice sheet is allowed to evolve freely, the Greenland ice topography in the ID-off configuration is fixed to its present-day value, because in the latter case no ice movement is possible.
To trace the equilibrium states of the climate-cryosphere system, we slowly change MF for the whole range corresponding to the glacial-interglacial transition between 126 and 115 kyr BP. We used here the method identical to that in Calov and Ganopolski (2005) except that we now vary MF by changing the precessional angle for fixed (large) eccentricity, while in Calov and Ganopolski (2005) we changed the eccentricity for two different precessional angles. As it was noted in Calov and Ganopolski (2005), a slow change in precessional angle leads to a rather similar stability diagram in MF space as slow variation of eccentricity. We fixed the eccentricity and the obliquity of the Earth's orbital parameters to their respective values at 115 kyr BP of e=0.04142 and ε=22.4 • , and the precessional angle (angle between vernal equinox and perihelion) was varied in negative direction from 90 • to −90 • . Such a variation in the precessional angle corresponds to a lowering in MF, which causes a change from warm to cold Northern Hemisphere summer conditions. The atmospheric CO 2 concentration is set to 271 ppm, which corresponds to its typical interglacial value prior to the onset of glacial inception. The mineral dust deposition rate is set to its present day value according to Mahowald et al. (1999).
The precessional angle was changed such that it resulted in a linear change in MF. In the ID-on configuration, MF decreases with a rate of 1.6×10 −4 W m −2 yr −1 outside of the vicinity of the bifurcation transition from interglacial to glacial state. This rate was additionally reduced by a factor of 1/20 in the vicinity of the bifurcation transition, because the climate-cryosphere system needs longer time to equilibrate near the bifurcation point. The model in the ID-on configuration has been run for about 1.5 million model years to trace the stability diagram, whilst for the ID-off experiment we used an order of magnitude larger rate of MF, because in absence of ice-sheet dynamics the model has a much shorter equilibration rate. Results of the stability analysis are shown in Fig. 1. Note that compared to Calov and Ganopolski (2005) we study here only that branch of the stability diagram, which corresponds to the trajectory from high to low MF. It has been shown in Calov and Ganopolski (2005) that with a gradual increase of MF the equilibrium trajectory does not follow the same path in MF space and that a pronounced hysteresis exists. However, here we are only concerned with the glacial inception and do not analyse the other branch of the stability diagram.
In addition to the method of slowly varying precessional angle, we performed a number of equilibrium simulations with fixed orbital parameters for certain values in MF ( Fig. 1). In the ID-on configuration, these equilibrium simulations were performed over 500 000 model years, while in the ID-off configuration the model was run for 5000 years, which in both cases is sufficient to attain equilibrium states with high accuracy. These equilibrium simulations were performed to assess the robustness of our stability diagram obtained by gradual changes in MF, because it is not clear a priori how slow MF should be varied to trace the equilibrium state of the climate-cryosphere system in the vicinity of bifurcation transitions. Near the glacial thresholds the chosen MF values of the equilibrium simulations differed only by 1 W m −2 . If of two neighbouring equilibrium sim-ulations one reached glacial state and the other one stayed in interglacial state, then the MF of the simulation which stayed in interglacial state was taken as the threshold value of MF for the respective model configuration (see the grey bars in Fig. 1). In Sect. 5, this technique is utilised to determine the dependence of the ID-on glacial threshold on atmospheric CO 2 concentration using CO 2 values different from the 271 ppm which applies here.
For the ID-on configuration, it can be seen that near the bifurcation the MF of the equilibrium simulations does not fully lie on the curve from the method with slowly varying parameter, which means that even such low rate of change of MF in the vicinity of the bifurcation was still not sufficiently low to trace the bifurcation point precisely. Although our bifurcation analysis is worked out in the MF phase space, our model computes the full orbital forcing in season and in latitude. The particular MF is chosen here, because our previous work showed that MF is a good proxy for orbital forcing in the context of interglacial to glacial changes (Calov et al., 2005a).
The glacial thresholds for the two model configurations (ice-sheet dynamics "on" and "off") are well separated in MF ( Fig. 1). For the ID-on model configuration, the critical value of MF corresponding to transition from interglacial to glacial state is S ID on =477 W m −2 , which is only slightly below the present-day value (about 480 W m −2 ). The glacial threshold (by which we mean the appearance of large areas with positive mass balance over the continents of the Northern Hemisphere) in the ID-off configuration is S ID off =442 W m −2 , far below S ID on and about to be as small as the minimum of the MF during last glacial inception (Fig. 5a). Obviously, the different positions of the glacial thresholds in MF space arise from the fewer feedback mechanism in the ID-off configuration compared to the ID-on configuration. Mainly, there are three types of mechanism which originate from the cryosphere: ice/snow-albedo feedback, surface-elevation mass-balance feedback and expansion of inland ice through ice flow (ice spreading), of which only albedo feedback operates in ID-off configuration. In addition, the area of positive surface mass balance in the ID-off simulation is notably smaller in comparison with the inland-ice area in the ID-on experiment for the same MF, which is again explained by missing feedbacks related to ice-sheet dynamics.

Transient response to instantaneous change in orbital forcing
Although useful for understanding of a nonlinear system, a stability diagram is not directly applicable to the transient response of the system. To investigate the time scale of transient response of the climate-cryosphere system to timedependent orbital forcing, we performed a set of experiments where we studied the model response (in ID-on configuration) to instantaneous changes of the orbital forcing. We selected two sets of orbital parameters. One set yields a moderate MF equal to 460 W m 2 , which represents cold (compared to present-day orbital configuration) boreal summer. Another set of orbital parameters results in a low MF equal to 430 W m −2 , which is due to extreme cold boreal summer conditions. The value of moderate MF is in between S ID off and S ID on , whilst the low MF is well below both glaciation thresholds. The moderate MF was obtained by prescribing an eccentricity of e=0.04142 and a precessional angle of ω=−32 • . For the low MF we used "cold" orbital configuration (i.e. precessional angle of ω=−90 • ) and, in addition, a high eccentricity of e=0.05405. As initial conditions we chose an equilibrium interglacial climate state (i.e. with Greenland ice sheet only) corresponding to present-day MF (i.e. 479.6 W m −2 ) and an interglacial CO 2 concentration of 271 ppm. Then we applied each of two orbital configurations described above and run the model for 100 000 years keeping CO 2 and orbital parameters fixed, which is equivalent to instantaneous change of MF from the present one to "medium" or "low" value. The time series of ice-sheet area, depicted in Fig. 2, expose an important difference in model response between moderate and low MF. While there is a fast response in inland-ice area with a 13 million km 2 increase during the first 100 years for low MF, there is no visible increase in inland-ice area for moderate MF (Fig. 2a). The fast response to low MF cannot be related to ice dynamics. A simulation without ice-sheet dynamics (ID-off) under the low MF shows about the same increase in area of positive mass balance as the corresponding simulation with ice dynamics (not shown in the figure). Therefore, snow-albedo feedback has to play a major role for the initially fast increase in inland-ice area. While the area and volume of the ice sheets increases very fast under low MF and approach their equilibrium values already after 20 000 years of integration, under moderate MF only rather small ice sheets are formed during this period of time and only after 20 000 years of integration the growth of the ice sheets accelerates considerably (Figs. 2 and 3). Note that the equilibrium state does not differ much between these two MF values.
The main differences in the initial phase of ice-sheet growth between moderate and low MF are illustrated with an ice-thickness plot (Fig. 3). While there are only rather small areas in the high north glaciated after 100 years in the simulation with moderate MF, a thin ice field with a thickness of about 10 m covers about half of the North American continent after the same model time under low MF (Fig. 3a, d). After 10 000 model years, there are small and relative thick ice caps in the high north for moderate MF, while there are about fully developed large ice sheets in North America and Northern Eurasia for low MF (Fig. 3b, e). At the end of the simulations, the ice thickness distribution for moderate and low MF is similar showing two big ice sheets in North America and Northern Eurasia (Fig. 3c, f).
In the following, we would like to further inspect the phases of ice build up. For moderate MF, one can detect three phases in ice-sheet build up: an initial ice-spreading phase (or spreading-only phase), when due to rather small sized ice caps, the ice-albedo feedback plays little role, a terminal relaxation phase and a transitional phase in between. For both, the initial und the terminal phase we developed a simple model, which we compare with our CLIMBER-2 simulations in the ID-on model configuration. The ice-spreading model (Appendix A) assumes a small quasi-equilibrium axisymmetric ice sheet, which grows under constant positive surface mass balance. The even simpler relaxation model (Appendix B) just assumes relaxation to an equilibrium inlandice volume with a time constant in an exponential function. A pronounced spreading phase during the initial phase of ice build-up can only be seen in the inland-ice volume curve with moderate MF, higher than S ID off , because for lower values of MF, the ice-albedo feedback plays an important role from the beginning. To test the robustness of our conceptual models, we performed two additional runs with MF equal 455 and 465 W m −2 . The spreading and relaxation phases of the three equilibrium simulations are shown in Fig. 4. It can be seen that in all three simulations the second derivative of the inland-ice-volume curves changes its sign from positive (concave shape) for the spreading-only phase to negative (convex shape) for the relaxation phase (Fig. 4). The concave shape of the volume curve is reproduced by the power law of Eq. (A9) and the convex shape by an exponential relation   Figure 4 shows that the ice spreading model is applicable only for rather small ice volumes, whilst the relaxation model describes quite satisfactorily the inland-ice-volume evolution after it reaches about half of its equilibrium value.

Transient response to realistic orbital forcing
In the following, we investigate the transient response of the climate-cryosphere system to realistic time-dependent orbital forcing (Berger, 1978) corresponding to the marine isotopic stages 1, 5 and 11 (MIS 1, 5 and 11). The simulations start at times of maximum MF, pass one minimum in MF and go on towards the next maximum in MF (Table 1, Fig. 5a). For the chosen time interval, we performed simulations with two different values in atmospheric CO 2 concentration. The CO 2 values were fixed to 280 and 220 ppm, respectively, throughout the entire simulations. In all experiments we used interglacial equilibrium state as initial conditions. Figure 5, displays time series of MF and inland-ice volume simulated under MIS 1, 5 and 11 orbital forcings. The MF during MIS 1 has a rather weak minimum, the MF minimum during MIS 5 is much stronger and the MIS 11 insolation minimum is approximately in between the MIS 1 and MIS 5 minima (Fig. 5a). The glacial thresholds for 220 and 280 ppm atmospheric CO 2 concentration shown in Fig. 5a are determined using the same technique as described in Sect. 3. It can be seen that these two glacial thresholds of S ID on (220)=486 W m −2 and S ID on (280)=476 W m −2 differ only by 10 W m −2 . Such a dependence of MF on atmospheric CO 2 concentration agrees well with results by Archer and Ganopolski (2005).
The transient simulations show pronounced glacial inception during MIS 5 for both atmospheric CO 2 concentrations which were applied here. Glacial inception is indicated by a strong increase in inland-ice area (Fig. 5b). In the transient MIS 5 simulation with 280 ppm atmospheric CO 2 concentration, glacial inception appears about 2500 years later and with about 40% less maximum inland-ice area compared to 220 ppm. The situation differs for the MIS 11 simulation: an appreciable increase of inland-ice area occurs for 220 ppm CO 2 concentration but not for 280 ppm. Finally, for MF corresponding to MIS 1 (including the future) extremely small changes in inland-ice area are simulated for both 220 and 280 ppm atmospheric CO 2 concentration, i.e., our model does not simulate glacial inception during MIS 1 even for a rather low CO 2 level. Figure 5a illustrates the complex relationship between crossing of glacial thresholds and the onset of glacial inception. During MIS 5 the MF drops far below the glacial thresholds for both CO 2 concentrations and there is always a strong glacial inception (Fig. 5b). But during MIS 1 and MIS 11 model response is different: considering the position of the insolation minima relative to the glacial thresholds, it is tempting to conclude that there should appear glacial inception during MIS 11 for CO 2 values of 280 ppm and that there would be glacial inception for MIS 1 for an atmospheric CO 2 concentration of 220 ppm. However, this conclusion is not supported by modelling results, which show that a small and short-term drop of MF below the glacial threshold is not sufficient to initiate glacial inception. Indeed, with a half period of about 10 000 years, the orbital precessional forcing is still too fast compared to the time scales of initial ice spreading discussed in the previous section. Our transient simulations indicate that MF has to drop below the glacial threshold by about 20 W m −2 to cause the onset of the Northern Hemisphere glaciation. This value emerges from analysis of the MIS-11 simulations: for the difference of 15 W m −2 between the glacial threshold for 280 ppm CO 2 and the minimum in the MIS-11 MF, our transient simulation does not show a pronounced glacial inception. But there is successful transient simulation of glacial inception for the difference of 25 W m −2 between the same MF minimum and the glacial threshold for 220 ppm CO 2 . From additional MIS-11 simulations with 240 and 260 ppm atmospheric CO 2 concentration (not shown in a figure), we found that a drop of about 20 W m −2 below the glacial threshold is needed to initiate glacial inception. Certainly, the value which MF has to fall below a glacial threshold to initiate glacial inception could be different for MIS 1, MIS 5 and MIS 11, because e.g. the decrease rate in MF varies for these marine isotopic stages. The thin solid lines indicate computed glacial thresholds for 220 and 280 ppm atmospheric CO 2 concentration. (b) Inland-ice area of the Northern Hemisphere except Greenland in 10 6 km 2 during MIS 1 for 220 ppm (thick red curve) and 280 ppm (thin red curve) atmospheric CO 2 concentration, during MIS 5 for 220 ppm (thick blue curve) and 280 ppm (thin blue curve) atmospheric CO 2 concentration as well as during MIS 11 for 220 ppm (thick orange curve) and 280 ppm (thin orange curve) atmospheric CO 2 concentration.

Asynchronous coupling of climate and ice-sheet components
Due to high computational costs it was not possible so far to perform transient simulations of glacial inception with complex Earth system models (AGCMs, AOGCMs). This problem, however, can be circumvented by using an acceleration technique.  This results in an acceleration A=P /t C =3, which implies that computation costs of climate component is reduced by factor three compared to the fully interactive run.
that typical time scales of the most computationally expensive components of complex Earth system models, the ocean and the atmosphere, are much shorter than the orbital time scales. This allows one to reduce run-time of the climate components of an Earth system model by an artificial accelerating of the external forcings, namely, orbital and GHGs. For example, acceleration by a factor of ten would imply that a transient run representing changes in orbital forcing and GHGs from 120 to 110 kyr BP would require only one thousand years of run-time of the climate component. Unlike the climate component, the ice-sheet module cannot be accelerated in the same manner, because the time scale of ice sheets is comparable with the orbital time scales. This, however, does not represent a problem, because current ice-sheet models are computationally relatively inexpensive compared to the AOGCMs and can be run without acceleration on the orbital time scale. The practical question is how much the climate component can be accelerated to produce still a realistic simulation of glacial inception. The CLIMBER-2 model, due to its sufficient degree of realism and low computational cost, is a useful tool to address this question. Below, we compare results of several accelerated runs with the reference run, where climate and ice-sheet components exchange information every year (Fig. 6a). An additional problem which arises for acceleration technique is related to large (and realistic) interannual climate variability simulated by AOGCMs. When the climate component is accelerated, for example, by factor of ten, the mass balance simulated in the climate component for a single year will be applied to the ice-sheet components over ten consecutive years, thus producing unrealistic high interdecadal variability in the mass balance. In a view of a strong nonlinearity of the ice-sheet component, it may cause additional biases in the simulation of ice-sheet dynamics. The most straightforward way to suppress this effect is to run the climate model for several consecutive years and use mass balance, which is averaged over this period of time, to drive the ice-sheet model. This will substantially reduce variability in the icesheet's forcing. If we denote the period of time which is used to reduce the temporal noise in climate model output by t C and the acceleration factor by A, then the time interval of information exchange between climate and ice-sheet components will be P =t C A (Fig. 6b). Therefore, this technique for a given acceleration factor will increase the interval of time of information exchange between climate and icesheet components. For example, for the acceleration factor A=10 and the averaging period of mass balance t C =10 yr, the exchange between ice-sheet and climate component will occur only once per hundred years of the ice-sheet model, which, again, can disturb the simulation of glacial inception. Since the CLIMBER-2 model does not produce realistic interannual variability, we cannot assess the effect of interannual variability on simulation of ice sheet growth in the accelerated runs. But we can assess the effect of reduction of coupling frequency associated with the averaging of climate model output over several years.
Here, we present a set of transient simulations of last glacial inception with different acceleration and time coupling parameters. In these simulations, the model was forced by atmospheric CO 2 concentration from Petit et al. (1999) and the Earth's orbital parameters according to Berger (1978). Equilibrium interglacial climate state serves as initial condition and the model was run from 128 to 100 kyr BP. Figure 7 compares simulated sea level during last glacial inception in experiments with different coupling parameters. Comparison of simulated sea level in the synchronously coupled run (P =1 yr, A=1, denoted in Fig. 7 as "run 1/1") with proxy data was already presented in Calov et al. (2005b) and is rather satisfactory. Simulated sea level is only slightly higher in run 30/3 (P =30 yr, A=3) compared to run 1/1. With a further increase of the time interval of information exchange in run 100/10 the deviation in sea level relative to the reference run is already clearly visible but qualitatively results are not very different. An information exchange interval of 1000 yr leads to extreme reduction in simulated sea level drop, which shows that such an information exchange interval is too long for realistic simulation of glacial inception. To quantify the quality or accuracy of the accelerated climate-cryosphere coupling, we performed a suite of experiments with a broad range of information exchange intervals and accelerations. As quality criterion, we chose the ratio of simulated inland-ice volume in accelerated runs to the reference run at 110 kyr BP. This ratio, measured in percentages, is shown in Fig. 8. The simulations show that time intervals of information exchange between 1 and 100 yr and accelerations inside the interval 1 to 10 result in a ratio between 85 and 100%. An increase of the information exchange interval beyond 300 and/or accelerating the climate module more than thirty fold lead to a quality ratio of less than 60%. The latter is not surprising, because such acceleration factors reduces the time scale of accelerated orbital forcing to the same order as the time scale of the components of the climate system. However, the fact that an acceleration by factor ten still enables to perform a relatively accurate simulation of glacial inception suggest that simulation of glacial inception with an AOGCM coupled to a ice-sheet model can be performed already with the current generation of supercomputers. Fig. 8. Isolines of ratio of ice volume to reference ice volume (run 1/1) at 110 kyr BP in percentages. The grey area marks settings when the climate-module-computation time becomes smaller than 1 year (one seasonal cycle).

Discussion
Our simulated glacial threshold in the MF space in the model configuration without ice-sheet dynamics corresponds to a much lower value of MF compared to that in the model experiment with ice-sheet dynamics. This has important implications for simulations of glacial inception with AGCMs or AOGCMs. Glacial inception in AGCMs/AOGCMs which do not include inland-ice dynamics can occur only for much lower MF than in models which do include cryosphere dynamics, because the ID-off threshold is considerably lower than the ID-on threshold (Fig. 1). In realistic transient simulations of glacial inception interactive ice-sheet dynamics have to be considered.
For constant MF values between S ID off and S ID on the growth of the ice sheet is characterized by three distinctive phases. Although results from both the simple ice-spreading model and the relaxation model agree only passably with those from CLIMBER-2, we think that these two simple models together reproduce first order effects of ice build up. The two simple models together and CLIMBER-2 exhibit a pronounced change in sign of the second derivative of the inland-ice volume. Further on, the increase of inland-ice volume of the simple spreading model strongly depends on the initial inland-ice volume, which, in turn, depends on the initial land area with positive mass balance. This explains why duration of the spreading phase is so sensitive to the value of MF (Fig. 4). If the MF is below the value of S ID off , no spreading phase exists, because snow-albedo feedback is strong enough to produce a large area of permanent snow cover over centennial time scale. In this case, the glaciation of the Northern Hemisphere proceeds much faster and large ice sheets can be formed on precessional time scales.
A threshold in our understanding is strictly tied to a system in equilibrium and, accordingly, our stability diagram ( Fig. 1) is determined with quasi-equilibrium simulations applying slowly varying forcing or, alternatively, with a set of pure equilibrium simulations. This choice is due to the fact that a unique determination of a threshold is impossible via transient simulations, because such a "threshold" will depend on the rate of change of the forcing. The sudden strong increases in ice area in the transient simulations seen in Fig. 5b (e.g. at about −4 kyr from MF min in the run through MIS 5 with 220 ppm CO 2 ) do not correspond to any glacial threshold. The time scale of initial spreading of the small ice caps in the high North can be rather large compared to the time scale of precession (Fig. 4). During the first some thousand years, the curves of inland-ice area in Fig. 5b have a similar characteristics as the ice-volume curves in Fig. 4, which are shown to indicate ice spreading. The strong increases in inland-ice area seen for the MIS 5 and MIS 11 simulations appear at those times when snow albedo feedback starts to contribute strongly to cooling. The timing of the strong increase in inland-ice area is a race between ice spreading, which gradually brings cooling in the system through albedo, and the change in MF, which will eventually start do increase after reached its minimum. It also depends on how deep and how fast the MF decreases. If MF does not reach deep enough (about 20 W m −2 , see Sect. 5) below the glacial threshold, there is no glacial inception in the sense of widespread glaciation as found in geological records.
Our simulations of glacial inceptions shed new light on a debate (Claussen et al., 2005;Crucifix et al., 2005) on the overdue Holocene glaciation hypothesis by Ruddiman (2003). As response, Ruddiman (2007) mentioned that Claussen et al. (2005) as well as Crucifix et al. (2005) did not apply his hypothesized drop in CH 4 in addition to the CO 2 decrease. Claussen et al. (2005) used only a CO 2 drop to 240 ppm, a value which was proposed by Ruddiman (2003). With the atmospheric CO 2 equivalent of 220 ppm, which we used in our simulations, we safely accounted for the radiative effect of the CH 4 drop from 700 to 450 ppb, which was proposed by Ruddiman. Such a decrease in methane corresponds to an additional radiative forcing of about 0.3 W m −2 , which follows from an empirical formula by Houghton et al. (1995) and usage of the efficacy of 1.4 by Hansen et al. (2005). This additionally radiative forcing of 0.3 W m −2 can be mimicked by a further drop in atmospheric CO 2 from 240 to 227 ppm. Still, even if we keep constant the 220 ppm CO 2 concentration during the whole Holocene (including next 10 000 yr), only a minute increase of the Northern Hemisphere ice area occurs even though MF drops below glaciation threshold for this CO 2 concentration over several thousand years. Therefore, during stage 1, unlike stages 5 and 11, the drop in CO 2 does not affect qualitatively model results and does not lead to the appreciable glaciation of the Northern Hemisphere. This is explained by the fact that due to a low eccentricity, the minimum of MF during MIS 1 is unusually weak and not sufficient to cause onset of glaciation.
By fixing eccentricity and obliquity under varying precessional angle, we opted for a rather artificial forcing in our stability analysis. Another possible procedure is to apply the real orbital forcing and, most important, slow down this forcing considerably. Again, slow change in the forcing is the key element in such an analysis.
Using the CLIMBER-2 model, we analysed the accuracy of simulation of the last glacial inception for different acceleration strategies. Under accuracy we understand here the degree of agreement between fully synchronous and accelerated runs. We found that this accuracy depends on both the acceleration factor and the time interval of information exchange. Our results demonstrate that for acceleration factors up to ten and exchange intervals of up to 100 yr (which implies averaging of mass balance over ten model years), the accuracy of simulation of glacial inception remains reasonably good, which indicates that the acceleration technique can reduce considerably computational costs of transient simulation with complex Earth system models.
It should be noted that in our simulations we do not include the carbon cycle. As it was shown by Brovkin et al. (2007), the time scale of response of the marine-carbon system to changing boundary conditions is several thousand years, i.e., it is longer than that of the physical components of the climate system. This fact can impose an additional constraint on the maximum possible acceleration of a comprehensive Earth system model, but this question is beyond the scope of this paper.

Conclusions
1. Two glacial thresholds are determined in different model setups of CLIMBER-2. The value of these thresholds depends also on the CO 2 concentration of the atmosphere (Fig. 5a). For 271 ppm atmospheric CO 2 concentration, the glacial threshold in the model configuration without ice-sheet dynamics corresponds to a very low MF value of 442 W m −2 , while the glacial threshold in the model configuration with ice-sheet dynamics relates to a MF value of 477 W m −2 , which is only slightly smaller than the present-day value in MF. Supposedly, the low threshold applies rather for AGCMs and AOGCMs where ice-sheet dynamics is not accounted for.
2. For moderate forcing (MF around 460 W m −2 , CO 2 concentration of the atmosphere 271 ppm), the CLIMBER-2 model with inclusion of ice-sheet dynamics exhibits rather long transient response times, which are related to initial spreading of small ice caps in the high north. The duration of this initial phase depends strongly on the value of MF. Because of the rather slow response of the climate-cryosphere system, a drop of MF of about 20 W m −2 below the glacial threshold is necessary to produce substantial glaciation during one precessional cycle. This implies that transient glacial inception, in the sense of glaciation of vast areas of the Northern Hemisphere, happens thousands of years later than the time of crossing the glacial threshold. We would like to stress that (glacial) thresholds correspond to a system in equilibrium and just crossing of the glacial threshold does not imply instantaneous glaciation of the Northern Hemisphere.
3. We demonstrated with transient simulations applying extreme scenarios in drop of CO 2 equivalent that glacial inception during the Holocene is all but impossible, because the MF does not drop far enough below the glacial threshold irrespectively of CO 2 concentration. In contrast, during MIS 5 and MIS 11 CLIMBER-2 is able to simulate transient glacial inception with reasonable values in CO 2 equivalent. Therefore, our results do not support Ruddiman's overdue Holocene glaciation hypothesis.
4. With CLIMBER-2 we investigate the consequences of asynchronous coupling of climate models with ice-sheet models. Our results are thought to be applicable to complex (AOGCM-based) Earth system models which include an ice-sheet component. In particular, an acceleration of climate component by up to the factor ten gives good agreement with the full synchronous simulation of glacial inception, whilst for higher acceleration factors the model underestimates considerably the growth of ice sheets during glacial inception. This result indicates that utilisation of asynchronous coupling in principle permit realistic simulation of glacial inception with complex Earth system models already with the existing computers performance.

Ice spreading in a quasi-equilibrium mathematical approach
Here, we would like to derive a simple relation for the dependence of inland-ice volume on time during the initial spreading phase of ice. We assume that small ice caps evolve under temporal constant average surface mass balance as quasiequilibrium profiles solely by ice flow. The inland-ice volume is calculated by assuming an axisymmetric ice sheet. Vialov (1958) and also Paterson (1994), page 243, give an expression for equilibrium profiles depending on the half span L and the summit thickness H as where n is the power exponent in Glen's flow law of ice. Further, we can relate the summit height with the ice sheets half span following Paterson (1994) applying the equation with the parameter k=3.4 (Paterson, 1994). The value for k is uncertain, what can be related to the flow properties of ice, e.g. its temperature dependence. Eq. (A2) is independent from the used power law and holds for any equilibrium profile. For a gross approximation of inland-ice volume, we assume an axisymmetric ice sheet. With Eqs. (A1) and (A2) we obtain by rotating the profile Eq. (A1) around the axis H (x=0) the integral Substitution of y= x L , dx=L dy leads to V = π kf L 5/2 (A3) with the abbreviation f = 2 1 0 y 1 − y n+1 n n 2(n+1) dy .
Geometrically, one can interpret f as a form factor. With higher exponents n the form factor becomes lower, because of the flatter edges of the profile (Fig. 9) www.clim-past.net/5/245/2009/ where ν=n/(n+1). The results for some values of n applying Eq. (A5) are summarised in Table 2. For a perfect plastic axisymmetric ice sheet (n→∞), Oerlemans (2003) found a similar expression for inland-ice volume (Eq. A3, f =8/15). Differences in the expressions are, because Oerlemans assumed an inclined bed, while we presume a flat ice-sheet base. Further on, we adopt f (n=3)≈0.66 in our estimation, because in Glen's flow law n=3 applies in a good approximation for most parts of an ice sheet. Such a value for n is implemented in several ice-sheet models. In particular, n=3 applies in the SICOPOLIS version which is utilised in our paper. Now, we have a relation between inland-ice volume and half span of an ice sheet (Eq. A3) together with values for all needed constants, which will be used in the desired relation between inland-ice volume and time.
Total inland-ice volume is governed by the equation with the inland-ice volume V , the inland-ice area S and the average surface mass balanceb. In accordance with the axisymmetry, the base of the ice-sheet area is assumed to be a circle S = π L 2 .
Combining Eqs. (A3), (A6) and (A7) lead to The differential Eq. (A8) can be integrated according where V 1 denotes the initial inland-ice volume at time t 1 . It should be mentioned that Eq. (A9) describes the initial ice spreading only. During this initial time one can assume approximately constant average surface mass balance in time.
It is clear that the strong power dependence with an exponent equal 5 in Eq. (A9) can lead to very large ice sheets if sufficient long time t has passed. In reality, mechanism like reduction of snowfall through the elevation desert effect limit the inland-ice volume if a large ice sheet develops. Certainly, Eq. (A9) is invalid for a mature ice sheet which can exist later in time.

Appendix B Equilibrium relaxation approximation
If an old ice sheet approaches equilibrium a relaxing exponential function is a good approximation. In the paper we use to demonstrate such a relaxation effect towards an equilibrium value V eq with the relaxation time τ starting with the inland-ice volume V 2 at time t 2 .