Greenland ice sheet contribution to sea level rise during the last interglacial period : a modelling study driven and constrained by ice core data

Abstract. As pointed out by the forth assessment report of the Intergovernmental Panel on Climate Change, IPCC-AR4 (Meehl et al., 2007), the contribution of the two major ice sheets, Antarctica and Greenland, to global sea level rise, is a subject of key importance for the scientific community. By the end of the next century, a 3–5 °C warming is expected in Greenland. Similar temperatures in this region were reached during the last interglacial (LIG) period, 130–115 ka BP, due to a change in orbital configuration rather than to an anthropogenic forcing. Ice core evidence suggests that the Greenland ice sheet (GIS) survived this warm period, but great uncertainties remain about the total Greenland ice reduction during the LIG. Here we perform long-term simulations of the GIS using an improved ice sheet model. Both the methodologies chosen to reconstruct palaeoclimate and to calibrate the model are strongly based on proxy data. We suggest a relatively low contribution to LIG sea level rise from Greenland melting, ranging from 0.7 to 1.5 m of sea level equivalent, contrasting with previous studies. Our results suggest an important contribution of the Antarctic ice sheet to the LIG highstand.


Introduction
Eustatic sea level reconstructions from sediment studies suggest that the last interglacial (LIG) sea level anomaly, relative to present day, stands for the highest in the last 200 ka BP (e.g.Vezina et al., 1999).A recent probabilistic assessment of sea level rise during the LIG based on local sea level indicators suggests a likely value greater than 6.7 m (Kopp et al., 2009).
Evidence about the extent and volume of the Greenland ice sheet (GIS) during the LIG period is relatively limited.This period corresponds to the deepest part of Greenland ice cores, where the signal is perturbed because of layer mixing, past surface melting and/or basal melting.Thus, the interpretation of the Greenland ice cores during this period is often difficult.However, LIG ice has been found at the six deep ice core drilling sites (GRIP, GISP 2, North GRIP, Camp Century, Dye 3 and the latest one, NEEM), suggesting a fairly limited ice reduction.Ice older than the LIG period has been found at DYE 3 in South Greenland, but its interpretation is still debated (Alley et al., 2010).Pollen and sediment studies can also provide valuable information.In particular, it is very likely that the southern GIS retreated further during the LIG than during the Holocene, as suggested by pollen studies (de Vernal and Hillaire-Marcel, 2008) as well as sediment studies (Colville et al., 2011).The southern GIS was however Table 1.Estimates of the GIS melting contribution to global sea level rise during the LIG period in previous ice sheet modelling studies.Two other estimates are also provided.
To date, several studies have been carried out to reconstruct LIG Greenland geometry with ice sheet models (ISMs).The reconstruction of past surface mass balance (SMB) in such studies is one of the major issues, as we have very few constraints on its natural variability on the millennial timescale.The most common approach has been to drive the model with a proxy for surface air temperature (we will refer to this formulation as the "index method").This is the approach followed by Cuffey and Marshall (2000); Huybrechts (2002); Tarasov and Peltier (2003), and Lhomme et al. (2005).Among these, considerable uncertainties about the estimates of GIS melting contribution to global sea level remain, ranging from 2.7 m to 5.5 m.Very few other studies have tackled the problem with a more physically based SMB.To reconstruct the GIS during the LIG period, Otto-Bliesner et al. (2006) use a one-way coupling of a general circulation model (GCM) with an ISM, however, without any feedback of the ISM on the GCM.Robinson et al. (2011) apply a regional energy-balance moisture orographic model, driven by a climate model of intermediate complexity, to force an ISM over the last glacial-interglacial cycle.The estimates of mean sea level rise contribution differ considerably in these two studies, with 2.2 to 3.4 m for Otto-Bliesner et al. (2006) and 0.4 to 4.4 m for Robinson et al. (2011).Table 1 lists the various previous estimates.
In the present study, we use the thermomechanically coupled ice sheet model GRISLI (GRenoble Ice Shelf and Land Ice model) in order to investigate scenarios of GIS reconstructions during the LIG period.This study is the first application of a hybrid model, mixing shallow ice and shallow shelf approximations (SIAs/SSAs), to reconstruct the LIG geometry of the GIS.We expect a hybrid model to reproduce the ice stream patterns better than SIA-only models.We also implemented the index method mentioned above.Compared with previous works, we however improved the formulation with a self-consistent index more representative of Greenland surface air temperatures.When working with proxy data, we do account for past surface elevation changes of the GIS by applying a time-dependent correction.We also improve on the classic index formulation used in previous GIS studies by introducing outputs from GCM simulations during the LIG.Another improvement was to re-implement a tracer transport model (Lhomme et al., 2005) to tie our reconstruction to ice core information.
In Sect.2, we first describe briefly the ISM used.The mass balance model formulation, and in particular the index method, is presented in depth.In this section, we also describe our model calibration method.Section 3 is focused on the LIG reconstruction.We examine the validity of our results with regard to proxy information and sensitivity experiments.

The GRISLI ice sheet model
GRISLI simulates ice sheet geometry and physical properties as a function of time.Given a specific climatic forcing, we compute the large-scale dynamic evolution.This ISM is three-dimensional and thermomechanically coupled.This model belongs to the hybrid type (see Kirchner et al., 2011b, for a recent review of ice sheet models).Thus, it accounts for the two major flow regimes observed in largescale ice sheets: slowly moving ice is computed following Hutter (1983) (shallow ice approximation, SIA), whereas fast ice streams and ice shelves are computed following the MacAyeal (1989) formulation (shallow shelf approximation, SSA).
SIA-only models do not explicitly represent the ice shelves, and in this type of model floating ice is either drastically calved or forced to ground.The presence of an ice shelf component in GRISLI gives a more realistic representation of the ice advance, including the advance onto the continental shelf during glacial times.This kind of model is thus particularly adapted for palaeoreconstruction of ice sheets, as it accounts for great variations of the grounding line position.
In addition, the SSA scheme is expected to provide better estimates of ice surface slopes at the margin of the ice sheet, whereas SIA models generally produce unrealistically high slopes.This error in modelled marginal slopes can have an impact on the extent of the ablation zone.
Several such palaeo-applications with GRISLI have already been done for the Antarctic ice sheet (Ritz et al., 2001;Philippon et al., 2006;Álvarez-Solas et al., 2011a) and for the Northern Hemisphere ice sheets (Peyaud et al., 2007;Álvarez-Solas et al., 2011b).The most recent version of the model used here was previously examined to assess its sensitivity to atmospheric forcing fields (Quiquet et al., 2012).
The technical characteristics of GRISLI are largely described in the previously mentioned studies; therefore, here we only describe the most relevant features.The distinction between the three types of flow is the following: -Ice shelves are based on a flotation criterion.Calving is assumed to happen when ice reaches a critical thickness.
-Ice streams (dragging ice shelves) strictly correspond to the location of bedrock valleys on a specified map.This map was already used in Quiquet et al. (2012) and is shown in the Supplement (Fig. S1).It remains constant over the whole transient simulation and represents locations where ice streams are allowed.Then, ice streams are activated when the temperature at the ice-bed interface reaches the melting point.SSA is then used as a sliding law, as proposed by Bueler and Brown (2009).We assume a linear viscous sediment type, with a basal drag proportional to basal velocity, with a coefficient β.
-Ice velocity is computed with the SIA equations only if the considered grid point is outside the prescribed ice stream map or if the temperature at the ice-bed interface is not warm enough.
The grounding line position is determined with a simple flotation criterion.
In the most recent version of the model, we reimplemented the tracer transport model formulation of Lhomme et al. (2005) (based on Clarke and Marshall, 2002).The particle advection problem is solved with a semi-Lagrangian formulation.Location and timing of deposition is computed on each grid point of the ISM.Depositional conditions, such as surface temperature, surface mass balance and surface elevation, are thus available for each vertical grid point within the ISM.A direct comparison with observed ice core profiles is then possible.Thus, as in Lhomme et al. (2005), we used the transport model to simulate the age-depth relationship at several ice core locations.Ice core records are compared with the simulations to define the optimal parameters of the ice sheet model.
GRISLI is run in this study on a 15-km Cartesian grid, with 21 vertical points within the ice.The bedrock data set is extracted from Amante and Eakins (2009).The initial ice thickness is provided by Bamber et al. (2001).We used the geothermal heat flux of Shapiro and Ritzwoller (2004).We locally modified this flux near ice core locations to take into account the values derived from measured temperature profiles in bore holes: a very high value of 135 mW m −2 at North GRIP (Dahl-Jensen et al., 2003) and a very low value of 20 mW m −2 at Dye 3 (Dahl-Jensen et al., 1998).This local modification is attenuated with the inverse square of the distance, within a fixed radius of 225 km.

The standard index method
We use an index method to reconstruct past surface mass balance (SMB).In this approach, present-day atmospheric temperature is perturbed by a spatially uniform anomaly of temperature deduced from proxy data.As mentioned earlier, this method has already been widely used for GIS palaeoreconstructions (Letréguilly et al., 1991;Ritz et al., 1997;Greve, 1997;Cuffey and Marshall, 2000;Huybrechts, 2002;Tarasov and Peltier, 2003;Lhomme et al., 2005) and within the EISMINT (European Ice Sheet Modelling INiTiative) project (Huybrechts et al., 1996).
The first reason why we chose this formulation is because it relies on a very small number of tunable parameters.Considering the large discrepancies in present-day SMB reconstructions between GCMs (Yoshimori and Abe-Ouchi, 2012), we can expect large uncertainties regarding the LIG SMB, because the climatic constraints are even weaker at this time.The index method can be easily tuned to reproduce present-day GIS topography and data derived from ice cores studies.The second reason is because this method relies strongly on proxy data (mainly δ 18 O), which represent the only constraints we have on past climatic conditions.
Present-day monthly near-surface air temperature is the one prescribed from the EISMINT experiments (Huybrechts et al., 1996) where temperature is a simple function of latitude and elevation.We have selected this temperature parameterisation, because it fits better the observed present surface temperature at the drilling sites, a condition necessary ulated accumulation rates using RCMs (RACMO and MAR) total precipitation rates, and precipitation rate, compared with measurements at GRIP (Andersen et al., 2006), North GRIP IP members, 2004), Camp Century (Robin, 1976), Dye 3 (Andersen et al., 2006), andNEEM andDahl-Jensen, 2007).Accumulation rates are expressed in metres of ice equivalent per year.29 Fig. 1.Simulated accumulation rates using RCM (RACMO and MAR) total precipitation rates, and combined precipitation rate, compared with measurements at GRIP (Andersen et al., 2006), North GRIP (North GRIP members, 2004), Camp Century (Robin, 1976), Dye 3 (Andersen et al., 2006), and NEEM (Buchardt and Dahl-Jensen, 2007).Accumulation rates are expressed in metres of ice equivalent per year.to reproduce the observed temperature profiles.Present-day total precipitation is deduced from outputs of two regional circulation models (RCMs).We used the RACMO2 (Ettema et al., 2009) and the MAR (Fettweis et al., 2011(Fettweis et al., ) models, averaged over the 1958(Fettweis et al., -2007 period. period.We developed a composite map of these two RCMs in order to have a good agreement between simulated present-day accumulation rates and measured accumulation rates at the five GIS deep ice core sites (Dye 3, GRIP, North GRIP, NEEM, Camp Century).This agreement is necessary to compare age-depth relationships.To do so, accumulation rates from MAR and RACMO have been compared with measurements at ice core locations.Where MAR exhibited a wet bias (DYE 3 and Camp Century), RACMO showed a good agreement, and where RACMO was too dry (GRIP, North GRIP and NEEM), MAR was close to the observations.An altitude and latitude weighting between these two precipitation fields yielded an overall better agreement (Fig. 1).
A time-dependent but spatially uniform perturbation, T palaeo s , is superimposed onto the present-day nearsurface air temperature field, T pd s | (x,y) , in order to evaluate the past near-surface air temperature field.Also, we take into account the feedback of surface elevation changes on temperature, T topo s | (x,y,t) , assuming a spatially uniform lapse rate.This lapse rate is assumed to have seasonal variations, being lower in summer than in winter.However, monthly values of this lapse rate remain constant over the transient simulation.
A sinusoidal annual cycle is constructed from the annual and July values given by Fausto et al. (2009).
The palaeo-near-surface-air-temperature field results in The T palaeo s | (t) perturbation is deduced from climate proxy data.The most common proxy used is the isotopic content in δ 18 O at GRIP or North GRIP.In previous ice sheet modelling studies, the δ 18 O proxy was completed by the deuterium record for ages older than the LIG period (e.g.Cuffey and Marshall, 2000;Huybrechts, 2002;Lhomme et al., 2005).Section 2.2.2 explains how we selected the proxy for temperature.
In any case, the δ 18 O proxy is converted into temperature anomaly with the following relation: (2)  2).The index used in this study exhibits a warming greater than +5 • C during the LIG.This value is much lower than that of Barker et al. (2011), which remains over +8 • C during 2 ka.
The α i coefficient, hereafter called the isotopic slope, is relatively unconstrained.The value of this slope is generally inferred while comparing simulated and measured temperature profiles at ice core drilling sites.This calibration is, however, dependent on the ice sheet model, because the simulated 3-D-velocity field impacts the simulated temperature profile.Also, the simulated temperature profile is strongly dependent on the past evolution of ice margins (e.g.Cuffey and Clow, 1997).
Temperature changes affect the precipitation rate, and the following correction is done, with P pd r | (x,y) , the present-day total precipitation rate, and P r | (t,x,y) , the past total precipitation rate: This formulation approximates the saturation pressure of water vapour (Charbit et al., 2002) and has been employed for past accumulation reconstructions from ice core layer counting (Johnsen et al., 1989;Dahl-Jensen et al., 1993).
In our approach, ablation is computed with a positive degree day (PDD) formulation (Reeh, 1991).PDDs are calculated with the monthly means of near-surface temperature.PDD coefficients follow Tarasov and Peltier (2002).We make a distinction between solid and liquid precipitation following Marsiat (1994).We take into account potential refreezing of melt water and rainfall in the firn layer, using a parameterisation adapted from Janssens and Huybrechts (2000).
We are aware that we are using a rather simplified approach and that some important processes might be missing, such as the LIG insolation anomaly compared with present day (Berg et al., 2011) and its role in the ablation rate.How-ever, even the most sophisticated RCMs disagree in simulating the present-day SMB of the ice sheet (e.g.Vernon et al., 2012).We prefer to adopt a simpler approach that is nonetheless strongly connected to proxy information.

A multi-proxy index for near-surface air temperature
The index method relies on temperature perturbations.The last interglacial period should of course be covered by this index, as well as several tens of thousands of years prior to this period in order to take into account the previous glacial stage and its effect on the ice sheet response.
The value of δ 18 O measured in Greenland ice cores (North GRIP members, 2004) is generally used as a proxy for temperature in this region.However, up to now, none of the Greenland ice cores shows a continuous measurement of δ 18 O reaching back to a period in time older than 123 ka BP.Alternatives have to be found to reconstruct past surface air temperature over this period of time.
Previous studies with a similar approach to ours have mainly used a composite index with a combination of GRIP δ 18 O and Vostok deuterium excess (Cuffey and Marshall, 2000;Huybrechts, 2002;Lhomme et al., 2005).However, ice cores show asynchronous climatic signals between the two hemispheres (e.g.Blunier et al., 1998).Using the Vostok core to construct a temperature history for Greenland is a drastic simplification which leads to large errors in timing as well as in amplitude.
In order to improve on this, we searched for proxies directly linked with the Northern Hemisphere temperature, or at least with the global temperature.
Considering its relatively high mixing ratio, methane is considered as a relatively robust indicator of past climate, in particular temperature, at the millennial timescale (Loulergue et al., 2008).We use here the EPICA DOME C methane concentration measurements (Lüthi et al., 2008;Loulergue et al., 2008).We scaled the methane amplitude to the North GRIP δ 18 O for the first 123 ka in order to have a reasonable temperature variability prior to the last interglacial period.To do this, we performed a linear regression to derive the conversion coefficient between the two variables.We found a slope of 0.027 ‰ ppbv −1 with a correlation coefficient of 0.70 (calculated over 721 points).
We acknowledge the fact that a direct calibration between Greenland temperature (via δ 18 O) and methane concentration is still a strong simplification.δ 18 O may be more representative of the local winter temperature (Denton et al., 2005), whereas methane is a proxy for global annual mean temperature.Moreover, the maximum methane concentration (around 128 ka BP) came several thousands of years earlier than the supposed maximum Greenland temperature (Masson-Delmotte et al., 2010).In order to have a better estimate of the LIG timing, we chose a sea surface temperature (SST) reconstruction from the ODP980 North Atlantic Fig. 4. Present day GIS geometry in the observations (Bamber et al., 2001) and simulated after a transient run with free surface evolution and with calibrated parameters.Ice thickness differences between simulation and observation are represented.Deep ice cores are located on the map.Contour spacing is 500 m.The ice sheet margin is represented by the red contour.The simulated present day ice sheet is generally too thick at the margins except in the northeastern region.The bedrock is far below sea level in this area, and the current ISM does not maintain this fjord englaciated.After calibration, the resulting simulated ice sheet presents a +0.8 m of sea level equivalent compared with observations.32 Fig. 4. Present-day GIS geometry in the observations (Bamber et al., 2001) and simulated after a transient run with free surface evolution and with calibrated parameters.Ice thickness differences between simulation and observation are represented.Deep ice cores are located on the map.Contour spacing is 500 m.The ice sheet margin is represented by the red contour.The simulated present-day ice sheet is generally too thick at the margins except in the northeastern region.The bedrock is far below sea level in this area, and the current ISM does not maintain this fjord glaciated.After calibration, the resulting simulated ice sheet presents a +0.8 m of sea level equivalent compared with observations.marine core (McManus et al., 1999;Oppo et al., 2006).
Here again, we scaled the SST amplitude to the δ 18 O record when the two signals overlapped.In this case, we found a slope of 1.337 ‰ • C −1 with a correlation coefficient of 0.86 (calculated over 61 points).
We took special care to have all our proxies on the same timescale.For that, we used the one from Lemieux-Dudon et al. ( 2010), as the ODP980 record was already placed on this common timescale by Masson-Delmotte et al. (2010).A comparison between these proxies around the LIG period is shown in Fig. 2.An additional improvement over previous work is that we corrected the measured North GRIP δ 18 O to take into account past surface elevation change.Isotopic values measured along the ice cores are indeed modulated by surface elevation changes, but the index must be representative of a temperature anomaly for a constant surface elevation (as written in Eq. 1).Both scaling of methane and SST to the North GRIP δ 18 O were done only after this surface elevation correction.
We blended the three proxies as follows: -we use the North GRIP δ 18 O for ages younger than 122.3 ka BP; -between 122.3 ka BP and 128.6 ka BP, SST-derived δ 18 O is used; and -for ages older than 128.6 ka BP, we used CH4-derived δ 18 O.
To avoid artificial abrupt changes in climate, the records were blended where they were close enough (at 122.3 and 128.6 ka BP).The transitions between the different proxies are shown in Fig. 2. The multi-proxy index obtained, expressed in terms of δ 18 O, is finally converted into temperature assuming a simple linear relationship (Eq.2).The resulting multi-proxy index is shown in Fig. 3.The temperature perturbation used in this study is available in the Supplement (Table S1).Barker et al. (2011) proposed a Greenland temperature anomaly over the last 800 thousand years based on a thermal bipolar seesaw model.However, a methane-driven sharp signal is observed in this work around 129 ka BP (Fig. 2).This signal seems to be lacking in other Greenland proxies (Masson-Delmotte et al., 2010).Thus, we did not use this reconstruction of temperature anomaly, but a comparison with our multi-proxy index is shown in Fig. 3.

Introduction of GCM simulations in the index method
The direct perturbation of present-day atmospheric fields with the previously described index does not account for changes in atmospheric circulation and its consequences on temperature and precipitation.To improve on this, we also used GCM snapshots run at 126 ka BP with the present-day topography of the GIS but under 126 ka BP orbital forcings.Pre-industrial greenhouse gases are prescribed.We selected two GCMs among the CMIP-3 models for their availability over the LIG period: -the IPSL-CM4 (Marti et al., 2010;Braconnot et al., 2008) and Salas-Mélia et al., 2005).
Several studies have shown that an anomaly method ("best" present-day climatology + perturbation deduced from GCMs) is preferable to a direct forcing with absolute fields from coarse-resolution GCMs (e.g.Quiquet et al., 2012).To account for this, we use a similar methodology to the one of Charbit et al. (2002) or Kirchner , are superimposed on top of present-day forcing fields: . (4) The obtained snapshot at 126 ka BP is then, here again, perturbed with an index.The previously prescribed index (Fig. 3) was designed to present a zero-anomaly at 0 ka BP ( T palaeo s | (t=0) in Eq. 1).We cannot use our temperature perturbation directly on top of these new snapshots as we would be double-counting the LIG warming (GCM 126 ka BP anomaly and T palaeo s | (t=126k) ).Thus, we apply a homogenous lowering of the previously described index by its own value at 126 ka BP (5 • C), in order to achieve a zero-anomaly at this time.
In the following, we will refer to the "no-anomaly experiment" for the direct application of the index on present-day forcing fields and respectively to CNRM and IPSL "anomaly experiment" for the two others.

Ice model calibration
To calibrate the ISM, we run transient experiments for the last 60 ka BP.We selected the parameters of the model to be as close as possible to the present-day state of the GIS.For that, we used the following constraints: -simulated geometry (ice covered area and ice volume) as close as possible to the one of Bamber et al. ( 2001) and -simulated surface velocity field close to the present-day observations of Joughin et al. (2010).
In addition to these large-scale features, we also included constraints from ice core drilling: -simulated temperature profile similar to the borehole measurements and -simulated age-depth relationship close to the GICC05 timescale (Rasmussen et al., 2006).
To score our simulations, we used some numerical estimates (present-day simulated volume and iced area, presentday basal temperature, Younger Dryas and Laschamp event depths) and also some qualitative estimates (high surface velocity area and temperature profile shape).The present-day simulated geometry of the GIS after calibration is compared with observations in Fig. 4.
As mentioned earlier, we owe the possibility of using an age-depth relationship as a constraint to the newly reimplemented tracer transport model.The simulated agedepth relationship at ice core locations is compared with the GICC05 timescale in Fig. 5.
The simulated geometry of the GIS is mainly driven by SMB and dynamical parameters of the ISM.In particular ice extent is governed by the choice of ablation coefficients of the PDD method and is modulated by the choice ulated age-depth relationship after calibration (plain line) and GICC05 (dashed lines) for four ores.The common GICC05 timescale is mainly available for GRIP and North GRIP for the d, and for Dye 3 for the early Holocene.The stars represent the approximate depths of the ryas (approx.11.5 ka BP) and the Laschamp event (approx.40.8 ka BP).We mainly improve rth GRIP dating with our calibration, as this record is used both as a forcing (age-δ 18 O) and as nt (age-depth).33 Fig. 5. Simulated age-depth relationship after calibration (plain line) and GICC05 (dashed lines) for four deep ice cores.The common GICC05 timescale is mainly available for GRIP and North GRIP for the long record, and for Dye 3 for the early Holocene.The stars represent the approximate depths of the Younger Dryas (approx.11.5 ka BP) and the Laschamp event (approx.40.8 ka BP).We mainly improve on the North GRIP dating with our calibration, as this record is used both as a forcing (age-δ 18 O) and as a constraint (age-depth).
of dynamical parameters (Ritz et al., 1997;Greve, 1997).On the other hand, the temperature profile and age-depth relationship mainly depend on the palaeoclimate perturbation parameters, such as the glacial-interglacial amplitude of temperature change, α i of Eq. ( 2), and precipitation ratio, γ in Eq. (3).
Our choice of parameters is summarised in Table 2.Note that we may differ from earlier work with a similar methodology (e.g.Clarke and Marshall, 2002;Tarasov and Peltier, 2003;Lhomme et al., 2005) because of the use of a different model (SIA/SSA here rather than SIA-only), the use of a different resolution and the temperature perturbation used (methane-SST rather than Antarctic deuterium).Furthermore, we have not calibrated only on the age-depth relationship but also on the temperature profile, and the simulated ice volume and extent.

Simulated topographies and sea level contribution
We performed transient simulations of the last 200 ka with the calibrated version of the model.Figure 6 presents the simulated ice volume of the GIS over the last glacialinterglacial cycle.The minimum ice volume is at around 121 ka BP in our reconstructions.The GIS geometry at this time is represented in Fig. 7 for the three experiments described previously (no-anomaly, CNRM and IPSL anomalies).The pattern of retreat is relatively similar in the three experiments: the northeastern and southwestern margins retreat more than the other regions.Ice is preserved at ice core locations in all three experiments.The South Dome is relatively stable, which is compatible with some geological evidence (Colville et al., 2011).Changes in elevation are however drastic in the no-anomaly experiment for the sites of Dye 3, Camp Century and NEEM.
It clearly appears that the no-anomaly experiment presents the most retreating GIS.It can be explained by the difference in surface temperature imposed during the LIG.At North GRIP, the two GCMs present a 126 ka BP annual mean nearsurface temperature very similar to the pre-industrial value.In contrast, the index presented in Fig. 3 exhibits a +5 • C increase at 126 ka BP in annual mean near-surface temperature anomaly.However, the GCMs exhibit a stronger LIG seasonal cycle compared with pre-industrial.Thus, the difference between the GCM temperatures and the index is however smaller when considering the summer temperature.The mean July temperature anomaly at North GRIP in the IPSL model is +5 • C, which is equivalent to the index.CNRM is, however, somewhat colder at North GRIP with only a +3 • C anomaly in July temperature.
We raise here one of the major questions concerning the index method: is the isotopic composition representative of the annual mean or of a season-specific temperature?Considering one or the other option will locally introduce some biases, depending on the geographical position and also the time period considered.
Our estimates of the contribution of GIS melting to mean sea level rise during the LIG period are as follows: -no-anomaly experiment of 2.90-3.71m; -CNRM anomaly experiment of 0.65-1.46m; and -IPSL anomaly experiment of 0.68-1.50m.
The lower bound of these estimates corresponds to the difference between the LIG minimum volume simulated by the ISM and the observed present-day volume.The upper bound represents the difference between the same LIG volume and the simulated present-day volume.Thus, the upper bound takes into account the ISM bias in simulating ice retreat because of an overestimation of the simulated present-day volume by 10-15 %.We wish to stress that the calibration was not done only on the volume, but also on the extent, the velocity field and on deep ice core drilling information in terms of temperature profile and age-depth relationship.We think that the resulting set of parameters is more robust than the one we would have obtained with a volume calibration alone.
The ranges of simulated volume for the two experiments with anomalies are very similar.Taking into account the changes in atmospheric circulation due to a change in orbital forcings seems to decrease the sensitivity of the GIS to the LIG warming.However, further studies should repeat these experiments with more GCMs in order to increase our confidence in this result.The different representations of physical and dynamic processes among GCMs may lead to different sensitivities to solar parameter-induced warming.The difference between the no-anomaly experiment and the two others may also reflect the inability of GCMs to reproduce the highly variable palaeoclimates (e.g.Masson-Delmotte et al., 2006).In addition, these GCM simulations do not take GHG variations into account, but simply assume fixed preindustrial conditions.They were also run under present-day GIS topography, whereas changes in the GIS topography during the LIG may affect the simulated climates.This assumption might introduce a cold bias.However, we carried out a drastic sensitivity study in which the CNRM model was run under 126 ka BP orbital forcing conditions, and where the whole Greenland ice sheet was removed.Results (not shown) indicate that the change in orbital forcings is the main driver even if the temperature signal caused by the modified ice sheet cannot be discarded.Further investigations should be performed in this direction.Finally, GCMs often suffer also from simplified snow representations.For example, the simple ageing of the snow pack in these models assumes no albedo changes for Eemian vs. pre-industrial, although these could have a significant impact leading to an increase in the LIG temperature anomaly (e.g.Punge et al., 2012).
The robustness of our results to the choice of parameters was briefly tested in sensitivity experiments.In particular, the accumulation ratio, γ , and the isotopic slope, α i , have an important effect on past climate reconstruction assumptions that are likely to affect the simulated ice sheet.Figure 8 shows that the two experiments using GCM anomalies are relatively robust.On the contrary, the precipitation ratio seems to greatly affect the results in the no-anomaly experiment.In this experiment, the value of 0.05 • C −1 for γ leads to an almost complete melting of the GIS.

Palaeo-data confrontation: is the warming realistic?
Our reconstructions depend mainly on our climatic assumptions.In particular, the chosen index governs the near-surface air temperature amplitude change over the GIS.Some reconstructions of the possible range of temperature change based on proxy estimates are however available.The estimated range of summer near-surface Fig. 9. Sensitivity of simulated GIS melting contribution to sea level rise during the LIG period to: isotopic slope which governs prescribed temperature anomaly amplitude, ranging from 0.30 ‰ • C −1 (lower than calibrated value and corresponding to higher values of temperature anomalies) to 0.42 ‰ • C −1 (higher than calibrated value and corresponding to lower values of temperature anomalies); precipitation ratio γ ranging from 0.05 • C to 0.09 • C, values lower than the calibrated ones, corresponding more to what has been used in previous studies; the use of Reeh (1991)  temperature in Greenland based on various proxies is 4-5 • C (CAPE Members, 2006).
Figure 9 shows the evolution of the July near-surface temperature at ice core locations.Note that, in our PDD formulation, melt is possible with surface temperatures greater than −5 • C. With this condition, the sites of Dye 3, Camp Century and even NEEM are affected by summer melting in all of our three experiments.If our estimates are correct, the NEEM ice core should show evidence of melt during the LIG.It is also worth noting that, even with melting events at NEEM, we still maintain a relatively extended Greenland ice sheet.
The July near-surface air temperature presented in Fig. 9 accounts for surface elevation changes.The no-anomaly experiment produces a change in surface elevation greater than one thousand metres at Camp Century, and results in a drastic July temperature increase (more than +12 • C compared with the present-day value).There is indeed no evidence in ice cores of such a great surface elevation change at Camp Century.We therefore consider our no-anomaly experiment to be probably not realistic.
For the two GCM anomaly experiments, the LIG warming compared with present day was higher at North GRIP than at Camp Century due to surface elevation changes.LIG surface elevation at North GRIP was slightly lower, whereas it was higher at Camp Century.At North GRIP, we have a +4 • C (CNRM) and +6 • C (IPSL) increase, whereas the increase is only +2 • C at Camp Century for both experiments.This range is however compatible with proxy estimates (CAPE Members, 2006).

Conclusions
This study provides estimations, strongly based on proxy information, of the possible GIS melting contribution to sea level rise during the LIG, using a hybrid ice sheet model.South Dome is relatively stable even in a warmer climate, while the northernmost ice core locations seem to be the ones at risk.We ran three different experiments.The no-anomaly experiment results in great variations of surface elevation in particular at the margin.Given that ice cores show no evidence of such drastic changes, we consider this reconstruction to be unrealistic.Keeping only the two experiments with GCM anomalies, we suggest a GIS melting contribution to sea level rise ranging from 0.7 to 1.5 m with a preferred estimate of approximately 1 m.This low value indicates a large contribution from the Antarctic ice sheet.
Even if we follow a methodology largely employed in the past (Cuffey and Marshall, 2000;Huybrechts, 2002;Tarasov and Peltier, 2003;Lhomme et al., 2005) to quantify the LIG Greenland melting, our conservative estimates disagree with these previous works.One reason could be the index used in this study.We consider our combination of North GRIP δ 18 O, North Atlantic SSTs, and Antarctic methane more reliable than the combined δ 18 O and Antarctic deuterium record used in the previous studies.Our temperature anomaly during the LIG period exhibits a relatively smooth pattern, which differs from the sharp signal obtained with the Antarctic deuterium record.But still, we acknowledge the limitations of using a proxy such as the δ 18 O to reconstruct large-scale Greenland temperature perturbations.Our reconstructions do not disagree with some other estimates, but are generally in the lower bound (Robinson et al., 2011;Colville et al., 2011;Stone et al., 2012).The pattern of the northern retreat of the GIS during the LIG is also similar to the one of Born and Nisancioglu (2012).
The reconstructions presented here are however tainted by the simplifications required to achieve them.We have shown that reconstructions are highly sensitive to changes in atmospheric circulation.One way to improve confidence in our reconstructions would be to include more temporal snapshots.Another way would be to increase the number of GCMs used.There are still only few GCM simulations of the LIG, but they are becoming available.One of the major simplifications of the index method is the use of a simple lapse rate to take into account the impact of surface elevation changes on temperature.To improve on that, further studies should include GCM simulations with various topographies of the GIS during the LIG.
Ideally, a direct coupling of a GCM with an ISM (at least a one-way coupling) would give a more realistic representation, but would also be strongly dependent on the GCM used, as mentioned in Quiquet et al. (2012).Another disadvantage of a direct coupling is the apparent inability of GCMs to reproduce the high variability observed in proxy data.Incorporating proxy constraints into GCMs is required (Masson-Delmotte et al., 2006) to obtain more accurate results.

Fig. 2 .Fig. 2 .
Fig. 2. Comparison between the 3 proxies used to generate the composite index around the LIG pe Blue: δ 18 O at North GRIP.The plain line represents the direct measurement and the dashed line re sents the value after the surface elevation changes correction.Red: SST reconstruction using ODP marine core, after scaling on the δ 18 O record, corrected from elevation changes.Green: Methane centration measured along the EPICA-DOME C ice core in Antarctica, after scaling on the δ 18 O rec corrected from elevation changes.Purple: Barker et al. (2011) synthetic Greenland δ 18 O based thermal bipolar seesaw model.Vertical bars stand for the transition from one proxy to another in construction of the multi-proxy index.SST reconstruction presents an important minimum at 130 ka When generating the index, transitions from methane to SST were done at 128.6 ka BP and 134 k to evaluate the impact of this collapse in the SST record.Both indeces produced the same LIG retre

Fig. 3 .Fig. 3 .
Fig. 3. Black: Multi-proxy index used in this study, expressed as an annual mean temperature anomaly relative to present day.Purple: Barker et al. (2011) synthetic Greenland δ 18 O, converted into an annual mean temperature anomaly with Eq. 2. The index used in this study exhibits a warming greater than +5 • C during the LIG.This value is much lower than Barker et al. (2011), which remains over +8 • C during 2 ka.

Fig. 6 .Fig. 6 .
Fig. 6.Simulated Greenland ice volume, expressed in m of sea level equivalent.For the 0-100 ka BP the atmospheric forcings are the same in the three experiments, resulting in a very similar volume evolution.

Fig. 7 .Fig. 7 .
Fig. 7. Simulated GIS at 121 ka BP for: (A) No-anomaly experiment; (B) CNRM anomaly experiment; and (C) IPSL anomaly experiment.Deep ice cores are located on the map.Contour spacing is 500 m.The ice sheet margin is represented by the red contour.

Fig. 8 .
Fig.9.Sensitivity of simulated GIS melting contribution to sea level rise during the LIG period to: isotopic slope which governs prescribed temperature anomaly amplitude, ranging from 0.30 ‰ • C −1 (lower than calibrated value and corresponding to higher values of temperature anomalies) to 0.42 ‰ • C −1 (higher than calibrated value and corresponding to lower values of temperature anomalies); precipitation ratio γ ranging from 0.05 • C to 0.09 • C, values lower than the calibrated ones, corresponding more to what has been used in previous studies; the use ofReeh (1991) ablation coefficients rather than Tarasov and Peltier (2002) ones.(A): no-anomaly experiment; (B): CNRM anomaly; (C): IPSL anomaly.The scales of volume variation are changed from (A) to (B) and (C).The two experiments with GCM snapshots seem to be relatively robust to the choice of parameters.The experiment presenting the greatest retreat is the most sensitive one.Note that we cannot produce a realistic present day ice sheet with the original Reeh (1991) ablation coefficients, as already mentioned inRitz et al. (1997) orGreve et al. (2011).

Fig. 8 .Fig. 9 .
Fig. 8. July temperature evolution during the LIG at the deep ice core sites.Plain line: no-anomaly experiment; Dashed line: CNRM anomaly; dotted line: IPSL anomaly.Anomaly experiments always suggest lower values for July temperature.Melt may potentially occur at Dye 3, Camp Century and NEEM sites, even in the anomaly method.

Table 2 .
Model parameters used in GRISLI for this study.