Deglacial ice sheet meltdown: orbital pacemaking and CO2 effects

One hundred thousand years of ice sheet buildup came to a rapid end ∼ 25–10 thousand years before present (ka BP), when ice sheets receded quickly and multi-proxy reconstructed global mean surface temperatures rose by ∼ 3– 5C. It still remains unresolved whether insolation changes due to variations of earth’s tilt and orbit were sufficient to terminate glacial conditions. Using a coupled three-dimensional climate–ice sheet model, we simulate the climate and Northern Hemisphere ice sheet evolution from 78 ka BP to 0 ka BP in good agreement with sea level and ice topography reconstructions. Based on this simulation and a series of deglacial sensitivity experiments with individually varying orbital parameters and prescribed CO 2, we find that enhanced calving led to a slowdown of ice sheet growth as early as ∼ 8 ka prior to the Last Glacial Maximum (LGM). The glacial termination was then initiated by enhanced ablation due to increasing obliquity and precession, in agreement with the Milankovitch theory. However, our results also support the notion that the ∼ 100 ppmv rise of atmospheric CO 2 after∼ 18 ka BP was a key contributor to the deglaciation. Without it, the presentday ice volume would be comparable to that of the LGM and global mean temperatures would be about 3 C lower than today. We further demonstrate that neither orbital forcing nor rising CO2 concentrations alone were sufficient to complete the deglaciation.

A number of recent transient modeling studies (Timm and Timmermann, 2007;Timmermann et al., 2009;Ganopolski and Roche, 2009;Roche et al., 2011;Smith and Gregory, 2012;He et al., 2013) quantified the effects of varying orbital parameters, GHG forcing, and ice sheet changes on the evolution of the atmosphere-ocean system during the last deglaciation.However, the ice sheet evolution in these numerical studies was prescribed rather than interactively computed.The ice sheet response to the simulated climate feedbacks, such as the elevation-desert effect (Yamagishi et al., 2005), ice albedo and stationary wave feedbacks (Roe and Lindzen, 2001;Abe-Ouchi et al., 2013), and ocean circulation changes (Gildor and Tziperman, 2000) were not represented.To account for the ice sheet-climate interactions, numerical models of varying complexity have been employed, including zero-dimensional conceptual models (Imbrie and Imbrie, 1980;Paillard, 1998), box models (Gildor and Tziperman, 2000), two-dimensional zonally averaged coupled climate-ice sheet models (Gallée et al., 1992), a threedimensional ice sheet model driven by climatic parameterizations derived from an atmosphere-ocean general circulation model (Abe-Ouchi et al., 2007, 2013), and a threedimensional ice sheet model coupled to an earth system model of intermediate complexity (Ganopolski et al., 2010;Ganopolski and Calov, 2011).To our knowledge, the transient effects of orbital parameter changes and carbon cycle feedbacks during the last deglaciation have not been previously quantified in a coupled ice sheet-climate model with three-dimensional ice sheet, atmosphere, and ocean components.
While the changes of earth's orbital and tilt parameters and the resulting shortwave radiative forcing at the top of the atmosphere are well understood, the carbon cycle feedbacks that led to the reconstructed atmospheric CO 2 changes are subject to ongoing research (Kohfeld and Ridgwell, 2009;Tagliabue et al., 2009;Brovkin et al., 2012;Chikamoto et al., 2012;Menviel et al., 2012).It was demonstrated, however, that the carbon cycle feedbacks played an important role during the last deglaciation and the Quaternary Period (last ∼ 2.6 million years), amplifying glacial-interglacial cycles (Gallée et al., 1992;Yoshimori et al., 2001;Ganopolski and Calov, 2011;Abe-Ouchi et al., 2013).Still, their role in triggering and promoting the climate change and ice sheet retreat during the last deglaciation has not been firmly established.
Here we estimate the transient effects of orbital and GHG changes during the last deglaciation on the global climate and the Laurentide, Eurasian, and Greenland ice sheets using the ice sheet model IcIES, coupled to the atmosphere, ocean and vegetation components of the intermediate complexity model LOVECLIM.The setup of this 3-D coupled ice sheetclimate model is described in Sect. 2. In Sect.3, the main results are presented, including a comparison to ice sheet and sea level reconstructions, an analysis of ice sheet mass balance changes, and sensitivity studies with respect to orbital and GHG changes.The paper concludes with a discussion of the main results, uncertainties, and possible caveats (Sect.4).

Model components
The numerical model used in this study is based on the Ice sheet model for Integrated Earth system Studies IcIES (Saito and Abe-Ouchi, 2004;Abe-Ouchi et al., 2007), bidirectionally coupled to the atmosphere-ocean-sea ice-land components of the intermediate complexity model LOVE-CLIM (Driesschaert et al., 2007) with LGM (Last Glacial Maximum) boundary conditions (Roche et al., 2007).The coupled model is driven by orbital parameter changes (Berger, 1978), and reconstructed atmospheric concentrations of CO 2 (Lüthi et al., 2008), CH 4 (Loulergue et al., 2008), and N 2 O (Schilt et al., 2010).
IcIES is set up to simulate land ice north of 30 • N on a 1 • by 1 • spherical grid.The shallow ice approximation is applied to compute the thickness and temperature evolution of the ice sheets.Ice shelves, which are the floating extensions of the grounded ice sheets, and the dynamics of the grounding line are not included.Note that therefore, the destabilization of grounded ice and acceleration of ice streams via the collapse of shelf ice, which is a process that was potentially important during the deglaciation (Alvarez-Solas et al., 2013), is not accounted for in this study.Ice sheet growth is restricted to a predefined geographic domain, which is kept fixed throughout the simulations (Fig. 1a).Ice sheets cannot grow outside of this domain, and they also cannot spread beyond the domain.The flux of ice through the boundary of the domain is treated as implicit passive calving.Note that ice sheets are allowed to grow in the present-day Hudson Bay, Barents Sea, and Kara Sea.However, ice sheets are not allowed to grow in the present-day Laptev, East Siberian, Chukchi, and Bering Sea.Offline experiments with IcIES (not shown) indicate that, if not prohibited, ice would build up also in the Laptev Sea during the LGM, inconsistent with reconstructions (Svendsen et al., 2004;Peltier, 2004).In addition to the implicit passive calving, a parameterization of active calving into the ocean and proglacial lakes is applied (Abe-Ouchi et al., 2013).Proglacial lakes are not simulated explicitly.Instead, the active calving, mimicked by additional surface melting at a rate of 10 m per year, occurs at grid cells that fulfill each of the three following conditions.First, the bedrock elevation at the grid point is below sea level.Second, the surface mass balance at the grid point is negative, corresponding to an ablation area.And third, at least one of the neighboring grid points fulfills the floating condition.Here, the floating condition for ice only depends on the ice thickness and bedrock elevation, since the sea level is held constant.Bedrock elevation is computed dynamically, assuming local isostatic rebound (mantle density of 4500 kg m −3 , timescale 5000 years).Abe-Ouchi et al. ( 2013) described IcIES sensitivity runs for the last four glacial cycles with respect to this active calving parameterization, mantle density and timescale for isostatic rebound.The surface mass balance is approximated by a positive degree day (PDD) scheme based on Reeh (1991), with PDD factors for water equivalent snow and ice melting of 3 mm day −1 K −1 and 8 mm day −1 K −1 , respectively.The standard deviation of the daily mean temperature from monthly means is assumed to be 5.5 • C and spatially homogeneous.The refreezing rate of melted snow is set to 60 % (for an intercomparison of different PDD models see Charbit et al., 2013).
The atmospheric component of LOVECLIM, ECBilt (Opsteegh et al., 1998), uses quasi-geostrophic equations with ageostrophic correction terms on three vertical levels and with a spectral truncation of T21, corresponding  to a horizontal grid-spacing of ∼ 5.6 • .Total cloud cover is prescribed according to satellite observations (ISCCP D2, Rossow et al., 1996).The effect of carbon dioxide (CO 2 ), methane (CH 4 ), and nitrous oxide (N 2 O) on the longwave radiation is computed with a radiation scheme that is linearized around present-day reference profiles (Schaeffer et al., 1999;Goosse et al., 2010).Here, to account for the low sensitivity of ECBilt to CO 2 variations, and to achieve a good match of the simulated and reconstructed ice sheet evolution during the last 78 ka, the effect of CO 2 deviations from the reference value of 356 ppmv on the longwave flux is scaled (Timm and Timmermann, 2007) by a factor of α = 3.The effects of changes of atmospheric methane and nitrous oxide concentrations on the longwave radiative transport are not modified.With the CO 2 scaling factor of 3, the present-day equilibrium climate sensitivity of IcIES-LOVECLIM to a CO 2 doubling amounts to about 4 • C. The importance and drawbacks of this scaling are discussed in Sect. 4. The ocean-sea ice component of LOVECLIM, CLIO (Goosse and Fichefet, 1999, Coupled Large-scale Ice Ocean), computes ocean currents, salinity, and temperature according to the primitive equations on two spherical subgrids covering the global oceans at a resolution of 3 • by 3 • horizontally, and 20 levels vertically; sea ice thermodynamics and advection are computed on the same horizontal grid.
The terrestrial biosphere component of LOVECLIM, VE-CODE (Brovkin et al., 1999, VEgetation COntinuous DEscription), assigns tree, grass, and desert fractions to land points not covered by ice.

Ice sheet-climate coupling
LOVECLIM and IcIES are coupled asynchronously.The faster-responding atmosphere-ocean-sea ice-vegetation components are subject to accelerated GHG and orbital forcings (Timm and Timmermann, 2007), by a factor of 20, while the land ice component has real time to adjust to the corresponding climate anomalies.Note that, because of this acceleration technique, the ocean in particular has less time in our model to adjust to forcing changes than it had in reality.This means for example that the cooling of the ocean during the deglaciation may be underestimated in our experiments compared to simulations without acceleration.
Surface temperature and precipitation climatologies are passed from ECBilt to IcIES every 1000 years, which is equivalent to 50 model years in LOVECLIM.These monthly mean climatologies are based on the entire 50-year-long LOVECLIM interval, covering the respective 1000 years of GHG and orbital forcing variations.IcIES then uses these climatologies to compute the surface mass balance of the respective 1000-year-long ice sheet simulation.Since the precipitation and temperature are computed on the lowerresolution (∼ 5.6 • ) ECBilt/VECODE grid, a downscaling method to the higher-resolution (1 • ) IcIES grid is required.Here, the surface temperature is bi-linearly interpolated onto the IcIES grid, and subsequently corrected for deviations of the higher-resolution surface height used in IcIES from the lower-resolution ECBilt surface height, assuming a spatially uniform lapse rate of 0.005 • C m −1 .Before this interpolation, the surface temperature is corrected for its present-day annual-mean bias.The bias is computed from the years 1961 to 1990 of a transient LOVECLIM simulation with prescribed GHG and orbital forcing variations, as well as prescribed present-day ice sheets, compared to 2 m air temperature and sea surface temperature observations for the same period (Jones and Moberg, 2003).Note that the same annual-mean bias correction is applied to all seasons.This means that the ice sheet model is driven with the temperature anomalies relative to Jones and Moberg (2003) with respect to the annual mean temperatures.The annual cycle of the simulated temperature is not altered by the bias correction.The surface temperature bias correction allows IcIES to build up more realistic ice sheets, in particular over North America, where the annual-mean surface temperature bias reaches up to 7 • C (Fig. 1b).The precipitation fields are simply bilinearly interpolated onto the higher-resolution IcIES grid; they are not corrected for their present-day biases.
The ECBilt/VECODE boundary conditions are updated based on the simulated ice sheet area and surface height at the end of each 1000-year-long IcIES interval.Each EC-Bilt/VECODE land grid cell is either defined as ice-covered or ice-free; fractional ice coverage is not applied.A cell is defined as ice covered if more than half of it is covered by ice at least 10 m thick at the end of the 1000-year interval using bi-linear interpolation from the finer IcIES grid.If an ECBilt/VECODE grid cell is ice covered, the surface background albedo is set to an ice albedo of 0.4.This low ice albedo is more typical for melting ice.During glacial periods, this low ice albedo hardly affects the surface albedo, because the ice is typically covered by snow, which has a much higher albedo.The snow albedo is set to 0.85 poleward of 72 • N/S, and 0.8 elsewhere.The surface albedo is set to the snow albedo if the surface is covered by at least 0.05 m of snow (water equivalent), and the albedo is linearly interpolated between the ice-and snow albedo for snow thicknesses between 0 and 0.05 m.The snow thickness in ECBilt/VECODE is limited to 10 m water equivalent.Any snow that exceeds this thickness is treated as runoff.The forest fraction in VECODE is set to 0 for ice-covered grid cells in order to avoid the effect of tree canopies on the surface albedo of ice-covered grid cells.Finally, the ECBilt surface elevation is updated based on the IcIES surface elevation at the end of each 1000-year interval.
The ocean model CLIO is only affected by ice sheet variations indirectly, via atmospheric changes.The precipitation that the land surface model receives, river catchment basements, ocean runoff points, and consequently ocean freshwater forcing are not affected by ice sheet variations in the present setup.We choose not to let the buildup and melting of ice sheets affect the ocean freshwater forcing (or ocean heat fluxes), because freshwater volume and freshwater flux cannot both be conserved in accelerated climate runs.Sea level and the ocean bathymetry are not adapted to ice sheet variations.The land-sea mask in ECBilt/VECODE and the ocean bathymetry in CLIO are prescribed according to LGM sea level, based on Roche et al. (2007).In particular, the Bering Strait is closed, and the parameterization of the water and sea ice transport through the Bering Strait in CLIO as described by Goosse et al. (1997) is not applied.

Spin-up of control simulation (CTR)
To achieve a realistic transient simulation of the Last Glacial Maximum (LGM), a two-step initialization approach is pursued.First, the model is initialized for the Eemian interglacial at 125 ka BP (Fig. 1c) using present-day climate and ice sheet initial conditions.However, in this transient simulation the glacial inception is not captured well; sizeable ice sheets only build up after 70 ka BP, in contrast to paleo-climate records, which show an early glacial inception around 115 ka BP.Moreover, it appears that the ice sheets do not reach the critical extent necessary to persist through the phase of large precession and obliquity at the start of marine isotope stage 3 (MIS 3; see Fig. 2c, d for orbital and CO 2 forcings).A possible reason for this late and weak inception is the low spatial resolution of our atmosphere component ECBilt (T21), which ignores small-scale topographic features that may have played an important role in the initial glacial ice sheet buildup (Abe-Ouchi and Blatter, 1993;Pollard and Thompson, 1997;Calov et al., 2005).The simulated ice volume in the first iteration at 63 ka BP is similar to the reconstructed global ice volume (sea level) at 78 ka BP.Hence, in a second iteration, the simulation is restarted at 78 ka BP using the ice sheet state from the initial run at 63 ka BP (Fig. 1c, blue lines).This second step leads to an improved glacial build up of ice during MIS 3, avoiding the unrealistic deglaciation that occurred in the initial iteration.The strong sensitivity of the glacial trajectory to the initial conditions indicates the presence of multiple equilibria for MIS 3 boundary conditions.The re-initialized run is used as the control run (CTR) in the present study.

Ice sheet evolution during the last 78 ka
Given large enough ice sheets prior to MIS 3, the evolution of the Northern Hemisphere ice sheets in CTR with prescribed time-varying GHG and orbital forcing is in good agreement with sea level proxies (Waelbroeck et al., 2002;Peltier, 2004;Yokoyama and Esat, 2011) and ice sheet reconstructions (Peltier, 2004) during the subsequent glacial ice sheet buildup, the LGM, the deglaciation, and the Holocene (Figs. 2a and 3).Furthermore, CTR features a relatively stable Greenland ice sheet across the glacial termination and into the Holocene, in accordance with observations (Fig. 3).
Comparing the deglacial evolution of the simulated ice sheets with the ICE-5G paleo-topographic reconstruction (Peltier, 2004, Figs. 2a and 3a, b), we find good agreement, in particular for the retreat of the Laurentide ice  sheet.On its deglacial retreat, the Laurentide ice sheet passes Lake Winnipeg and the Great Slave Lake in northwestern Canada around 10 ka BP.During this progression the Laurentide ice sheet undergoes a saddle collapse (also called "zipper-effect"), a separation of the Cordilleran and Labrador/Keewatin ice sheets at around 11 ka BP (Fig. 3c), consistent with previous modeling results (Ganopolski et al., 2010;Gregoire et al., 2012;Abe-Ouchi et al., 2013).For the LGM, the simulation and reconstruction of the Laurentide ice sheet differ over Alaska and the Aleutian Islands, and the Laurentide ice sheet dome is further eastward in our model.Compared to ICE-5G, the Greenland and Laurentide ice sheets in our model are thicker towards the margins, which leads to an overestimation of the LGM Northern Hemisphere ice volume in CTR by about 20 m s.l.e. (Fig. 2a).
The simulated Eurasian LGM ice sheet is smaller than the reconstruction and lacks ice cover over the Baltic Sea, the Barents and Laptev seas, central Europe and the British Isles.Note that the prescribed ice sheet domain (Fig. 1a) prevented the Eurasian ice sheet from spreading across the North Sea and into the British Isles.This unrealistic feature will be corrected in future experiments.The large Laurentide ice sheet in our simulation, overcompensating the too-small Eurasian ice sheets, leads to an LGM Northern Hemisphere ice volume of 140 m s.l.e., which is about 20 m s.l.e.greater than in ICE-5G.The large simulated Northern Hemisphere LGM ice volume is still consistent with the reconstructed global sea level lowstands.For example, Yokoyama and Esat (2011) estimated an LGM global mean sea level about 145 m lower than at present.But, assuming that the Antarctic ice sheet contained about 10-20 m s.l.e. more ice during the LGM compared to today, even this relatively high estimate of 145 m by Yokoyama and Esat indicates that the Northern Hemisphere ice volume in our simulation is overestimated by 5-15 m s.l.e.

Orbital and CO 2 effects during the deglaciation
To elucidate the individual roles of atmospheric CO 2 changes and orbital forcing during the deglaciation, two transient sensitivity experiments are performed using LGM (21 ka BP) initial conditions from CTR.In the first experiment, the atmospheric GHG concentrations are kept constant at LGM values, and only the orbital forcing is varied from 21 ka BP to today (Fig. 2a, blue line).Initially, the correspond- ing Northern Hemisphere ice volume decreases until about 9 ka BP, with a similar onset of the deglaciation compared to CTR.However, in response to the switch from a warm to a cold orbit with decreasing obliquity and summer insolation, the ice sheets start to recover.This experiment documents that direct orbital forcing alone was not sufficient to complete the deglaciation.Hence, the observed increase of CO 2 between 18 ka BP to about 9 ka BP was an essential element for the ice sheet retreat during the last glacial termination.If atmospheric CO 2 concentrations had remained at LGM values, large parts of North America would still be covered by thick ice sheets today (Fig. 3d).The second experiment, with fixed LGM orbital forcing and time-varying GHG concentrations (Fig. 2a, orange line), shows a 30-40 % deglacial ice sheet retreat, which is delayed by about 3 ka compared to CTR as a result of the late increase of atmospheric CO 2 concentrations starting around 18 ka BP.The fact that greenhouse gas changes alone are not sufficient for the deglaciation illustrates the importance of orbital forcing both for the timing and the amplitude of the last glacial termination.These sensitivity experiments make a compelling case for the im-portance of orbital and greenhouse gas forcing synergies, as well as for the presence of nonlinearities in the coupled ice sheet-climate system that amplify the response to the combined forcing relative to the sum of the individual responses.A third experiment, starting at 40 ka BP with varying GHGs and obliquity but fixed precession, shows relatively small deviations from CTR (Fig. 2a, purple line).This illustrates that, in our model, obliquity changes play a larger role for the timing of the LGM and the deglaciation than precession variations.

Mass balance analysis
The evolution of the Northern Hemisphere ice sheet volume is determined by the time-integrated positive and negative ice sheet mass balance terms, namely accumulation, ablation, basal melting, and calving (Fig. 2b).From 78 ka BP until after the deglaciation at 10 ka BP in CTR, the total accumulation integrated over the ice sheets is proportional to the ice sheet area (dashed line in Fig. 2b).In other words, the mean accumulation rate over the ice sheets remains almost (zipper effect), d, as simulated for present (0 ka BP), assuming that atmospheric greenhouse gases (GHGs) remained constant at their LGM values (21 ka BP), e, at present in CTR.f, Ice thickness at present in CTR relative to observations (Bamber et al., 2001;Layberry and Bamber, 2001).Black contours show observed present-day coastlines for reference.constant.This supports the notion (Gallée et al., 1992;Yoshimori et al., 2001) that the growth and decay of the Northern Hemisphere ice sheets during this period are caused by variations of melting and calving, and not by variations of accumulation.
At the beginning of CTR, from 78 ka BP to about 60 ka BP, ablation, basal melting, and ice sheet calving are smaller than the accumulation rate, leading to a continuous growth of the ice volume from an initial value of about 40 m to about 90 m sea level equivalent (m s.l.e.).Around 60 ka BP the combined effects of a high precession index, which leads to shorter Northern Hemisphere summers and intensified peak insolation, increasing obliquity (Fig. 2d), and rising CO 2 concentrations (Fig. 2c) lead to a short reversal of the net mass balance and thus decreasing ice volume at the onset of marine isotope state 3 (MIS 3; 60-24 ka BP).In the subsequent period (50-32 ka BP), obliquity and precessional effects appear to almost cancel each other, resulting in a near-neutral mass balance and hence slow ice sheet buildup.This lull is further promoted by the fact that atmospheric greenhouse gas concentrations stay relatively constant during MIS 3. At 32 ka BP, ice sheet growth resumes.The Eurasian ice sheet starts to expand after it had mostly disappeared in our model simulation.During this time, the southern margin of the Laurentide ice sheet migrates southward, and the western margin into Alaska towards the Bering Strait.The larger ice sheet area leads to more accumulation (Fig. 2b).Slightly reduced calving and surface ablation rates further contribute to the larger net ice sheet growth (dV /dt).The ice sheet area and mean accumulation rate keep rising until they peak at the LGM around 19 ka BP in our model.However, despite the positive accumulation trend, dV /dt already reaches its maximum around 28 ka BP, and subsequently decreases due to enhanced calving.Around 21 ka BP, ablation at the surface of the ice sheets intensifies, leading to a further reduction of dV /dt.At 19 ka BP, dV /dt eventually changes sign, which defines the LGM in our model and the beginning of the deglaciation.Note that at the LGM, the surface mass balance is still positive in most areas of the NH ice sheets.Ablation zones are restricted to small areas at the ice sheet margins (Fig. 4a), and the implicit passive calving mechanism (Sect.2.1) balances most of the accumulation.As described in Sect.3.2, the rapidly magnified ablation that leads to the onset of the deglaciation is mostly caused by increasing obliquity, and is later enhanced by rising atmospheric CO 2 .
During the subsequent Holocene period, a very stable Greenland ice sheet remains, for which accumulation, calving and ablation terms on the order of 0.05 Sv cancel each other almost completely out (Fig. 2b).However, the Greenland ice sheet in CTR at present is too thick compared to observations (Bamber et al., 2001;Layberry and Bamber, 2001, Fig. 3f).Compared to high-resolution model results by Ettema et al. (2009, Fig. 4b, c), the surface mass balance in CTR at present is overestimated, in particular in northeast and central Greenland.The observed surface mass balance gradient over southern Greenland from high in the east to negative in the west is not resolved.

Global temperature changes during the deglaciation
The simulated global mean temperature rise during the deglaciation can be explained in large part by rising CO 2 (Fig. 5a).In CTR, which includes both orbital forcing and GHG changes, the global mean temperature rise from 19 to 9 ka BP amounts to about 4 • C, which compares well with recent paleo-climate reconstructions (Shakun et al., 2012;Marcott et al., 2013).Note that, as mentioned in Shakun et al. (2012), the proxy data that form the basis for the reconstructed global mean temperature change are spatially biased towards ocean margins.Hence, the reconstructed amplitude of the global mean temperature change must be regarded with caution.Orbital forcing alone and its effect on surface albedo (causing a global reduction of 0.01, Fig. 5b), only account for about 0.5 • C, whereas the GHG changes alone are responsible for ∼ 2. , for the control run with both orbital and GHG forcing (black, CTR), only GHG forcing (orange, fixed 21 ka BP orbit), and only orbital forcing (blue, fixed 21 ka BP GHGs).Gray shading in panel a shows the reconstructed global mean surface temperature anomaly (Shakun et al., 2012;Marcott et al., 2013).6. a, Northern hemisphere ice sheet volume in meters sea level equivalent (msle) for CTR (black), and for fixed LGM (21 ka BP) greenhouse gas concentrations and orbital parameters (blue line); gray shading indicates sea level reconstructions (Waelbroeck et al., 2002).b, Bedrock height and ice sheet thickness at -90 ka BP (see Fig. 3 for color scales).c, Accumulation (blue area), ablation (light gray shading), basal melting and calving (dark gray shading), net mass balance (thick black line), and ice sheet area for comparison (dashed line) scaled as in Fig. 2b.(Shakun et al., 2012;Marcott et al., 2013).
ing temperature rise of 0.8 • C between 19-9 ka BP is due to synergy effects, characterizing the nonlinear response of the coupled climate-ice sheet system to combinations of external forcings.

Non-equilibrium LGM
Acknowledging the fact that ice sheets are integrators of mass imbalances, it is evident that they are typically not in equilibrium with external forcings on orbital timescales.The transient nature of the LGM ice sheet-climate system is further illustrated by branching off another sensitivity experiment from the LGM state of CTR with fixed LGM orbital and CO 2 forcings (Fig. 6a).After about 90 000 years of ice sheet integration time, the Northern Hemisphere ice sheets reach a near-equilibrium volume of about 300 m s.l.e.and, in contrast to the transient LGM ice sheets in CTR, cover vast areas of northern Asia/Siberia (Fig. 6b).At near-equilibrium, the large ice sheet melting and calving terms of together more than 40 m s.l.e. per thousand years (almost 0.5 Sv) show little variability, and are compensated by accumulation (Fig. 6c).The simulated North American LGM ice sheet in CTR has a similar extent to that in the fixed LGM boundary condition sensitivity experiment, suggesting that the Laurentide ice sheet under LGM conditions in CTR may have been close to its long-term steady state.These results are consistent with previous results by Abe-Ouchi et al. (2013).They illustrated in their Fig.2b that, during glacial cycles, in insolation/ice volume phase space, the transiently simulated Laurentide ice sheet oscillates much closer to the respective equilibrium solutions than the Eurasian ice sheet.

Discussion and conclusions
Supporting the astronomical theory of ice ages, the presented model results demonstrate that the initial timing of the last glacial termination was determined by orbitally induced insolation changes triggering rapid Northern Hemisphere ice sheet retreat shortly after the LGM.This finding is in agreement with recent sea level estimates (Clark et al., 2012), which indicate early meltwater releases around 19 ka BP.However, our model results also illustrate that orbital forcing alone was not sufficient, and that the ∼ 100 ppmv CO 2 rise needs to be taken into account to explain the fast and full retreat of the North American and Eurasian ice sheets during the deglaciation.This supports previous model studies, which have identified GHG variations as an important contributor to Pleistocene glacial cycles (Gallée et al., 1992;Yoshimori et al., 2001;Ganopolski and Calov, 2011;Abe-Ouchi et al., 2013).And the results emphasize how important it is to identify the mechanisms responsible for the orbitalscale modulation of atmospheric CO 2 concentrations (Kohfeld and Ridgwell, 2009;Tagliabue et al., 2009;Chikamoto et al., 2012;Menviel et al., 2012;Brovkin et al., 2012).
In the present study, the ice sheet surface mass balance is computed using a positive degree day (PDD) scheme (Sect.2.1).PDD schemes are empirical relations designed to compute the melt potential solely from temperature, and the PDD scheme used here was derived from present-day observations in Greenland (Reeh, 1991).PDD schemes are not physically based, and do not account explicitly for shortwave radiation changes.It is therefore questionable whether PDD schemes are applicable under variable orbital forcing (van de Wal, 1996;Robinson et al., 2010).Since the surface mass balance is crucial for the simulation of glacial cycles, we plan to test our results with an energy-balance-based surface mass balance scheme, replacing the PDD scheme, but this is beyond the scope of the present study.
As described in Sect.2.1, the effect of CO 2 variations on the longwave radiative transport in the atmosphere is scaled up by a factor α=3.The scaling is necessary in our setup to match the simulated evolution of the NH ice volume to the  (Waelbroeck et al., 2002).b, Bedrock height and ice sheet thickness at -90 ka BP (see Fig. 3 for color scales).c, Accumulation (blue area), ablation (light gray shading), basal melting and calving (dark gray shading), net mass balance (thick black line), and ice sheet area for comparison (dashed line) scaled as in Fig. 2b.imulated Northern Hemisphere land ice volume to the parameter α, which scales the effect of CO 2 on the longosphere (Section 2.1).The default value as used in this study is α=3 (CTR).Gray shading indicates sea level iniferal oxygen isotope ratios (Waelbroeck et al., 2002).
Figure 7. Sensitivity of the simulated Northern Hemisphere land ice volume to the parameter α, which scales the effect of CO 2 on the longwave radiation in the atmosphere (Sect.2.1).The default value as used in this study is α = 3 (CTR).Gray shading indicates sea level reconstructions from benthic foraminifera oxygen isotope records (Waelbroeck et al., 2002).
reconstructed sea level change during the last 78 ka (Fig. 7).Sensitivity runs with α = 2.5, α = 3, and α = 3.5 illustrate that the choice of α has a large effect on the results, and that α = 3 yields the best fit with sea level reconstructions.The high sensitivity of our results to the scaling indicates that a realistic glacial cycle can only be achieved within a small range of climate sensitivities.While this scaling is artificial, it is motivated by the fact that the CO 2 sensitivity of LOVECLIM is in the lower range compared to more complex general circulation models (Timm and Timmermann, 2007), possibly due to missing feedback processes.For example the convective cloud feedback (Abbot and Tziperman, 2009), which potentially amplifies the orbital and CO 2 forcing at high lati-tudes, is neglected, because clouds are prescribed according to satellite observations (Sect.2.1).Future simulations including an empirically derived interactive cloud parameterization (Sriver et al., 2014) are planned.
The simulated 3 ka lag of the onset of the deglaciation in the sensitivity run with only greenhouse gas forcing compared to CTR suggests that the deglacial CO 2 rise affected the Northern Hemisphere ice sheets only 3 ka after the orbital forcing initiated the deglaciation.However, this quantitative result has to be regarded with caution, because the atmosphere and ocean model components are run with orbital and GHG forcing accelerated by a factor of 20.For each climate-ice sheet coupling interval, the atmosphere and ocean only have 50 years to adjust to the 1000-year orbital and GHG changes before a new climatology is passed to the ice sheet model.This means that the 3 ka lag is equivalent to only 3 coupling time steps or 150 atmosphere-ocean years.Because of this acceleration, the simulated warming and ice sheet response to the rising CO 2 during the deglaciation may be underestimated, and slower than in reality.Uncertainties with respect to the dating of the deglacial CO 2 increase, which are on the order of 1 ka (Schwander et al., 2001;Parrenin et al., 2007;Pedro et al., 2012;Parrenin et al., 2013), add to the uncertainty of the estimated lag between orbital and CO 2 effects during the deglaciation.For the period after 22 ka BP, the composite CO 2 record used here (Lüthi et al., 2008) was based on measured CO 2 concentrations in gas bubbles enclosed in the Antarctic Dome C ice core (Monnin et al., 2001), and the EDC3 gas age chronology (Parrenin et al., 2007).
Paleoclimate data and model studies indicated that the atmosphere during the LGM transported more dust than at present (Mahowald et al., 1999, and references therein).Dust, while it is in the atmosphere, can increase the planetary albedo, which causes a cooling.Dust deposition on snow and www.clim-past.net/10/1567/2014/ice reduces their albedo, which causes a warming.None of these processes are accounted for in the present study.During the LGM, the atmospheric dust was estimated to cause an additional global cooling of ∼ 0.5 to 1.4 • C (Calov et al., 2005;Schneider von Deimling et al., 2006;Mahowald et al., 2006).Hence, on the one hand, the reduction of atmospheric dust concentrations during the deglaciation may have been a positive feedback process.On the other hand, the reduction of dust deposition rates on snow and ice during the deglaciation may have been a negative feedback process.Previous studies with the coupled ice sheet-climate model CLIMBER-2, including a dust model, indicated that aeolian dust feedbacks played an important role during glacial cycles (Ganopolski et al., 2010;Ganopolski and Calov, 2011).
Ice sheet meltwater in the form of river runoff or icebergs has long been identified as a potential cause for ocean circulation changes (Ruddiman and McIntyre, 1981).Reconstructions of the deglacial drainage chronology and drainage locations (Tarasov and Peltier, 2006), and atmosphere-ocean model studies (Mikolajewicz et al., 1997;Roche et al., 2009;Liu et al., 2009;Okazaki et al., 2010;Menviel et al., 2011) suggested that meltwater from the Laurentide ice sheet led to reduced North Atlantic Deep Water formation after the LGM.The associated reduction of the northward heat transport in the Atlantic may have first led to a slowdown of the deglaciation.However, the smaller heat transport in the Atlantic may have been compensated by an intensified heat transport in the Pacific (Okazaki et al., 2010).And the subsequent recovery of the Atlantic meridional overturning circulation (Liu et al., 2009) may have accelerated the deglaciation again, for example via the destabilization of grounded ice due to the collapse of shelf ice (Alvarez-Solas et al., 2013).In the accelerated experiments presented here, the ice sheets and the oceans do not exchange freshwater (Sect.2.2), and the deep ocean circulation remains qualitatively similar to the presentday circulation, with strong North Atlantic Deep Water formation (not shown).Future unaccelerated runs that include ice shelves and an interactive hydrological cycle are required to study thermohaline circulation changes and their role during the deglaciation.
Fig.1.Model setup and initialization.a, Prescribed domain in which ice sheet growth is possible (green area), compared to the ocean areas in which no ice sheet growth is possible, despite present-day water depths below 500 m (gray area, present-day topography interpolated from ETOPO5, 1988).b, LOVECLIM annual mean surface air temperature (SAT) bias computed for the period 1961 to 1990 relative to a climatology based on observations(Jones and Moberg, 2003).c, Simulated Northern Hemisphere (NH) ice sheet volume in meters sea level equivalent (msle) for the initial run starting in the Eemian (125 ka BP, blue curve), and the control run in the present study (CTR, black curve), compared to a reconstruction of global sea level based on benthic foraminifera oxygen isotope records(Waelbroeck et al., 2002, light gray shading indicates confidence interval +/-13 msle, which accounts for regression and dating uncertainties).Note that a sizeable inception in the initial run only sets in around 70 ka BP.The blue arrow illustrates the re-initialization, CTR is started at 78 ka BP with initial conditions from the initial run at 63 ka BP.

Fig. 2 .
Fig. 2.Ice sheet evolution, mass balance, and forcing.a, Simulated Northern Hemisphere (NH) ice sheet volume in meters sea level equivalent (msle) for a run with orbital and GHG forcing (control run CTR, thick black), only GHG forcing (orange, fixed 21 ka BP orbit), only orbital forcing (blue, fixed 21 ka BP GHGs), and fixed precession (purple, eccentricity and longitude of perihelion fixed at 40 ka BP, obliquity and greenhouse gases vary as in CTR).Reconstructed NH ice volume(Peltier, 2004, ICE-5G, teal dashed line), and reconstructed global sea level changes from Barbados coral record(Peltier and Fairbanks, 2006, vertical gray lines), from benthic foraminifera oxygen isotope records(Waelbroeck et al., 2002, light gray shading indicates confidence interval), and from coral terraces at Huon Peninsula, Papua New Guinea(Yokoyama and Esat, 2011, gray circles).b, Ice volume rate of change due to accumulation (blue area), ablation (light gray area), basal melting and calving (dark gray area), and their balance (net, thick black line) from CTR. Dashed gray line indicates the simulated NH ice sheet area scaled by max(accumulation)/max(area).c, CO2 forcing(Lüthi et al., 2008).d, Precessional index (black line, high values indicate NH summers closer to the sun), and obliquity(Berger, 1978, dashed blue line).The dotted gray vertical lines indicate boundaries between marine isotope stages (MIS).

Figure 2 .
Figure 2. Ice sheet evolution, mass balance, and forcing.(a) Simulated Northern Hemisphere (NH) ice sheet volume in meters sea level equivalent (m s.l.e.) for a run with orbital and GHG forcing (control run CTR, thick black), only GHG forcing (orange, fixed 21 ka BP orbit), only orbital forcing (blue, fixed 21 ka BP GHGs), and fixed precession (purple, eccentricity and longitude of perihelion fixed at 40 ka BP, obliquity and greenhouse gases vary as in CTR).Reconstructed NH ice volume(Peltier, 2004, ICE-5G, teal dashed line), and reconstructed global sea level changes from Barbados coral records(Peltier and Fairbanks, 2006, vertical gray lines), from benthic foraminifera oxygen isotope records(Waelbroeck et al., 2002, light gray shading indicates confidence interval), and from coral terraces at Huon Peninsula, Papua New Guinea(Yokoyama and Esat, 2011, gray circles).(b) Ice volume rate of change due to accumulation (blue area), ablation (light gray area), basal melting and calving (dark gray area), and their balance (net, thick black line) from CTR. Dashed gray line indicates the simulated NH ice sheet area scaled by max(accumulation)/max(area).(c) CO 2 forcing(Lüthi et al., 2008).(d) Precession index (black line, high values indicate NH summers closer to the sun), and obliquity(Berger, 1978, dashed blue line).The dotted gray vertical lines indicate boundaries between marine isotope stages (MIS).

Figure 3 .
Figure 3. Bedrock height and ice sheet thickness (a) at the Last Glacial Maximum (LGM) as simulated with IcIES-LOVECLIM, (b) at theLGM according to ICE-5G reconstruction(Peltier, 2004), (c) at 11 ka BP, which is when the Cordilleran and Labrador ice sheets separate in CTR (zipper effect), (d) as simulated for present (0 ka BP), assuming that atmospheric greenhouse gases (GHGs) remained constant at their LGM values (21 ka BP), (e) at present in CTR.(f) Ice thickness at present in CTR relative to observations(Bamber et al., 2001; Layberry  and Bamber, 2001).Black contours show observed present-day coastlines for reference. Fig.

Figure 4 .
Figure 4. Surface mass balance (a) in CTR at LGM (19 ka BP), (b) in CTR at present, and (c) interpolated onto the IcIES grid from a high-resolution model simulation by Ettema et al. (2009).

Fig. 5 .
Fig.5.Evolution of global mean annual mean surface temperature anomaly relative to 21 ka BP (a), and surface albedo (b), for the control run with both orbital and GHG forcing (black, CTR), only GHG forcing (orange, fixed 21 ka BP orbit), and only orbital forcing (blue, fixed 21 ka BP GHGs).Gray shading in panel a shows the reconstructed global mean surface temperature anomaly(Shakun et al., 2012;Marcott et al., 2013).
Fig.6.a, Northern hemisphere ice sheet volume in meters sea level equivalent (msle) for CTR (black), and for fixed LGM (21 ka BP) greenhouse gas concentrations and orbital parameters (blue line); gray shading indicates sea level reconstructions(Waelbroeck et al., 2002).b, Bedrock height and ice sheet thickness at -90 ka BP (see Fig.3for color scales).c, Accumulation (blue area), ablation (light gray shading), basal melting and calving (dark gray shading), net mass balance (thick black line), and ice sheet area for comparison (dashed line) scaled as in Fig.2b.

Figure 5 .
Figure 5. Evolution of global mean annual mean surface temperature anomaly relative to 21 ka BP (a), and surface albedo (b) for the control run with both orbital and GHG forcing (black, CTR), only GHG forcing (orange, fixed 21 ka BP orbit), and only orbital forcing (blue, fixed 21 ka BP GHGs).Gray shading in panel (a) shows global mean surface temperature anomaly reconstructions(Shakun et al., 2012;Marcott et al., 2013).
Fig.6.a, Northern hemisphere ice sheet volume in meters sea level equivalent (msle) for CTR (black), and for fixed LGM (21 ka BP) greenhouse gas concentrations and orbital parameters (blue line); gray shading indicates sea level reconstructions(Waelbroeck et al., 2002).b, Bedrock height and ice sheet thickness at -90 ka BP (see Fig.3for color scales).c, Accumulation (blue area), ablation (light gray shading), basal melting and calving (dark gray shading), net mass balance (thick black line), and ice sheet area for comparison (dashed line) scaled as in Fig.2b.

Figure 6 .
Figure 6.(a)Northern Hemisphere ice sheet volume in meters sea level equivalent (m s.l.e.) for CTR (black), and for fixed LGM (21 ka BP) greenhouse gas concentrations and orbital parameters (blue line); gray shading indicates sea level reconstructions(Waelbroeck et al., 2002).(b) Bedrock height and ice sheet thickness at −90 ka BP (see Fig.3for color scales).(c) Accumulation (blue area), ablation (light gray shading), basal melting and calving (dark gray shading), net mass balance (thick black line), and ice sheet area for comparison (dashed line) scaled as in Fig.2b.
(Waelbroeck et al., 2002zation.a,Prescribeddomain in which ice sheet growth is possible (green area), compared to the ocean areas in which no ice sheet growth is possible, despite present-day water depths below 500 m (gray area, present-day topography interpolated from ETOPO5, 1988).b,LOVECLIMannualmeansurface air temperature (SAT) bias computed for the period 1961 to 1990 relative to a climatology based on observations(Jones and Moberg, 2003).c,SimulatedNorthern Hemisphere (NH) ice sheet volume in meters sea level equivalent (msle) for the initial run starting in the Eemian (125 ka BP, blue curve), and the control run in the present study (CTR, black curve), compared to a reconstruction of global sea level based on benthic foraminifera oxygen isotope records(Waelbroeck et al., 2002, light gray shading indicates confidence interval +/-13 msle, which accounts for regression and dating uncertainties).Note that a sizeable inception in the initial run only sets in around 70 ka BP.The blue arrow illustrates the re-initialization, CTR is started at 78 ka BP with initial conditions from the initial run at 63 ka BP.