Interactive comment on “ Coupled regional climate – ice sheet simulation shows limited Greenland ice loss during the Eemian ”

Helsen et al. present in this paper the first modeling exercise of the last interglaciation (LIG) Greenland Ice Sheet (GIS) that is forced by GCM climate downscaled through a regional climate model with more precise modeling of the surface mass balance. They find that the GIS contributed ∼2.1 m to LIG sea level rise (range 1.2-3.5 m). Given the >6 m highstand, an Antarctic contribution is also required.


Introduction
Eemian summer temperatures over most Arctic land reached peak values of 4-5 K higher than during the pre-industrial era (CAPE-Last Interglacial Project Members, 2006;Overpeck et al., 2006).Regionally, Eemian temperature anomalies were lower or even negative, in particular sea surface temperature (SST) over the Labrador Sea and Nordic Seas (Bauch et al., 1999(Bauch et al., , 2012;;Winsor et al., 2012).Reconstructions of peak eustatic global sea level rise (SLR) vary from 4 to 9 m above the present-day level (Jansen et al., 2007;Kopp et al., 2009;Dutton and Lambeck, 2012).However, the contribution of the Greenland Ice Sheet (GrIS) to the Eemian sea level highstand is uncertain.Interpretation of ice core records is controversial (Alley et al., 2010); some studies suggest an almost complete deglaciation of southern Greenland and a substantially reduced central dome (Koerner, 1989;Koerner and Fisher, 2002); later studies claim the survival of the southern dome, based on the presence of Eemian ice at Dye-3 (NGRIP members, 2004).Moreover, the inclusion of biological material in the basal silty ice of the Dye-3 ice core dating back to > 400 kyr BP suggests that the ice sheet has not retreated from this location since then (Willerslev et al., 2007).In addition, only minor ice thickness changes in central/northern Greenland are deduced from ice cores in this area (NGRIP members, 2004), a conclusion that is supported by the recent interpretation of the NEEM ice core, pointing to a minimum elevation of only 130 ± 300 m below the present-day surface (NEEM community members, 2013).
Ice sheet modelling experiments result in a broad spectrum of Eemian ice loss estimates, ranging from a sea level equivalent of 0.4 to 5.5 m, with a tendency towards higher mass loss estimates (Cuffey and Marshall, 2000;Huybrechts, 2002;Tarasov and Peltier, 2003;Lhomme et al., 2005;Otto-Bliesner et al., 2006;Robinson et al., 2011;Fyke et al., Fig. 1.Schematic outline of the asynchronous model coupling.An Eemian run of the ECHO-G global climate model is used as lateral boundaries for the regional atmospheric climate model, which in turn computes SMB, T s and refreezing that are used as forcing for the ice sheet model.Every 1.5 kyr, the changed ice sheet topography and extent are updated in the regional atmospheric climate model, after which new fields of SMB, T s and refreezing are calculated, etc. 2011; Born and Nisancioglu, 2012;Quiquet et al., 2013;Stone et al., 2013).The outcome of modelling experiments depends strongly on the method chosen to calculate the surface mass balance (SMB), being the sum of snowfall, sublimation and meltwater run-off (R).Often, a single temperature chronology is used derived from ice core records.This temperature is then used to estimate changes in accumulation and R, both based on lapse-rate considerations (Cuffey and Marshall, 2000;Tarasov and Peltier, 2003;Lhomme et al., 2005;Quiquet et al., 2013).This approach has the advantage of being computationally cheap, but strongly simplifies the complex nature of SMB changes.Alternatively, forcing fields from climate models can be used, such as what has been done e.g. for the Eemian with time-slice (one-way) coupled climate-ice-sheet model experiments, using a positive degree day (PDD) approach to calculate melt (Otto-Bliesner et al., 2006;Born and Nisancioglu, 2012;Stone et al., 2013).
A physically more realistic approach to estimate the SMB is to use a regional climate model coupled to a multi-layer snow and ice melt model with prognostic albedo calculations.These models simulate modern GrIS SMB very well (Fettweis et al., 2008;Ettema et al., 2009), and can be considered a major improvement over coarse-gridded GCMs.This approach has also been successfully applied for past climate states (Van de Berg et al., 2011), using a fixed ice sheet topography.Here we two-way-couple the regional atmospheric climate model RACMO2 ( Van de Berg et al., 2011) to the thermomechanical ice sheet model ANICE (Helsen et al., 2012), using the extensively validated GCM ECHO-G as lateral forcing (Kaspar et al., 2005) (Sect.2), to conduct a 15 kyr transient ice-sheet-climate simulation of the GrIS over the Eemian interglacial period.This approach offers two major improvements with respect to earlier Eemian GrIS simulations.Firstly, we take advantage of the improved performance of regional climate models over GCMs, and secondly, the two-way coupling ensures a physically realistic representation of the feedbacks between the continuously changing ice sheet topography and the climate forcing, yielding more confidence in the SMB reconstruction, and hence in the ice volume calculation over the Eemian.

Coupled model set-up
The computational costs of high-resolution (∼ 20 km) atmospheric climate models limit maximum simulation periods to decades.Since we simulate the GrIS evolution over 15 kyr, one-to-one coupling of an ice sheet model and a climate model is not feasible at present for this type of study.As a viable alternative, we use an asynchronous two-waycoupling strategy, involving the 3-D thermomechanical ice sheet model ANICE (e.g.Van de Wal, 1999;Bintanja and Van de Wal, 2008;Helsen et al., 2012) and the regional atmospheric climate model RACMO2 (e.g.Ettema et al., 2009;Van de Berg et al., 2011).Every 1.5 kyr, fields of surface topography and ice mask are fed from the ice sheet model to the regional climate model.We compute a new climatology using appropriate boundary conditions (see below) and subsequently transfer modelled SMB, surface temperature and amount of refrozen meltwater to the ice sheet model.The latter is then run for another 1.5 kyr, and so forth (Fig. 1).Changes are only applied to the ice cover in Greenland; the ice in the Canadian Archipelago is left unchanged in the regional climate model.This coupling was performed 10 times, allowing for an ice-sheet-climate reconstruction from 130 to 115 kyr BP.
During ice sheet model simulations (in between model couplings), we continuously correct for the SMB-elevation feedback using the SMB gradient method (Helsen et al., 2012).We recalculate the coefficients of SMB gradients for each Eemian climatology using all RACMO2 grid points within a minimum search radius of 150 km.This radius is enlarged until a minimum number of H s -SMB pairs (n = 37) is found, following Helsen et al. (2012).
Our modelling strategy involves multiple transfers of data products from the regional climate model to the ice sheet model and vice versa.As RACMO2 is run on a geographically based coordinate system, and ANICE on a (different) rectangular x-y grid, back and forth mapping and interpolating is inevitable; this is carried out using OBLIMAP, a coupler involving oblique stereographic mapping and a centre of projection at 40 • W, 72 • N (Reerink et al., 2010).

The Regional Atmospheric Climate Model
The Regional Atmospheric Climate Model version 2.1 (RACMO2) of the Royal Netherlands Meteorological Institute (KNMI), adjusted at UU/IMAU for use over ice sheets, is used to calculate all Eemian climatologies used in this study (Ettema et al., 2009;Van de Berg et al., 2011).RACMO2 is run on an 18 km × 18 km grid, with a domain size of 206 × 264 grid points (Fig. 1).It is forced at its boundaries by 6-hourly prognostic fields of temperature, humidity, wind, surface pressure, sea surface temperature and sea ice fraction from the ECHO-G global climate model (Kaspar et al., 2005(Kaspar et al., , 2007)).A detailed description of this set-up is given in Van de Berg et al. (2011).The performance of ECHO-G for the Eemian climate has been extensively validated with proxy records (Kaspar et al., 2005).The interior of the RAMCO2/GR domain is allowed to evolve freely.Each Eemian time slice climatology is based on the last 25 yr of a 30 yr simulation in order to avoid spin-up effects.Retreating ice on land is replaced by tundra, with an albedo of 0.18.Solar incoming radiation is adjusted for each Eemian time slice simulation to account for the influence of orbital forcing (Table 1; Fig. 2, Berger and Loutre, 1991).For each Eemian time slice experiment, the ECHO-G forcing fields are adjusted to account for a changing ambient climate during the simulated period.The changing boundary conditions are based on two ECHO-G climatologies: Eemian optimum (EO) and glacial inception (GI) (Kaspar et al., 2007).We used an interpolated climatology as forcing for RACMO2.Differences in monthly means between GI and EO are calculated for each forcing field (X).Subsequently, a fraction (f ) of these differences is added to the EO fields for each time slice experiment in order to account for the changing boundary conditions through the Eemian: (1) The value of f is set proportional to JJA (June-July-August) insolation at 60 • N, with f = 1 at 115 kyr BP and f = 0 at 125 kyr BP (Fig. 2).This approach allows for a gradual transition of the climate forcing from the Eemian optimum to glacial inception.We used f = 0 for time slices between 130 and 125 kyr BP (hence maximum ambient air temperature), since we only wish to use the scaling procedure for interpolation between the two climatologies, not for extrapolation.This is justified by palaeoclimatic evidence for already high global temperatures (Winograd et al., 1997;CAPE-Last Interglacial Project Members, 2006;EPICA community members, 2006) and greenhouse gases (Petit et al., 1999;Siegenthaler et al., 2005) during this time.Note that no other than modern ice cover is present outside Greenland in our climate model simulations.In particular (remnants of) the Laurentide Ice Sheet have probably had influence on atmosphere and ocean circulation in the early Eemian, but this is currently difficult to quantify.

The ice sheet model
We use the 3-D thermomechanical model ANICE ( Van de Wal, 1999;Bintanja and Van de Wal, 2008;Graversen et al., 2011;Helsen et al., 2012) based on the shallow ice approximation (SIA, Hutter, 1983), and explicitly accounting for the temperature-dependent stiffness of the ice.Ice temperature is calculated based on the 3-D advection, diffusion, friction, geothermal heat flux (G) at the bottom and annual surface temperature (T s ) adjusted for the effect of refreezing (Helsen et al., 2012).The vertical dimension is scaled with the local ice thickness, and consists of 15 layers with increasing resolution near the bed in order to better resolve the local large gradient in ice velocity.For areas where basal temperatures reach the pressure melting point, we allow the ice to slide over its bed by using a Weertman-type sliding law (Weertman, 1964) corrected for the effect of subglacial water pressure (Bindschadler, 1983).Formation of ice shelves is not allowed; as soon as the ice starts to float, it breaks off, Table 1.Climate model forcing parameters.Orbital parameters are calculated for the midpoint of each time slice (Berger and Loutre, 1991).Greenhouse gas concentrations between 130 and 125 kyr BP follow the ECHO-G settings (Kaspar et al., 2007).The concentrations between 125 and 115 kyr BP follow the observed concentrations in the Vostok ice core (Petit et al., 1999;Sowers, 2001).ANICE describes sheet flow (SIA), which is the dominant flow type over the major part of the GrIS.Although it is known that fast-flowing outlet glacier dynamics are not well captured by this approximation, we argue that the errors introduced by the omission of this type of flow are minor, since the ice sheet quickly retreated onto land during our Eemian simulation, strongly reducing the influence of marine-terminating outlet glacier dynamics.
Although Ellesmere Island is partly within the domain of the ice sheet model, and hence ice sheet changes in this region are also simulated, we do not include these in the computation of our ice volume changes.

Initialisation
The size and shape of the GrIS at the start of the Eemian are important but poorly constrained initial conditions for this simulation.The response of the ice sheet to Eemian warming will partly depend on the temperature history before the Eemian, due to the temperature dependence of the ice viscosity (Fig. 3).For a realistic simulation of the GrIS response to Eemian warming, the influence of the preceding glacial cycle on size, shape and ice temperatures should therefore be taken into account.To this end, we performed an experiment that simulates GrIS evolution through the penultimate glacial cycle, from ∼ 237 to ∼ 130 kyr BP.
Ice core records from Greenland do not cover the penultimate glacial cycle, so we need to rely on an ice-corederived climate record from Antarctica. Figure 4a shows that temperature histories for Greenland (GRIP, Johnsen et al., 2001) and Antarctica (Vostok, Petit et al., 1999) differ considerably.We show the Vostok record for the period before 105 kyr BP, since the GRIP record is disturbed due to ice-flow irregularities prior to this date (NGRIP members, 2004).Most notable differences are a much smaller temperature amplitude and less high-frequency variability in the Vostok record.Therefore, both records are converted to a glacial index GI t (Fig. 4b), such that GI = 1 represents Last Glacial Maximum (LGM) conditions and GI = 0 corresponds to present-day conditions, following Greve (2005).The black line represents the GI t signal after applying a Gaussian filter with a 2 kyr filter width.The red line shows the original variability in the GRIP record, which is much higher.The inset in Fig. 4b shows a scatter plot of the standard deviation in GI over a 3 kyr moving average as a function of the mean GI, and illustrates that glacial conditions are characterised by larger variability than during Holocene interglacial conditions.It is likely that the penultimate glacial cycle in Greenland was also characterised by such high-frequency variability.Using the linear regression in this scatter plot (blue line, σ GI = 0.05 + 0.24 GI) and a random-number generator we added high-frequency variability to the 3 kyr mean GI curve (Fig. 4c), which can then be regarded as a potential temperature history over Greenland for the period prior to the Eemian.This temperature history was used as forcing for the ice sheet model prior to 130 kyr BP.We use the SMB gradient method (Helsen et al., 2012) to calculate a direct adjustment of the SMB as a function of temperature perturbations.The present-day SMB field (Ettema et al., 2009) is used as input for this method, and spatial SMB gradients within this field are calculated.Following Helsen et al. (2012), the temperature forcing as derived above is translated into a climatic elevation change ( H climate ) using an atmospheric lapse rate γ atm = −7.4K km −1 .This value is obtained from a regression of elevation and RACMO2 surface temperature.Due to the spatially variable SMB gradients, a single temperature perturbation is translated into a spatially distributed SMB adjustment.
The resulting simulated ice volume over the last two glacial cycles is shown in Fig. 4d.Maximum ice volume (3.70 × 10 15 m 3 ) occurs at 136 kyr BP, following growth during the penultimate ice age.Ice sheet retreat occurs during the glacial termination, as well as interior thickening as a consequence of increased accumulation.For example, surface elevation at NEEM increased by 140 m between 135 and 130 kyr BP (not shown).We start our coupled regional climate-ice-sheet model experiments at 130 kyr BP, with an ice sheet volume of 3.60 × 10 15 m 3 (indicated by the red cross in Fig. 4d), 24 % larger than today.This is also the state of the ice sheet at the start of the Supplement Movie.

Eemian GrIS changes
The immediate response to the warm Eemian forcing (Fig. 2) is a decrease in ice volume (Fig. 5).The ice discharge, D, remains significant as long as large parts of the ice sheet are marine terminating, but quickly diminishes when the ice sheet retreats onto land.From then on, the ice loss is dominated by a decreasing SMB, mainly driven by increased summertime meltwater run-off.The dashed lines in Fig. 5a represent the ice sheet evolution when simulations are continued without additional coupling with the regional climate model.Acceleration of the mass loss occurs after each coupling interval until the insolation maximum at 125 kyr BP (Figs. 2  and 5).Ablation rates along the southwestern margin (Fig. 6) are enhanced by the growing area of dark tundra that is exposed when the ice sheet retreats.Maximum local ablation values reach 3500 kg m −2 yr −1 .
A persistent feature throughout the Eemian is the local SMB maximum at 50 • W, 72 • N (Fig. 6a).This is caused by high bedrock elevation, forming a promontory in the ice sheet, enhancing snowfall.As such, it acts as a topographic pinning point for the ice sheet, as wet areas react less sensitively to warming than dry areas.Frequent snowfall events keep the surface albedo high, and more melt is needed before the winter snow is removed in the ablation season, thus delaying exposure and melt of the darker ice.This locally prevents a retreat of the ice sheet and as a result this promontory demarcates the northern limit of rapid ice sheet retreat.
Another consequence of topographically forced precipitation in the west is the early re-advance of the southwestern margin at 122.5 kyr BP, when the ice sheet still retreats along the dry northern margin (Supplement Movie).GrIS volume increase starts again at 121 kyr BP, following decreasing summertime insolation and Northern Hemisphere cooling (Fig. 2).
Although retreat of the northern margin occurs, ice loss in this area is limited in comparison to other recent model studies of Eemian GrIS evolution (Robinson et al., 2011;Fyke et al., 2011;Born and Nisancioglu, 2012;Stone et al., 2013).A strong retreat of the northern margin does not agree with available ice core evidence (see Sect. 5).According to our results, 64 % of the total volume loss occurred south of 70 • N by strong retreat from the southwestern ice margin (Fig. 6b  and c).Nevertheless, the southern dome survived, and remained connected to the central dome due to persistent high accumulation (SMB>0) along the southeastern ice margin, comparable to the present-day pattern.As a result, the remaining ice volume in south Greenland is substantially larger than in several other studies (Cuffey and Marshall, 2000;Huybrechts, 2002;Tarasov and Peltier, 2003;Lhomme et al., 2005;Otto-Bliesner et al., 2006, see Sect. 5).
In our simulation, this area of high accumulation along the southeastern margin maintains significant ice discharge into the Irminger Sea.This mass flux remains relatively constant throughout the Eemian (blue dotted line in Fig. 5b), and accounts for most (∼ 90 %) of the calving during the period of minimum ice sheet extent.Except for this southeastern section, the ice sheet loses contact with the ocean in the first few kyr of ice sheet retreat, after which the magnitude and pattern of ice loss is mainly determined by the SMB.During this latter period, details of outlet glacier dynamics are likely to be of lesser importance.As a result, the uncertainty in SMB outweighs other sources of uncertainty, such as that in ice discharge by marine-terminating outlet glaciers (Sect.4.3 and Fig. 7).

Ice sheet contribution to sea level
Figure 5b shows the total mass balance, MB, which is the sum of the surface mass balance and the ice discharge (MB = SMB − D).During the retreat phase (128-122 kyr BP), mass loss of the order of ∼ 180 Gt yr −1 is sustained for several millennia (Fig. 5b).This translates to an equivalent sea level rise contribution of 0.48 m kyr −1 .Although surface lowering in the ablation zone increases the surface mass loss, this effect is counteracted by the fact that the lowest parts of the ablation area, where run-off rates are highest, disappear when the ice margin retreats.This limits the magnitude of the mass loss that can be sustained, similar to what happens to retreating valley glaciers.
These feedbacks between changing ice sheet extent and its climate forcing also explain the difference with an earlier study using the same regional climate model ( Van de Berg et al., 2011) but assuming an unchanged present-day value of ice discharge and a fixed topography, resulting in a larger GrIS mass loss equivalent to 1.1 m SLR kyr −1 .Eemian SMB attains a minimum value of ∼ 60 ± 76 Gt yr −1 at 124 kyr BP in our simulation (Fig. 5b), which agrees closely with Van de Berg et al. (2011), although spatial differences exist due to the use of different ice sheet topography and extent.
The minimum ice sheet volume of 2.35 × 10 15 m 3 is reached at 121 kyr BP, and equals a +2.1 m contribution to eustatic sea level.This value has been obtained by comparing simulated ice volume at 121 kyr BP to the simulated steadystate ice volume for the present-day climate using the same regional climate-ice-sheet modelling approach (Helsen et al., 2012).This means that a zero contribution to sea level rise is Clim.Past, 9, 1773-1788, 2013 www.clim-past.net/9/1773/2013/assumed for an ice volume of 3.2 × 10 15 m 3 , which is 10 % larger than the observed present-day ice volume.This excess ice in both the simulated steady-state ice sheet (Helsen et al., 2012) and the Eemian ice sheet is predominantly found along the eastern margin (Fig. 6d) due to suboptimal representation of narrow outlet glacier dynamics.Since there is currently no ice in this area, we can safely assume that this excess ice was also absent during the minimum Eemian extent.Shown in green in Fig. 5b is the freshwater flux (FWF), which is the sum of the ice discharge and the ice sheet meltwater run-off (FWF = D + R), and therefore larger in magnitude than MB.A maximum FWF of ∼ 800 Gt yr −1 (0.028 Sv) occurs at the start of the coupled simulation, but FWF decreases due to the reducing influence of D, as the ice sheet margin retreats increasingly onto land.A relatively stable FWF of ∼ 600 Gt yr −1 (0.021 Sv) prevails between 128 and 124 kyr BP, but note that the FWF from the tundra is not included here, which likely is an increasing number due to the expanding tundra area.This simulated Eemian value is lower than the present-day FWF from Greenland (Bamber et al., 2012), which can be explained by the much smaller role of D in our simulation.The maximum Eemian FWF contribution from ice sheet run-off (350 ± 76 Gt yr −1 ) is of the same order of magnitude as the modern-day value.

Sensitivity experiments
We explored the uncertainty of our results by performing eight sensitivity experiments (Table 3 and Fig. 7).These experiments could not be performed in the same coupled setup as our control simulation, since this would require an additional ten 30 yr RACMO2 climate runs for each sensitivity experiment, which is computationally not feasible.We therefore conducted offline ice sheet model sensitivity experiments using the 10 control Eemian SMB climatologies, and calculated deviations of this SMB using the SMB gradient method (Helsen et al., 2012).This allows for the ice topography to evolve freely throughout the simulated period, while SMB adjustments are made as a function of the difference with the surface elevation in the control run.Note that this method is also used in the control simulation, in between the couplings with the climate model.

Start of the
One of the uncertainties is the start of Eemian warm period, for which 130 ± 1 kyr BP is a commonly used date (Jansen et al., 2007).Starting the simulations 1 kyr earlier/later obviously leads to an earlier/later onset of ice sheet retreat, but no additional nonlinear effects are introduced; that is, the rate of ice sheet retreat remains unaffected (Fig. 7a).This effect results in a 0.3 m uncertainty in the contribution to sea level at the time of minimum ice sheet extent.

Sliding
A moderate uncertainty is introduced by the parameterisation of sliding over the bed.We doubled and halved the sliding velocity by adjusting the sliding coefficient, resulting in a −0.4 to +0.6 m SLR uncertainly (Table 3 and Fig. 7a).It should be noted that the sliding coefficient is usually used as a tuning parameter.Testing the full model sensitivity for a change in this parameter would require a new tuning and initialisation for each new parameter setting.This would considerably reduce the sensitivity of our results.The experiment conducted here represents a sudden increase/decrease of sliding velocity at 130 kyr BP, providing an upper estimate of the uncertainty in sliding behaviour.

Marine basal melt
An extra mass loss term is introduced by allowing melt at the marine interface of the ice sheet.We tested two values of this melt rate, based on realistic modern-day ice front retreat observations, i.e. 100 and 1000 m yr −1 .The ice sheet model uses a fixed grid, so we cannot directly apply these values to grounding line retreat.Alternatively, these retreat rates are assumed to operate on the vertical ice front that is in contact with the ocean.The equivalent mass loss occurring at the iceocean-surface interface (in kg) is transferred to ice thinning of the grounded ice point bordering the ocean.This eventually results in ungrounding and thus calving in our model.Obviously this mechanism only affects marineterminating glaciers.As most of the ice sheet perimeter has retreated onto land at ∼ 128 kyr BP, the influence of marine basal melt ceases after the first 2 kyr of the simulation.The experiment with a marine basal melt of 100 m yr −1 did not significantly change the result of the control simulation.Using 1000 m yr −1 melt rate at the ice-ocean interface resulted in an additional ice volume loss of 0.3 m SLR.

SMB
Constraining the uncertainty resulting from the SMB forcing in our study is not straightforward.Van de Berg et al. (2011) report an uncertainty of 76 Gt yr −1 for their Eemian GrIS SMB reconstruction.We essentially use the same SMB product, although the SMB is based on a different ice sheet topography for each time slice.Since we calculate SMB adjustments based on a deviation between the simulated ice sheet elevation and the elevation (H s ) as used in the subsequent climate model simulations, we are able to introduce an additional SMB offset by shifting our H s -SMB functions over a certain H s .Tests using the present-day ice sheet topography also used in Van de Berg et al., 2011), revealed that a 100 m shift of H s results in a SMB of 76 Gt yr −1 .The sensitivity runs for SMB (Fig. 7b) are obtained using this ± 100 m offset in H s .This SMB uncertainty (contributing an additional +1.2 m SLR in the minimum SMB scenario) dominates the total uncertainty estimate (Fig. 7b), indicating that the simulated Eemian minimum is mainly determined by the accuracy of the SMB calculations.
The 76 Gt yr −1 uncertainty in SMB represents the uncertainty in the model results for RACMO2.On top of that, there is an uncertainty from the choice of the GCM used as lateral forcing fields for RACMO2.Here we chose ECHO-G, but a choice for a different GCM would result in different results.Next to that, the lack of (remnants of) the Laurentide Ice Sheet in ECHO-G may lead to an overestimation of early Eemian warming.These sources of uncertainty are hard to quantify, since it would need a suite of experiments using different lateral forcing fields, which has not been performed for Eemian conditions yet.Nevertheless, the Eemian climate run from ECHO-G is extensively validated against proxy data on temperature, which shows a ∼ 2 K warming in Europe (Kaspar et al., 2005).This is also in line with an independent reconstruction based on pollen data found in marine sediments in the Bay of Biscay (Sánchez Goñi et al., 2012).We can further compare simulated July temperatures on Baffin Island with proxy records as this is an area within the RACMO2 domain.We find 3 to 4 K higher July temperatures here, which is close to, but slightly lower than, the 4 to 5 K higher summertime temperatures as reconstructed from a compilation of palaeoclimatic evidence (CAPE-Last Interglacial Project Members, 2006), and also in agreement with proxies for lake water temperatures (Axford et al., 2011).The Eemian optimum climate from RACMO2 does show warming, but not uniform over the entire domain; cooling also occurs locally.For example, near-surface temperatures over the Labrador Sea show a distinct negative anomaly (down to 4 K) (van de Berg et al., 2013).Evidence for local oceanic cooling in this region is also found in Eirik Drift sediments and by palaeoclimate modelling (Winsor et al., 2012;Sánchez Goñi et al., 2012).Overall, the comparison with the available palaeoclimatic evidence yields confidence in the quality of the Eemian climate from ECHO-G.

Total uncertainty
We realise that the total model uncertainty might be larger than captured by the above-mentioned components.However, palaeoclimatic evidence from ice and marine cores can also be used to constrain the total uncertainty of our results (see Sect. 5).This comparison refutes a larger retreat of the southern part of the GrIS than the minimum SMB scenario, while larger retreat of the northern part of the GrIS is not in line with ice cores in this area (see Sect. 5).We therefore argue that the uncertainty in the SMB, and henceforth in the total GrIS contribution to SLR, is adequately estimated.
The uncertainty estimates in the different time series in Figs. 5, 7 and 8 reflect the RMS difference between the individual sensitivity experiments and the control run, assuming that the different contributions are independent.Note that the uncertainty in total mass balance (black line in Fig. 5b) is very small.This is due to the strong dependency of calving on SMB; experiments with high (low) SMB result in larger (smaller) ice sheets that have more (less) contact with the ocean, and hence a larger (smaller) calving flux.Thus, total ice sheet mass balance features much less variability in the different experiments than its individual components SMB and calving flux.

Comparison to other ice sheet simulations
In comparison with previous simulations of Eemian GrIS retreat, we have employed a different approach to force our ice sheet model.As the SMB forcing strongly influences GrIS response, it is not surprising that our results differ from those of earlier studies.Here we attempt to identify some key characteristics of our Eemian ice-sheet-climate simulation in order to clarify these differences.
As shown above, high accumulation along the southeastern margin prevented disintegration of the southern dome, notwithstanding substantial retreat of the southwestern ice margin.This contrasts with earlier work that resulted in higher (2.7-5.5 m) GrIS Eemian SLR estimates (Cuffey and Marshall, 2000;Huybrechts, 2002;Tarasov and Peltier, 2003;Lhomme et al., 2005).The (near) absence of ice in southern Greenland in these previous studies is likely caused by underestimated accumulation in the southeast, as these studies used the accumulation map from Ohmura and Reeh (1991), which is drier in this region.More recent studies, all using different GCM products as forcing fields, resulted in more ice cover in southern Greenland, and different magnitudes of ice sheet retreat in the north (Robinson et al., 2011;Fyke et al., 2011;Born and Nisancioglu, 2012;Quiquet et al., 2013;Stone et al., 2013).The northern region is vulnerable because of the low accumulation.Small changes in SMB result in large differences of ice sheet volume response.One  members, 2004) are shown as dashed grey lines; solid grey lines and shading represent results from the NEEM ice core record with its uncertainty (NEEM community members, 2013), and the lower (brown) dashed line in bottom panel denotes bedrock.Note that due to their close proximity, we do not make a distinction between GRIP and GISP2 in our results.
of the major issues concerns the calculation of the ablation, which is often based on a PDD model approach.Van de Berg et al. (2011) showed that it is problematic to correctly reproduce north-south gradients of ablation in different climate states with a PDD model.Therefore, Stone et al. (2013) suggested that a more physically based representation of insolation-driven ablation in the Eemian should result in a stronger northern ice sheet retreat than their own results.However, with a regional climate model as forcing, we find the opposite, i.e. less northern retreat.This discrepancy can probably be attributed to cooler conditions over northern Greenland in our regional climate model, compared to their set of GCM simulations.The latter might be due to a priori assumptions of their Eemian ice sheet topography, lacking ice cover in northern Greenland (Stone et al., 2010(Stone et al., , 2013)).
In general, our coupled regional climate-ice-sheet modelling approach offers important advantages over previous simulations.Firstly, the high resolution (18 km) of the regional climate model results in an explicit, physical and highresolution representation of the precipitation pattern.Moreover, the resolution is very similar to the resolution of the ice sheet model, which circumvents the need for downscaling.Next to that, the ablation is calculated using a physically based multi-layer snow model and energy balance calculations, with explicit albedo-melt feedback, circumventing the use of temperature-melt relationships such as a PDD model.Lastly, the two-way-coupling strategy allows for feedbacks between topography and atmosphere.Ideally, the coupling interval would be smaller than 1.5 kyr, but this was not possible due to computational costs.

Comparison to marine records
Well-studied marine sediment sequences have been retrieved from a site close to southern Greenland in Eirik Drift.This site provides evidence for a reduced but continuous presence of ice in southern Greenland during the Eemian (Carlson et al., 2008;Colville et al., 2011), which is in agreement with our model results.Also the timing of our simulated ice sheet minimum extent at 121 kyr BP matches with Eirik Drift sediments (Carlson et al., 2008;Colville et al., 2011).Next to that, De Vernal and Hillaire-Marcel (2008) use pollen records from this site to provide evidence of an extensive increase in vegetation in southern Greenland.This also supports our reconstruction of the exposure of a large tundra area in southwest Greenland.
Marine proxies from the same Eirik Drift site point to reduced Labrador Sea subsurface temperatures and reduced salinity between 130 and 122 kyr BP (Winsor et al., 2012).They attribute this to reduced Labrador Sea convection, due to increased FWF from Greenland, although the influence of sea ice anomalies are not ruled out.These findings are supported by climate models of intermediate complexity, pointing to a locally weakened deep convection in the Irminger and Labrador seas, as a result of surface water freshening during the Eemian (Cottet-Puinel et al., 2004;Bakker et al., 2012;Sánchez Goñi et al., 2012).Without additional FWF forcing, this temperature anomaly is also present in the combined ECHO-G-RACMO2 Eemian climatology with the present-day ice topography; cooling over the Labrador Sea occurs due to a local sea ice expansion and less entrainment of warm Atlantic water (van de Berg et al., 2013).
Our simulated FWF (Fig. 5b) is probably not large enough to significantly impact overturning circulation, although Bakker et al. (2012) suggest that a FWF of 0.039 Sv can locally weaken deep convection in the Irminger and Labrador seas.Although simulated annual mean FWF is lower than this value, summer peak values exceed this threshold easily.However, our modelling approach does not include oceanographic feedbacks, which precludes a robust analysis of the impact of ice loss on ocean circulation, and vice versa.

Comparison to ice core records
Figure 8 shows simulated Eemian surface height and temperature changes for five deep drilling locations together with reconstructions from ice core interpretation.The NEEM ice core is the only record from Greenland that covers the complete Eemian (up to an age of 128.5 kyr BP), which is shown as the dark-grey lines in Fig. 8.The other ice cores are assumed to contain Eemian ice, but no complete reconstruction of the Eemian climate can be retrieved due to disturbed chronologies (NGRIP members, 2004).Nevertheless, a consistent 3 ‰ higher isotopic composition in the (presumably) Eemian segments of the Camp Century, NGRIP, GRIP and GISP2 ice cores suggests only a minor drop (∼ 100 m) in surface elevation of central and northern Greenland during the Eemian (NGRIP members, 2004).This lower limit of ice sheet thinning is depicted by the grey dashed lines in the left panels of Fig. 8.This value agrees with the NEEM reconstruction, which points to a minimum elevation in the Eemian of 130 ± 300 m below the present-day level (NEEM community members, 2013).
The simulated elevation at NGRIP remained constant, while GRIP, NEEM and Camp Century experienced limited reductions in surface height (∼ 140, 250 and 300 m, respectively).If we compare these values with ice core reconstructions, the simulated change in surface elevation compares well for NEEM, only the timing of the minimum is 1.3 kyr later in our simulation.Also, at NGRIP our simulation remains within the maximum estimate of 100 m elevation decrease.Slightly larger thinning is found for Camp Century and GRIP, but remains within the uncertainty estimate for these locations.The larger (downward) uncertainty range demonstrates the sensitivity of surface elevation to mass loss rates in the case of stronger warming scenarios (see below).
The deep drilling sites are currently located on the ice divide.However, this was likely not the case during the entire Eemian; the ice divide moved southward at NEEM and Camp Century, whereas it moved northward from GRIP.NGRIP has been in the vicinity of the divide throughout the Eemian, resulting in a relatively stable surface elevation.The timing of the minimum surface elevation at the different locations is not synchronous, as a result of the out-of-phase response of the southwestern and northern parts of the GrIS.
For Dye-3 a maximum elevation drop during the Eemian of 500 m has been proposed (NGRIP members, 2004).Our results deviate significantly from this; in our reconstruction for Dye-3, ice-free conditions start at 123 kyr BP, but the distance to the ice margin never exceeds 2 grid points (∼ 57 km).Although Dye-3 is ice-free during ∼ 5 kyr in most of the sensitivity experiments (Fig. 7), the most conservative experiment retains a continuous, shallow ice cover at Dye-3.Our results therefore do not exclude that localised ice cover persisted at Dye-3, which would explain the preservation of ancient (> 400 kyr BP) biological material in the basal ice (Willerslev et al., 2007).Otherwise, our results compare favourably with reconstructions of elevation changes based on ice cores from central and north Greenland (NGRIP members, 2004;NEEM community members, 2013), which suggest that Eemian ice thickness of the north and central ice sheet was close to its present value, whereas ice thickness in the south was substantially reduced.
A different picture emerges when we compare our results to ice-core-derived temperature reconstructions (right panels of Fig. 8).NEEM community members (2013) deduce a maximum warming of 8 ± 4 K above the mean for the last millennium for NEEM, whereas NGRIP members (2004) reconstructed a +5 K Eemian optimum temperature for central and northern Greenland.Our simulated temperatures remain well below these temperature reconstructions (Fig. 8).The ice core record concerns mean annual temperature, which implies that summertime warming must have been even larger, considering the increased summertime insolation in Greenland during the Eemian.This would induce even larger melt rates.
An important factor in the temperature reconstruction from NEEM community members ( 2013) is that corrections are made for elevation changes such that temperature on a constant elevation (the present-day surface) is reported.Especially in the early Eemian, ice-core-derived elevation is above this value (and above our simulated elevation), which leads to higher reconstructed temperatures.Our simulated temperatures are actual surface temperatures, and hence affected by elevation changes.Temperature differences are still substantial between 125 and 117 kyr BP, but for this period our simulated temperatures overlap within the NEEM uncertainty estimate.In spite of these large differences in temperature, we do find extensive occurrence of surface melt at the NEEM location, in agreement with the ice core evidence (NEEM community members, 2013).
An Eemian warming of 8 K at NEEM would imply an additional 6 K warming at 124 kyr BP in our simulation.At that time, NEEM elevation is closest to the local equilibrium line altitude (ELA = 1582 m), though still 850 m higher.Taking γ atm = 7.4 K km −1 , this implies that 6 K additional warming would bring the ELA to approximately the local elevation.Although our control simulation does not come close to these conditions, a short duration of such strong warming cannot be ruled out.Experiments with such warming (not shown) reveal that the additional mass loss downstream of NEEM would be very large, and the associated surface drawdown (∼ 0.3 m yr −1 ) would lead to ablation conditions at NEEM within two centuries.This in turn is in disagreement with ice core evidence, which suggests only minor surface height changes and continuous accumulation conditions throughout the Eemian.
To summarise, the comparison of elevation and temperature changes reveals a profound disagreement.If Greenland temperatures in our simulation had been higher, this would not only improve the agreement with ice-core-derived temperatures but also induce significantly larger ice sheet retreat (and even conditions with negative SMB at NEEM), deteriorating the comparison with proxy elevation changes.It is difficult to reconcile the ice core interpretation of both a relatively stable ice sheet elevation in combination with a strong (5-8 K) Eemian warming over several millennia, certainly considering that summertime warming must have been even greater.
There are numerous complicating factors when a comparison is made between simulation results and proxy data: for example, ice-core-derived reconstructions of elevation and temperature rely on transfer functions involving assumptions on the relation between gas content and elevation, and the isotopic composition and temperature, respectively; a difference exists between elevation of the accumulation site and current ice core location; seasonality effects have to be taken into account.The presence of Eemian ice in ice cores from Camp Century, NEEM, NGRIP, GRIP and GISP2 that accumulated at altitudes similar to present day provides an important constraint on the estimate of total uncertainty estimate in our results; our minimum SMB scenario results in very strong surface lowering compared to the palaeo-evidence, which appears to preclude an even larger uncertainty in the downward (warm) direction.

Implications for sea level reconstructions
Our sensitivity experiments combined with the available geological evidence from ice (NGRIP members, 2004;Willerslev et al., 2007;NEEM community members, 2013) and marine cores (Carlson et al., 2008;De Vernal and Hillaire-Marcel, 2008;Colville et al., 2011) strongly suggest that the GrIS contribution to the Eemian sea level highstand was between 1.2 and 3.5 m, with a best estimate of 2.1 m, since this is the solution of the coupled model (Fig. 5).The timing of the minimum ice volume at 121 kyr BP also compares well with marine proxies, which suggest reduced input of meltwater from Greenland after ∼ 120 kyr BP (Carlson et al., 2008;Colville et al., 2011).Kopp et al. (2009) show an earlier eustatic sea level maximum, at 124 kyr BP, followed by a secondary peak at 118 kyr BP.This contrasts with Dutton and Lambeck (2012), who infer a later eustatic sea level maximum, but mainly stress that the current limited availability of far-field data and dating uncertainties make the reconstruction of a eustatic Eemian sea level curve uncertain.We find a single minimum in GrIS volume, which is an obvious result from the gradual transition between the Eemian optimum and glacial inception in our forcing fields, and the lack of possible oceanic feedbacks.Our results therefore also lack any indication of millennial-scale sea level oscillations, as found in several coral stratigraphies (Rohling et al., 2008;Thompson et al., 2011).These may also register variations in Antarctic ice volume, which is not considered here.
Freshwater sources other than the GrIS are needed to explain the global sea level maximum during the Eemian (Kopp et al., 2009;Dutton and Lambeck, 2012); a difference of 3.1 m remains between our maximum estimate (3.5 m) of the GrIS contribution and the minimum SLR of 6.6 m (95% probability) (Kopp et al., 2009).The thermosteric plus mountain glacier and ice cap contribution probably did not exceed 1 m (Radić and Hock, 2011;McKay et al., 2011), implying that the Antarctic Ice Sheet must have contributed at least 2.1 m.If we adopt our minimum GrIS estimate (1.2 m) and the maximum estimated sea level highstand of 9.4 m (33 % probability) (Kopp et al., 2009), this leaves the possibility of an Antarctic contribution of 7.2 m, implying a contribution from both the West and East Antarctic ice sheets (Bamber et al., 2009;Pollard and DeConto, 2009).
From our sensitivity experiments we infer a maximum rate of GrIS contribution to SLR between 0.34 and 0.66 m kyr −1 , with a best estimate of 0.48 m kyr −1 .This is only 9 % of the total rate of SLR of ≥ 5.6 m kyr −1 at 127 kyr BP.This value (Kopp et al., 2009) includes a contribution from the Laurentide Ice Sheet, as it is reconstructed for global sea levels between 0 and −10 m.After 126 kyr BP the rate of SLR decreases to ≥ 2.1 m kyr −1 (95 % probability) for sea levels ≥ 0 m (R. Kopp, personal communication, 2012); the GrIS contribution to this Eemian SLR totalled ≤ 24 %.Hence, also in terms of rate of change, our results require a significant Antarctic contribution to explain the estimates of the eustatic Eemian sea level curve.

Summary and conclusions
Based on a two-way-coupled regional climate-ice-sheet modelling approach, we conclude that GrIS retreat during the Eemian was dominated by SMB changes, and contributed 1.2-3.5 m to the maximum sea level highstand, with a best estimate of 2.1 m.This estimate broadly agrees with palaeoice-sheet reconstructions from marine records (Carlson et al., 2008;De Vernal and Hillaire-Marcel, 2008;Colville et al., 2011) and ice cores (NGRIP members, 2004;Willerslev et al., 2007;NEEM community members, 2013), although a substantial disagreement with ice-core-derived temperatures remains to be explained.The limited sea level contribution of the GrIS implies that a significant contribution of the Antarctic Ice Sheet is needed to explain the total Eemian sea level highstand.

Fig. 2 .
Fig. 2. The scaling factor (green) used to adjust the ECHO-G Eemian fields at the RACMO2 lateral boundaries towards glacial inception values are set proportional to the JJA (June-July-August) insolation at 60 • N (blue).Dashed green line indicates values of the scaling factor if it had also been coupled to insolation in the period before 125 kyr BP.

Fig. 3 .Fig. 4 .Fig. 5 .
Fig. 3.The dependence of ice sheet sensitivity to their temperature history: (a) schematic temperature forcings and (b) resulting ice sheet volume evolution.Dashed lines in (b) are the ice sheet volume response of the different forcings from (a).Each of these simulation is continued until the ice volume of the steady-state ice sheet (grey dashed line) is reached.Results in (b) are horizontally shifted such that they have identical ice volume at t = 0.After t = 0, all simulations are forced with the same Eemian climatology from Van de Berg et al. (2011) for 10 kyr.

Fig. 6 .
Fig. 6.SMB, minimum extent and ice thickness anomaly.(a) SMB pattern during Eemian optimum conditions, (b) simulated minimum GrIS at 121 kyr BP, (c) difference in ice thickness compared to present-day simulation and (d) difference in ice thickness compared to present-day measurements.An animation of GrIS evolution from 130 to 115 kyr BP is available online.

Fig. 7 .
Fig. 7. Simulated GrIS volume change and equivalent sea level contribution through the Eemian in different sensitivity experiments: (a) experiments using a different start of the Eemian warming and a perturbed sliding velocity, (b) experiments with additional mass loss due to ice-ocean contact and perturbations in the SMB forcing.Circles (squares) indicate the timing ice margin retreat (re-advance) at Dye-3.Complete retreat from the Dye-3 location did not occur in the SMB min experiment, when the ice margin remained exactly on that location during the minimum extent.Background shading represents the uncertainty estimate based on the RMS differences of the sensitivity experiments.

Fig. 8 .
Fig.8.Surface elevation and temperature at deep drilling sites.Time series of surface height change (left panels) and mean annual surface temperature (right panels) at five deep drilling locations, relative to present day.Coloured lines and shading represents results from our reference run with uncertainty.Estimates of minimum elevation and maximum air temperature as reconstructed from ice cores(NGRIP  members, 2004) are shown as dashed grey lines; solid grey lines and shading represent results from the NEEM ice core record with its uncertainty (NEEM community members, 2013), and the lower (brown) dashed line in bottom panel denotes bedrock.Note that due to their close proximity, we do not make a distinction between GRIP and GISP2 in our results.

Table 2 .
Ice sheet model parameter values.

Table 3 .
Summary of sensitivity experiments.Each experiment was carried out with one perturbed parameter.Resulting minimum ice volume and associated SLR are given.