Greenland Ice Sheet Model Parameters Constrained Using Simulations of the Eemian Interglacial

Using a new approach to force an ice sheet model, we performed an ensemble of simulations of the Greenland Ice Sheet evolution during the last two glacial cycles, with emphasis on the Eemian Interglacial. This ensemble was generated by perturbing four key parameters in the coupled regional climate-ice sheet model and by introducing additional uncertainty in the prescribed " background " climate change. The sensitivity of the surface melt model to climate change was determined to be the dominant driver of ice sheet instability, as reflected by simulated ice sheet loss during the Eemian Interglacial period. To eliminate unrealis-tic parameter combinations, constraints from present-day and paleo information were applied. The constraints include (i) the diagnosed present-day surface mass balance partition between surface melting and ice discharge at the margin, (ii) the modeled present-day elevation at GRIP; and (iii) the mod-eled elevation reduction at GRIP during the Eemian. Using these three constraints, a total of 360 simulations with 90 different model realizations were filtered down to 46 simulations and 20 model realizations considered valid. The pa-leo constraint eliminated more sensitive melt parameter values , in agreement with the surface mass balance partition assumption. The constrained simulations resulted in a range of Eemian ice loss of 0.4–4.4 m sea level equivalent, with a more likely range of about 3.7–4.4 m sea level if the GRIP δ 18 O isotope record can be considered an accurate proxy for the precipitation-weighted annual mean temperatures.


Introduction
Prediction of the future response of the Greenland Ice Sheet (GIS) to global warming is of great practical importance since the GIS can contribute up to 7 m to global sea level rise.On short (centennial) time scales, the response of GIS is primarily controlled by changes in surface mass balance (which can now be modeled relatively accurately) and by changes in fast flow (which is still poorly understood).On millennial time scales, when the GIS can lose a considerable portion of its volume, the situation is additionally complicated by a number of climate-ice sheet feedbacks.Explicit simulation of all of these processes requires the use of fullycoupled, high resolution Earth system models, which is still impractical on this time scale.In addition, information about the present-day GIS does not provide sufficiently strong constraints on its long-term evolution.Therefore, a study of past climate changes, and the response of the ice sheet to these changes, could help improve and better constrain ice sheets models.
The Eemian Interglacial (ca.130-115 ka BP) was characterized by high maximum summer insolation and warmer conditions in the high latitudes of both hemispheres.Paleo data suggest that sea level was higher than today by 4-6 m (Overpeck et al., 2006), or even as much as 6-8 m (Kopp et al., 2009).These numbers suggest a considerable contribution from both the Greenland and Antarctic ice sheets.At the same time, the presence of Eemian and older ice in several ice cores indicates that a large portion of the GIS survived the Eemian Interglacial.Moreover, gas concentrations in Eemian ice from the Greenland summit can be interpreted to show that the height of the summit was not much lower during the Eemian compared to the present day (Raynaud et al., 1997).These data can potentially provide useful constraints for ice sheet models.
Published by Copernicus Publications on behalf of the European Geosciences Union.
A. Robinson et al.: Greenland ice sheet model parameters constrained using the Eemian A number of attempts to simulate the response of the GIS to Eemian climate conditions have been made with different models and approaches (Letréguilly et al., 1991;Cuffey and Marshall, 2000;Huybrechts, 2002;Tarasov andPeltier, 2002, 2003;Greve, 2005;Lhomme et al., 2005;Otto-Bliesner et al., 2006).These works reveal large uncertainties in the simulated contribution of the GIS to the sea level highstand during the Eemian Interglacial, ranging from almost no contribution to over 5 m sea level.The main problem with simulating the ice sheet's evolution during the Eemian Interglacial is the definition of surface boundary conditions for the ice sheet model.At a minimum, the ice sheet models require prescription of the seasonal variations of temperature and precipitation when using the standard positive degree day approach.However, not only do the Greenland ice cores lack continuous records of temperature and precipitation through the entire Eemian Interglacial, they also do not extend any further back in time.Simulation of the GIS response to Eemian climate conditions also requires data from the previous glacial cycle to properly initialize the ice sheet.In addition, even the existing record only provides information about precipitation-weighted temperatures (annual mean temperatures weighted more heavily for months with higher precipitation), while models require annual mean and, even more importantly, summer temperatures.Several attempts have been made to construct the temporal evolution of climatic forcing by combining Greenland climate reconstructions for the last glacial cycle and Antarctic ice core records for the penultimate glacial cycle (Cuffey and Marshall, 2000;Huybrechts, 2002;Tarasov and Peltier, 2003;Greve, 2005).Although the Greenland and Antarctic temperature records do bear a certain similarity, this is still a rather crude approach, since both the magnitude and temporal dynamics of temperature changes in the high latitudes of the northern and southern hemispheres differed considerably.Moreover, the most crucial characteristic of summer temperature anomalies at the margins are not well constrained by the existing ice core records.In particular, the δ 18 O content of Eemian ice may have been affected by both changes in surface elevation and the seasonality of precipitation (Jouzel et al., 1997).This approach also does not explicitly account for the additional effect of insolation changes, which are of comparable importance to temperature changes on the surface melting of ice for this period (Robinson et al., 2010).A coupled climate model (with an additional downscaling procedure) can provide all necessary climate information needed to force an ice sheet model, but running coupled GCMs through the whole glacial cycle, or even just for the interglacial period, is still computationally too expensive.Until now, only a time-slice simulation has been used to force a model of the GIS in this way (Otto-Bliesner et al., 2006).However, the GIS has a multimillennial time-scale response to climate change and the climate did not remain constant during the interglacial.Such an approach also does not solve the problem with the initialization of the GIS at the onset of the Eemian Interglacial.
Earth system models of intermediate complexity (EMICs, Claussen et al., 2002) are useful tools for the study of the evolution of ice sheets on millennial to orbital time scales, since some of them (like CLIMBER-2, Petoukhov et al., 2000) are computationally efficient enough to be run through the glacial cycles (e.g., Bonelli et al., 2009;Ganopolski et al., 2010) and the only required boundary conditions are readily available information: Earth's orbital parameters and the atmospheric concentration of greenhouse gases.The problem is, however, that this type of model has a very low spatial resolution and it can neither accurately simulate the regional climate over the GIS, nor the feedbacks associated with the melting of the ice sheets.To overcome this problem, we have developed a regional energy and moisture balance model (REMBO, Robinson et al., 2010), which produces a physically-based downscaling of the climate over the GIS using a present-day climatology, and anomalies simulated by a coarse resolution climate model.Unlike most previous studies, we have replaced the traditionally used positive degree day (PDD) melt scheme with a surface melt scheme that explicitly accounts for changes in insolation.The latter is crucial for simulation of the Eemian Interglacial when summer insolation was much higher than at present.
The aim of this study is to investigate whether existing paleo data from the Eemian Interglacial, in conjunction with present-day data about the GIS and its mass balance, provide sufficient constraints on the choice of model parameters and the GIS contribution to the sea level high-stand during the Eemian Interglacial.To this end, we performed an ensemble of simulations of the GIS over the last two glacial cycles by varying two ice sheet model parameters, one melt scheme parameter and one regional climate model parameter (four model parameters in total).In addition, we varied the magnitude of the Eemian background warming.Present-day observations and paleoclimate data have been used to select the subset of model versions that is consistent with empirical constraints.

Model setup and experimental design
The coupled REMBO-ITM-SICOPOLIS model used in this study is identical to that described by Robinson et al. (2010), except for a slight modification to the surface albedo parameterization (see Appendix A).A brief description of the main components is provided below.REMBO (Regional Energy-Moisture Balance Orographic model) produces daily climatological fields of temperature and precipitation for Greenland with 100 km resolution.Monthly ECMWF Reanalysis data (ERA-40) of 2-m temperature and relative humidity (Uppala et al., 2005), averaged over the period 1958-2001 and linearly interpolated in time to produce daily fields, are used as boundary conditions around Greenland.In a previous analysis, Robinson et al. (2010) show that REMBO is able to capture the present-day seasonal cycle of temperatures over Greenland and that it improves the fit with observations as compared to conventional bilinear parameterizations.While the precipitation field mismatches observations in several local areas, the large-scale field is correct in scale and sufficient enough for long timescale simulations.This lends confidence to the idea that present-day climatic forcing is fairly accurate, and that the energy reaching the surface determined by the model would not be far from reality.
The daily temperature and precipitation output from REMBO is used to determine the surface mass balance of the ice sheet and the evolution of surface albedo via the surface mass balance model ITM (Insolation-Temperature Melt).ITM represents an alternative to the commonly used PDD (Positive Degree Day) approach (Braithwaite and Olesen, 1989;Reeh, 1991), in that both are computationally efficient and semi-empirical.However, ITM has a more physical basis that provides two advantages: it explicitly accounts for short-wave radiation (while PDD does implicitly, via different melt coefficients for snow and ice); and, since it incorporates a simple snowpack model, it has a memory which the standard PDD scheme does not.The ITM model is tuned to use climatological fields of temperature, precipitation and insolation at the top of the atmosphere and is driven by the elevation-corrected output of REMBO.ITM simulates surface mass balance and mean annual ice sheet surface temperature on the same grid used by the ice sheet model.
The ice sheet model used for these simulations is SICOPOLIS, version 2.9 (Greve, 1997a,b).It is a 3-D thermomechanical ice sheet model based on the shallow ice approximation (SIA) and it is comparable to other SIA models (e.g., Huybrechts et al., 1991).SICOPOLIS has been used in numerous studies of ice sheet evolution both in the past and future (e.g., Greve, 2000Greve, , 2005;;Calov et al., 2002Calov et al., , 2005)).It is run on a 20 km resolution grid and is forced from above by the annual surface mass balance and surface temperature fields obtained from REMBO.At the base of the ice sheet, the geothermal heat flux is prescribed, and the lithosphere deforms locally with a time lag of 3 ka.
In this study, we use the results of transient simulations of the last two glacial cycles performed with the intermediate complexity model CLIMBER-2 to drive REMBO.The model and experimental setup is identical to that described by Ganopolski et al. (2010), except that here we only apply the results from the last two glacial cycles.CLIMBER-2 (Petoukhov et al., 2000) is a coarse-resolution EMIC that has been shown to produce a realistic present-day climate, climate sensitivity (Ganopolski et al., 2001) and changes in climate over the last glacial cycle (Ganopolski et al., 2010).The CLIMBER-2 simulation was forced by variations in Earth's orbital parameters, computed following Berger (1978) and by a prescribed time-series of the equivalent atmospheric CO 2 concentration.The resulting global temperature anomaly time series agrees well with paleoclimate reconstructions, in that, for example, the Greenland annual temperature at the Last Glacial Maximum (LGM) was approximately 20 • C cooler than today and the Eemian Interglacial was ca. 2 • C warmer than today.
To force the transient simulations performed here, we used the same insolation at the top of the atmosphere (July insolation at 65 • N is shown in Fig. 1a as an illustration) and equivalent CO 2 concentration (Fig. 1b) as in the CLIMBER-2 run.The boundary temperature field of REMBO was modified by adding the spatially uniform, regional monthly temperature anomaly computed by CLIMBER-2 to the presentday climatological fields (the anomaly time-series is shown in Fig. 1d).We assumed no changes (compared to the present day) in relative humidity at the Greenland borders.Analysis of the CLIMBER-2 output shows that relative humidity around Greenland changed by a maximum of 5% during glacial cycles, while during the Eemian, it differed negligibly from that of the present day.The anomalies from the CLIMBER-2 simulation provide us with self-consistent, regionally relevant time-series and enable us to avoid problems involved with relying on discontinuous oxygen isotope records.
A. Robinson et al.: Greenland ice sheet model parameters constrained using the Eemian Changes in sea level were prescribed using an appropriately scaled SPECMAP time-series (Imbrie et al., 1984), shown in Fig. 1c.The ice sheet extent was limited to areas above the contemporaneous sea level.We found that the inclusion of sea level changes allowed more ice to grow during the glacial periods, but it played only a small role in the evolution of the ice sheet in warmer periods.Nonetheless, we included sea-level adjustments to the margin to allow for its effect on total volume.Note that sea level changes modeled by CLIMBER-2 agree favorably with the SPECMAP reconstruction and, in principle, could instead have been used to drive the ice sheet model without any appreciable differences.
All simulations began at 250 ka BP and were run until the present.The ice sheet model was initialized with the presentday GIS topography from Bamber et al. ( 2001) and the standard set of model parameters as initial conditions.Since the model was run for more than 100 ka before the Eemian Interglacial, there should be no dependence on the choice of initial conditions.

Modern and paleo empirical constraints
Assuming that the climate obtained from REMBO is fairly well represented but other processes (such as fast flow in the ice sheet model) are ignored, we believe that there are then at least two aspects of the ice sheet that are possible to model with reasonable certainty and that can provide information about the sensitivity of various model parameters: the diagnosed, present-day partition of GIS mass balance between surface melting and ice discharge into the ocean (i.e., calving), and the elevation (and elevation changes) of the summit of the GIS.Arguably, both of these characteristics are less affected by the lack of fast ice streams, ocean interaction and low resolution than the spatial extent and the volume of the GIS.
In addition, we would like to explore the possibility of using available paleoclimate information from the Eemian Interglacial to further reduce the range of possible combinations of model parameters.This leads to the three constraints used in this study, described below.

Present-day surface mass balance partition
The first constraint is the estimate of the present-day partition between GIS surface mass balance (SMB) and ice discharge into the ocean for the fixed topography of present-day Greenland.As we do not explicitly model calving (ice that flows into the ocean is simply considered to have calved), we use the ratio of diagnosed surface mass balance to the total precipitation over the GIS.Since the SMB is positive and the present-day ice sheet is thought to be in equilibrium, this mass is eventually lost to ice discharge, while the remaining mass is lost via surface runoff.The proposed ratio thus indicates the present-day balance between runoff and ice discharge.The benefit of using this ratio is that (1) it eliminates biases in the absolute values that result from different ice sheet surface area masks, (2) allows an analysis of the reasonable sensitivity range of the surface mass balance model without incorporating additional uncertainty from the ice sheet model or paleoclimate forcing and (3) keeps the focus on the partition between surface mass balance and ice discharge.We believe the last consideration is especially important for properly estimating the correct sensitivity of the ice sheet to long-term climatic forcing.We also found that the values of the surface mass balance partition for the modeled present-day ice sheet topography are rather similar to those obtained using the prescribed actual topography.We have used the latter values to eliminate additional biases related to the ice sheet model, but this criterion appears to be robust in this sense.Furthermore, it should be noted that such a criterion implicitly assumes that the distribution of melt and precipitation is largely correct.REMBO has been validated against present-day observations and compares well to results from regional climate models (Robinson et al., 2010).
For the present-day GIS, several estimates from regional climate models (RCMs) show that positive surface mass balance should account for about 50% (ranging from 48-63%) of the total incoming precipitation (Box et al., 2006;Fettweis, 2007;Ettema et al., 2009), while the remaining incoming mass is lost via runoff.These results do not necessarily encompass the range of all possible values, and furthermore our approach does not account for all physical processes at the surface (e.g., blowing snow or evaporation).Thus, we allow for additional uncertainty by choosing the range 40-65%.

Present-day GRIP elevation
The present-day summit elevation at GRIP (73 • N, 38 • W, 3230 m) is a robust feature of the ice sheet that provides a useful constraint for our ensemble.Because it is located in the middle of the ice sheet at essentially the thickest location, changes in this elevation reflect large-scale changes in the ice sheet surface mass balance, rather than highly dynamic changes that occur near the margins (Alley et al., 2010).Raynaud et al. (1997) indicated that even with a margin position varying in the range of 100-200 km, no more than 100 m elevation difference is modeled.A similar sensitivity is also reflected in our own simulations.Although biases in simulated accumulation at the summit can affect the modeled elevation, the precipitation field obtained from REMBO agrees best with observations for the large-scale field over the interior of the ice sheet.Indeed with several combinations of model parameters, it is possible to model the summit location and elevation correctly.Thus, as a second constraint, we assume that the present-day GRIP elevation should be obtained to within ±100 m.

Eemian GRIP elevation
The ice core drilled at the Greenland Ice Sheet summit (GRIP, 1993) provides important information that can be used to further constrain the paleo simulations.A reliable time reference for this ice core cannot be determined for ice older than ca. 100 ka, which is likely due to stratigraphic folding of the ice (Alley et al., 1995).However, the total gas content in the ice core for the Eemian period is comparable to that of the Holocene, with slightly higher values also present (Raynaud et al., 1997).This indicates that the minimum elevation of the GRIP location during the Eemian was very likely no more than a few hundred meters less than that of the present day, but there is also approximately 350 m of uncertainty in the measurements (Cuffey and Marshall, 2000).Thus, an elevation loss greater than 300-600 m would be difficult to explain.In this study, we apply the constraint that the elevation at GRIP was at most 400 m below present day, but we also test how choosing different values for this constraint affects the selection of valid ensemble members.
Using the three constraints outlined above, we are able to determine which model simulations are more likely to be realistic, and more important for the future stability of the ice sheet, which model parameter combinations should be considered valid.

Additional possible paleoclimate constraints
It could also be possible to further constrain the ensemble using information from other ice core locations, such as Camp Century (77 • N, 61 • W, 1890 m) or DYE-3 (65 • N, 43 • W, 2490 m).However, our model's resolution, lack of fast processes and spatially-constant boundary temperature anomaly means that we have less confidence in our ability to accurately capture the behavior of the ice sheet in these locations.This is complicated by the fact that little is known with certainty about regions closer to the margin.For example, DYE-3 may or may not have remained ice-covered during the Eemian Interglacial (Jansen et al., 2007;Alley et al., 2010).Eemian ice layers do exist at the base, however older ice has not been found.Meanwhile, δ 18 O anomalies indicate an elevation difference of perhaps 500 m (Johnsen et al., 2001;NGRIP, 2004), and DNA evidence from the silty layers beneath the ice sheet indicate that plant growth only occurred much earlier, perhaps during MIS11 (Willerslev et al., 2007).The latter does not rule out the possibility of ice-free, permafrost conditions at this location, as opposed to only a reduced-thickness ice sheet.In light of the existing controversy and our model's poor representation of these locations, we do not consider data from them as hard constraints.Nonetheless, we will comment on the results from different model versions for these locations.1 provides a list of the parameters and the values used in this study.

Geothermal heat flux
The geothermal heat flux (the lower boundary condition of the ice sheet model in the bedrock) was set to a constant value everywhere.Little is known about the exact values underneath the ice sheet, but the large-scale field is likely to be fairly uniform and recent estimates put the value near 50 mW m −2 , with some variation around the margins (Shapiro and Ritzwoller, 2004;Vinther et al., 2009).It cannot be ruled out, however, that this value is significantly larger in some localized areas below the GIS (e.g., Fahnestock et al., 2001).To account for the uncertainty in the broad sense, we chose geothermal heat flux values of 50, 60 and 70 mW m −2 .The geothermal heat flux is one factor that determines how much of the ice sheet base is temperate (at the pressure melting point) and thus is sliding at the base.

Sliding coefficient
The sliding law used in this study, is a Weertman-type equation with a third-order dependence on the gradient of elevation, where v s is the ice sliding velocity, H is the ice thickness, ρ is the density of ice, g is the gravitational force and h is the surface elevation.purposes, the constant γ s was chosen as the uncertain parameter varied in this study.Increasing this parameter has the effect of increasing the ice velocity at the base of regions of temperate ice.This should occur mostly around the margins and it acts to adjust the slope of the ice sheet.Higher values of sliding also generally decrease the total volume of ice.

Melt model parameter, c
The surface melt of snow and ice has been calculated in ITM from a simple energy balance equation (van den Berg et al., 2008), where the potential melt rate, M s , is assumed to largely derive from two main contributions.In the first term, representing incoming short-wave radiation at the surface, S is the incoming solar radiation at the top of the atmosphere, τ a is the transmissivity of the atmosphere and α s is the surface albedo.
The transmissivity is assumed to be a linear function of elevation (Robinson et al., 2010) and the surface albedo is calculated following the method described in the Appendix.The second term c + λ T is a linear parameterization of the longwave radiation and turbulent heat flux, where T is the daily mean temperature and λ and c are constants.The choice of λ derives from the empirical value of long-wave radiation absorbed by snow and ice, which corresponds to a melt rate of 3 mm w.e. per degree (analogous to the choice of degree day factor in the PDD approach, see Reeh, 1991).The remaining term, c, is then assumed to be a free constant parameter, the value of which can vary widely depending on the domain in question (van den Berg et al., 2008) and the albedo parameterization used.Robinson et al. (2010) determined that for modeling present-day conditions over Greenland, c = −55 W m −2 gives the best partition of surface mass balance components.When performing equilibrium simulations fully coupled to an ice sheet model, however, it was found that several values of c could still produce a reasonable present-day ice sheet.Therefore, in this study we chose to incorporate melt model uncertainty into the ensemble via the parameter c, allowing a range of −45 to −65 W m −2 .
We further considered that the parameterization of surface albedo (as described by Robinson et al., 2010) is also imperfect and is likely the source of great uncertainty.However, several simulations combining various values of the melt parameter c with changes to the surface albedo scheme show that the first order changes to surface albedo can be captured simply by changing the value of c.

Atmospheric moisture diffusion coefficient
Uncertainty in the climatic forcing likely plays the key role in determining both the past and future stability range of the GIS.The parameters used in REMBO to determine temperature are tuned to match present-day observations, with both annual mean temperatures and the seasonal cycle well reproduced (Robinson et al., 2010).Therefore, no uncertainty was considered in the temperature equation itself; rather, uncertainty in the temperature forcing is represented via the paleo factor (see below).To introduce uncertainty in the representation of precipitation by REMBO, we chose to vary the strength of the overall precipitation field.This was achieved via the moisture diffusion constant, κ Q .By decreasing this value, less moisture at the boundaries was able to diffuse inward, decreasing the total amount of precipitation.Although this does not improve the local deficiencies of the model, it does allow adjustment of the large-scale field (especially over the interior of the ice sheet), which plays a role in determining the total volume and the maximum elevation of the ice sheet.

Paleo climate forcing
When performing paleo simulations using an anomaly approach, great uncertainty can exist, depending on the choice of anomaly forcing.In REMBO, changes in precipitation are controlled by changes in temperature and elevation, while relative humidity is prescribed.This means that aside from the prescribed variations in CO 2 equivalent radiative forcing, the only anomaly forcing needed as input is for the temperature field.We chose to prescribe a spatially-uniform monthly temperature anomaly at the boundaries of Greenland.The anomaly temperature time-series was determined in the glacial cycle simulations of CLIMBER-2 by Ganopolski et al. (2010).We obtained temperatures by averaging values from the two CLIMBER-2 grid cells that encompass Greenland for each month and interpolating these into daily values for input to REMBO.
Our approach does not account for possible regional differences in temperature changes around Greenland, which would introduce additional uncertainties into the analysis.GCM simulations of the Eemian disagree about both the pattern and magnitude of warming.In the absence of more information and to be consistent with previous modeling studies, we therefore applied a spatially constant anomaly to the boundary temperature field.
The maximum Eemian summer temperature anomaly in the "Greenland" grid cells simulated in CLIMBER-2 is just under 2 • C (Fig. 1d).Due to its coarse spatial resolution, it is likely that CLIMBER-2 underestimates the amount of local warming over the GIS at that time.The CAPE project (CAPE, 2006) compiled a map of Eemian summer temperature anomalies obtained from marine sediment cores for several locations in the Arctic region.Most data show that warming was anywhere from 2-5 • C, compared to pre-industrial values.Several coupled GCMs used to performed Eemian climate simulations give the range of 2.5-5 • C for summer temperature anomalies over the GIS as well (Masson-Delmotte et al., 2010).To account for potentially higher Eemian warming than simulated by CLIMBER-2, we chose to apply a scaling factor to the positive temperature anomalies.Negative anomalies were not modified in any way, since the negative temperature anomalies during glacial periods are well simulated and also much less important.This approach is intended to provide a reasonable adjustment to the climate anomalies to account for uncertainty.As shown by Cuffey and Marshall (2000), it is not so much the duration of the warm period that determines the mass loss during the Eemian, but rather the maximum temperature anomaly.Therefore, we applied factor values of 1, 1.5, 2 and 2.5, resulting in maximum prescribed Eemian summer warming around Greenland of 1.7, 2.5, 3.4 and 4.3 • C, respectively.The resulting Eemian temperature anomaly time series for summer are shown in Fig. 2.

Transient simulations of the GIS
Permuting the above five parameters using the values in Table 1 produced 360 simulations, which correspond to 90 independent model versions (since the paleo factor only modifies the boundary forcing and, thus, does not produce additional model versions).Figure 3 shows the temporal evolution of the ice sheet computed in the ensemble of model simulations over the last two glacial cycles.During glacial periods, different model versions produce rather similar results, which is not surprising since under cold climate conditions, the GIS occupies the whole land area and surface melt is essentially absent.However, it appears that any reasonable combination of parameter values (based on present-day tuning) can result in dramatically different evolution histories for the GIS.With some model versions, the GIS melts entirely during the Eemian and with others, minimal changes occur for the same period (see Fig. 4).Similarly, simulated precipitation-weighted mean annual temperatures over the summit remain very close in all model versions over the glacial period, but during the Eemian, they differ by up to 15 • C, which is much more than the differences imposed by using different temperature anomalies around Greenland (i.e., different paleoclimate factors).This large range in simulated temperatures is explained by both changes in surface elevation and changes in surface albedo (in the model versions where a substantial portion of GIS melts away during the Eemian).The annual precipitation at GRIP (Fig. 3d) shows changes largely consistent with other estimates (e.g., Tarasov and Peltier, 2003).At the last glacial maximum, precipitation was about four times less than at present -mainly a result of the much lower glacial temperature.During the Eemian, the total precipitation at GRIP was strongly dependent on the temperature anomaly and the changes in topography.
In Q geo = 50 mW m −2 ).For the present day, all simulations have a similar distribution of ice, covering almost the entire land area.The total present-day volume is simulated to be 2% smaller than the present to up to 35% larger, mainly due to the additional ice at the margin.The interior distribution of ice in all simulations reflects the present day reasonably well, in terms of elevation and surface slope.By contrast, the minimum volume simulated for the Eemian Interglacial differs drastically between the five model versions.Depending on the melt parameter c, the GIS can experience minimal changes or melt almost completely.The age at which the minimum volume is reached is related to the total amount of ice melted and can vary between 124 ka BP and 121 ka BP.When more volume is lost, this state is maintained longer and the minimum is reached later in time.
Figure 5a-e shows the precipitation-weighted annual temperature distribution for the same simulations when the Eemian minimum GIS volume has been reached.Over the ocean, the temperature anomaly is prescribed, however the inland distribution is largely affected by the distribution of ice.In the cases where a large portion of ice disappeared, anomalies of up to 16 • C can persist long after the boundary warming has decreased.This temperature anomaly is largely due to elevation changes, but changes in surface albedo also affect it.
The precipitation field is also strongly affected by topographical changes, as shown in Fig. 5f-j.In REMBO, precipitation is assumed to increase for regions with strong surface gradients.For simulations in which the ice sheet retreats, precipitation near the margins of Greenland is reduced, but it increases near the margins of the ice sheet.The precipitation anomalies closely follow the Clausius-Clapeyron relation with respect to temperature, for a GIS geometry similar to today.However, when the GIS elevation and/or extent is considerably different (as was possible during the Eemian), orographic effects also become very important -and are accounted for by REMBO.Our approach cannot account for large-scale changes in atmospheric circulation (e.g., storm tracks) outside of Greenland, but this is not different from other approaches so far.Even accounting for these large changes in precipitation, overall changes in the surface mass balance of GIS during the Eemian were primarily controlled by temperature and insolation changes.
It is interesting that model versions which produce more realistic present-day GIS distributions (less ice volume and area, e.g., c −45 W m −2 ) simulate an almost complete disappearance of the GIS during the Eemian, even under very modest warming anomalies (f p = 1).This is unrealistic, but in line with arguments presented by Robinson et al. (2010).Due to the lack of fast processes and/or model resolution, too much surface melt is required to model the GIS extent close to the observed one.This violates the observed partition between surface melt and ice discharge, and shifts the model much closer to an instable threshold than it should be in reality.Therefore the "realistic" (in the traditional meaning of this term) simulations of the modern GIS provide no constraints on the magnitude of Eemian melting.However, as we discuss below, other constraints indeed help to narrow the range of valid model parameters.

Constraining the model parameters
Figure 6 is a schematic representation of the constraints as they apply to all simulations in the ensemble.Using this plot, we are able to identify which constraints are responsible for rejecting different simulations.For example, the present-day surface mass balance partition serves to eliminate the most and the least sensitive melt parameter values (−45 and −65 W m −2 ), whereas the Eemian summit constraint generally only eliminates the more sensitive melt parameter values.This is discussed in more detail below.
It is important to distinguish between the set of simulations and the set of different models in the ensemble.This distinction arises because each model version produced four simulations, corresponding to the modified Eemian climate (via the paleo factor, f p ).Thus, a simulation is considered valid if it does not violate any constraints.A model version is considered valid if at least one of its simulations has been considered valid.
Figure 7a shows that the dependence of the modern surface mass balance partition on the melt parameter c is essentially linear.In our ensemble, only the moisture diffusion Fig. 6.Schematic table of constraints applied to the paleo simulations of the Greenland Ice Sheet.Each color band contains 72 simulations, corresponding to the choice of melt parameter c, where lighter shades correspond to more sensitive secondary parameter combinations for added distinction (this color scheme is used for all plots).The first row shows all simulations.The middle rows correspond to specific constraints -so if a simulation is consistent with the constraint, it is plotted in that row (otherwise, white regions indicate rejected simulations).The last row shows the simulations that were consistent with all constraints and thus are considered valid.Figure 7b shows the present-day GRIP elevation versus geothermal heat flux (Q geo ).Of the five parameters varied, the geothermal heat flux plays the largest role in determining the elevation at GRIP.In fact, with a value of 60 mW m −2 , already only a couple of simulations are able to obtain a present-day GRIP elevation greater than 3200 m.The more likely valid choice for Q geo would thus be 50 mW m −2 , which conforms to recent estimates (e.g., Shapiro and Ritzwoller, 2004;Vinther et al., 2009).NGRIP (75 • N, 42 • W, 2920 m), which lies only a few hundred kilometers away along the ice divide, shows a similar relationship to geothermal heat flux, although some more simulations using 60 mW m −2 are able to produce the right elevation.For comparison, we also looked at the present-day elevations at Camp Century and DYE-3; however, for reasonable simulations, no systematic relationship between Q geo and elevation could be found.The third row in Fig. 6 shows that this constraint filters out simulations with high values of geothermal heat flux.
Figure 7c shows the difference between the modeled minimum Eemian and present-day elevation at GRIP.The melt model sensitivity is the strongest factor that determines the Eemian GRIP elevation reduction, followed by the sliding coefficient.For many simulations, all ice is lost at GRIP and these are clearly too sensitive.In fact, none of the model versions with c = −45 W m −2 and only a few with c = −50 W m −2 satisfy this constraint.Therefore, we are able to exclude essentially the same subset of the most sensitive model versions as excluded by the constraint on the presentday SMB partition.This does not affect models with a small response of the GIS to Eemian warming.Data from δ 18 O records indicate that the annual mean temperature at GRIP (or, more precisely, the precipitation-weighted mean annual temperature) was 4-6 • C warmer than at present (Johnsen et al., 2001).Using this information would additionally reduce the number of valid simulations (primarily by only permitting runs with high paleoclimate forcing and excluding model versions with a low sensitivity to Eemian warming).It is very encouraging that both paleoclimate constraints and the modern constraint on the SMB partition were consistent in the simulations that were eliminated.However, even for the reduced range of valid model runs, a wide range of possible GIS responses to Eemian warming remains.

Discussion
From 360 simulations, the above three constraints reduce the ensemble to just 46 valid simulations and 20 valid combinations of model parameters.In Fig. 4, showing the transient evolution of the GIS, the lighter colored lines indicate "invalid" simulations.The melt parameter c and the magnitude of boundary warming essentially determine the amount of ice lost in the Eemian, while the other parameter choices play only minor roles.The applied constraint on the GRIP elevation reduction is able to eliminate the most sensitive model versions and reduce the valid range of Eemian melt.
Figure 8 shows the various margin positions of the valid simulations when the Eemian minimum volume was reached.All simulations agree in that at least a significant northern dome centered around GRIP and a smaller southern dome must have remained ice covered through the last interglacial.Depending on the sensitivity of the model and the boundary warming, these two domes could have been separated by an ice-free region, they could have been joined or the ice sheet could have maintained coverage over all of Greenland.Most simulations show ice coverage of the four northern ice core locations (GRIP, NGRIP, NEEM and Camp Century).At DYE-3, several configurations are possible, indicating that the region is very sensitive to model parameters and climatic forcing.
In terms of overall volume loss, Fig. 9 shows the fraction of the ice sheet that melted during the Eemian versus the maximum warming experienced at GRIP for the same period 6 0 °N 7 0 °80 °-75°-60°-45°-30°-15°E Fig. 8. Eemian minimum ice margin contours for the 46 accepted simulations.The colors correspond to the choice of parameters (as in Fig. 6).The black diamonds show the locations of ice core drill sites (from South to North: DYE-3, GRIP, NGRIP, NEEM and Camp Century).
(invalid results are shown by the much lighter crosses).For the ensemble of valid simulations, the range of ice lost during the Eemian was 5-60% (or 0.4-4.4m sea level), relative to the present-day modeled values.The cluster of valid points can be filtered corresponding to each choice of the paleo factor, f p , which was used to increase the Eemian temperature anomaly at the boundaries.As f p increases from 1.0 (1.7 • C) to 2.5 (4.3 • C), the maximum ice loss increased from about 45% to 60%, and the minimum estimate increased from 5% to 50%.Furthermore, the ice core record at GRIP implies that the annual (precipitation weighted) temperatures anomaly reached 5 • C (Johnsen et al., 2001).The only simulations that present an anomaly at GRIP close to 5 • C are those that had a lower GRIP elevation and lost more volume.Therefore, it is more likely that the GIS volume lost during the Eemian was closer to 50-60% (or 3.7-4.4m sea level.)Previous estimates of ice loss from Greenland during the Eemian are generally consistent with the results presented here, even though other mass balance schemes were applied.In at least three cases, the annual PDD approach was applied and a wide range of climatic forcing was tested: Cuffey and Marshall (2000) provided a plausible range of 4.0-5.5 m sea level and Huybrechts (2002) provided a similar maximum contribution.Tarasov and Peltier (2003) estimated a range of 2.7-4.5 m sea level using a similar approach to that of Cuffey and Marshall (2000), but they included additional constraints based on δ 18 O tracers at the ice core locations.In a more recent study, using an ice sheet model and climate time slices from a Global Circulation Model, Otto-Bliesner et al. ( 2006) estimated a range of 2.0-3.5 m sea level.All estimates fall within a similar range, and highlight the fact that without more accurate information about the amount of Eemian warming at the margin of the ice sheet, it may be difficult to provide a narrower estimate.
For other locations on the ice sheet, the constrained model ensemble produced various results.Most simulations considered valid by our three chosen constraints maintained ice at DYE-3 throughout the Eemian Interglacial.The range of elevation loss was quite large, however, ranging from almost no change to a 1300 m decrease for the most extreme cases (or total loss of DYE-3 ice in 6 cases).If the DYE-3 elevation did decrease by 1300 m and, correspondingly, the GRIP elevation decreased by the maximum allowed amount of 400 m, the relative elevation change between the two locations would be near 1000 m.This would imply an additional 6 • C of Eemian warming at DYE-3, which would not be consistent with δ 18 O record here.Nonetheless, it is worth mentioning that in our model it is possible to melt a significant portion of the GIS while maintaining reducedthickness ice at DYE-3.
The maximum limit of volume loss during the Eemian is mainly determined by the constraint applied to the elevation loss at GRIP during this period.We chose to apply a realistic constraint of 400 m here.However, the uncertainty in this value is large.To demonstrate how the choice of this constraint could affect the results, we reanalyzed the results using different limits, ranging from 0-1500 m of elevation loss at GRIP.As can be seen in Fig. 10, as the GRIP elevation was allowed to decrease further, the maximum Eemian contribution to sea level rise increased non-linearly.For a realistic range of 300-600 m of elevation loss, the uncertainty in the results is only 0.4 m in either direction from our original choice of 400 m, which is rather small compared to the overall range of uncertainty.This range is also consistent with a 5 • C temperature anomaly at GRIP, while an elevation drop of much more than 600 m would imply a much higher temperature anomaly.

Conclusions
Simulations of the evolution of the Greenland Ice Sheet have been performed for the last 250 ka BP using a coupled regional climate-ice sheet model that is physically based and applicable for a wide range of climatic and topographic change scenarios.Several key model parameters were perturbed to produce an ensemble of model versions, which were then constrained using information about the ice sheet for the Eemian Interglacial and the present day.The stability of the GIS is predominantly determined by the surface mass balance, and particularly, by the sensitivity of the melt equation to external climatic changes.Other parameters, such as the geothermal heat flux and sliding coefficient, play appreciable but less significant roles in determining the past evolution and present geometry of the simulated GIS (as found by others, e.g., Ritz et al., 1997).We have shown that even modest variations of a single model parameter lead to large uncertainties in the simulated volume of the Eemian GIS.
Combined information about the present day and the Eemian helps to reduce the range of valid parameter combinations and model simulations considerably.The modern and paleo constraints produce consistent and, in some respects, mutually complimentary, limitations of the model parameters.Using both modern and paleo constraints together reduces the number of valid model versions to 22% of the initial subset (20 out of 90).Without the use of additional information about the range of Eemian temperature changes, the estimate of the contribution of the GIS to Eemian sea level rise is rather uncertain.From our study, an acceptable range of Eemian melt is 5-60% mass loss (0.4-4.4 m sea level).However, when additional constraints on the boundary warming are considered, the likely range narrows to 25-60% mass loss (1.9-4.4 m sea level).The highest values in this range are the most likely (50-60% mass loss; 3.7-4.4m sea level), given the estimate of up to 5 • C warming at GRIP.These numbers are, of course, only rough estimates, given the poor representation of the ice sheet at the margins.But assuming that the Eemian high stand was 6-8 m above present, this estimate for the Greenland melt still requires a considerable contribution to sea level rise from the Antarctic Ice Sheet.
None of our model versions produce a sufficiently realistic present-day GIS in terms of volume and spatial extent to choose a best set of parameter values.The most realistic simulations of the modern GIS (i.e., with less volume and surface area) were obtained in the experiments that produced completely unrealistic simulations of the Eemian GIS (almost complete melting).At the same time, many model versions that satisfied the Eemian paleo constraint were also consistent with the modern constraint on the GIS mass balance partition.This would support the idea that in view of the imperfectness of existing ice sheet models, the latter constraint (criterion) is more appropriate for the selection of model parameters for past and future simulations of the GIS -at least, on the millennial time scale.
Finally, in spite of the limitations of the model used and remaining uncertainties, our work indicates that using past and present constraints together can help to rule out both too sensitive and too insensitive model versions.
albedo is calculated as in the original formulation to provide an estimate for the planetary albedo used in the energymoisture balance equations of REMBO (since the albedo of land would play a role here).
Furthermore, we modified the critical snow depth in Eq. (A1) to depend on the type of vegetation that would be present.To do so, we calculated the available positive degree days based on temperature and converted this to a land type.
We tested the effect of changing the albedo parameterization via the critical snow depth and the minimum ground albedo, along with variations to the ITM parameter c.These three parameters are interrelated and changing one can offset the effects of the other.However, to simulate realistic glacial cycles (in that the ice sheet regrows completely during cold periods), it was necessary to increase the minimum ground albedo (as described above) to that of an ice-covered surface (0.4).In the annual average, the model produces comparable results to those presented by Robinson et al. (2010).

Fig. 1 .
Fig. 1.Examples of forcing used for simulations of the last two glacial cycles: (a) maximum solar insolation at 65 • N; (b) equivalent CO 2 concentration; (c) global sea level anomaly from SPECMAP; (d) Summer (June-July-August averaged) temperature anomaly prescribed at the boundaries of Greenland.

Fig. 3 .
Fig. 3. Temporal evolution of (a) the volume and (b) the surface area of the GIS, as well as (c) the precipitation-weighted annual temperature anomaly at GRIP and (d) the annual snowfall at GRIP, as simulated by the ensemble of model versions.Lighter lines correspond to "invalid" simulations.Each color band is determined by the melt parameter c, as explained in the caption of Fig. 6.

Fig. 4 .
Fig. 4. GIS distribution for different values of the melt parameter c, as simulated for (a-e) the present and (f-j) the Eemian minimum volume.The black diamonds show the locations of ice core drill sites (from South to North: DYE-3, GRIP, NGRIP, NEEM and Camp Century).
Fig. 5. (a-e) Precipitation-weighted annual temperature anomaly and (f-j) the change in annual precipitation relative to present day, at the Eemian minimum extent for each of the five simulations in Fig. 4.

Fig. 7 .
Fig. 7. Ensemble results for the three constraints versus the most influential parameter in each case: (a) diagnosed SMB partition for the present-day GIS, compared to results from regional climate models (PolarMM5: Box et al., 2006; MAR: Fettweis et al., 2007; RACMO2/GR: Ettema et al., 2009); (b) simulated present-day GRIP elevation versus the prescribed geothermal heat flux at the base of the ice sheet;(c) maximum reduction in the modeled GRIP elevation for the Eemian Interglacial relative to the modeled present-day elevation.Darker circles and lighter crosses correspond to "valid" and "invalid" simulations, respectively, and the grey lines show the outer limits used for each constraint.
Fig. 9. (a)Maximum percent of Eemian ice loss versus the maximum precipitation-weighted temperature anomaly experienced at GRIP during the Eemian.The percentage of loss is also converted into sea-level equivalent, assuming 100% correspondence to the presentday volume of the GIS.Darker circles and lighter crosses correspond to "valid" and "invalid" simulations, respectively.The valid range corresponding to each choice of the paleo factor, f p , is indicated by the dark lines.The temporal evolution during the Eemian of (b) the Greenland contribution to sea level and (c) the precipitation-weighted annual temperature at GRIP is shown for valid simulations.

Fig. 10 .
Fig.10.Sensitivity of the range of Eemian volume loss to the allowed Eemian GRIP elevation change.Two circles are shown for each constraint level, representing the minimum and maximum simulations.The color of the circle corresponds to the choice of parameter c of the simulation that determines the extreme (as described in Fig.6).No simulation had less than 150 m of elevation reduction during the Eemian.The grey band shows the likely range of acceptable constraint values and the vertical black line is the constraint value chosen as default in this study (400 m).

Table 1 .
All parameter values used to generate the ensemble of model versions.
large-scale precipitation).In addition, we considered possible uncertainties in temperature anomalies simulated by CLIMBER-2.Namely, we changed the magnitude of warming around Greenland during the Eemian Interglacial.These areas of uncertainty are discussed in detail below, and Table Robinson et al.:Greenland ice sheet model parameters constrained using the Eemian www.clim-past.net/7/381/2011/Clim.Past, 7, 381-396, 2011 A.