Sensitivity of atmospheric forcing on Northern Hemisphere ice sheets during the last glacial-interglacial cycle using output from PMIP3

Abstract. We use the three-dimensional Parallel Ice Sheet Model (PISM) to simulate Northern Hemisphere ice sheets evolution through the last glacial-interglacial cycle. The simulation is driven by the NGRIP δ18O index combined with climate forcing at two time slices, the Last Glacial Maximum (LGM) and present day (PD). In order to investigate the sensitivity of the ice sheets to the atmospheric forcing, atmospheric output from nine climate models from the Paleoclimate Modeling Intercomparison Project Phase III (PMIP3) are used to force the ice sheet model with the same set-up. The results show large diversity in simulated ice sheets between different models. By comparing the atmospheric forcing, we found that summer surface air temperature pattern resembles the ice sheet extent pattern at the LGM, which shows great sensitivity to summer surface air temperature. This implies that careful constrains on climate output is essential for simulating reliable glacial-interglacial Northern Hemisphere ice sheets. The ablation process is of vital importance for high-latitude Northern Hemisphere ice sheets. Besides, the absent nonlinear interactions between ice sheet and atmosphere and ocean, which have different signals regionally, also contribute to the mismatches between simulated ice sheets and geological evidences. Hence, we highlight the needs for coupling an ice sheet model to GCM to take into account these missing processes.


Introduction
During the late Pleistocene, climate went through vast changes known as glacial-interglacial cycle with a periodicity of around 100,000 years (100 kyr), accompanied by large extent of Northern Hemisphere ice sheets advancing and retreating.Milankovitch theory, which is that insolation change resulting from orbital changes (eccentricity, obliquity and precession), is widely accepted as the main driving and modulating factor for glacial-interglacial cycles (Milankovitch, 1941;Hays et al., 1976).However, orbital forcing alone can not explain the strong 100 kyr cycle of Northern Hemsiphere ice sheets, which have larger amplitude, slower build up and faster retreat process.This indicates that internal climatic feedbacks acting as nonlinear amplifiers are also of vital importance (Imbrie et al., 1993;Huybers and Wunsch, 2005;Lisiecki, 2010;Huybers, 2011).
The interactions between ice sheets and other climate systems can amplify, pace or drive global climate change on different time scales (Clark et al., 1999).Such processes include the surface-albedo feedback, ocean circulation feedback and Glacial Isostatic Adjustment (GIA).For Northern Hemisphere ice sheets with a large proportion of land based ice, the atmospheric forcing plays an important role to the mass balance variation.Snow accumulation contributes to positive mass balance while surface ablation, basal melting and calving lead to negative mass balance.At the same time, ice sheets change albedo and surface topography, which in return can influence the heat transport and the atmospheric circulation, and result in changing temperature and precipitation (Abe-Ouchi et al., 2007;Roe and Lindzen, 2001;Liakka et al., 2012;Löfverström et al., 2014Löfverström et al., , 2015)).Accurate reconstruction of Northern Hemisphere ice sheets for both extent and geometry is essential for fully understanding the physical mechanisms of ice sheets themselves and their influence on global climate change.
Generally there are two ways of ice sheet reconstruction: geologically based and dynamically based.By using radiocarbondating, geomorphological features, relative sea level observations or other types of geological data, Northern Hemisphere ice sheet dynamics are well established since the Last Glacial Maximum (see Fig. 1; Dyke, 2004;Hughes et al., 2016;Margold et al., 2015;Gowan et al., 2016) but are poorly constrained for pre-LGM, with relatively vague outline at a wide range of times (Kleman et al., 2010;Svendsen et al., 2004).The geometry, volume and exact timing of ice sheet evolution are hard to be inferred from the geological record.An alternative way to fully reconstruct ice sheets is the numerical ice sheet modeling, which requires proper boundary conditions and essential ice sheet feedback mechanisms.
Coupling an ice sheet model to a sophisticated General Circulation Model (GCM) is an effective way to produce numerical ice sheet reconstructions.It is computationally expensive to do this due to the long time period of simulations.Alternatively, an approach of coupling ice sheet models to Earth system Models of Intermediate Complexity (EMICs, Claussen et al., 2002) has been used (e.g., Charbit et al., 2005;Ganopolski et al., 2010;Bauer and Ganopolski, 2017), while the missing processes in EMICs may also have large effects on regional ice sheet distributions.Another approach for long-term ice sheet simulations is by using time slice simulations at the LGM or present day, and interpolating between them by using the 18 O time series record from ice cores, which is called the "glacial index" method (e.g., Greve et al., 1999;Marshall et al., 2000Marshall et al., , 2002;;Charbit et al., 2002Charbit et al., , 2007;;Rodgers et al., 2004).The interactions between ice sheet and other climate components can not be investigated by using this method.
Fully understanding of the climatic effect on ice sheets and the feedbacks between ice sheets and climate is the key for not only simulating ice sheet evolution correctly but also for understanding the physical mechanisms behind it.In this study, we focus on the atmospheric effect on Northern Hemisphere ice sheets during the last glacial-interglacial cycle.The model set-up is presented in Sect. 2. By using the "glacial index" method, firstly we simulate Northern Hemisphere ice sheets evolution with one GCM output.We tuned the model to match the past sea level record at the LGM.Comparison between the model results and geological evidence is conducted to investigate the ability of GCM to reproduce reasonable ice sheet configurations (Sect. 4.1).Then we test the uncertainties linked to the atmospheric forcing.By using the same ice sheet model set-up but different climate forcing from different GCMs (The Paleoclimate Modeling Intercomparison Project Phase III, PMIP3 models, Sect.3), we examine the sensitivity of ice sheets to the forcing (Sect.4.2).

Model set-up
We use the Parallel Ice Sheet Model (PISM, version 0.7.3), which is an open source, three-dimensional thermo-mechanically coupled shallow ice sheet model (the PISM authors, 2016;Aschwanden et al., 2012;Bueler and Brown, 2009;Winkelmann et al., 2011)  The model run starts from the Last Interglacial, 122.9 kyr Before Present (BP) (see Sect. 2.2), and ends at 0 kyr BP.The initial condition of ice sheet and topography configuration is from the present day state.

Ice dynamics
The stress balance is computed by using a combination of shallow ice approximation (SIA) and shallow shelf approximation (SSA).SIA and SSA are dominant in grounded and floating ice respectively.Other regions such as ice streams, with significant internal ice deformation and basal sliding are solved by combing the velocity solution of the two approximations (Winkelmann et al., 2011;Bueler and Brown, 2009).
A non-dimensional enhancement factor for both SIA and SSA is added to the flow law.Following the recommendation from Cuffey and Paterson (2010), we set the SIA enhancement factor E SIA to 5. The values of the other parameters used in our simulation and related references are summarized in Table 1.Surface gradients are computed by finite differences to determine the driving stress.The conservation of energy is solved within the ice, the subglacial layer and a layer of thermal bedrock.A geothermal heat flux input is included at the lower boundary (Aschwanden et al., 2012;the PISM authors, 2016).

The subglacial dynamics
The basal resistance to ice flow is computed based on the hypothesis that a layer of till is underlain the ice sheet (Clarke, 2005;Bueler and Brown, 2009).A pseudo-plastic sliding law is used for determining sliding: where ⌧ b is the basal shear stress, ⌧ c is the yield stress, u is the sliding velocity and u threshold is a parameter called threshold velocity and q is the pseudo-plastic sliding exponent (Table 1).
The Mohr-Coulomb criterion (Cuffey and Paterson, 2010) is used to determine the yield stress, which is related to the till material property (the till friction angle ) and the effective pressure on the saturated till N till .
= 30 is a typical value from lab experiments (Cuffey and Paterson, 2010, p. 268).The till cohesion value, c 0 , is set to 0 (Schoof, 2006).The effective pressure N till is determined by the following parameterization: This is based on laboratory experiments with till extracted from an ice stream in Antarctica (Tulaczyk et al., 2000a).P o is the ice overburden pressure, W till is the effective thickness of water in the till, W max till is the maximum amount of water in the till, e 0 is the till reference void ratio, C c is the till compressibility coefficient and is the effective fraction overburden pressure in the till.The water in the base is not conservative in our simulations, it is stored in the till up to the maximum thickness (W max till ).For the water exceeding that thickness, it is drained off instantaneously (Tulaczyk et al., 2000b;Bueler and Brown, 2009).
The ice sheet evolution is accompanied by the deformation of the solid earth (GIA).We use the Lingle and Clark method (Lingle and Clark, 1985;Bueler et al., 2007) to calculate it.

The ice shelf dynamics
The model for marine portions of ice sheets is from the PIK (Potsdam Institute for Climate Impact Research) component of PISM (Albrecht et al., 2011;Levermann et al., 2012;Winkelmann et al., 2011).The flotation criterion for determining whether the ice floating or grounded is combined with a temporal evolution mask, with an input of relative sea level change.The ice shelf calving mechanism is controlled by two schemes.The calving rate is a function of the horizontal strain rates.The ice shelf is automatically removed when the ice thickness is thinner than the threshold (H calthres ) of 200 m.
For the boundary condition at the base of the ice shelf, we follow the setup from Martin et al. (2011).The ice at the basal boundary is at pressure melting temperature.The mass flux from ice shelf to ocean is related to the heat flux Q heat between ocean and ice (Martin et al., 2011, Eq. 3-5), which is: ⇢ oce is density of ocean water.c poce is specific heat capacity of ocean mixed layer.T is the thermal exchange velocity.T f is a virtual temperature which represent the ocean water freezing temperature at different depths.The ocean temperature, T oce , is set to a constant value of -1.7 C (Beckmann and Goosse, 2003).F melt is a dimensionless model parameter and we set to 1 ⇥ 10 2 to increase the melt rate.

The surface mass balance
The surface mass balance is computed by the positive degree-day scheme (Reeh, 1989).Monthly mean temperature and precipitation are used as inputs.For accumulation, precipitation when temperature is below 0 C is considered to be snow, and temperature above 2 C is considered to be rain.With the temperature in between, the percentage of snow and rain is linearly interpreted.For ablation, the positive degree-day (PDD) is calculated (Calov and Greve, 2005, Eq. 6): A "white noise" variation is added onto the seasonal cycle temperature to account for synoptic variability.The standard deviation is 5 K.For the other model parameters, see Table 1.

Climate inputs
The simulation is driven by a "glacial index" I(t) combined with the climatic monthly near-surface air temperature (T mon ) and precipitation (P mon ) fields at two time slices, the LGM (21 kyr BP) and present day (PD).The "glacial index" is derived from the North Greenland Ice Core Project (NGRIP) 18 O 50 year-average records (Andersen et al., 2004, 122.95 kyr BP -0 BP, Fig. 2).The value of I is 1 for LGM (I LGM = 1, t = 21kyrBP ) and 0 for PD (I P D = 0, t = 0), which are the full glacial condition and interglacial condition respectively, and linear interpolated in between: The present day surface air temperature fields are from NCEP/NCAR Reanalysis long term monthly mean datasets (1981-2010, 2.5 ⇥ 2.5 , Fig. 3c-d) (Kalnay et al., 1996).Precipitation fields are from GPCP long term monthly mean datasets (1981-2010, 2.5 ⇥ 2.5 , Fig. 4c-d) (Adler et al., 2003).For LGM, the temperature and precipitation is from the Earth System The data is bilinearly interpolated onto the model grids.
According to the GCM simulation, the surface air temperature (SAT) at the LGM is colder than PD globally expect for some areas in the southern part of Eurasian continent (Fig. 3e-f).Due to higher elevation and albedo, the SAT where ice sheets existed at the LGM can be lower by up to -36 C from PD both in winter and summer.For precipitation (Fig. 4), there was less precipitation in the eastern part of North America and Scandinavia in winter, while there was more precipitation in the northwestern part of North America and southwest of Great Britain.In summer, there is less precipitation in northern North America and more precipitation in southern North America at the LGM than PD.
Paleoclimate fields are based on the linear relationship between PD and LGM: This approach is similar to the one previously used in Greve et al. (1999) and Marshall et al. (2000).Equation ( 9) is used to avoid negative precipitation.Equation ( 10) is the precipitation correction due to surface elevation (H) change, based on the exponential relationship between water vapor saturation pressure and temperature in the upper atmosphere. is a tunning parameter (unit: km 1 ) used for matching the far field sea level record at the LGM, which is set to = 0.75.A slight modification is made for present day precipitation to eliminate the error caused by the precipitation elevation correction: The scalar relative sea level reconstruction we use for the land-sea mask is taken from the eastern Mediterranean (Rohling et al., 2014), see Fig. 2b.

Basal inputs
The basal topography data we use is ETOPO1, which integrates land topography and ocean bathymetry from numerous global and regional datasets, and has a resolution of 1 arc-minute (Amante and Eakins, 2009).The basal geothermal heat flux data is from Davies (2013), with 2 ⇥ 2 equal area grid, which is a combination of data measurements and correction methods (Fig. S1).All the basal input data is bilinearly interpolated onto the 20 km model grids.

PMIPmodel comparison experiment set-up
The model set-up described before is based on the atmospheric forcing from one GCM output (COSMOS-AWI).While GCMs vary between each other and contain model deficiencies.In the experiment we focus on the forcing uncertainties related to GCM uncertainties.In order to test how sensitive simulated ice sheets are to climate forcing, we ran simulations by using the same ice sheet model parameters described in Sect.2, with the only difference that the atmospheric forcing at the LGM is from different GCM output.
In winter, generally all the models have similar large-scale surface air temperature patterns at the LGM, but differ regionally (Figs.S2, S5).In GISS-E2-R, northern North America is less cold and northern Siberia is much colder.In summer, the temperature differences between different model results are even more extensive (Fig. 9, Fig. S6).A large warmer area over the Comparing with the present day pattern, there is less precipitation in eastern North America and western coastal areas at the LGM in winter (Fig. 10).What should also be noticed is that more precipitation in the Canadian Arctic region and western Laurentide region are simulated (Fig. 10).In summer, precipitation is less at the LGM than present day in the whole inland area of northern North America (Fig. S7).We have also compared the individual ice volume changes in the Greenland, Eurasian and North American ice sheets (Fig. 5b).The sea level contribution of different ice sheets to the total sea level drop are different during different marine isotope stages.For MIS 5e, the Greenland Ice Sheet (green line) is the main contributor to sea level fall, and after that remained relatively stable with 10 m SLE until the deglaciation.From MIS 5d, North American ice sheets (blue line) started to build up and make up to 2/3 of total SLE during that stage.The Eurasian ice sheets (black line) build up is restricted before MIS 4 with a SLE of less than 10 m, then peaked at 60 kyr BP with 26 m SLE, and is coincident with the timing of North American ice sheets (which had 64 m SLE).Compared to North America, the amplitude of ice volume change is smaller during Dansgaard-Oeschger events.The maximum sea level drop of Eurasian ice sheets was during MIS 2 with 30 m, and around 80 m for North American ice sheets.At 15 kyr BP, North American ice sheets were slightly larger than at 20 kyr BP, while the Eurasian ice sheets were smaller than before.After that, both the North American ice sheets and Eurasian ice sheets retreated rapidly.The North American ice sheets grew significantly during the Younger Dryas with around 6 m sea level drop while the Eurasian ice sheets did not advance significantly.

Spatial distribution of ice sheets
The spatial distribution of simulated ice sheets through the last glacial-interglacial cycle is also investigated.Snapshots of the spatial distribution at different key periods are shown in Fig. 6.Consistent with the SLE time series shown in Fig. 5, the extent of Greenland Ice Sheet was smaller than present day during MIS 5e (Fig. 6).As the temperature index goes down, ice sheet that was around 500 m thick started to build up along the northeast coast of North America, in the region of Baffin Island and Labrador-Quebec sector.It advanced westward into the Interior Plains during MIS 5d, when Cordilleran Ice Sheet and the Scandinavian Ice Sheet also started to build up with flat, thin ice.During MIS 5c, the ice sheets retreated again with bare ice remaining on Baffin Island and Ellesmere Island, then re-advancing.Compared to MIS 5d, the ice sheets were much thicker in Baffin Island and the Labrador-Quebec sector with around 1500 m.The Cordilleran Ice Sheet also grew notably during MIS 5b.Retreated happened during MIS 5a, with a remaining ice sheet on the northeast coast of North America and ice caps in Barents-Kara area.
During MIS 4, the ice sheets extended far further south in both North America and Eurasia.For North America, the ice sheets built up in Labrador-Quebec, Keewatin, the Great Lakes and Cordilleran areas separately, leaving the Hudson Bay Lowlands, the western part of Hudson Bay and south of Keewatin almost ice free.For Eurasia, the Barents-Kara Ice Sheet, the Scandinavian Ice Sheet and the British-Irish Ice Sheet built up, while the Scandinavian Ice Sheet and Barents-Kara Ice Sheet separated.
During MIS 3, the total ice volume increased gradually, accompanied with fluctuations due to the Dansgaard-Oeschger events.
During this period, the Scandinavian Ice Sheet and Barents-Kara Ice Sheet merged, western Laurentide Ice Sheet and eastern Laurentide Ice Sheet merged, and the Cordilleran Ice Sheet and the Laurentide Ice Sheet merged.At around 21 kyr BP, the Northern Hemisphere ice sheets reached their maximum extent.
For the pre-LGM ice sheet history, we compared our results with some geological evidence of ice flow directions.Our results show some consistency with the geological constrains for North American ice sheets.According to the review by Kleman et al. (2010), an ice divide close to the Labrador coast in the Quebec sector existed earliest at MIS 5b or 5d during the last glacial-  et al., 2004).In other words, the Eurasian ice sheets advanced further southwest from MIS 4 to the LGM.In our simulation, the Barents-Kara Ice Sheet did not build up prior to MIS 4 and there was no change in ice sheet extent through time.This is likely because the large-scale North American ice sheet build-up changed the atmospheric stationary waves, which in favor of the growth of southwest Eurasian ice sheets (Roe and Lindzen, 2001;Liakka et al., 2012;Löfverström et al., 2014).To simulate this would require feedbacks between ice sheet model and atmosphere model.
The ice sheet configuration during the LGM is relatively well known.According to Margold et al. (2015), there were three major domes of Laurentide Ice Sheet: Labrador, Keewatin and Foxe, which can also be observed in our simulation.The North American ice sheets extended southward up to 40 N with ice sheet thickness up to 3000 m.The interior of Alaska was ice free during the LGM.For Eurasia, the ice sheet covered the Barents-Kara Sea, the Scandinavia and extended southwest to the British-Irish area.
For the deglacial period, most of the Northern Hemisphere ice sheets started to retreat at around 15 kyr BP, while the British-Irish Ice sheet retreated earlier at 16.5 kyr BP.By around 13 kyr BP, the total ice volume decreased to half of its maximum volume (Fig. 5), with ice covered regions persisting on Hudson Bay and the Canadian Shield, the center of Cordilleran region, most of Barents-Kara Sea and part of Scandinavia.By 9 kyr BP, all the ice sheets completely retreated except the Greenland Ice Sheet.The modeled ice sheet retreat process was too fast compared to the geological evidence of the deglaciation (Dyke, 2004;Hughes et al., 2016), which is probably because the "glacial index" method we use is from Greenland, while the signal in other region might include other feedbacks.This discrepancy is likely due to the nonlinear processes between ice sheet, atmosphere and ocean are of vital importance during abrupt climate changes and the termination.

Sensitivity of atmospheric forcing from PMIP3 experiment
Figure 7 shows the SLE time series with the same set-up used in Sec.4.1 but using climate forcing from different PMIP3 models.Most of the models can result in ice sheets that produce a maximum sea level fall between 100 m to 150 m, except CNRM-CM5 and MRI-CGCM3, with a maximum sea level fall of 16 m and 49 m respectively.For GISS-E2-R, the total ice volume is relatively large during cold stages, and smaller during warm stages compared to the other models.(Fig. 7, green line).However, there is no evidence of ice sheet buildup in Siberia during MIS2 (Niessen et al., 2013).For the other models, the general patterns are similar, except for CCSM4 and FGOALS-g2 with extra ice sheets on the East Siberian Sea, Laptev Sea and Chukchi Sea.
In order to investigate why the results show such a diversity of configurations with the only difference being atmospheric forcing, and how sensitive the ice sheet build up process is to the forcing, we compared the resultant surface air temperature and precipitation patterns (Fig. 9, S2, S3, S4) and their anomalies (Fig. 10, S5, S6, S7).We found that the ice sheet extent pattern during the LGM resembles the summer surface air temperature (SAT) pattern, where the SAT is less than around -5 C at that time (Fig. 9).In CNRM-CM5, the areas were in the western coast of North America, Keewatin region, Labrador, southern Greenland and Scandinavia mountains.The same pattern is also observed in other models, with much colder area in Siberia in GISS-E2-R, relatively warm results in MRI-CGCM3 and relatively cold summer SAT along the Arctic coast of the Eurasian continent in CCSM4 and FGOALS-g2.The summer surface air temperature seems to have the dominative influence on the ice sheet extent with a threshold of around -5 C.This is probably because summer temperature controls the surface ablation of ice sheets, which is a dominant surface process for Northern Hemisphere ice sheets in the summer.In winter, accumulation is a more prominent process than ablation because of the relatively low temperatures.More ablation may cause a negative surface mass balance, which then influence the ice sheet build up process.
From our simulation, the multi-domed pattern at the LGM can be observed in almost all the model results (Fig. 8).According to the present day precipitation pattern (Fig. 4c-d), there was more precipitation along the coast of North America and Europe, while in the middle of the continents was relatively dry, especially Keewatin region.So how the ice sheet dome existed in these region?Comparing with the temperature and precipitation forcing pattern, we found that in all the models, there was more precipitation in winter in Keewatin at the LGM than present day (Fig. 10), which resulted in accumulation in that region.The change of precipitation pattern could strongly influence the geometry change of the ice sheet.
Here, we have chosen one ice sheet model and varied the forcing in a similar way as in Dolan et al. (2015) for the Pliocene.
The main focus is on testing the sensitivity of an ice sheet model to atmospheric input fields.Varying the ice sheet parameters but with constant forcing is out the scope of the present paper.

Conclusions
The "glacial index" method is used for simulating the Northern Hemisphere ice sheets through the last glacial-interglacial cycle.The sensitivity of ice sheets to atmospheric forcing is examined by using the output from PMIP3 GCMs (Braconnot et al., 2012;Meinshausen et al., 2011).The most important findings are as follows: 1.The simulated ice sheets and geological evidence of ice sheet extent shows some agreements.During glacial inception, the ice sheets first built up along the coast of Quebec-Labrador sector.The growth of eastern Laurentide Ice Sheet was earlier than the western Laurentide Ice Sheet during the build-up stage (Kleman et al., 2010).For the LGM, the ice sheet extent resembles the geological reconstruction quite well, with the ice sheet extent extending southward to 40 N, maximum ice sheet thickness up to 3000 m, ice free Alaska region and a British-Irish Ice sheet (Dyke, 2004;Hughes et al., 2016).The Northern Hemisphere 2. There are also some disagreements of the results with geological constraints.The most notable one is that the retreat process during deglaciation was too fast.By 9 kyr BP, all of the ice sheets except the Greenland Ice Sheet retreated.According to the geological reconstruction, there were still an Innuitian Ice Sheet and significant ice covered areas in Hudson Bay and the Canadian Shield at this time.The ice sheet evolution pattern for Eurasian ice sheets during the period pre-LGM also differed.
The Eurasian ice sheets advanced southwest from MIS 4 to the LGM according to the geological evidence (Svendsen et al., 2004), while this phenomenon is not observed in our simulation.
3. The climatic output from nine PMIP3 GCMs are used to force the ice sheet model, which resulted in diversity in simulated ice sheet configurations.This shows the great sensitivity of glacial-interglacial Northern Hemisphere ice sheets to the atmospheric forcing.We found that the ice sheet extent most resembles the summer surface air temperature pattern.This is because summer surface air temperature is related to ablation, which contributes to a large portion of surface mass loss and plays a dominant role in the buildup of Northern Hemisphere ice sheets.Meanwhile, precipitation related to ice sheet accumulation is a secondary control factor.The change of precipitation pattern could strongly influence ice sheet geometry.
This study clearly demonstrates the great sensitivity of Northern Hemisphere ice sheets to the summer surface air temperature.The implications are that additional constrains on climate output, including model deficiencies themselves, should be considered carefully for simulating glacial-interglacial Northern Hemisphere ice sheets.Part of the mismatch between the simulated ice sheets and geological evidence is probably because of the lack of nonlinear interactions between ice sheet, atmosphere and ocean, which have different signals regionally, and are essential for simulating the ice sheets evolution.Fully coupled systems are available until now for EMICs (e.g., Ganopolski et al., 2010), while the missing processes in EMICs may also have large effects on regional ice sheet distributions.
Surface ablation is important for Northern Hemisphere ice sheets evolution, so careful treatment of ablation processes is essential for a reliable simulation results.For future studies, an alternative ablation scheme to PDD, surface energy balance, will be used for checking the influence of surface ablation.In addition, coupling between ice sheet model and a GCM is a logical next step for considering the feedbacks between ice sheet and other components of climate systems.(Peltier, 2009;Argus and Peltier, 2010;Peltier and Drummond, 2008;Peltier et al., 2015), MOCA from Lev Tarasov (Tarasov andPeltier, 2003, 2002;Tarasov et al., 2012) and ANU from Kurt Lambeck (Lambeck and Johnston, 1998;Lambeck et al., 2003;Lambeck and Chappell, 2001;Lambeck et al., 2002).Fennoscandia and British-Irish Ice Sheets (Gowan et al., 2016;Dyke, 2004;Hughes et al., 2016).The locations of the three domes of Laurentide Ice Sheet: Labrador-Quebec, Keewatin and Foxe.The areas mentioned in this study include the Hudson Bay (HB), the Great Lakes (GL), Baffin Island (BI), Ellesmere Island (EI), Taimyr Peninsula (TP), Laptev Sea (LS), East Siberian Sea (ESS) and Chukchi Sea (CS).The yellow area is the Interior Plains, the pink area is the Canadian Shield and the purple area is the Scandinavia Mountains (SM).

Model
21 (Sect.2.1).We implemented an atmospheric module into the model with a new climate forcing scheme, based on PISM's extensible atmosphere and ocean coupling feature.The experiments require climate inputs including atmospheric forcing and relative sea level change (Sect.2.2), and basal inputs including topography and geothermal heat flux (Sect.2.3).The spatial domain is on a 600 ⇥ 600 polar stereographic grid with a resolution of 20 km, covering the Northern Hemisphere extending to 21.15 N. The vertical spacing is uneven (ranging from 6.344 m to 43.656 m, maximum height is 5000 m), with 201 levels above the bedrock.The bedrock deformation computation domain has a vertical thickness of 2000 m with 21 levels.
Model COSMOS (ECHAM5-JSBACH-MPIOM) with T31 resolution fromZhang et al. (2013)  at Alfred Wegener Institute Helmholtz Center for Polar and Marine Research (COSMOS-AWI), in which a steady climate state at the LGM is simulated (Fig.3a-b and 4a-b).External forcing and boundary conditions are set according to the PMIP3 protocol (see details in Sect.3).
Clim.Past Discuss., https://doi.org/10.5194/cp-2017-105Manuscript under review for journal Clim.Past Discussion started: 30 August 2017 c Author(s) 2017.CC BY 4.0 License.northern Eurasian continent, northern North America and northern Greenland is in CNRM-CM5.The simulated summer state in MRI-CGCM3 is also warmer than most of the other models.In GISS-E2-R, northern Siberia is much colder than the other models.For precipitation, the model results are relatively close to each other (Fig. 10, Figs.S3, S4, S7).

4
Figure5ashows the comparison between sea level equivalent (SLE) calculated from modeled ice sheet volume (red) and the relative sea level reconstruction fromRohling et al. (2014, blue).Generally, the modeled ice volume resembles a sawtooth curve indicating slower ice sheet build up and faster retreat processes.Comparing with the NGRIP climate forcing time series (Fig.2), the timing of the ice sheet volume decreases are consistent with the timing of the abrupt increase signal.This indicates the abrupt warming in Greenland resulted in negative surface mass balance.Since we started from the present day Greenland Ice Sheet configuration, the total ice volume slightly decreased because of higher temperatures during the last interglacial.The total ice volume reached its relative high value at around 109 kyr BP (Marine Isotope Stage 5d, MIS 5d) with a 30 m equivalent sea level fall, then went through a rapid decrease.After that, it reached another two extreme values at around 91 kyr BP and 86 kyr BP (MIS 5b) with 40-60 m of sea level change.The ice sheets reached its second peak in volume at 60 kyr BP (MIS 4) with up to 90 m SLE.The amplitude of the SLE variability is not as large as the reconstruction time series, especially during the glacial inception.This is probably because we use a PDD scheme for calculating the surface mass balance, which may have underestimated the contribution of short-wave and long-wave radiation or albedo change.Between 60 kyr BP and 25 kyr BP, the ice sheets went through high frequency advance and retreat cycles with millennial-scale sea level variation, which is likely related to the Dansgaard-Oeschger signals in NGRIP climate forcing time series.The total ice sheet volume reached its maximum at around 24 kyr BP with approximately 120 m SLE.This lasted for almost 10 kyr until 15 kyr BP, after which it underwent fast retreat to half of its maximum volume at 13.5 kyr BP.Except for a slight increase of ice volume corresponding with the Younger Dryas at 11.7 kyr BP, the total ice volume continued decreasing until 9 kyr BP, then became stable with 6-7 m SLE remaining.

Clim.
Past Discuss., https://doi.org/10.5194/cp-2017-105Manuscript under review for journal Clim.Past Discussion started: 30 August 2017 c Author(s) 2017.CC BY 4.0 License.interglacial cycle.This may indicate the location of glacial inception for North American ice sheets started around the northeast coast.The ice domes in the Keewatin and Labrador sectors were probably dynamically independent for most of the time before the LGM.The Labrador dome expanded southward earlier than the Keewatin sector at around MIS 4. For Eurasian ice sheets, geological evidence indicates that the Barents-Kara Ice Sheet extended further east to the Taimyr Peninsula prior to the LGM, and Barents-Kara Ice Sheet got smaller while Scandinavian Ice Sheet got bigger during each successive glaciation (Svendsen

Figure 8
Figure8shows the ice thickness at the LGM (21 kyr BP) from different models.The models exhibit large differences both in ice sheet thickness and extent from one simulation to the other.Consistent with the SLE time series, the ice sheets from CNRM-CM5 almost never built up during the whole glacial-interglacial cycle, with only small ice sheets in the western coast of North America, Keewatin region, Labrador, southern Greenland, and Scandinavia mountains.MRI-CGCM3 shows similar results with more ice covering Hudson Bay, Greenland and Barents-Kara sea.Compared with results from COSMOS-AWI, there is less ice in North America in GISS-E2-R, but ice build up in Siberia, which contributed to the sea level during the LGM Clim.Past Discuss., https://doi.org/10.5194/cp-2017-105Manuscript under review for journal Clim.Past Discussion started: 30 August 2017 c Author(s) 2017.CC BY 4.0 License.icesheets contribute about 120 m SLE, with North American ice sheets around 80 m, Eurasian ice sheets 30 m, and Greenland ice sheet 10 m at the LGM.A multi-domed Laurentide Ice Sheet was observed in our simulation consistent with observations(Margold et al., 2015).

Figure 2 .
Figure 2. (a) The NGRIP ice core 18 O record (Andersen et al., 2004) and the corresponding value of the "glacial index".(b) The relative sea level change (Rohling et al., 2014, dark blue) with 1 error bars (light blue).

Figure 3 .
Figure 3. Winter (DJF) and summer (JJA) surface air temperature for Last Glacial Maximum (LGM, a-b) from COSMOS-AWI (Zhang et al., 2013) and Present Day (PD, c-d) from NCEP Reanalysis data(Kalnay et al., 1996), and the difference between LGM and PD surface air temperature (LGM minus PD, e-f).

Figure 5 .
Figure 5. (a) Modelled sea level equivalent (SLE) of the Northern Hemisphere ice sheets (red) using COSMOS-AWI compared with the relative sea level change (blue, Rohling et al., 2014) with 1 error bars.(b) Separated sea level equivalent (SLE) of Greenland Ice Sheet (green), Eurasian ice sheets (black), North American ice sheets (blue) and Northern Hemisphere ice sheets (red) through the last glacialinterglacial cycle.

Figure 6 .
Figure 6.Modelled ice thickness (m) evolution through the last glacial-interglacial cycle at typical climate stages.The simulation is forced by the climatology monthly mean surface air temperature and precipitation from COSMOS-AWI.

Figure 7 .
Figure 7. Modelled sea level equivalent (SLE) of Northern Hemisphere ice sheets change through the last glacial-interglacial cycle using the PMIP3 model output.

Figure 8 .
Figure 8. Modelled ice thickness at the Last Glacial Maximum (LGM, 21 kyr BP) using the PMIP3 model output.

Figure 9 .
Figure 9.The surface air temperature (SAT) at the Last Glacial Maximum (LGM) in summer (JJA) in different models participating in PMIP3.

Figure 10 .
Figure 10.The Precipitation (Precip) difference between Last Glacial Maximum (LGM) and Present Day (PD) in winter (DJF) in different models participating in PMIP3 (LGM minus PD).

Table 1 .
Typical parameter values used for glacial cycle simulations.