Investigating the sensitivity of numerical model simulations of the modern state of the Greenland ice-sheet and its future response to climate change

Ice thickness and bedrock topography are essential boundary conditions for numerical modelling of the evolution of the Greenland ice-sheet (GrIS). The datasets currently in use by the majority of GrIS modelling studies are over two decades old and based on data collected from the 1970s and 80s. We use a newer, high-resolution Digital Elevation Model of the GrIS and new temperature and precipitation forcings to drive the Glimmer ice-sheet model offline under steady state, present day climatic conditions. Comparisons are made of ice-sheet geometry between these new datasets and older ones used in the EISMINT-3 exercise. We find that changing to the newer bedrock and ice thickness makes the greatest difference to Greenland ice volume and ice surface extent. When all boundary conditions and forcings are simultaneously changed to the newer datasets the ice-sheet is 33 larger in volume compared with observation and 17 larger than that modelled by EISMINT-3. We performed a tuning exercise to improve the modelled present day ice-sheet. Several solutions were chosen in order to represent improvement in different aspects of the GrIS geometry: ice thickness, ice volume and ice surface extent. We applied these new parameter sets for Glimmer to several future climate scenarios where atmospheric CO2 concentration was elevated to 400, 560 and 1120 ppmv (compared with 280 ppmv in the control) using a fully coupled General Circulation Model. Collapse of the ice-sheet was found to occur between 400 and 560 ppmv, a threshold substantially lower than previously modelled using the standard EISMINT-3 setup. This work highlights the need to assess carefully boundary conditions and forcings required by ice-sheet models, particularly in terms of the abstractions required for large-scale ice-sheet models, and the implications that these can have on predictions of ice-sheet geometry under past and future climate scenarios. Â© 2010 Author(s).


Introduction
Complete melting of the Greenland ice-sheet (GrIS) would raise sea level by as much as 7.3 m (Bamber et al., 2001), and could be associated with other major climatic effects such as changes in the thermohaline circulation and oceanic heat transport due to enhanced freshwater fluxes (Fichefet et al., 2003). Estimates of the GrIS's contribution to sea level change during the period 1993 to 2003 range between +0.14 to +0.28 mm yr −1 (IPCC, 2007), although recent estimates suggest as much as +0.75 mm yr −1 for ) linked with significant recent increases in GrIS melt, runoff and mass loss Rignot et al., 2008). Recent model projections suggest that the GrIS could be eliminated within a few millennia for global warming between 1.9 to 4.6 • C relative to pre-industrial temperatures (Gregory and Huybrechts, 2006). These projections are based on a numerical model which does not include a representation of fastflowing outlet glaciers. These glaciers have been observed to undergo dynamic changes in recent years, resulting in faster ice flow and consequent ice loss (Howat et al., 2007;Joughin et al., 2004;Luckman et al., 2006;Rignot et al., 2008;Rignot and Kanagaratnam, 2006), meaning that the model probably underestimates the rate of mass-loss from the GrIS.
The majority of recent modelling studies of the GrIS use the data assembled for the EISMINT (European Ice-sheet Table 1. List of default parameters and physical constants used in the model. Those highlighted in bold are varied in the tuning experiments (for a complete set see Rutt et al., 2009).

Symbol Value
Units Description ρ i 910 kg m −3 Density of ice g 9.81 m s −2 Acceleration due to gravity a 1.733 × 10 3 Pa −3 s −1 Material constant for T * ≥ 263 K a 3.613 × 10 −13 Pa −3 s −1 Material constant for T * < 263 K Q 139 × 10 3 J mol −1 Activation energy for creep for T * ≥ 263 K Q 60 × 10 3 J mol −1 Activation energy for creep for T * < 263 K R 8.314 J mol −1 K −1 Universal gas constant α i 8 mm water d −1 • C −1 Positive degree day factor of ice α s 3 mm water d −1 • C −1 Positive degree day factor of snow L G −6.227 • C km −1 Atmospheric temperature lapse rate n 3 -Flow law exponent f 3 -Flow enhancement factor G −0.05 W m −2 Uniform geothermal heat flux Modelling INiTiative) model intercomparison project as a present day representation of the GrIS. Because the description of the data is included in the report from the 3rd EIS-MINT workshop (Huybrechts, 1997), we refer to them here as the EISMINT-3 data. The data consist of a Digital Elevation Model (DEM) of ice thickness and bedrock elevation, and parameterised temperature and precipitation fields, onto which climate anomalies are typically superimposed (e.g. Driesschaert et al., 2007;Greve, 2000;Huybrechts and de Wolde, 1999;Ridley et al., 2005;Lunt et al., 2008Lunt et al., , 2009) . The high-resolution bedrock and ice thickness used in EISMINT-3 are nearly two decades old and are based on data collated during the 1970s and 1980s. More recent and accurate datasets for the boundary conditions of bedrock topography and ice thickness (Bamber et al., 2001) as well as temperature  and precipitation (ECMWF, 2006) forcings are now available. Differences in these datasets could have considerable impacts on the modelled evolution of the GrIS and hence the resulting ice-sheet volume and geometry, for simulations of past, modern and future climates.
In this paper, we use the Glimmer ice-sheet model (Rutt et al., 2009) to investigate and compare the impact on the modelled steady-state ice-sheet of two sets of boundary conditions: those used in the EISMINT-3 exercise and the more recent and up-to-date datasets. Furthermore, we perform a tuning exercise with respect to the most recent datasets in order to determine the values of various ice-sheet model parameters which give the best fit between modelled and observed geometry for present day conditions. Finally, we use the results from the tuning exercise to assess the impact of different parameter combinations on future warming scenarios with atmospheric CO 2 held at 400 ppmv, 560 ppmv and 1120 ppmv (compared with 280 ppmv in the control) where the ice-sheet model is driven offline using output from a fully-coupled General Circulation Model (GCM). Most recent sensitivity studies have only used one set of ice-sheet model parameters (e.g. ablation coefficients) for simulations of future ice-sheet evolution (e.g. Alley et al., 2005;Driesschaert et al., 2007;Mikolajewicz et al., 2007;Ridley et al., 2005). Our results highlight the need to use a range of ice model parameter sets in order to assess their impact on future ice-sheet climate scenarios.

Model description
We use the 3-D thermomechanical ice-sheet model Glimmer version 1.0.4 (Rutt et al., 2009). Although not the most recent version of the model, we use this version for consistency with our previous work (e.g. Lunt et al., 2008Lunt et al., , 2009. The core of the model is based on the ice-sheet model described by Payne (1999). All physical constants and parameters discussed in this section are given in Table 1. Here we describe the parts of the model which pertain to the model parameters which we tune in the subsequent sections. A full description of the model can be found in Rutt et al. (2009).
The ice thickness (H ) evolution is driven by the mass conservation equation where u is the horizontal velocity andū is the horizontal velocity averaged over the ice thickness, B is the surface mass balance rate and S is the basal melt rate. Equation (1)  The ice dynamics are represented with the widely-used shallow-ice approximation, which assumes ice deformation occurs as shear strain only, so that u(z) = u(b)−2(ρ i g) n |∇s| n−1 ∇s z b A T * s−z n dz , (2) where s is the ice-sheet surface altitude, b is the bedrock altitude, g is the acceleration due to gravity, ρ i is the ice density, u(b) is the sliding velocity, x and y the horizontal coordinates and z the vertical coordinate, positive upward. Equation (2) implicitly uses the non-linear viscous flow law (Glen's flow law) to relate deformation rate and stress. The two parameters are the exponent, n, and the ice flow law parameter, A(T * ), which is an empirical parameter where T * is the absolute temperature corrected for the dependence of the melting point on pressure. This parameter follows the Arrhenius relationship where a is a temperature-independent material constant, Q is the activation energy and R is the universal gas constant. In Eq.
(3) f is the flow enhancement factor, a tuneable factor, which can be used to change the flow law parameter, and, hence, change the ice flow velocity. The flow enhancement factor accounts for ice impurities and development of anisotropic ice fabrics, effects not represented by separate parameters in the model. The model is formulated on a Cartesian x − y grid, and takes as input the surface mass balance and mean nearsurface air temperature at each time step. In the present work, the ice dynamics time step is one year. To simulate the surface mass balance, we use the Positive Degree Day (PDD) scheme described by Reeh (1991). The basis of the PDD method is the assumption that the melt, m, that takes place at the surface of the ice-sheet is proportional to the timeintegrated temperature above freezing point, known as the positive degree day where T (t) is the near-surface air temperature and α is the PDD factor. Two PDD factors which describe the rate of melting are used, one each for snow (α s ) and ice (α i ), to take account of the different albedo and density of these materials. The integral in Eq. (4) is calculated on the assumption of a sinusoidal annual variation in temperature, and takes as input the mean annual temperature and half-range. Diurnal and other variability is taken into account using a stochastic approach. This variability is assumed to have a normal distribution with a standard deviation of 5 • C. The use of PDD mass balance models is well-established in coupled atmosphereice-sheet modelling studies of both paleoclimate (e.g. De-Conto and Pollard, 2003;Lunt et al., 2008) and future climate (e.g. Ridley et al., 2005;Mikolajewicz et al., 2007).
All precipitation is assumed to be potentially available for accumulation within the Glimmer annual PDD scheme. The following possibilities are taken into account when considering the total annual ablation. Melting snow is allowed to refreeze to become superimposed ice up to a fraction, w, of the original snow depth. When the ability of the snow to hold meltwater is exceeded but the potential snow ablation is less than the total amount of precipitation (amount of snow available), run-off can occur. If the potential snow ablation is greater than precipitation, snow will melt first, and then ice, such that the total ablation is equivalent to the sum of snow melt (total precipitation minus the amount of meltwater held in refreezing) and the sum of ice melt (calculated by deducting from the total number of degree days from the number of degree days need to melt all snow fall and converted to ice melt). Therefore, the net annual mass balance is the difference between the total annual precipitation and the total annual ablation. Glimmer also includes a representation of the isostatic response of the lithosphere, which is assumed to behave elastically, based on the model of Lambeck and Nakiboglu (1980). The timescale for this response is 3000 years. In all model runs described below, the isostasy model is initialised on the assumption that the present day bedrock depression is in equilibrium with the ice-sheet load. Although this assumption may not be entirely valid, any rates of change will not have a significant influence for present day geometry (Huybrechts and de Wolde, 1999).
Geothermal heat flux (G) can be supplied to the model as a constant or a spatially varying field (both of which are explored in Sect. 5.2), and a thermal bedrock model (Ritz, 1987) takes the thermal evolution of the uppermost bedrock layer into account where initial conditions for the temperature field are found by applying the geothermal heat flux to an initial surface temperature.
The forcing data (temperature and precipitation) are transformed onto the ice model grid using bilinear interpolation. In the case of the near-surface air temperature field (T a ), a vertical lapse-rate correction is used to take account of the difference between the high-resolution (20 km in this case) surface topography seen within Glimmer (s G ), and that represented by the forcing data (s) (in this case a latitude longitude 1 • by 1 • grid or approximately 111 km resolution), such that 400 E. J. Stone et al.: Greenland ice-sheet and its future response to climate change 3 The datasets

EISMINT-3 intercomparison experimental design
In order to evaluate the consistency in predictions between different ice-sheet models, the EISMINT validation exercise was set up (Huybrechts and Payne, 1996). EISMINT-3 (Huybrechts, 1997) was the final part of this exercise which involved modelling changes in ice mass given a climate scenario for a number of different ice-sheet models with prescribed parameters and climate forcings (van der Veen and Payne, 2004). This included the evolution of GrIS mass changes under (a) steady-state present climate conditions, (b) a transient climate such as the last climatic cycle based on GRIP ice core data and (c) future greenhouse warming. By modelling present day steady-state conditions, it is possible to test the validity of the reconstructions that the models produce, by comparing the model predictions with observations of the present day ice-sheet. In the EISMINT-3 standard, the initial condition of bedrock and surface elevation was compiled by Letreguilly et al. (1991) on a 20 km Cartesian grid. The precipitation forcing is from Ohmura and Reeh (1991) and the temperature forcing is given by the following simplified parameterisations (Huybrechts and de Wolde, 1999;Ritz et al. 1997) which were themselves based on observed surface temperature data (Ohmura, 1987) where H surf is the surface elevation (m), is the geographical latitude (in degrees and positive), T ann is the mean annual temperature, T s is the summer temperature (both in • C), and L a = −7.992, L s = −6.277 are annual and summer atmospheric lapse rates respectively (in • C km −1 ).

Recent boundary conditions/forcings
New and more accurate bedrock and surface elevation datasets are now available with significant differences in ice volume (∼ 4% increase) and ice thickness (factor of 10) around the margins compared with the Letreguilly dataset (Bamber et al., 2001). This new dataset utilises improvements in the boundary conditions of surface elevation. Ice thicknesses were derived from combining data collected in the 1970s with new data obtained from an ice penetrating radar system from 1993 to 1999. The bedrock topography was subsequently derived from a DEM of the ice-sheet and surrounding rocky outcrops. As such the DEM is produced from a combination of satellite remote sensing and cartographic datasets. In contrast, the Letreguilly dataset is based on cartographic maps for ice-free regions and radio echoing sounding for determination of ice thickness. No satellitederived products were used. The Bamber dataset has the advantage of significantly more sources of accurate data and better coverage. The Bamber dataset is on a 5 km resolution grid; for the purposes of the present work, it was interpolated onto a 20 km resolution grid, generated by pointwise averaging on the same projection. Henceforth, we will refer to the EISMINT-3 bedrock and ice thickness dataset as the "Letreguilly" dataset and the more recent dataset as the "Bamber" dataset. The precipitation data used in EISMINT-3 (Ohmura and Reeh, 1991) is based purely on precipitation measurements from meteorological stations (35) and pits and cores in the interior of the ice-sheet. Not only is this based on a small number of data locations but the accuracy of measurements is also a matter of contention. Catch efficiency, particularly for solid precipitation, by gauges is somewhat reduced by turbulent winds along with the potential for snow to be blown out of gauges (Yang, 1999). Measurement error may reach 100% during the winter months, when accumulation is most important for mass balance (Serreze et al., 2005). We make use of precipitation data derived from ERA-40 reanalysis from 1979(ECMWF, 2006) on a regular latitude-longitude 1 • by 1 • resolution grid. ERA-40 reanalysis is produced using a data assimilation technique which consists of a number of analysis steps (Uppala et al., 2005). Background information is produced from a short-range forecast and combined with observations for this same period of the forecast to produce an "analysis". Statistically-based estimates of errors are used for the synthesis of background forecast and observation. Each forecast is initialised from the most recent previous analysis step. Observations do not consist of all meteorological variables but the analysis is complete in terms of the variables chosen. As such, variables can be produced from analysis (e.g. temperature) while others are purely based on forecast and are, therefore, not constrained by observations (Uppala et al., 2005). In ERA-40, precipitation is one such variable produced by the forecast rather than by the analysis in the ECMWF model. However, it has been shown to be reasonable for Greenland (Serreze et al., 2005). Validation against Danish Meteorological Institute (DMI) coastal stations results in a 36% mean excess for ERA-40 (Hanna and Valdes, 2001), although the inaccuracies in gauge measurements mean that this should be treated with some caution. In terms of other reanalysis products available, comparison studies have shown ERA-40 to be superior to NCEP/NCAR datasets in terms of smaller biases, ability to capture largescale patterns of precipitation and its depiction of interannual variability, deeming ERA-40 a more suitable choice (Bromwich et al., 1998;Hanna et al., 2006;Serreze et al., 2005;Serreze and Hurst, 2000).
The near-surface air temperature forcing used in the EISMINT-3 exercise is based on a parameterisation of surface temperature compiled by Ohmura (1987), which has a latitudinal and altitude dependency (see Eqs. 6 and 7). Two lapse rate values are used: the mean annual lapse rate and a summer lapse rate. Currently, lapse rate in Glimmer is not temporally or regionally varying so the summer Table 2. Summary of sensitivity experiments performed by changing each boundary condition/forcing individually from that used in the EISMINT-3 exercise to the more recent datasets. Also included is an experiment where all boundary conditions/forcings are updated together. For details on the EISMINT-3 datasets see Sect. 3.1.
Bedrock and ice thickness Precipitation Temperature Bamber et al. (2001) EISMINT-3 EISMINT-3 EISMINT-3 ERA-40 EISMINT-3 EISMINT-3 EISMINT-3 Hanna et al. (2005) Bamber et al. (2001 ERA-40 Hanna et al. (2005) lapse rate is used since this is when the ablation process is strongest. The parameterisations were constructed to fit data from 49 meteorological stations. Instead, we use, to be consistent with precipitation, near-surface air temperature data derived from ERA-40 "corrected" 2-m near-surface air temperatures . The temperatures were corrected based on their derived surface lapse rates and differences between the ECMWF orography and a DEM derived from the Ekholm (1996) Greenland grid . Reasonable agreement exists between these model-derived temperatures and observations at the DMI station locations and GC-Net stations . We use bilinear interpolation to transform the high-resolution dataset from its Cartesian 5 km resolution grid onto a 1 • by 1 • latitude longitude grid. Since the dataset only covers the regions where there is ice, the temperature parameterisation used in EISMINT-3 temperature is used in the ice-free regions of Greenland in conjunction with the Ekholm orography. This means that the sensitivity to temperature is specifically a sensitivity to the near-surface air temperature of the ice-sheet and not the ice-free regions.

Sensitivity to boundary conditions and forcings
In order to test the sensitivity of the ice-sheet model to the various forcing inputs and boundary conditions, we performed a set of steady-state experiments shown in Table 2, initialised from the present day geometry of the ice-sheet. The model is run for 50 000 years in order to reach equilibrium. The configuration of the ice-sheet model is kept at that of EISMINT-3 with standard parameter values as shown in Table 1. For each simulation in the set, one forcing/boundary condition is changed to the most recent dataset, keeping all others at that used in EISMINT-3. An additional experiment is performed where all the forcings and boundary conditions are changed to the most recent. Figure 1 shows the evolution of ice volume and ice surface extent with time for EISMINT-3 and the four sensitivity experiments.  Bamber et al. (2001) and Letreguilly et al. (1991) are also shown for comparison.

Bedrock and ice thickness
The quality of the bedrock topography is important in icesheet models since it largely determines the ice thickness at regional scales. This is because topography influences where the build up of snow and ice can occur and, therefore, is a major control on the threshold of ice-sheet initiation. Furthermore, topography influences the convergence and divergence 1 Figure 2. a) The ratio of the difference between ice thickness of the Bamber dataset 2 and ice thickness of the Letreguilly dataset (z bamber -z letreguilly /z letreguilly ) expressed as a 3 percentage. The regions of largest relative difference occur around the margins with 4 good agreement between the datasets in the ice-sheet interior. b) The ratio of the 5 difference in initial bedrock topography of Bamber dataset and the topography of 6 Letreguilly expressed as a percentage. Again the largest differences occur around the 7 margins of Greenland and also in the central region where the bedrock is below sea 8 level. c) The ratio of the difference in relaxed bedrock topography after the removal of 9 ice and isostatic equilibrium has been reached expressed as a percentage. The 10 resultant orography shows the relative difference around the margins of up to 500%, 11 with the Bamber orography significantly higher. z bamber − z letreguilly /z letreguilly expressed as a percentage. The regions of largest relative difference occur around the margins with good agreement between the datasets in the ice-sheet interior. (b) The ratio of the difference in initial bedrock topography of the Bamber dataset and the topography of Letreguilly expressed as a percentage. Again the largest differences occur around the margins of Greenland and also in the central region where the bedrock is below sea level. (c) The ratio of the difference in relaxed bedrock topography after the removal of ice and isostatic equilibrium has been reached expressed as a percentage. The resultant orography shows the relative difference around the margins of up to 500%, with the Bamber orography significantly higher.
of ice flow such that flow into lowland basins and valleys from surrounding higher relief regions will result in faster build up of ice compared with flow from an isolated upland region into a lower basin (Payne and Sugden, 1990). As a result, the topography influences the stress, velocity and thermal regimes of the ice-sheet (van der Veen and Payne, 2004).
At the outset there are differences in ice thickness and bedrock topography between the two bedrock and icethickness datasets (see Fig. 2a and b). The bedrock topography around the margins is consistently higher for the Bamber dataset compared with the Letreguilly dataset, with the ice thickness difference up to a factor of 10 to 20 thicker. When simulated to steady-state, the Bamber bedrock and ice thickness datasets results in significantly (13.7%) greater ice volume and 11.5% larger ice surface extent compared with the Letreguilly dataset. Ice extends further to the northern and western margins of Greenland with a higher central dome. The initial higher elevation of the ice-free bedrock of the Bamber dataset provides favourable conditions for ice growth where temperatures are low enough for mass balance to become positive. In these regions ice velocities are low compared with other marginal regions, allowing the ice-sheet to build-up with minimal ice loss. The basal temperatures are also lower than when the Letreguilly dataset is used, resulting in marginally lower velocities for ice flow. This arises because the ice in the Bamber dataset is thicker at the beginning of the simulation. The increase in ice volume and surface extent, however, can be attributed predominately to a stronger temperature-elevation feedback mechanism for the Bamber dataset.

Precipitation
Changing the precipitation forcing, from that of Ohmura and Reeh (as in EISMINT-3) to ERA-40, results in an increase in equilibrium ice-sheet surface extent of 2.1%. However, there is almost no effect on the ice-sheet volume. This can be explained by the fact that all precipitation that falls is assumed to fall as snow in the annual PDD scheme. Since the temperature forcing has no effect on the amount of snow, it is the quantity and distribution of precipitation that results in the difference in ice surface extent. Figure 3 shows that the annual precipitation is up to two times greater on the eastern and western margins of Greenland for ERA-40 compared with Ohmura and Reeh (1991). The accumulation rate is greatest in south-east Greenland for both precipitation datasets but with the high values extending further north along the eastern margin for ERA-40. The extra precipitation falling over the western and eastern margins coupled with a positive temperature-elevation feedback results in growth and extension of the ice-sheet into previously ice-free regions. However, the precipitation falling over central and north Greenland is three times less for ERA-40, resulting in less accumulation in the interior and lower maximum altitude of the ice sheet. These opposing effects result in similar ice-sheet volumes. However, Hanna et al. (2006) show that ERA-40 is ∼50% too "dry" in the central northern parts of Greenland, as validated using ice-core data. Furthermore, it seems increasingly likely that both the Ohmura and Reeh (1991) and ERA-40 precipitation datasets underestimate precipitation and accumulation in south-east Greenland, where recent regional climate model results suggest much higher than previously observed precipitation rates Burgess et al., 2010). 50 1 Figure 3. Change in precipitation over Greenland between EISMINT-3 (Ohmura and 2 Reeh, 1991) and ERA-40 re-analysis (Uppala et al. 2005) expressed as a ratio of 3 EISMINT-3:ERA-40. Annual near-surface temperature (in °C) contours also shown 4 (Ohmura, 1987).  EISMINT-3 (Ohmura and Reeh, 1991) and ERA-40 re-analysis (Uppala et al., 2005) expressed as a ratio of EISMINT-3:ERA-40. Annual near-surface temperature (in • C) contours are also shown (Ohmura, 1987).

Temperature
Changing the temperature forcing to the modified Hanna dataset results in a similar ice volume (1.6% larger) compared with EISMINT-3 and an almost identical ice-sheet extent. Figures 4 and 5 show the temperature distribution and the surface mass balance respectively at the beginning and end of the experiments for EISMINT-3 temperature and the Hanna modified temperature datasets. As expected, at the beginning of the simulation temperatures around the margins of the GrIS are similar (same datasets) but the Hanna ERA-40 corrected temperatures over the ice-sheet are several degrees lower ( Fig. 4a,b) . By the end of the simulations, temperatures over much of Greenland have become lower as a result of the positive temperature-elevation feedback (Fig. 4c,d) resulting in an increase in positive net mass balance in southern Greenland (see Fig. 5c,d). However, the regions around the margins remain ice-free as a result of continued ablation with a net negative mass balance. The model is particularly sensitive to the temperature forcing around the margins of the ice-sheet, where temperatures are at zero or above and so close to ablation as opposed to those in the interior where the primary mass balance change is from accumulation . It is, therefore, important that marginal temperatures close to where the net mass balance becomes negative are resolved accurately in order to model the ablation process and the resulting geometry of the GrIS.  Table 3 summarises the results of changing bedrock and ice thickness, precipitation and temperature independently from EISMINT-3 to the newer datasets. Bedrock and ice thickness result in the largest ice volume and ice surface extent change while changing precipitation and temperature have a significantly smaller effect on the ice volume.

Summary
Updating all the boundary conditions and forcings together results in a modelled GrIS ice volume 33% larger than observed (Bamber et al., 2001) and 17% larger than EISMINT-3. The system is effectively linear since adding together the difference between the EISMINT-3 case and the individual response of the ice-sheet to each forcing/boundary condition results in a modelled GrIS very similar to when all forcings/boundary condition are varied together (see Fig. 1). This is the case for ice volume (1.7% smaller) and ice surface extent (0.1% smaller).
These results show that when using alternative boundary conditions and forcings Glimmer gives a poorer representation of the modern ice-sheet compared with observations. It is likely that some of the internal ice-sheet model parameters were tuned to work with the boundary conditions used  in EISMINT-3. In order to produce a reasonable best fit between modelled and observed geometry we tune a number of ice model parameters to work with the new datasets.

Tuning methodology
Several parameters in large-scale ice-sheet modelling are still poorly constrained, resulting in highly variable ice-sheet volume and extent depending on the values prescribed in the model (Ritz et al., 1997). This necessitates tuning of the icesheet model with the recent datasets in order to determine the optimal ice-sheet for steady-state conditions (i.e. closest geometry to reality). Previous work (e.g. Ritz et al., 1997) has looked at the sensitivity of ice-sheet volume and extent to a number of parameters, including flow enhancement factor (f ) in the flow law (see Eq. 3), the sliding coefficient, the geothermal heat flux (G) and the coefficients (PDD factors) of the ablation parameterisation for ice (α i ) and snow (α s ) (see Eq. 4). In addition, Hebeler et al. (2008a) also looked at the effect on ice volume and extent of the Fennoscandian ice-sheet during the Last Glacial Maximum from uncertainty in model parameters (e.g. lapse rate in addition to those mentioned above) and climate forcing by performing a parametric uncertainty analysis using Glimmer, and found a variation of 65% in equilibrium ice sheet extent due to uncertainty in the parameters used in the ice sheet model and up to 6.6% due to uncertainty in topographic input.
The most common methodology in glaciological modelling sensitivity studies is to vary one parameter at a time within a prescribed range while holding all others constant (e.g. van de Wal and Oerlemans, 1994;Essery and Etchevers, 2004;Fabre et al., 1995;Huybrechts and de Wolde, 1999;Pattyn, 2003;Ritz et al., 1997). We build on the methodology used in this previous work by using the statistical method of Latin-Hypercube Sampling (LHS) (an efficient variant of the Monte Carlo approach) which generates a distribution of plausible parameter sets within a prescribed set of ranges (McKay et al., 1979). It uses a stratified-random procedure where values are sampled from the prescribed distribution of each variable. The cumulative distribution of each variable is divided into N equiprobable intervals and a value selected randomly from each interval. The N values obtained for each variable are paired randomly with the other variables. The method assumes that the variables are independent of one another (which is the case here) and ensures a full coverage of the range of each variable. LHS has been used in a number of applied scientific disciplines including analysing uncertainty in vegetation dynamics (Wramneby et al., 2008), rainfall models for climate assessment (Murphy et al., 2006) and climate/ocean models (Edwards and Marsh, 2005;Schneider von Deimling et al., 2006). However, it has yet to be used in large-scale ice-sheet modelling. The advantage of this methodology is that it is an efficient method to test the response of the ice-sheet to many different combinations of parameters by ensuring sufficient coverage of the parameter space without having to test all possible model combinations (which would be extremely computationally expensive). In this way, by varying more than one parameter at a time (as for any multivariate sampling method) it also allows the influence of each parameter on the outcome of the model simulations to be assessed while taking interactions with other parameters into account.
We investigate not only the result of uncertainty in the following parameters, but also which combination gives the optimal fit to the present day GrIS. The geometry of the GrIS is controlled by the flow of ice from the ice divide in the interior towards the coastal regions due to internal deformation where at relatively low altitudes, typically < ∼ 2000 m, ice mass is lost by melting according to the PDD scheme. Ice mass can also be lost by basal melt and/or the process of basal sliding which can increase the flow of ice to regions of ablation at the edge of the ice-sheet. Since basal sliding is not included in these simulations, this process will not be considered but the likely impact of this missing process is highlighted in the discussion section. We choose the following parameters to tune since they fundamentally affect the processes described in Sect. 2. Firstly, the flow rate of ice can be tuned with the flow enhancement factor, f (see Eq. 3), to simulate ice flow reasonably accurately. Secondly, the surface mass balance can be tuned using the PDD factors and vertical lapse rate. The melting of ice at low altitudes is determined by ablation, which in this study is calculated according to the annual PDD scheme. Since this uses an empirical relationship, we choose to vary the PDD factors for ice (α i ) and snow (α s ) within the ranges obtained through measurement studies (see below), and, therefore, influence the amount of melting that can occur in the ablation zones. These parameters will not, however, alter the position of these zones. This instead can be achieved by varying the vertical atmospheric lapse rate (L G ), which can influence the regions where ablation has the potential to occur. Thirdly, ice loss by basal melt without sliding can be achieved by varying the geothermal heat flux (G), which can raise the basal ice layer temperature to its pressure melting point.
LHS requires a maximum and minimum bound for each tuneable parameter to be defined. Here we discuss the bounds we have selected for each value, shown in Table 4.
The range for the flow enhancement factor for this study is between 1 and 5. According to Dahl-Jensen and Gundestrup (1987), borehole measurements from the Dye-3 ice core give a mean enhancement factor of around 3 with a maximum value of 4.5 and a minimum value of around 1 for ice deposited during the Wisconsin. This is the range used by Ritz et al. (1997) and Hebeler et al. (2008a) for their sensitivity studies. Values within this range have also been used in other work (e.g. Fabre et al., 1995;Greve and Hutter, 1995;Huybrechts et al., 1991;Letreguilly et al., 1991). Table 4. List of five parameters varied according to the ranges determined from the literature. The parameters α i , α s , G and f are similar to those used in Ritz et al. (1997).

Parameter
Minimum value Maximum value Positive degree day factor for snow, 3 5 α s (mm d −1 • C −1 ) Positive degree day factor for ice, The global average geothermal heat flux (oceans and continents) is estimated at 87 mW m −2 (Banks, 2008). Since it is difficult to measure geothermal heat flux beneath the ice directly, many studies (e.g. Calov and Hutter, 1996;Huybrechts and de Wolde, 1999;Ritz et al., 1997) assume that the average value for Pre-Cambrian Shields (Greenland bedrock) is ∼ 42 mW m −2 (Lee, 1970) although a value of 50 mW m −2 is used in EISMINT-3, and values as high as 65 mW m −2 have also been used (Greve, 2000). In terms of more recent measurements inferred from ice cores, the lowest recorded heat flux over Greenland is 38.7 mW m −2 from Dye-3 (Dahl-Jensen and Johnsen, 1986). The average value for continents is 61 mW m −2 (Lee, 1970). Although values as high as 140 mW m −2 have been measured at NGRIP (Buchardt and Dahl-Jensen, 2007;NGRIP, 2004) and values as low as 20 mW m −2 modelled (Greve, 2005), we use the range between 38 and 61 mW m −2 for the geothermal heat flux over the whole of Greenland. This is similar to the ranges used by previous sensitivity studies (Greve and Hutter, 1995;Ritz et al., 1997). We also investigate the effect of a spatially varying geothermal heat flux over Greenland (Shapiro and Ritzwoller, 2004) with all other parameters set at the default EISMINT-3 values. We compare this with the standard setup where the geothermal heat flux is 50 mW m −2 over Greenland.
Ice and snow ablation is related to near-surface air temperature by the PDD factor, which represents a simplification of processes that describe the energy balance of the glacier and overlying boundary layer. The implausibility of using one universal factor being valid for all of Greenland presents a challenge. The standard value used for ice by many modellers is 8 mm d −1 • C −1 (e.g. Huybrechts andde Wolde, 1999, Ritz et al. 1997). However, Braithwaite (1995) concluded that PDD factors for ice are generally larger than the standard value and could be as high as 20 mm d −1 • C −1 . The PDD factor for snow has also been estimated to range between 3 and 5 mm d −1 • C −1 with a standard value of 3 used by most modelling studies (Braithwaite, 1995). Modelling of PDD factors using a regional climate model in southern Greenland found ranges for the ice PDD factor between 8 and www.the-cryosphere.net/4/397/2010/ The Cryosphere, 4, 397-417, 2010 40 mm d −1 • C −1 and the snow PDD factor between 3 and 15 mm d −1 • C −1 (Lefebre et al., 2002). Other GrIS modelling studies have used higher PDD factors than the standard (e.g. Greve, 2000;Vizcaíno et al., 2008). We use a range for the ice PDD factor between 8 mm d −1 • C −1 and 20 mm d −1 • C −1 and a range for the snow PDD factor between 3 mm d −1 • C −1 and 5 mm d −1 • C −1 . The near-surface atmospheric lapse rate varies both spatially and temporally over Greenland. Lapse rate is known to vary significantly throughout the year due in part to changes in moisture content of the atmosphere. Observations from automatic weather stations indicate a mean annual lapse rate along the surface slope of −7.1 • C km −1 with seasonally varying lapse rates varying between −4.0 • C km −1 (in summer) and −10.0 • C km −1 (in winter) (Steffen and Box, 2001). Relationships derived from ERA-40 reanalysis data also yield less negative summer lapse rates of −4.3 • C km −1 at the margins and a more negative annual lapse rate of −8.2 • C km −1 for the bulk of the GrIS . Since Glimmer only uses one value for lapse rate we vary it between −4 and −8.2 • C km −1 which corresponds to the seasonal variation in lapse rate. This also encompasses the range used in the EISMINT-3 standard experiment for annual and summer lapse rate given in Eqs. (6) and (7).

Sensitivity to tuning parameters
We generate 250 plausible parameter sets using LHS and run the ice-sheet model for 50 000 years under a steady-state present day climate. Figure 6 shows the distribution of the 250 experiments with each experiment represented by a circle for three of the five tuneable parameters and the other two represented by size and colour of the circle.
In order to analyse the 250 experiments' ice-sheet geometries, four diagnostics are chosen and analysed using two skill scores. Three of these diagnostics are ice surface extent, total ice volume and maximum ice thickness. Their ability to replicate observation is described by the absolute error skill score, where zero is a perfect match. In addition, the Normalised Root Mean Square Error (NRMSE) skill score for ice thickness is used to measure the spatial fit of ice thickness over the model domain. Again, zero would describe a perfect match between modelled ice thicknesses and observed. We calculate the diagnostics with respect to the DEM derived by Bamber et al. (2001), interpolated to 20 km resolution. Figure 7 summarises the sensitivity of maximum ice thickness error, ice surface extent error and ice volume error to the five tuneable parameters.
Maximum ice thickness and ice volume are dependent on the flow law enhancement factor since faster flow will result in a thinner (and hence smaller) ice-sheet as a result of lowering the ice viscosity. An error of approximately +10% to −10% for maximum ice thickness occurs between enhancement factors 1 and 5 respectively with an optimum maximum ice thickness occurring between enhancement factors 2.5 In three dimensions geothermal heat flux (G), PDD factor for snow (α s ) and atmospheric vertical lapse rate (L G ) are shown. In addition, for each experiment the PDD factor for ice (α i ) is shown in terms of the colour-scale and the enhancement flow factor (f ) in terms of the size of circle. and 3. In contrast, the optimum enhancement factor is not reached for ice volume within the limits of the range (1 to 5) investigated. However, the enhancement flow factor has little effect on the ice surface extent due to opposing feedbacks. Faster flow will result in an increase in the flux of ice towards the ice-sheet margins. However, as the surface lowers as a result of this faster flow the ablation zone will increase at the margins leading to loss of ice. This result is similar to that found by Ritz et al. (1997) and Hebeler et al. (2008a), in terms of ice volume and maximum ice thickness. However, Hebeler et al. (2008a) found no increase in ice surface extent of their modelled region, comparable to results shown here. In contrast, Ritz et al. (1997) found an initial slight increase in ice surface extent. It is possible that this arises due to the different topography and climate configurations used as hypothesised by Hebeler et al. (2008a).
There is low sensitivity of all three diagnostics to variation in the geothermal heat flux. Since this influences basal temperatures of the ice-sheet it affects the fluidity of the ice and the flow, as well as any basal melt. Ice velocity also depends on the geothermal heat flux via the basal melt rates and in turn determines the rate of sliding of the ice-sheet. This basal sliding is predicted to occur only when the basal temperature is equal to the pressure melting point of ice. However, the original EISMINT-3 experiment did not include basal sliding and in order for a clean comparison basal sliding has also been switched off in this suite of experiments. At the icesheet margins, the basal temperature is already at the melting point and, therefore, the geothermal heat flux is not expected to influence greatly the ice volume or ice surface extent. It is, therefore, more important in the central parts of the ice-sheet where it could influence the flow of ice and affect the ice volume and maximum ice thickness via basal melt. Although basal temperatures in the interior are close to this threshold for all cases, even those with the highest geothermal heat flux, are not significant enough to cause basal melting in central parts of Greenland. As a result the geothermal heat flux parameter is unlikely to have become more important if basal sliding had been included in this suite of simulations. This is because the implication of sliding concerns the outer parts of the ice-sheet where the ice base is at melting point for all geothermal heat flux values investigated. A similar result was found by Hebeler et al. (2008a) for the Fennoscandian icesheet where very low mean annual atmospheric temperatures resulted in very low ice temperatures. As a consequence, the influence of geothermal heat flux on the thermal regime of the ice-sheet was minimal.
We also performed an experiment where the geothermal heat flux was spatially varying over Greenland (Shapiro and Ritzwoller, 2004) with all other parameters set at the default values. This was compared with the standard setup where the geothermal heat flux was uniform over Greenland. The differences are minimal with ice volume reduced by 0.3%, the ice surface extent increased by 0.4% and the maximum ice thickness reduced by 0.1%. Since basal sliding is switched off, the only effect this could have is on the basal melt and temperature of the ice at the base affecting the flow by changing the viscosity of ice.
Several parameters influence the near-surface air temperature in the EISMINT-3 experiment, including latitudinal dependency, seasonal variation and atmospheric lapse rate. Due to the PDD formulation of mass balance, these factors also directly affect ablation and ice-sheet evolution. Equilibrium ice surface extent increases with an increase in negative lapse rate (Fig. 7). A similar relationship holds for ice volume but is less pronounced. This is because a less negative lapse rate results in relatively higher near-surface air temperatures at high altitude, thereby expanding the area available for ablation. The least negative lapse rates result in the least error but are not typical of the annual lapse rate of −6.5 to −8 • C km −1 used in several studies (e.g. Ridley et al., 2005;Huybrechts and de Wolde, 1999;Vizcaíno et al., 2008). However, those that use −8 • C km −1 also include a summer lapse rate. Since Glimmer only utilises one lapse rate and since the majority of melting is assumed to occur during the spring/summer months a summer lapse rate is justified as the input lapse rate correction in the model. Maximum ice thickness is completely insensitive to lapse rate. This arises because at the ice divide, where the ice thickness is highest, temperatures are already significantly below zero. Any lapse rate correction will not influence the surface mass balance greatly.
Maximum ice thickness is also insensitive to the PDD factors for ice and snow. This is because no ablation occurs in the central part of the GrIS. However, the ice surface extent is strongly affected, decreasing with increasing PDD factors. Ice volume is also sensitive to the PDD factors but less pronounced than ice surface extent. Although varying these parameters has an effect on melting rates it does not alter the position of the ablation zones. Similar results were found by both Ritz et al. (1997) and Hebeler et al. (2008a).
The results of these sensitivity experiments show which parameters control different aspects of the geometry of the GrIS. Ice surface extent is fundamentally dependent on those parameters which control ablation (PDD factors and lapse rate) while maximum ice thickness and ice volume is controlled by parameters affecting ice flow (flow enhancement factor). All three diagnostics are insensitive to variation in the geothermal heat flux. From this suite of experiments it is possible to select one or more parameter sets which reproduce the present day GrIS with a good fit.

Selecting the optimal parameter sets
In order to select an optimal set of parameters which produce the best fit for present day ice-sheet geometry, the 250 sensitivity experiments were ranked according to each of the four diagnostics. Figure 8 shows ranking for the three absolute error skill scores on the left-hand axis and the ranking for NRMSE for ice thickness on the right-hand axis. First note that the percentage error is consistently smaller for maximum ice thickness compared with ice volume and ice surface extent.
We independently select a subset from the best-performing experiments for each diagnostic in order to assess the effect that different parameters sets could have on GrIS modelling experiments for past and future ice-sheet evolution experiments. By having parameter sets which represent different aspects of the geometry of the ice-sheet, some idea of the uncertainty in ice-sheet evolution can be obtained: for example, future warming events. One possible way to select a subset is to arbitrarily choose an ensemble size, and then choose an equal number from each diagnostics' skill score. Here we use an alternative methodology which selects the best performing experiments by identifying a step change in gradient in the best ranked experiments, as demonstrated in the insets of Fig. 8. This removes any need for an arbitrary choice and also excludes any experiments which are significantly worse but selected because an equal number from each diagnostic is required. Two experiments have been chosen according to ice volume error, two according to ice surface extent error and two according to maximum ice thickness error. The three experiments according to NRMSE for ice thickness are the same as two selected for ice volume and one selected  (1) to the closest agreement with observation (251). Observations are taken from Bamber et al. (2001) on the 20 km resolution grid. The larger symbols represent the rank position of the standard EISMINT-3 experiment. The inset graphs show the optimal experiments zoomed in for ranking from 230 to 251 for (i) maximum ice thickness, (ii) ice volume, (iii) NRMSE for ice thickness and (iv) ice surface extent. The y-scale for each inset is independent for each diagnostic in order to see the change in gradient more clearly. Filled circles/diamonds represent the optimal parameter sets for reproducing the modern day GrIS. according to ice surface extent. This provides six possible parameter sets which could be used to model the GrIS more accurately in terms of different aspects of its geometry. Figure 9 and Table 5 shows the six experiments selected and their corresponding parameter values. Figure 10 shows how well the six chosen parameter sets compare for the different diagnostic skill scores. A full unit circle would represent the experiment that out-performs all other experiments for all diagnostic skill scores. Likewise, an empty segment shows the experiment which performed worst of all experiments for that diagnostic. By comparing this measure of skill score between all 250 experiments (see Fig. 10a) one out of the six chosen parameter sets performs better than average for all diagnostics (experiment 165). Those selected according to ice volume and NRMSE for ice thickness perform significantly better than average for all diagnostics apart from maximum ice thickness (experiments 10 and 233), while those selected according to maximum ice thickness (experiments 67 and 240) perform slightly below or about average for the other diagnostics. Finally the experiment selected according to ice surface extent (experiment 99) performs better than average for all Table 5. Tuned parameter values for the six optimal experiments chosen according to diagnostic skill score identified by their experiment number.  diagnostics excluding maximum ice thickness. Figure 10b shows how well each chosen experiment compares with the other selected experiments. One will perform the worst and one the best for each diagnostic. The experiments chosen according to maximum ice thickness perform poorly for all other diagnostics, while those chosen according to ice volume and NRMSE for ice thickness perform worst for maximum ice thickness. The experiment chosen according to ice surface extent only, performs poorly for all other diagnostics while the one chosen according to ice surface extent and NRMSE ice thickness performs better than average for all diagnostics compared with the other five experiments.
Finally, the geometry of the GrIS is shown in Fig. 11 for all six tuned parameter sets and is compared with the Bamber observation (Fig. 11a). Four adequately represent the limited extent of the ice-sheet in the north and west (Fig. 11b,d,e,f) but the shape of the ice-sheet in the interior is somewhat different. However, the experiments chosen according to maximum ice thickness (Fig. 11c,g) overestimate the extent of the ice-sheet in the west and the north but represent the maximum ice thickness in the interior adequately.

Sensitivity of the Greenland ice-sheet to tuned parameter sets under future warming scenarios
In order to assess how the results from tuning affect a perturbed GrIS climate from pre-industrial, we investigate the evolution of the GrIS under differing warming scenarios. This work builds on the future warming experiments described in Lunt et al. (2009). In that study, under otherwise pre-industrial boundary conditions, CO 2 concentrations were perturbed from pre-industrial (280 ppmv Table  3 5 and Figure 10 Table 5 and Fig. 10 (experiment numbers 10, 67, 99, 165, 233 and 240 respectively). 2000). These simulations were run for a time integration of 400 model years. In addition, a future warming experiment where pre-industrial CO 2 is quadrupled to 1120 ppmv was performed. However, in order to reach equilibrium a longer time integration (665 model years) was required using a version of the GCM, HadCM3L, with a lower-resolution (2.5 • × 3.75 • compared with 1.25 • × 1.25 • for HadCM3) ocean. The ice-sheet model set-up in Lunt et al. (2009) used EISMINT-3 but with ERA-40 reanalysis reference climatology for precipitation. Anomaly coupling is used to force the ice-sheet model offline. The tuneable parameters are the same as the defaults in Table 1 but with a lapse rate at −7 • C km −1 . We also use ERA-40 precipitation for the reference climatology but where this work differs is the use of new near-surface air temperature (modified Hanna temperature) and bedrock/ice thickness (Bamber dataset) datasets, and of course the tuned parameter sets. Figure 12 shows the resultant configurations of the ice-sheet for the three warming scenarios. Figure 12a shows the results from Lunt et al. (2009) for comparison with the results using the optimal tuned parameter sets.
The original methodology with a 400 ppmv climate results in a similar ice-sheet to modern (reduced less than 2% of the modern ice-sheet). In contrast, our results using the six optimal tuned parameter sets with the more recent boundary conditions and forcings ( Fig. 12b-g) give a range of different ice-sheet configurations under a 400 ppmv climate. Although not completely collapsed, the 400 ppmv ice-sheets for Fig. 12b, d-f are somewhat reduced in the north of the island, with a reduction in ice volume compared with the modern day ice-sheet volume ranging between 20 to 23%. However, the scenario in Fig. 12c shows almost complete collapse at 400 ppmv with a reduction in ice volume of 82% while the scenario in Fig. 12g shows only a 5% reduction in ice volume. The main difference in parameter values between Fig. 12c and the other five experiments is the atmospheric lapse rate which is at least 2 • C km −1 more negative than any of the other lapse rates chosen. During icesheet retreat a more negative lapse rate will act to warm the region further and cause more surface melt than a less negative lapse rate via the temperature-elevation feedback mechanism. A warmer climate compared with pre-industrial results in increased melting during summer months. In all Fig. 12. Ice-sheet geometry for future warming scenarios (400 ppmv, 560 ppmv and 1120 ppmv CO 2 ) for (a) standard EISMINT-3 setup as shown in Lunt et al. (2009) and (b) to (g) the selected parameter sets from tuning (experiment numbers 10, 67, 99, 165, 233 and 240 respectively). See Table 5 for the tuned parameter sets corresponding to these particular experiments. cases a "tipping point" is reached whereby the temperatureelevation feedback results in ablation increasing relative to accumulation as the ice-sheet lowers and the temperature increases. This however in the case of Fig. 12c, is re-enforced by having a more negative lapse rate value resulting in rapid loss of the ice-sheet with only the highest eastern regions of the island occupied by ice. However, the other experiment selected according to maximum ice thickness (Fig. 12g) shows almost no loss of mass under a 400 ppmv climate. Although the flow enhancement factors are similar the lower PDD factors and less negative lapse rate result in less melt and no collapse of the ice-sheet.
Under a 560 ppmv climate, the GrIS is markedly reduced compared with modern with a reduction in ice-sheet volume ranging from 52 to 86%. This is not the case for the set-up used in Lunt et al. (2009) where only 7% of ice mass was lost compared with modern.
The further warming associated with quadrupling CO 2 concentrations results in almost complete elimination of the GrIS in all cases (loss of ice volume ranging from 85 to 92%). This result agrees with Lunt et al. (2009), where the ice-sheet is also shown to almost completely disappear apart from ice in the southern tip of the island and the high altitude eastern regions.
For the standard EISMINT-3 setup, results indicate a critical threshold for GrIS collapse somewhere between 560 ppmv and 1120 ppmv. However, the new parameter sets indicate a critical threshold for the GrIS becoming unstable somewhere between 400 and 560 ppmv in the majority of the simulations. There is also another possible threshold between pre-industrial (280 ppmv) and 400 ppmv where ice is lost in the north for four out of the six simulations and almost complete collapse of the ice-sheet for one of the remaining two experiments.
Comparisons can also be made with similar studies using different GCMs and/or ice-sheet models. For instance, Ridley et al. (2005) showed the ice-sheet collapsed to 7% of its original volume under a quadrupled CO 2 climate. The extra ice mass in our simulations (1 to 8% extra) can partly be accounted for by the ice present in southern Greenland which is absent in Ridley et al. (2005). This is likely due to the absence of the ice-albedo feedback between climate and ice-sheet, which is included in their simulations by interactive coupling of the GCM to the ice-sheet model. Interestingly the study of Mikolajewicz et al. (2007) shows that under a 560 ppmv climate using a fully coupled climate icesheet model the GrIS could result in significant melting in the long-term (simulation only carried out for 600 years). Furthermore, Alley et al. (2005) showed that under a doubled CO 2 climate the GrIS would eventually almost completely disappear.

Discussion and conclusions
In this section we discuss the sources of uncertainty and the missing processes in the experimental design and the influence this has on the conclusions drawn.
Firstly, several other potential temperature datasets over Greenland exist to force the Glimmer ice-sheet model. A new parameterisation based on more up-to-date automatic weather station data, for instance, is now available with a similar form to Eqs. (6) and (7) (Fausto et al., 2009). However, for this work we chose the ERA-40 derived temperature product for consistency with the new precipitation dataset. Furthermore, datasets also exist in terms of satellite products. For satellite datasets, temperature data are available from the Advanced Very High Resolution Radiometer (AVHRR) Polar Pathfinder (APP) from 1982-2004 which is collated twice a day at the local solar times of 14:00 and 04:00. Although the data is initially on a 5 km resolution it is subsampled at 25 km pixels. The APP-x product includes allsky surface temperature with the cloudy-sky surface temperatures calculated using an empirical relationship between clear-sky surface temperature, wind speed, and solar zenith angle (daytime). However, this only applies to surface temperatures over sea-ice and not land. Therefore, temperatures over Greenland are based only on data from clear-sky retrieval with temperatures in cloudy regions interpolated from clear-sky areas. Although useful for comparing with present day surface temperatures from climate models, this dataset is not suitable to directly force an ice-sheet model over Greenland because (a) the largest uncertainties are likely to be over Greenland (J. R. Key, personal communication, 2010), (b) no associated orography exists which is used to downscale from the resolution of the forcing data onto the high-resolution of the ice-sheet model and (c) sensitivity studies using Glimmer indicate that the APP-x temperatures were significantly too low, in observed ice-free regions such as western Greenland, (by up to 12 • C in western Greenland compared with EISMINT-3 temperatures which have at least been derived from surface observation) to reproduce a reasonable modern day ice-sheet without tuning ice-sheet model parameters beyond uncertainty ranges. This could, in part, be due to the satellite recording ice surface temperatures rather than air temperature. Furthermore, clear-sky retrievals errors are predominantly due to uncertainties in cloud detection (Key et al., 1997) particularly during the night. The low temperatures, bright surface and high elevation make remote sensing over Greenland particularly difficult in terms of accurate cloud detection.
Secondly, in contrast with many studies, we spin up the model from present day initial conditions without taking the climate history into account. Since the GrIS is still affected by past climatic change this assumption must be justified. The main method used to spin up the ice-sheet model over several climatic cycles has caveats of its own. It uses a temperature forcing derived from a smoothed ice core record and has been used in several studies (e.g. Huybrechts and de Wolde, 1999;Ridley et al., 2005;Vizcaíno et al., 2008). However, uncertainty exists in the functions used to derive a reliable temperature record and subsequent accumulation record from an oxygen isotopic record although new and more sophisticated methods are being developed (Cuffey and Marshall, 2000;Lhomme et al., 2005). The effect of ice flow processes on deeper parts of ice cores also makes them somewhat unreliable and extending beyond the last interglacial is somewhat unrealistic (Grootes et al., 1993;Johnsen et al., 1997). For these reasons we only initiate the ice-sheet model from the present day initial conditions, which we can be certain are relatively accurate.
Thirdly, the process of basal sliding was not included in the experimental design, which has implications for the amount of ice mass lost dynamically. An increase in the ice velocity, by incorporating the sliding velocity (see Eq. 2), would result in more ice transferred from the accumulation zone to the ablation zone and, therefore, reduce the volume of the ice-sheet. Inclusion of this missing process could result in lower PDD factors than those obtained in the tuning exercise presented here. Indeed, the study by Parizek and Alley (2004) showed an increase in GrIS sensitivity to various warming scenarios due to surface meltwater lubrication of flow. Recent modelling developments have also investigated the potential positive feedbacks from including basal sliding on the inland migration of fast-flowing glaciers increasing the drawdown of the ice-sheet interior (e.g. Price et al., 2008). Currently, Glimmer has a simplified representation of basal sliding and the basal hydrology. Furthermore, there is no representation of the sediment deformation. The presence of unconsolidated sediments alters the hydrological system by incorporating melt water until saturation is reached. This reduces the yield stress of the material substantially and deformation of the basal till by the overlying ice load inducing glacier motion. However studies have mainly focussed on the local scale of ice streams rather than the continental scale of ice-sheets (Tulaczyk et al., 2000;Sayag and Tziperman, 2008).
Fourthly, current large-scale ice-sheet models lack higherorder physics, and although able to simulate slow moving ice dynamics adequately, they are not yet able to represent the dynamics of fast-moving ice streams. Recent work has indicated that current net mass loss from the GrIS is roughly equally partitioned between surface mass balance changes and changes in dynamics . Development of ice-sheet models in these areas is currently being researched with improvements to ice dynamics (e.g. Pattyn, 2003;Soucek and Martinec, 2008), and inclusion of an accurate representation of the fast ice streams and ice shelves (Pattyn et al., 2006;Schoof, 2006Schoof, , 2007. Recent observations of glaciers in Greenland have documented rapid changes in marginal regions of the ice-sheet with increased flow velocities observed on Jakobshavn Glacier (Joughin et al., 2004) and on other glaciers (e.g. Howat et al., 2007;Rignot and Kanagaratnam, 2006). The inclusion of these fast flowing ice streams in ice-sheet models could lead to larger dynamical changes in the ice-sheet than currently predicted by models at least on relatively short timescales of hundreds of years. Incorporation of these fast flow features in the ice-sheet model could also result in lower PDD factors from tuning. Furthermore, if these dynamical changes are marine-driven then for long-term future ice-sheet predictions, once the ice-streams are no longer marine terminating, the dynamical changes will cease.
It has also been shown that processes at the ice margin have a strong influence on the surface extent of the ice-sheet but are poorly accounted for with a coarse grid of 20 km resolution. The use of energy-balance/snow pack models (EBSM) to predict surface mass balance (e.g. Bougamont et al., 2007) as opposed to the PDD approach has been shown to give contrasting results under a four times CO 2 climate with the PDD scheme significantly more sensitive to a warming climate generating runoff rates almost twice as large compared with an EBSM. However, some aspects of these results are not undisputed (P. Huybrechts, personal communication, 2009). The ablation zone on Greenland varies from only 1 km wide along the southeast coast and up to 150 km wide along the southwest coastline and, therefore, requires a very high horizontal resolution if ablation is not to be over or underestimated in the model (van den Broeke, 2008). Future development of the EBSM approach using a finer grid of 5 km resolution could result in a marked improvement for modelling ablation processes. It would also be highly beneficial to downscale to a 1 × 1 km resolution using a PDD approach (e.g. Janssens and Huybrechts, 2000) and the highresolution Greenland DEMs now available (e.g. Bamber et al. 2001).
Fifthly, the grid on which the ice-sheet dynamics are solved could influence the model outcome. An alternative to the finite difference modelling approach used here could be to instead implement the finite element modelling method. This has the advantage that the element size can be reduced in areas of high gradient and increased in areas of low gradient. Furthermore, the model can conform to irregular boundaries that are awkward to model with rectangular elements used in the finite differences technique. Currently this methodology is used over smaller domains such as individual glaciers (e.g. Zwinger et al., 2007) or within flow line models of icesheets (e.g. Parizek, 2005).
Finally, overcoming the abstraction required for largescale ice-sheet models, in order to keep computing demands to a minimum while ensuring spatial variability at the sub-scale level is captured, presents a challenge. However, subgrid parameterisation for the calculation of ablation/accumulation has been shown to be effective in compensating for dependencies on scale while incurring only a small additional computational cost (Hebeler and Purves, 2008b;Marshall and Clark, 1999).
We evaluate the sensitivity to boundary conditions and climate forcings in the context of modelling the evolution of the GrIS under present day, steady-state conditions and show the geometry and size of the ice-sheet is highly sensitive to the initial condition of bedrock and ice thickness. An icesheet volume 13.7% larger than that produced with the Letreguilly dataset results with the new and improved Bamber dataset. Overall, our study indicates that using the more recent datasets for forcings and boundary conditions with the standard set of model parameters (Table 1) give a poorer representation of the modern ice-sheet, with an ice-sheet volume 33% larger than observation. The results further show that topography and its inherent uncertainty has a significant effect on ice-sheet geometry obtained from large-scale models of considerable abstraction such as Glimmer. Therefore, the use of more realistic topography and climate data on an original resolution significantly higher than that used in Glimmer may not be entirely suitable for current large-scale ice-sheet modelling.
Several parameters are not well-constrained in large-scale ice-sheet modelling and can influence ice-sheet volume and extent. We performed a sensitivity/tuning study in order to assess the importance of certain parameters on the geometry and size of the GrIS. The method of LHS was used in order to efficiently vary more than one parameter at a time to obtain a best fit between modelled and observed geometry. The maximum ice thickness and ice volume were shown to depend on the factors affecting ice flow. In this case increasing the flow enhancement factor makes the ice flow faster which lowers the height of the ice dome. The ice surface extent is predominantly dependent on the PDD factors and the atmospheric lapse rate. Although geothermal flux can affect ice flow since it acts to melt the ice, which is a prerequisite for basal sliding, this had little effect on the simulations presented here because basal sliding was switched off.
By selecting "best fit" experiments according to different skill score diagnostics a range of parameter sets can be used for assessing the uncertainty in ice-sheet modelling experiments by analysing the resultant geometries. The sets of parameters that give the best fit to the present measured icesheet are somewhat different from the standard set most commonly used by ice-sheet modelling studies. Higher PDD factors than the standard (10.2 to 19.9 mm d −1 • C −1 for α i and 3.6 to 4.8 mm d −1 • C −1 for α s ) are required in all cases in order to account for both ablation and calving processes at the margin. The lack of basal sliding in these simulations means that these higher PDD factors are likely partially compensating for this missing process. Furthermore, less negative atmospheric lapse rates (five out of the six tuned parameter sets ranged between −4.1 and −6.0 • C km −1 ) are generally needed to produce a good fit in terms of volume by reducing the growth of the ice-sheet.
The parameters varied using LHS are strictly independent in a mathematical sense. However, it is possible that the values chosen could have similar and opposite effects on accurately predicting the present day GrIS geometry. For example, high PDD factors in combination with low lapse rates could simulate a good representation of the GrIS. In our conclusions we do not attempt to make a probabilistic interpretation of the results such that certain combinations are more likely than others in producing an accurate representation of the ice-sheet.
The optimal parameter sets chosen to best represent the modern day GrIS were used to assess their effect on the evolution of the ice-sheet under future warming scenarios. We obtained a different threshold for ice-sheet collapse, occurring somewhere between 400 ppmv and 560 ppmv compared with previous work which suggested a threshold between 560 and 1120 ppmv (Lunt et al., 2009) when using the same models. Differences in ice-sheet geometry and volume also occur between the optimal parameter sets. Although all icesheets modelled for present day showed complete glaciation of Greenland, one particular parameter set (Table 5, experiment 67) showed complete collapse at 400 ppmv. We show under perturbed climates from present day the evolution of the GrIS behaves differently for the parameter sets tuned in the model. This work suggests that, if possible, tuning exercises should be applied to the GrIS under several different climatologies. Since observations are required for comparison this is somewhat restrictive. However, examples of alternative climates to the present day could be the last deglaciation or the Last Glacial Maximum, for which there exist some data on ice-sheet extent.
We have shown that future predictions of the GrIS are highly sensitive to a number of factors relating to the physical basis of the ice-sheet model. Most current models neither have a robust representation of the fast flowing processes, nor are the parameters which influence the ice physics tightly constrained. As a result future development of the ice-sheet model to improve the representation of these processes may lead to different behaviour under warm climate conditions. The lack of higher-order physics, low resolution, absence of basal sliding and subglacial hydrology and highly parameterised surface mass balance, inevitably means that the tuning presented here compensates for these absent processes in order to replicate as closely as possible the present day GrIS. As a result, future predictions of the GrIS should be evaluated with some caution in the context of these sensitivities and deficiencies of the ice-sheet model.