Interactive comment on “Investigating the evolution of major Northern Hemisphere ice sheets during the last glacial-interglacial cycle” by S. Bonelli et al

parameterisations on the coarse CLIMBER grid. Therefore, the parameterisations cannot directly affect the surface mass balance of the ice sheets and might have a minor impact only. The authors of the reviewed paper should demonstrate how much their dust and snow aging parameterisations influence the results. How strong do these parameterisations affect mass balance and ice volume? Further, a similar dust parameterisation following Warren and Wiscombe (1980) was already described and applied to last glacial inception by Calov et al. (2005). Calov et al. (2005) introduced the dust parameterisation in their high resolution surface energy balance module. Their paper should be cited at the place where the dust parameterisation is introduced. The albedo is computed on several vertical layers; therefore, the impact of dust parameterization on the albedo is accounted for at the topography corresponding to that of the ISM. As a consequence, in our model, contrary to the Dr. Calov’s statement, this parameterization directly affects snow mass balance. To better assess its effect, we have performed a test where no dust is accounted for. In this simulation, the NH ice sheets do not completely retreat during the deglaciation period and an ice sheet is produced over Siberia, reaching a volume of 35 x 1015 m3 at the LGM. This effect has been discussed in the revised version of the manuscript (see section 5, discussion). Compared to the previous study by Calov et al. (2005b), the impact of dust on the first phase of the glacial cycle is minor, probably because in the early phase of the glacial inception the atmospheric CO2 concentration is close to its pre-industrial level and dust weight is therefore negligible. For further discussion, please see point 7, answers to comments by reviewer #1. The paper by Calov et al. (2005) is now cited where dust parameterization is described.

1 Introduction Milankovitch's theory (1941) states that the succession of glacial-interglacial cycles is primarily driven by changes in the seasonal distribution of insolation induced by the vari-Correspondence to: S. Bonelli ( stefano.bonelli@lsce.ipsl.fr)ations of the Earth's orbital parameters.Nevertheless, the duration and the intensity of cold and warm periods is also strongly influenced by a series of nonlinear amplification mechanisms involving the atmospheric CO 2 concentration (Berger et al., 1999;Shackleton, 2000), changes in the thermohaline circulation (Adkins et al., 1997;Stocker, 2000;Khodri et al., 2001;Rahmstorf, 2002;Kageyama et al., 2006) and in the vegetation cover (DeNoblet et al., 1996;Kageyama et al., 2004;Claussen et al., 2006), as well as a set of feedbacks due to ice sheets.These include the icealbedo feedback, the elevation effect or the ice sheet influence on atmospheric and oceanic circulation (Gallée et al., 1992;Clark and Pollard 1999;Kageyama et al., 2004;Calov et al., 2005a,b;Tarasov and Peltier, 2006;Abe-Ouchi et al., 2007).
The development of major ice sheets does not follow a linear behaviour through the last climatic cycle, but proceeds by alternating phases of expansion and regression (Dyke and Prest, 1987;Boulton and Clark, 1990;Clark et al., 1993;Dyke et al., 2002;Svendsen et al., 2004).The average eustatic sea level directly reflects the evolution of grounded ice volume.According to most reconstructions inferred from benthic records and coral reefs, it dropped by 120 to 140 m around 21 kyr BP at the Last Glacial Maximum (LGM) as a direct consequence of continental ice build-up (Camoin et al., 2001;Lambeck and Chappell, 2001;Yokoyama et al., 2001;Waelbroeck et al., 2002;Siddall et al., 2003;Bintanja et al., 2005).The corresponding extent of the ice cover during the last glacial cycle is provided by a set of geological and geomorphological reconstructions (Boulton and Clark, 1990;Dyke et al., 2002;Svendsen et al., 2004).These works suggest that at the LGM large ice sheets covered extended parts of North America and Eurasia at intermediate and high latitudes in addition to Greenland and Antarctica.However, uncertainties still remain about the shape, volume and thickness Published by Copernicus Publications on behalf of the European Geosciences Union.
S. Bonelli et al.: Investigating the evolution of major Northern Hemisphere ice sheets of these former ice sheets.Refined models of glacial isostatic adjustments, constrained by relative sea-level observations (Yokoyama et al., 2001;Lambeck et al., 2001;Lambeck et al., 2002;Milne et al., 2002) and by geological reconstructions of the ice margins (Peltier, 1994(Peltier, , 2004;;Lambeck et al., 2006) can produce maps of the ice extent and thickness.They are also useful to study the internal viscoelastic structure of the Earth.Nevertheless, being constrained by relative sea-level data, the reconstructed ice thickness is often under-constrained in regions in which no data are available.In addition, these models cannot provide a unique solution to reconstruct the temporal evolution of the ice sheets.
Further knowledge of the Earth's glacial and climatic history can be obtained by the numerical modeling of ice sheets.Several studies have been performed to simulate the extent and volume of the Northern Hemisphere (NH) ice sheets during glacial-interglacial cycles or during specific periods such as the last glacial inception or the last deglaciation.To achieve this goal, various approaches are possible.A classical approach consists in forcing an ice sheet model (ISM) with paleoclimatic data (e.g.Siegert and Dowdeswell, 2004) or general circulation model (GCM) outputs (Verbitsky and Oglesby, 1992;Fabre et al., 1998;Yamagishi et al., 2005).The aim of these studies is to examine whether the simulated ice sheets are in equilibrium with the forcing climate.Some of these studies underline the critical role of snow mass balance in determining the topography and the evolution of modelled ice sheets (Fabre et al., 1998;Yamagishi et al., 2005).This latter point is critical, since models sometimes need to be corrected to better account for precipitation and, hence, snow accumulation (Fabre et al., 1998).A possibility to examine the ice volume evolution through time consists in forcing ISMs with GCM time-slices experiments interpolated through time using a climatic index of the surface temperature obtained for example from isotopic records from ice cores (Marshall and Clark, 2002;Zweck and Huybrechts, 2005;Charbit et al., 2007).This approach allows a detailed study of the response of the ISM to climate forcings (Zweck and Huybrechts, 2005;Charbit et al., 2007).It successfully provides consistent reconstructions of NH ice sheets topography, but the ice evolution does not feedback on climate.An alternative to this experimental set-up consists in the full coupling between simplified climate models driven by insolation and atmospheric CO 2 concentration and ISMs.Following this strategy, the evolution of the NH ice volume over the last glacial-interglacial cycle was successfully simulated with a simplified 2-D ISM under the insolation and atmospheric CO 2 forcings (Gallée et al., 1992;Berger et al., 1999;Loutre and Berger, 2000;Crucifix et al., 2001), whereas Tarasov andPeltier (1997, 1999) reproduced the 100-kyr cycle using an energy balance model coupled with a vertically integrated ISM (Tarasov and Peltier, 1997) or with a 3D-ISM (Tarasov and Peltier, 1999).These works show the role of changes in the orbital parameters as triggering mechanisms for glacial inception, as well as the critical effect of atmospheric CO 2 levels in determining a glacial history consistent with geomorphological data.However, the results may be influenced by oversimplified processes in the climate models or in the ISMs.
In this paper we increase the degree of process representation by performing transient simulations with a climate model of intermediate complexity (atmosphereocean-vegetation), CLIMBER 2.3 (Petoukhov et al., 2000;Ganopolski et al., 2001), coupled with a 3-D thermomechanical ISM for the NH ice sheets, GREMLINS (Ritz et al., 1997).CLIMBER 2.3 incorporates more physical mechanisms than energy balance models and is able to represent, even though rather simplistically, the major features of the atmosphere-vegetation-ocean system (Petoukhov et al., 2000).On the other hand, 3-D thermo-mechanical ISMs such as GREMLINS produce a more realistic topography of simulated ice complexes than 2-D models; they also explicitly account for the full 3-D temperature field in the ice sheet, the basal melting and the ice dynamics (Ritz et al., 1997).The fast computational time of the coupled model makes it appropriate to run multimillennia transient simulations, which are still too computationally expensive to be performed with more refined GCMs.An earlier version of this model was successfully used to simulate the last glacial inception (Kageyama et al., 2004) and the last deglaciation (Charbit et al., 2005).In particular, Kageyama et al. (2004) showed the importance of ice-climate feedbacks and vegetation changes for the last glacial inception, while Charbit et al. (2005) studied the critical effect of atmospheric CO 2 concentration in modulating the timing of the Fennoscandian ice sheet deglaciation (Charbit et al., 2005).The present work enlarges the perspective of these previous studies to the simulation of the whole climatic cycle through a transient 126-kyr long run.
The aim of this paper is first to evaluate the model performances in simulating the last glacial-interglacial cycle and secondly to study the responses of the two major NH ice sheets (Laurentide and Fennoscandia) to the atmospheric CO 2 concentration and insolation forcings, as well as their sensitivity to these parameters.

The CLIMBER climate model
The CLIMBER 2.3 model is a climate model of intermediate complexity describing the atmosphere, ocean, sea ice, land surface processes and terrestrial vegetation cover (Petoukhov et al., 2000;Ganopolski et al., 2001).The atmosphere and the land surface scheme modules have a spatial resolution of 51 • in longitude and 10 • in latitude.The atmospheric module is a 2.5 dimension statistical-dynamical model.This means that the prognostic variables depend on the latitude and longitude but the variations of temperature and humidity with height are computed using a fixed law.It only resolves large scale processes, the heat and moisture transports due to synoptic mid-latitude weather systems being parameterized.Land cover is divided into six different surface types (i.e.forests, grasslands, deserts, open ocean, sea-ice and glaciers) and each grid cell may contain the different surface types.The dynamical terrestrial vegetation model, VECODE, (Brovkin et al., 1997) provides the distribution of trees, grass and deserts over each continental cell as a function of climate.Each of these three surface types may be or not covered by snow.The ocean model is composed of three zonal 2-D latitude-depth (2.5 • ×20 uneven vertical layers) basins (Atlantic, Indian and Pacific), including the Arctic and the Austral Oceans.The mean zonal effect of the longitudinal transport is taken into account through imposed advective transport at key latitudes (e.g.Drake Passage, Austral Ocean).The ocean module also includes a thermodynamic sea-ice model which predicts, for each grid cell, the sea ice fraction and thickness with a simple treatment of diffusion and advection of sea ice.Despite its coarse spatial resolution and the simplified representation of climate processes, CLIMBER has successfully been used for various paleo-climate studies and is therefore appropriate to investigate climatic contexts very different from the present-day one (Ganopolski et al., 1998;Kubatzki et al., 2000;Brovkin et al., 2002).The other advantage of CLIMBER compared to other EMICs is that it incorporates a full description of the hydrological cycle (i.e.humidity transport, precipitation and evaporation are explicitly computed) and is therefore well suited to be coupled with an ice-sheet model.However, some processes having a huge impact on the snow albedo and hence on the ice-sheet surface mass balance are poorly or not represented in the model.As an example, in the original version, the albedo of glaciers was set to that of fresh snow (0.65 to 0.95 depending on the wave length and on the cloudiness) and the snow aging process was not accounted for, although Gallée et al. (1992) have shown that it was critical to properly reproduce a full glacial-interglacial cycle.Moreover, ice-core records show that dust concentration was much larger during glacial periods than during interglacials (Petit et al., 1999).A recent study (Krinner et al., 2006) has suggested that periods of strong dust deposition in Northern Asia prevent the formation of a perennial snow cover in this region, through a drastic snow albedo reduction in response to the dust blackening effect.Therefore, to improve the surface mass balance calculation, we have developed a simple parameterization of the effect of dust deposition on snow albedo.Snow albedo is computed on several vertical layers and the impact of dust parameterization on the albedo is therefore accounted for at the topography corresponding to that of the ISM.We assume that at a given time, the dust deposition Dep (i, n, t) at each location of latitude i and longitude n is expressed as a linear combination of the LGM and the Present-Day depositions taken from Mahowald et al. (1999): where Dep LGM (i, n) and Dep PD (i, n) represent the deposition at the LGM and the present-day periods, with ω(t) calibrated as a function of the atmospheric CO 2 variations: The dust parameterization implemented here is similar to the one adopted by Calov et al. (2005a,b) in the coupled CLIMBER-SICOPOLIS model, with the exception that we compute the dust weight ω(t) as a function of the atmospheric CO 2 concentration, instead of the simulated ice volume (Calov et al., 2005a,b).This enables us to keep the dust content independent from the ice volume itself (i.e.dust is considered here as an external forcing).Dust concentration is assumed to be a function of dust deposition and precipitation PRC(i, n, t) (i.e.dry deposition is ignored): This means that the larger the precipitation, the lower the dust concentration.The snow albedo over glacier areas is then computed as follows: where α 0 is the albedo of fresh snow, and the reduction δ(c) of snow albedo due to dust deposition is directly dependent on dust concentration.This dependency has been determined using the results from a model designed for the calculation of snow albedo (Warren and Wiscombe, 1980).
The correction δ(c) ranges between 0 (for dust concentrations below 45 ppm) to 0.4 (for dust concentrations greater than 275 ppm).This new parameterization considerably improves the model performances under glacial conditions (not shown).Similarly to Calov et al. (2005b), the major effect of dust parameterization on the simulated ice evolution is to produce a lower NH ice volume due to the decreased albedo effect.Nevertheless, for the early phase of glacial inception, its impact on the simulated ice volume is less than in Calov et al. (2005b) because the atmospheric CO 2 concentration is close to its pre-industrial level and dust weight is therefore negligible.The effect of the dust parameterization becomes more important after 70 ky BP, since it prevents the formation of a large ice sheet over Siberia, in agreement with the results of previous studies (Krinner et al., 2006).

The GREMLINS ice-sheet model
GREMLINS is a three dimensional thermo-mechanical icesheet model which predicts the evolution of the geometry (extension and thickness) of the ice and the coupled temperature and velocity fields (Ritz et al., 1997) Tarasov andPeltier (1997, 1999), and by the rigidity of the lithosphere.The temperature field is computed both in the ice and in the bedrock by solving a time-dependent heat equation.Changes in the ice thickness with time are computed from a continuity equation and are a function of the ice flow, the surface mass balance and the basal melting.The ice flow results both from internal ice deformation and basal sliding.It is calculated with the zero-order shallow ice approximation (Hutter, 1983;Ritz et al., 1997).Since the ice-shelf dynamics is not included in the model, the ice is lost by calving by setting the ice thickness to zero when ice begins to float.The ice can advance onto the continental shelf provided that the ice thickness is 75% or higher than the basin depth.This implies that ice can advance over shallow seas only and, therefore, that the calving rate may be underestimated.
The surface mass balance is the sum of accumulation and ablation.The accumulation term is computed by CLIMBER 2.3 (see below) and the ablation term is computed using the positive-degree-day (PDD) method, which is based upon an empirical relationship between air temperatures and melt processes.In the present study, this method is used exactly as described in Reeh (1991).It accounts for the possibility of reaching positive temperatures, and hence melting, even if the mean daily temperature is below the melting point.The standard deviation of the daily temperature is assumed to be 5 • C from 126 to 21 kyr BP.The basic assumption of this method is that ablation is linearly dependent on the sum of positive-degree-days over the entire year.The ablation is computed by assigning different ablation rates for snow and ice to account for their albedo difference.Following Reeh (1991), we assume that 60% of snow melting refreezes and produces superimposed ice and warming due to the phase change.The degree-day factors for snow and ice are respectively set to 3 and 8 mm day 1 C 1 , following values suggested by Braithwaite (1995) from observations in central-west Greenland.The use of these PDD parameters allows to satisfactorily simulate the present-day Greenland ice sheet (Ritz et al., 1997) and are the most commonly used in the ice-sheet modeling community (e.g.Ritz, et al., 1997;Huybrechts and de Wolde, 1999;Greve, 2000;Marshall, et al., 2000;Abe-Ouchi, et al., 2007), with the exception of Tarasov and Peltier (2002), who used PDD factors dependent on the mean summer temperature.In reality, even if there is a clear correlation between the surface air temperatures and the ablation, there is no evidence of a single universal value for the degree-day factors.The melting rate of ice not only depends on mean temperature, but also on albedo and turbulence, and may be higher than 8 mm day 1 C 1 in cold regions, where the ablation may be underestimated up to 10% (Braithwaite, 1995).Nevertheless, in the absence of strong constraints linking climatic parameters and ablation, we still adopt the standard formulation of the positivedegree-day method.

The coupling procedure
The main challenge we are faced with in coupling CLIMBER 2.3 with GREMLINS is related to the major difference in their spatial resolutions.Therefore, to compute the surface mass balance of the ice sheets, the mean annual and summer surface air temperatures and the annual snowfall computed by CLIMBER must be distributed onto the finer GREMLINS grids through an appropriate downscaling procedure.
In this work, the procedure for computing the variations of these three variables with altitude has been modified compared to previous studies (Kageyama et al., 2004;Charbit et al., 2005).We simply compute the temperature profile for each grid box and each surface type by using the surface air temperature and free atmospheric lapse rate computed by CLIMBER, as described in Petoukhov et al. (2000).For the land ice surface type, an inversion effect is now added to this parameterization to compute the lapse rate above the ice sheets, following the results obtained by the regional model MAR on present-day Greenland (Fettweis et al., 2007).
where T CLIMBER represents the atmospheric temperature computed by CLIMBER at the altitude of a given GREM-LINS cell.This inversion is defined as an additional cooling T INV : where nday is the day of the year.The obtained downscaled temperatures compare favourably with the MAR climatology of Fettweis et al. (2007) and, for the LGM climate, to the PMIP results from GCM experiments.
The precipitation profile is obtained as a function of the altitude z taking into account the depletion in humidity with the decreasing temperatures obtained through the procedure described above: where z CLIMBER is the altitude of the CLIMBER grid point, and q SAT (T ) the saturated water vapor pressure for the temperature T (z) at the altitude z.
In our simulations, the ISM is called every 20 years and is run also during 20 years with the same climate forcing; in turn, the new ice-sheet geometry (altitude and surface), as well as the new land and/or ocean fraction computed by GREMLINS are then averaged on the CLIMBER grids through an aggregation procedure and provide new boundary conditions for the climate model.The freshwater flux coming from ice melting is released to the CLIMBER oceans and the freshwater feedback on climate is therefore accounted for.Meltwater is distributed in the oceans according to a runoff mask based on present-day topography.Conversely, if the ice sheets grow, the equivalent liquid water is subtracted from the runoff of the oceans.This is a new implementation compared to previous studies from Kageyama et al. (2004) and Charbit et al. (2005).The coastline is interactively computed at each time step on the basis of the grounded ice volume and the isostatic rebound.
The resulting coupled model has a fast computational time and, therefore, can be used to perform long-term simulations and to investigate the long-term behavior of the NH ice sheets through the last glacial-interglacial cycle.Our approach is a compromise between the need for a comprehensive climate model fully coupled to a refined 3-D ISM and the necessity of performing long-term transient simulations.Moreover, these experiments enable us to assess the feedbacks between climate and major ice sheets: we observe in a simplified way the impact of climate changes on simulated ice-sheets and, conversely, that of ice-induced feedbacks on computed climatic fields.

Experimental set-up
In the present work we simulate a complete climatic cycle, from the last interglacial period (126 kyr BP) until present time.In the baseline experiment (BSL), the model is forced by the variations of orbital parameters (Berger, 1978) and the atmospheric CO 2 concentration from the Vostok ice-core data (Petit et al., 1999).
Since the extent of the Greenland ice sheet during the Eemian (isotopic stage 5e, MIS5e) is not well constrained, the initial state of the ice sheet at 126 kyr BP (i.e. initial topography) is the same as the present-day one.However, the Eemian was warmer than the pre-industrial time in the Arctic and the Greenland ice sheet was probably smaller (Chapman et al., 2000;Marshall and Cuffey, 2000;Kukla et al., 2002;Otto-Bliesner et al., 2006).The initial climatic state is provided by a 5 kyr integration of CLIMBER alone under the 126 kyr BP insolation and atmospheric CO 2 concentration (274 ppm), with the present-day topography of Greenland given as a boundary condition.Both insolation and the atmospheric CO 2 concentration vary all along the period under study (Fig. 1 (Berger, 1978); (b) atmospheric CO 2 concentration (ppm) (Petit et al., 1999).
A critical problem to simulate a full glacial-interglacial cycle, and especially the inception of the ice sheets, is to obtain realistic precipitation patterns.In our model, the precipitation from a given CLIMBER grid box uniformly falls on the corresponding GREMLINS grid points.Therefore, due to the difference of resolution between CLIMBER and GREM-LINS, the simulated precipitation patterns often suffer from shortcomings, in particular over high altitude regions where most precipitations actually fall on the wind exposed slope of the first mountain range encountered by the air parcels.The comparison between the simulated present-day precipitations (PRC CLIMBER ) with those provided by the CRU climatology (PRC CRU ) shows that the climate model tends to significantly underestimate the precipitation over the areas extending from the eastern part of the Scandinavian Alps to the Barents-Kara sea region.Conversely, over Beringia precipitations are overestimated.To account for these discrepancies, we apply to the simulated CLIMBER precipitation a corrective factor PRC CRU /PRC CLIMBER (set to 2.75 for Scandinavia and 0.77 for Beringia) and we assume that this correction remains unchanged throughout the entire glacial-interglacial cycle. www.clim-past.net/5/329/2009/

Results
The objective of the baseline experiment is to examine whether the coupled model is able to simulate a realistic evolution of major NH ice sheets during the last glacialinterglacial cycle.In particular, we are interested in capturing the main features of glacial inceptions, expansions or retreats shown by the available records in terms of timing, geography and amplitude of sea-level variations.Our approach focuses on the large-scale processes and on the phase relationships between the external forcings and the periods of growth and retreat of major ice sheets.Due to the absence of representation of rapid ice flow dynamics in the ISM, the issue of rapid climate variability is not addressed in the present work.

Simulated ice sheet topography
Figure 2 shows the spatial distribution of the simulated ice sheets during the first phase of glacial inception, respectively at 122, 115 and 110 kyr BP.
At 122 kyr BP (Fig. 2) the model produces small ice caps in northwestern Canada, in the northern Rocky Mountains and over both the Ellesmere Island and the eastern part of the Baffin Island in the Canadian Archipelago.These simulated ice caps rapidly evolve into a much larger ice sheet at 115 kyr BP.Indeed, these areas represent the first simulated nucleation zones; immediately before ice inception, they are characterized by negative summer average surface temperatures (Fig. 3a).
At 115 kyr BP the simulated summer temperature is negative over most of Greenland, the Canadian Archipelago and the Rocky Mountains close to the Canada-Alaska border (Fig. 3), where we simulate a positive mass balance of more than 0.5 m/y (not shown).Continental ice spreads over most of Alaska's mainland, the northern Rocky Mountains and reaches the northwestern borders of Hudson-Bay.Its simulated thickness exceeds 2000 m in its central parts.We also simulate areas of negative summer temperature over Fennoscandia.A relatively large ice sheet appears between 60 • N and 70 • N in the Scandinavian Alps.In some places its thickness is higher than 2000 m; no continental ice is simulated in the Baltic region south of present-day Finland.
The total NH ice volume further increases at 110 kyr BP despite the fact that summer insolation is larger than at 115 kyr BP.At this period, the North American ice sheet now thickens to more than 3000 m in central Alaska, and occupies a large portion of the 60-70 • N latitudinal belt west of Hudson Bay (Fig. 2).Most of the Canadian Archipelago is icecovered, but the isolated island ice sheets do not coalesce into a larger ice complex, which is likely due to a crude representation of the ice advance over marine areas.At this time, the model also produces a thicker ice sheet over Eurasia, but its extent is significantly smaller than that of the North American ice sheet.Continental ice is mainly located over the Scandinavian Peninsula and does not spread south of 60 • N.
The timing of glacial onset is consistent with the reconstructed early phase of glaciation, as inferred from sea-level data (Camoin et al. 2001;Lambeck and Chappell 2001;Waelbroeck et al., 2002;Siddall et al., 2003;Bintanja et al., 2005), geomorphological observations (Boulton and Clark, 1990;Svendsen et al., 2004) and previous modelling studies (Tarasov and Peltier 1997;Berger et al., 1998;Kageyama et al., 2004;Calov et al., 2005a,b;Kubatzki et al., 2006).Nevertheless, observational data (Andrews and Barry, 1978) indicate that the regions of ice-sheet inception in North America were those bordering the eastern coast, and in particular Baffin Island, the Quebec-Labrador sector and the northeastern Keewatin region.The model correctly captures the glacial onset on Baffin and Ellesmere Islands, but does not produce ice accumulation over Quebec.The underestimation of ice volume in this region is correlated with surface air temperatures, as shown in Fig. 3. Indeed, at 123 kyr BP, before glacial inception, summer surface air temperatures over the Quebec-Labrador sector are mostly positive, between 5 and 15 • C.This may be linked to strong ocean circulation, since the simulated North Atlantic Deep Water intensity ranges between 23 and 25 Sv during this early period (not shown).Conversely, the model tends to overestimate ice growth on northwestern Rocky Mountains.Paleoenvironmental records indicate that, at the early phase of glacial onset, this region was likely as warm as today, or even warmer, and substantially ice-free (Clark et al., 1993).Indeed, the Cordilleran ice sheet does not seem to have appeared before the late isotopic stage 4 or 5 (i.e.∼75 kyr BP), when the Laurentide ice sheet (LIS) was high enough to displace the jet stream causing precipitations to fall over the Rocky Mountains (Roe and Lindzen, 2001).This is in disagreement with our results, because the coupled model simulates low summer temperatures in some regions of the Rocky Mountains at 123 kyr BP (down to −5 • C just before ice inception), triggering early ice sheet nucleation.Furthermore, precipitation is overestimated in the interior of the Rockies because, in the model, they do not represent a strong barrier to prevent precipitation from falling over this region.As a result, the simulated mean annual precipitation is greater than 1000 mm/yr at 123 kyr BP.This is also valid for present-day conditions, since the model overestimates precipitation by ∼150% compared to the CRU climatology in central areas of the Rockies.Glacial inception in the northwestern Rocky Mountains is not unusual in model studies and has been found by Kageyama et al. (2004) with a previous version of this model, and by Charbit et al. (2007) in GCM-driven simulations.The model also produces a glacial onset over Alaska, in disagreement with observations (Boulton and Clark, 1990).This is a common problem for coupled climate-ISM simulations (Deblonde and Peltier, 1991;Peltier and Marshall, 1995).Previous studies conducted with GCMs suggest that a more realistic modelling of vegetation over Beringia has a strong impact on local albedo and atmospheric circulation patterns (Wyputta and McAvaney., 2001).Refined vegetation indirectly leads to circulation changes that generate a warm anomaly over Alaska and Eastern Siberia, which, in turn, reduces the accumulation of grounded ice.Furthermore, previous works argue that the lack of a simulated Kuroshio Current may be also responsible for colder climatic conditions in these regions, thus fostering snow accumulation (Peltier and Marshall, 1995).Simulated ice inception and growth over Alaska and Beringia in the present study results from the combination of both coarse climate model resolution and from these regional and important features of the climate system, poorly represented in our coupled model.However, by using the CLIMBER-SICOPOLIS model, Calov et al. (2005a,b) suggest another possible explanation: they performed a transient simulation of the first phase of the glacial inception, from 126 kyr BP to 100 kyr BP, with an approach similar to the one adopted here, with the exception of the dust parameterization and the atmospheric CO 2 forcing inferred from Barnola et al. (1987).Their results do not show any ice build-up over North-western Alaska, whereas ice is produced in Eastern Siberia.On the contrary, when no dust is accounted for, Alaska is completely icecovered at 115 kyr BP (Calov et al., 2005b).At 110 kyr BP, geomorphologial reconstructions suggest the presence of an ice dome in the Quebec-Labrador region, and a coalesced ice sheet over the Canadian Archipelago (Boulton and Clark, 1990).The coupled CLIMBER-GREMLINS model does not produce any Quebec-Labrador dome due to warm North Atlantic (and Hudson Bay) SSTs simulated by the CLIMBER model, affecting the climate of nearby regions: at 50 • N the simulated mean summer SSTs range between 10 and 12 • C in northwestern Atlantic, whereas they range between 8 and 12 • C on Quebec's mainland at the same latitude.Moreover, the crude representation of the ice shelves does not favor the production of a single massive ice sheet over the Canadian Archipelago and the advance of ice over Hudson Bay, contrary to what is found in Calov et al. (2005a,b).The formation of an ice cover in this region in Calov et al. (2005a,b) may be due to a different parameterization of ice calving, which could aim at crudely reproducing the ice shelves.In the same way, our model does not simulate ice over the Barents-Kara sea region, as suggested by Svendsen et al. ( 2004), whereas continental ice over Scandinavia is properly reproduced.In addition to the different parameterization implemented to account for the impact of dust S. Bonelli et al.: Investigating the evolution of major Northern Hemisphere ice sheets on the snow/ice albedo and possibly for the representation of floating ice, most of the differences between Calov et al. (2005a,b) and the present study may come from the differences in the coupling procedures.We account for the inversion phenomenon and we adopted the PDD method for the computation of snow mass balance, whereas they adopted a surface energy mass balance interface scheme.Moreover, the two ice sheet models have a very similar physics.The major difference is that SICOPOLIS includes polythermal treatment in the heat equation.This gives the possibility of having temperate ice at the base of the ice sheet which may affect ice viscosity and deformation.This method is not included in GREMLINS but we believe that the impact of this feature is weak at the time of glacial inception compared to the uncertainty on sliding parameters.Observational data show that the 110 kyr BP period coincides with the maximum ice extent after the first glacial inception and is characterized by a major contribution from the LIS.A smaller ice sheet exists in the mountainous region of the Scandinavian Peninsula, but its contribution to NH ice volume is minor (Kleman et al., 1997).This feature clearly appears in our model results.
The spatial distribution of the simulated ice sheets from 110 kyr BP to present-day is plotted on Fig. 4 for six different key periods: 100, 75, 60, 30, 21 and 0 kyr BP.These snapshots refer to different stadial or interstadial states and correspond to pronounced phases of glacial expansion or retreat.
Between 110 and 100 kyr BP (Figs. 2 and 4), the ice sheet covering Alaska and the northwestern part of Canada increases, mainly due to the elevation of the ice thickness.Conversely, the southern part of Baffin Island is ice-free and no ice is simulated in this region at latitudes lower than 70 • N. At this time, observational data suggest that the LIS was partly melted (Andrews et al., 1983).The model does not produce any large ice cap in the Fennoscandian region.Only a small mountain ice sheet survives in high altitude regions of the Scandinavian Alps, but its extent is considerably smaller than the one simulated at 110 kyr BP, in agreement with the reconstructions from Mangerud et al. (2004).The simulation correctly produces the melting of the Fennoscandian ice sheet (FIS), as well as that of the southern part of Baffin Island.Nevertheless, ice thickness increases over northwestern Canada and Alaska due to the snow/ice albedo and elevation effects.The North American ice sheet is huge enough to sustain itself under increasing insolation (Fig. 2).
At 75 kyr BP, the simulated ice sheet covers Alaska's mainland and spreads further south, down to 50 • N (Fig. 4); a massive ice sheet is produced over the Canadian Northwestern Territories, reaching a thickness of more than 4000 m.Ice also occupies most of the islands of the Canadian Archipelago, including the southern part of Baffin Island.As suggested by observations (Kleman et al., 1997), no ice sheet builds up in non-mountainous areas of the Fennoscandian Peninsula.Indeed, consistently with reconstructions, the extent of the simulated FIS at 75 kyr BP (MIS5a) is smaller than the one corresponding to the 110 kyr BP period (MIS5d) and comparable to the 100 kyr BP snapshot (Kleman et al., 1997;Svendsen et al., 2004).
At 60 kyr BP, the model produces a large ice sheet west of Hudson Bay (Fig. 4).
This complex develops as far as ∼50 • N, spreading over most of the North American continent.The ice volume covering the Canadian Archipelago and Alaska has also slightly increased.Continental ice builds up in the Quebec-Labrador region, reaching a thickness of more than 2000 m in central areas and forming the second largest ice sheet in North America mainland.The model simulates the presence of a Keewatin ice dome in agreement with observational data (Boulton and Clark, 1990), but underestimates the Laurentide ice complex: the Labrador dome is smaller than that inferred from reconstructions and there is no ice over Hudson Bay.
The simulation of an ice sheet over the Scandinavian region is consistent with geomorphological reconstructions  (Svendsen et al., 2004).However, in contradiction with geological data, the Barents-Kara sea region remains ice-free (Mangerud et al., 2002;Svendsen et al., 2004).As previously mentioned, this discrepancy is likely related to the fact that the ice shelves are only crudely represented in the ISM.
Geomorphological data show that the 60 kyr BP period is characterized by the presence of large, persistent ice sheets over both North America and Eurasia (Boulton and Clark, 1990;Mangerud et al., 2002;Svendsen et al., 2004).With the exception of the above-mentioned regional discrepancies, the model reasonably produces this major feature.Indeed, the MIS5/MIS4 transition (∼75 kyr BP) outlines a pronounced switch of the climate-cryosphere system towards cold conditions, as inferred from both temperature reconstructions from Greenland ice cores (Jouzel et al., 2007) and sea-level reconstructions (Waelbroeck et al., 2002).The coupled model captures this transition and simulates the presence of continental ice sheets that do not substantially melt until the post-LGM period.Thus, the MIS5/MIS4 transition may be considered as the beginning of the true glacial phase.Glacial conditions are now widely spread and larger ice-induced feedbacks contribute to a further cooling of the system.Between 60 kyr and 30 kyr BP, the Eurasian ice volume does not significantly change, whereas the ice volume over the Quebec-Labrador sector and the Keewatin region has increased.
After 30 kyr BP, the simulated ice volume increases again until ∼21 kyr BP (LGM) (Fig. 4).The system switches towards more pronounced glacial conditions, similarly to the cold shift observed during the MIS5/MIS4 transition.At the LGM, the Labrador and Quebec regions are ice-covered and the LIS advances south of 50 • N. Consistently with reconstructions, the simulated ice masses completely cover the St. Lawrence Gulf, as well as the Canadian regions bordering the Atlantic Ocean (Dyke et al., 2002).However, geological reconstructions indicate the existence of two domes centered over the Keewatin region and the Quebec-Labrador plateau.The absence of a bi-domed ice sheet may be due to the ISM, which does not account for sediment deforma-tion and for the ice fast flow resulting from water-saturated sediments.According to Calov et al. (2002) and to Tarasov and Peltier (2004), this is a prerequisite to properly simulate a multi-domed ice surface topography.Another discrepancy with observations (Dyke et al., 2002) lies in the fact that a huge ice sheet is simulated over Beringia.
At the LGM, the coupled model produces negative summer temperatures over most of high latitude regions in the NH, with values as low as −25 • C over central parts of the Canadian plateau (Fig. 5).Over Alaska, summer temperatures are also negative, and range between −15 and −20 • C.These temperatures are largely affected by the albedo and altitude effects induced by the simulated ice sheets.The FIS spreads over the whole Baltic region and the British Islands, reaching a thickness of more than 3000 m in its central dome.The model also simulates the junction between the British Isles and the southern part of Scandinavia.However, it still underestimates the ice cover in the Barents-Kara sea region (Svendsen et al., 2004).Interestingly, a small mountain ice cap develops in the Alps, but its volume is negligible compared to major ice sheets.Iceland and the Faroe Islands are also completely ice covered.
The timing of the LGM is consistent with sea-level reconstructions (Waelbroeck et al., 2002;Siddall et al., 2003;Bintanja et al., 2005).Despite some discrepancies with observational data, often due to identified missing processes in the ISM or to the coarse spatial resolution of CLIMBER, the model satisfactorily reproduces the main features of LGM ice topography and simulates large ice sheets over both Eurasia and North America, in agreement with geomorphological investigations (Boulton and Clark, 1990;Mangerud et al., 2002;Dyke et al., 2002;Svendsen et al., 2004).
At 0 kyr (Fig. 4), the model produces the complete retreat of continental ice over most of the North American continent, as well as over Eurasia.The only exceptions are small residual ice sheets in northern Alaska and in mountainous areas of the Scandinavian Peninsula.With the exception of these regions, of Greenland and areas around the Canadian Archipelago, the simulated average summer surface air temperatures at the end of the transient simulation are positive.These temperatures are mostly consistent with observations, except over Alaska, where the model produces lower summer temperatures (Fig. 5).Low summer temperatures over Alaska, compared to other areas located at similar latitudes, are also simulated at 123 kyr BP, before the early phase of glacial onset (Fig. 3).The coupled model tends to underestimate summer temperatures in this region.

Ice-equivalent sea level change
The atmospheric CO 2 concentration and the summer insolation, as well as the corresponding simulated ice-equivalent sea level change due to major NH ice sheets, are plotted in Fig. 6.The computed ice-equivalent sea level is also directly compared to the relative sea-level curve from Waelbroeck et al. (2002) inferred from benthic isotopic records (Fig. 6b).It is important to note that the simulated sea level only refers to the contribution of ice changes in the NH and does not account for the Antarctic contribution (Philippon et al., 2006), whereas the curve from Waelbroeck et al. (2002) represents a global signal.
Sea-level change is directly dependent on ice sheet evolution.The ice-equivalent sea level inferred from our model simulation is computed as follows: We assume that the oceanic area is constant and equal to the present-day one since it does not significantly change throughout the last glacial-interglacial cycle.
A first drop of sea level is observed between 122 kyr and 110 kyr BP and is correlated with a sharp decrease of summer insolation at 65 • N.This timing is consistent with the sea level reconstruction from benthic records (Waelbroeck et al., 2002).At 110 kyr BP, the simulated seal level has decreased by 28.5 m, that is about 15 m higher than the one obtained by Waelbroeck et al. (2002).This is consistent with the fact that no ice is simulated over the Labrador-Quebec sector and the Barents-Kara region (Fig. 4), in contradiction with the observations (discussion in Sect.3.2.1).The presence of an ice sheet covering a large part of Alaska (Fig. 4) increases the simulated sea level drop, but does not compensate the lack of ice over these regions.However, as previously mentioned, the sea level reconstruction (Fig. 6b) corresponds to a global signal (Waelbroeck et al., 2002), whereas the sea level simulated in our experiment only accounts for the contribution of NH ice sheets.A potential contribution from the Antarctic ice sheet is thus not considered.Indeed, Duplessy et al. (2007Duplessy et al. ( , 2008) ) suggest that the West Antarctic ice sheet was partly destabilized during the last interglacial period and that its further expansion likely contributed to the sea level drop recorded after the Eemian.
Between 110 and 100 kyr BP, the ice-equivalent sea level evolution simulated by the coupled model slightly increases by 6.5 m (Fig. 6c).This corresponds to the disappearance of the FIS (Fig. 4).During this period a rapid increase of summer insolation at 65 • N is observed and reaches one of the maxima of the last 126-kyr (Figs. 1 and 6a), whereas the atmospheric CO 2 undergoes a first significant decrease from 250 to 225 ppm.A new sea level stand is observed at 90 kyr BP (−38.5 m), then followed by a slight increase until ∼75 kyr BP.
During this first period of the glacial cycle (126−75 kyr BP), the evolution of the simulated ice sheets appears therefore to be correlated with the summer insolation signal, suggesting that this is the major factor triggering glacial expansion or retreat of the ice sheets.
A significant development of continental ice sheets is simulated after 75 kyr BP, causing a further sea level decrease (−64 m at 65 kyr BP).This drop is correlated with a significant decrease of summer insolation at 65 • N, combined to low atmospheric CO 2 concentrations (190 ppm at 65 kyr BP).The decrease of both forcing factors makes the system switch towards colder conditions, characterized by a rapid sea-level drop (∼24 m in less than 15 kyr in our simulations) and by the existence of long-lasting ice sheets over both Eurasia and North America (Fig. 4), which corresponds to the actual beginning of the full NH glaciation (e.g.Jouzel et al., 2007).The newly re-formed FIS persists even under increasing summer insolation after ∼70 kyr BP (Fig. 4), which proves that the atmospheric CO 2 concentration is low enough to counterbalance the insolation radiative effect.As shown for the first phase of glacial inception, the FIS may build up in presence of low summer insolation and relatively high atmospheric CO 2 , but the associated positive feedbacks due to increased albedo and surface elevation are not efficient enough to prevent it from melting under unfavorable orbital parameters (i.e.increasing insolation).This sequence of glaciation events shows that the decrease of summer insolation alone is sufficient to explain the expansion of grounded ice over North America and Eurasia.Nevertheless, the FIS needs low atmospheric CO 2 concentration to survive to subsequent increases of insolation.This explains why the simulated FIS can sustain itself only after ∼75 kyr BP and until the LGM (Fig. 7).
The sea level drop between 60 and 30 kyr BP is less pronounced than that occurring between 75 and 65 kyr BP.During this period the simulated sea level evolution is smoother than the one inferred from benthic foraminifera records (Waelbroeck et al., 2002) and no abrupt expansion or regression of the ice volume is reproduced.The simulated ice sheets are large enough to counterbalance, via the ice albedo and elevation effects, the relatively high values of insolation prevailing during the 50-30 kyr BP period.Moreover, the maintenance of the ice sheets is also supported by rather low CO 2 levels (200 to 220 ppm).The coupled model is not sensitive enough to decreasing insolation when summer values remain higher than 420 W/m 2 and the simulated ice volume remains stable; some improvements are therefore needed to reproduce rapid ice volume variations during this period.Adding a more refined representation of the ice shelves and of basal hydrology processes in the ISM could help produce results better fitting to data.Around 30 kyr BP we observe a further pronounced expansion of main NH ice sheets in response to a new decrease of both insolation and atmospheric CO 2 concentration.This expansion is associated to a sea level decrease of ∼39 m in less than 10 kyr.At the LGM, the simulated ice equivalent sea level stand is −121.5 m (Fig. 6c).The timing of the maximum sea level drop is consistent with reconstructions (Waelbroeck et al., 2002;Siddall et al., www.clim-past.net/5/329/2009/Clim.Past, 5, 329-345, 2009Past, 5, 329-345, 2003)).The simulated LGM ice volumes of the LIS and the FIS amount to 41.8×10 15 m 3 and 5.3×10 15 m 3 , respectively (Table 2).These results are in agreement with the findings from Lambeck et al. (2000); on the contrary, the simulated LIS is larger than the one described in the ICE-5G reconstruction (Peltier, 2004), whereas the simulated FIS volume is lower.At the LGM, the Antarctic contribution to sea-level change computed by Philippon et al. (2006) with the same climate model coupled to a 3-D Antarctic ISM (Ritz et al., 2001) ranges between 9.5 and 17 m.Therefore, we can estimate the simulated ice-equivalent sea-level at the LGM to range between −131 m and −138.5 m.These values lie within the upper bound of the most commonly accepted reconstructions, ranging between ∼−110 and −143 m at the LGM (Yokoyama et al., 2001;Lambeck and Chappell, 2001;Camoin et al., 2001;Waelbroeck et al., 2002;Siddall et al., 2003;Bintanja et al., 2005).Compared to previous modeling works (Tarasov and Peltier, 1997;Zweck and Huybrechts, 2005), the simulated LGM sea level is ∼10 m lower due to the presence of a large ice sheet covering Alaska/Beringia, whereas this estimation is consistent with the findings by Abe-Ouchi et al. (2007).

Experimental set-up
To better investigate the effect of both atmospheric CO 2 concentration and insolation on major NH ice sheets, we perform a set of sensitivity tests.These experiments are conducted as previously described for the baseline (BSL) experiment, but under constant insolation or constant atmospheric CO 2 concentration (Table 1).For all the experiments, the dust parameterization is implemented as described in Sect.2.1 for the BSL run, and it is indexed on the CO 2 concentration inferred from Petit et al. (1999).The INSO-126 experiment is run under constant 126-kyr BP insolation and pCO 2 evolution inferred from the Vostok ice core (Petit et al., 1999).In this test, summer insolation at 65 • N is equal to 452 W/m 2 and is therefore particularly unfavorable for glacial onset.The sensitivity tests CO2-280, CO2-250, CO2-235 and CO2-200 are performed by driving the coupled CLIMBER-GREMLINS model with BSL-insolation (Berger, 1978) and constant CO 2 concentration fixed to 280, 250, 235 and 200 ppm, respectively.These tests are designed to investigate a potential CO 2 threshold effect on the evolution of major NH ice sheets and their different sensitivity to external forcings.

Results
The temporal evolution of the simulated ice volumes is plotted in Fig. 7 for all the transient sensitivity experiments performed in this study.These results are compared with those obtained in the baseline (BSL) experiment.In the CO2-280 test (Fig. 7), the simulated LIS ice volume during the early phase of the glacial inception (126−110 kyr BP) does not significantly differ from the one produced in the BSL experiment (Fig. 7).This confirms that the early glacial build-up is primarily driven by decreasing summer insolation (Figs. 1 and 6a) since the model still produces ice development over North America and Eurasia for a CO 2 level of 280 ppm, typical of interglacial periods.Conversely, after 110 kyr BP the CO2-280 ice volume remains significantly lower than the BSL one.A similar result is obtained in the CO2-250 experiment, although the ice volume is slightly larger than in CO2-280, but still lower than in the baseline experiment (28.2×10 15 m 3 vs.52.5×10 15m 3 for the NH ice sheets at the LGM, Table 2).
When the atmospheric CO 2 further decreases to 235 ppm (CO2-235), half-way between typical glacial and interglacial values (Petit et al., 1999), the LIS build-up is faster and slightly larger than in BSL until 65 kyr BP, but remains lower than in BSL during the 65−21 kyr BP interval.For this latter period the 65 • N summer insolation is higher than 420 W/m 2 until ∼30 kyr BP.This value is too high to trigger any significant ice accumulation when the atmospheric CO 2 concentration is fixed to 235 ppm (or higher), and the growth of the ice sheets stops and stabilizes between 65 and 30 kyr BP in tests CO2-280, CO2-250 and CO2-235 (Fig. 7).Throughout this period, the atmospheric CO 2 concentration in BSL actually varies between 190 and 220 ppm, thus contributing to a further cooling of the climate system.This cooling adds up to the positive feedback effects due to existing ice sheets and results in a larger LGM ice volume in BSL.On the contrary, for constant CO 2 concentrations higher than 235 ppm these feedbacks only prevent the LIS from melting, but they are not sufficient to trigger a further glacial expansion until the pronounced decrease of summer insolation observed after 30 kyr BP.
In the CO2-200 experiment, the CO 2 concentration is low enough to produce a large glacial inception over North America and Eurasia.At 110 kyr BP, the North American and Fennoscandian ice volumes are 49% and 100% larger than those simulated in BSL.At the LGM the simulated LIS and FIS ice volumes are respectively 53.5×10 15 m 3 and 13.6×10 15 m 3 , compared to 41.8×10 15 m 3 and 5.3×10 15 m 3 in the BSL experiment.
Although the behavior of both ice sheets appears to be quite similar, the excess of ice simulated over Fennoscandia under a 200 ppm CO 2 level, compared to the BSL test, is considerably larger than the ice excess obtained over North America (157% for FIS, 28% for LIS at the LGM).Moreover, a stable, long-lasting FIS is only formed after ∼75 kyr BP in BSL and in the CO2-200 experiments.
These results show that a threshold behavior appears in our simulations: when the atmospheric CO 2 decreases from 235 ppm to 200 ppm, the effect of orbital configurations favorable for glacial advance (Berger, 1978) is amplified.The lower CO 2 values also provide a further cooling, necessary to prevent the FIS from melting during the 65−30 kyr BP interval and resulting in the build-up of a large European ice complex at the LGM.Conversely, higher CO 2 concentrations prevent from significant ice growth over Fennoscandia at the LGM in experiments CO2-235, CO2-250 and CO2-280.Therefore, the Fennoscandian ice sheet appears to be more sensitive to atmospheric CO 2 concentration than the North American ice sheet.This confirms our conclusions obtained from the analysis of the relationship between external forcings and the evolution of ice volume in the baseline experiment (Sect.3).A possible explanation may be represented by the influence of the Atlantic Ocean on regional climate, which contributes to mitigate the temperature of nearby regions such as Europe and the Quebec-Labrador sector (Fig. 3c for the BSL run).Therefore, to trigger an extensive, long-lasting glaciation, a further cooling due lower CO 2 concentrations needs to be added to the orbital forcing.Conversely, with the exception of the Quebec-Labrador sector, the American ice sheet is less influenced by the presence of the Atlantic Ocean and is characterized by a more "continental" climate: in central areas, the effect of the solar forcing is not buffered by ocean mitigation, and surface temperatures are more strongly dependent on insolation changes.

Discussion: comparison with previous modeling studies
The approach adopted here is similar to the one proposed in previous studies for the last glacial inception (Kageyama et al., 2004) and (Calov et al., 2005a,b).Here, we produce a higher NH ice volume compared to Kageyama et al. (2004), as well as the early formation of the FIS, in better agreement with data reconstruction.These differences are due to the improved coupling procedure between CLIMBER and GREMLINS (Sect.2.3); in particular, the parameterization of the inversion phenomenon favors the build-up of NH ice sheets during the first phase of glacial inception (not shown).
Both differences in the coupling procedure and in the dust parameterization (Sect.2.1) likely explain the differences in simulated ice geometry with Calov et al. (2005a,b), as discussed in Sect.3.2.1.Regarding the dust parameterization, its major effect on the simulated ice evolution is to produce a lower NH ice volume due to the decreased albedo effect.Nevertheless, for the early phase of glacial inception, its impact on the simulated ice volume is less than in Calov et al. (2005b) because the atmospheric CO 2 concentration is close to its pre-industrial level and dust weight is therefore negligible.Indeed, at 110 kyr BP, our model produces a total NH ice volume of 14.5×10 15 m 3 , vs. 14.7×10 15 m 3 when dust is not accounted for (i.e. the simulated NH ice volume is 1.5% larger).The effect of the dust parameterization becomes more important after 70 kyr BP, since it prevents the formation of a large ice sheet over Siberia (35×10 15 m 3 at the LGM), in agreement with the results of previous works (Krinner et al., 2006).
Our study is based on transient, fully coupled climate-3D ISM simulations, which are now run throughout the full glacial-interglacial cycle.This approach extends the time period covered by the analysis to the 100-kyr ice age cycle, similarly to the previous works by Gallée et al. (1992), who was the first to investigate this time period with a coupled climate-ISM, Berger et al. (1998Berger et al. ( , 1999) ) and Tarasov andPeltier (1997, 1999).Nevertheless, our modelling tool is different, since Gallée et al. (1992) and Berger et al. (1998Berger et al. ( , 1999) ) reproduced past glacial-interglacial transition with a climate model of intermediate complexity coupled with a 2-D ISM (LLN 2D-NH), whereas Tarasov and Peltier (1999) used an energy balance model coupled with a vertically integrated ISM (1997) or with a 3-D-ISM (1999).The coupled CLIMBER-GREMLINS model enables us to include major climate-ice feedbacks in the analysis, and to describe the main characteristics of the long-term variations of the climate system over the last 126 kyr.Our results are generally consistent with these previous studies, thus confirming their robustness.In particular, Tarasov and Peltier (1997) pinpoint the stronger climate sensitivity of the Eurasian ice sheet.Furthermore, they show that the radiative impact of varying atmospheric CO 2 concentration is critical to achieve adequate glaciation.Indeed, as highlighted by our sensitivity www.clim-past.net/5/329/2009/Clim.Past, 5, 329-345, 2009 tests and by previous studies (Berger et al., 1998(Berger et al., , 1999)), the atmospheric CO 2 concentration affects the intensity and regional distribution of glaciation events, whereas ice inception is triggered by changes in summer insolation.The coupled model does not produce any glacial onset in the INSO-126 sensitivity test, in agreement with Loutre and Berger (2000) and Calov et al. (2005b).
The nonlinearity associated with internal climate feedbacks, driven by insolation and atmospheric CO 2 concentration, can describe the dynamics of the 100-kyr ice-age cycle (Tarasov and Peltier, 1997).Indeed, the effect of feedback processes between ice sheets and atmosphere is crucial (Abe-Ouchi et al., 2007).Furthermore, using the GREMLINS ISM forced by GCMs outputs, Charbit et al. (2007) showed the importance of the climate-ice sheet coupling the final ice evolution and the strong sensitivity of the ice sheets to the climate forcing.All these features need to be accounted for to reproduce a full glacial-interglacial cycle.

Summary and conclusions
The CLIMBER-GREMLINS model reproduces the main phases of advance and retreat of Northern Hemisphere ice sheets during the last glacial cycle, although the amplitude of these variations is less pronounced than those based on global sea level reconstructions.At the LGM, the simulated ice volume is 52.5×10 15 m 3 and the spatial distribution of both the American and Eurasian ice complexes is in reasonable agreement with observations, with the exception of the marine parts of these ice sheets.
We investigate the responses of the two major NH ice sheets to atmospheric CO 2 concentration and insolation, as well as their sensitivity to these forcings.Our simulations suggest that the Laurentide ice sheet has a different dynamics of development compared to the Fennoscandian ice sheet, since they do not respond with the same sensitivity to atmospheric CO 2 concentration.We simulate the early build-up of continental ice over North America as soon as the summer insolation decreases at 122 kyr BP.The evolution of this ice complex is largely driven by the insolation forcing throughout the whole 126-kyr cycle.Conversely, the model produces glacial inception on the Fennoscandian Peninsula in the early phases of glaciation, but this early ice sheet is not stable enough to survive to following periods of increasing summer insolation.Indeed, both weak summer insolation and low atmospheric CO 2 concentration are necessary to trigger a long-lasting glaciation over Eurasia.Consistently with reconstructions (Svendsen et al., 2004), a long-lasting FIS only appears after 75 kyr BP in correspondence of the MIS5/MIS4 cold transition.
As proposed by Milankovitch (Milankovitch, 1941), changes in summer insolation trigger the ice sheets inception.Nevertheless, variations in atmospheric CO 2 concentration, as well as climate-ice sheet feedbacks, are important mecha-nisms affecting the duration and intensity of cold/warm periods and the extent and altitude of major NH ice sheets.Sensitivity tests conducted under standard insolation and for various constant CO 2 levels demonstrate that the FIS is more sensitive to the atmospheric CO 2 concentration than the LIS.In our simulations, the build-up of a stable FIS responds to a threshold behavior to CO 2 concentration and its development is drastically different when the atmospheric CO 2 decreases from 235 to 200 ppm.
This work also underlines the limits of the coupled model.These are either due to the CLIMBER coarse spatial resolution, unable to capture some important regional features of the climate system (i.e.influence of the Rockies on the precipitation over the American continent and the existence of the Kuroshio Current in North Pacific), or to identified missing processes in the ISM, such as the absence of rapid ice flow dynamics and sediment deformation.The addition of these mechanisms in the ISM should improve the spatial distribution of simulated ice over marine regions such as Hudson Bay or the Barents-Kara sea sector.
); the lowest value of summer insolation (JJA) at 65 • N occurs at 113 kyr BP (395 W/m 2 ) and the highest ones are reached at 125 kyr BP and 101 kyr BP (>450 W/m 2 ).

Fig. 3 .
Fig. 3. Simulated average summer surface air temperatures (JJA) at different key periods (in • C): 123 kyr BP (left panel), 115 kyr BP (middle panel) and 70 kyr BP (right panel).These periods represent the starting phase of pronounced cold episodes characteristics of the last climatic cycle.

Fig. 5 .
Fig. 5. Simulated average summer surface air temperatures (JJA) (in • C) at the LGM (left panel) and at the end of the simulation (0 kyr, middle panel).The right panel shows the present-day average summer surface air temperatures inferred from observations (CRU climatology).

Table 1 .
Overview of model 126 kyr-long transient experiments.

Table 2 .
Simulated LGM volumes for all major NH ice sheets.Ice volumes are expressed in [10 15 m 3 ].