Quantification of the Greenland ice sheet contribution to Last Interglacial sea level rise. Climate of the Past , 9 , 621-639.

. During the Last Interglacial period ( ∼ 130–115 thousand years ago) the Arctic climate was warmer than today, and global mean sea level was probably more than 6.6 m higher. However, there are large discrepancies in the estimated contributions to this sea level change from various sources (the Greenland and Antarctic ice sheets and smaller ice caps). Here, we determine probabilistically the likely contribution of Greenland ice sheet melt to Last Interglacial sea level rise, taking into account ice sheet model parametric uncertainty. We perform an ensemble of 500 Glimmer ice sheet model simulations forced with climatologies from the climate model HadCM3, and constrain the results with palaeodata from Greenland ice cores. Our results suggest a 90 % probability that Greenland ice melt contributed at least 0.6 m, but less than 10 % probability that it exceeded 3.5 m, a value which is lower than several recent estimates. Many of these previous estimates, however, did not include a full general circulation climate model that can capture atmospheric circulation and precipitation changes in response to changes in insolation forcing and orographic height. Our combined modelling and palaeodata approach suggests that the Greenland ice sheet is less sensitive to orbital forcing than previously thought, and it implicates Antarctic melt as providing a substantial contribution to Last Interglacial sea level rise. Future work should assess additional uncertainty due to inclusion of basal sliding and the direct effect of insolation on surface melt. In addition, the effect of uncertainty arising from climate model structural design should be taken into account by performing a multi-climate-model comparison.


Introduction
Past time periods provide important case studies for evaluating the performance of earth system models, because model results can be compared with geological records. In particular, warm climates of the past are useful because they can provide an analogue for possible future warming. The Last Interglacial (LIG) provides such a case study as globally averaged sea level was thought to be several metres higher than today, and high-latitude temperatures warmer. Estimates of maximum sea level increase, derived from sedimentary deposits and coral sequences, typically range from 4 to 6 m (Rostami et al., 2000;Muhs et al., 2002). A recent sea level data synthesis shows that sea level likely exceeded 8 m higher than today with the highstand extremely likely (95 % probability) greater than 6.6 m (Kopp et al., 2009), consistent with less glacial ice on Earth during the LIG. The likely contributors to the sea level rise are ice losses from the Greenland and Antarctic ice sheets along with high-latitude Arctic ice fields, such as those in the Canadian Arctic, together with thermal expansion of seawater.
Further evidence from proxy data located in the Arctic and European regions suggests the LIG climate featured temperatures, at least regionally, several degrees warmer than today CAPE members, 2006). This estimated temperature increase is supported by climate model simulations indicating summer Arctic warming up to 5 • C relative to modern day, with the greatest warming over Eurasia and in the Baffin Island/Greenland region (Montoya et al., 2000;Kaspar et al., 2005;Otto-Bliesner et al., 2006). Palaeopollen, macrofossil and soil records suggest the expansion of boreal forests northwards into areas now occupied by tundra in Russia, Siberia and Alaska during peak LIG warmth (Muhs et al., 2001;Kienast et al., 2008). On Greenland itself, ice core measurements from the Summit region (NorthGRIP, GRIP and GISP2 ice cores) indicate ice was present during the LIG NorthGRIP, 2004;Suwa et al., 2006), with the surface elevation no more than a few hundred metres lower than present day, based on the total gas content of the ice (Raynaud et al., 1997). In addition, basal ice from the northwestern ice core, Camp Century, has been proposed to be of LIG age and ice from the bottom section of a core from the Renland peninsula in eastern Greenland is dated older than 130 thousand years ago (ka) (Johnsen et al., 2001) although there is uncertainty in the dating of these two ice cores (Alley et al., 2010). New results from the NEEM ice core project, located in northern Greenland, indicate basal ice in this location is of LIG age (NEEM community members, 2013). The dating of basal ice at Dye-3 in southern Greenland, however, remains highly uncertain (e.g. Koerner and Fischer, 2002;NorthGRIP, 2004).
Estimates of the Greenland ice sheet (GrIS) contribution to sea level rise during the LIG range from 0.4 to 5.5 m based on a wide range of modelling techniques. These include palaeothermometry from ice cores coupled with thermodynamical ice sheet models (Letréguilly et al., 1991;Ritz et al., 1997;Cuffey and Marshall, 2000;Huybrechts, 2002;Greve, 2005) with similar methodological studies also constraining their results by matching model-predicted isotopic stratigraphy from ice cores with data (Tarasov and Peltier, 2003;Lhomme et al., 2005). Another method uses coupled climate -ice sheet models of varying complexity (Otto-Bliesner et al., 2006;Fyke et al., 2011;Robinson et al., 2011;Born and Nisancioglu, 2012) to predict LIG GrIS geometry and sea level contribution.
Here we assess the contribution of Greenland ice loss to global sea level rise, derived from simulations of the LIG global climate and evolution of the GrIS from 130 to 120 ka, using the general circulation model (GCM) HadCM3 pseudo-coupled to the ice sheet model Glimmer. We use an efficient offline coupling methodology to account for ice sheet-climate interactions (DeConto and Pollard, 2003) and estimate the range in GrIS contribution to LIG sea level change by considering ice sheet model uncertainty. As a result, this method aims to better understand the GrIS response to a warmer-than-present climate, which is critical for the assessment of future climate change.

Model description and experimental design
Here, we outline the models used to estimate the GrIS contribution to LIG sea level change. Due to computational expense we have developed a method to pseudo-couple our climate model to an ice sheet model which takes into account the effect of the albedo feedback mechanism without the need to run fully coupled (two-way) climate-ice sheet simulations. We begin by describing the main features of the climate and ice sheet models, followed by a detailed description of the experimental design and this coupling methodology. An ensemble is performed to take into account parametric uncertainty in the ice sheet model with the aim of estimating a range in GrIS contribution to LIG sea level. We use palaeodata to disregard simulations which do not satisfy these robust palaeodata ice sheet constraints (see Sect. 3.2). Finally, from the ensemble a probability density distribution of maximum sea level contribution from the GrIS to LIG sea level rise is constructed (see Sect. 3.2.2 for details of the probabilistic method).

The climate model
The GCM simulations described in this paper are carried out using the UK Met Office coupled atmosphere-ocean GCM, HadCM3, version 4.5 (Gordon et al., 2000), which has been used in the third and fourth IPCC assessment reports. The atmosphere component of HadCM3 is a global grid-point hydrostatic primitive equation model, with a horizontal gridspacing of 2.5 • (latitude) by 3.75 • (longitude) and 19 levels in the vertical with a time step of 30 min. The performance of the atmosphere component has been shown to agree well with observations (Pope et al., 2000). The land surface scheme (MOSES 2.1) includes representation of the freezing and melting of soil moisture and the evaporation process. Within this land surface scheme ice sheets are prescribed and are fixed.
The resolution of the ocean model is 1.25 • by 1.25 • with 20 levels in the vertical. The ocean model uses the mixing scheme of Gent and McWilliams (1990) with no explicit horizontal tracer diffusion. The sea ice model uses a simple thermodynamic scheme and contains parameterisations of ice concentration (Hibler, 1979) and ice drift and leads (Cattle and Crossley, 1995). For simulations of the presentday climate, the ocean model has been shown to simulate sea surface temperatures in good agreement with modern observations, without the need for flux corrections (Gregory and Mitchell, 1997).

Model description
We also use the three-dimensional thermomechanical ice sheet model Glimmer version 1.0.4 (Payne, 1999;Rutt et al., 2009), which is forced with monthly temperature and precipitation from HadCM3. The horizontal resolution of the model is 20 km with 11 vertical layers. The ice dynamics are represented with the widely used shallow-ice approximation (SIA) approach, which neglects longitudinal stresses in the ice sheet. This simplification is appropriate for ice masses that are thin compared with their horizontal extent.  (Stone et al., 2010). Also included are the mean and standard deviation for each parameter used in Eq. (11).

Standard Parameter
Range Mean (µ) deviation (σ ) Positive degree day factor for snow, α s (mm water d −1 • C −1 ) 3.0 to 5.0 4.0 ± 1.2 Positive degree day factor for ice, α i (mm water d −1 • C −1 ) 8.0 to 20.0 14.0 ± 6.9 Enhancement flow factor, f 1.0 to 5.0 The principle advantage of using the SIA for modelling the GrIS on palaeotimescales is that it is computationally cheap, allowing large, multi-millennial ensembles to be easily performed. Although the method is accurate for the interior of a large ice sheet such as the GrIS, this is not the case at the margins. Here, streams of fast-flowing ice and coupling to ice shelves complicate the ice dynamics such that the SIA is unable to capture the currently observed changes in ice sheet geometry and velocity, which occur on short timescales. The lack of higher-order physics has resulted in the majority of ice sheet models overestimating the present-day ice sheet volume and extent (e.g. Ritz et al., 1997;Stone et al., 2010;Greve et al., 2011;Robinson et al., 2011). The surface mass balance is simulated using the positive degree day (PDD) approach described by Reeh (1991). The basis of the PDD method is the assumption that the melt that takes place at the surface of the ice sheet is proportional to the time-integrated temperature above the freezing point, known as the positive degree day. Two PDD factors are used, one each for snow and ice, to take account of the different albedos and densities of these materials. The use of PDD mass balance models is well-established in coupled atmosphereice sheet palaeoclimate modelling studies (DeConto and Pollard, 2003;Lunt et al., 2008Lunt et al., , 2009. Limitations of using the PDD scheme are discussed in Sect. 4. Glimmer also includes a representation of the isostatic response of the lithosphere to a change in ice mass. The response of the lithosphere is assumed to behave elastically, based on the isostasy model of Lambeck and Nakiboglu (1980). The forcing data from HadCM3 are transformed onto the ice model grid using an adaption of bilinear interpolation termed "mean-preserving interpolation", which preserves the mean of the global grid box on the Glimmer grid and ensures that precipitation is conserved in the atmosphere-ice sheet coupling. More details of the interpolation scheme are given in Sect. 6.2 of Rutt et al. (2009). In the case of the surface air temperature field, a spatially homogenous vertical lapse rate correction is used to take account of the difference between the high-resolution topography seen within Glimmer and that represented within HadCM3. The use of a lapse rate correction to better represent the local temperature is established in previous work (e.g. Pollard and Thompson, 1997;Vizcaíno et al., 2008).
One limitation of the ice sheet model simulations presented here is that they do not include the process of basal sliding which has implications for the amount of ice mass lost dynamically. An increase in the ice velocity, by incorporating the basal sliding velocity, would result in more ice transferred from the accumulation zone to the ablation zone and, therefore, would likely reduce the volume of the ice sheet under a warm climate. Inclusion of this missing process could lead to a smaller GrIS during the LIG. 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. However, previous studies Robinson et al., 2011) have shown that, although the sliding coefficient parameter affects GrIS geometry, it is less significant compared with other parameters in determining the past evolution and present geometry of the modelled GrIS.

Generation of an ice sheet model ensemble
For the baseline climate to which the GCM temperature and precipitation anomalies are applied we use those described in Stone et al. (2010). The temperature climatology is derived from ERA-40 observations (Hanna et al., 2005) and precipitation also from ERA-40 reanalysis (Uppala et al., 2005). We use the Greenland bedrock topography of Bamber et al. (2001) on a 20-km-resolution grid.
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 . Previous work (Stone et al., 2010) investigated the sensitivity of ice sheet evolution for the modern GrIS to five tuneable parameters which affect the ice sheet dynamics and surface mass balance. These parameters are the PDD factors for ice and snow, near-surface lapse rate, flow enhancing factor and the geothermal heat flux (see Table 1).
Here we generate an ensemble of 500 simulations using the statistical method of Latin Hypercube Sampling (LHS) in order to efficiently sample the five-dimensional parameter space. This method generates a distribution of plausible parameter sets within a prescribed set of ranges (McKay et al.,  1979). Parameter values are sampled from the prescribed distribution of each variable and paired randomly with the other variables, assuming that the variables are independent of one another, ensuring the parameter space is as evenly covered as possible (which is the case here). The LHS distribution is given in Fig. 1. For more details on parameter choices refer to Stone et al. (2010).

Experimental design and coupling methodology
Computationally, it is not yet feasible to run HadCM3 fully coupled (two-way) with Glimmer for timescales of thousands of years, such as through the LIG. A methodology is developed based on that of Deconto and Pollard (2003) in order to account for a transient climate which evolves as the ice sheet volume evolves, whilst minimising computational expense. It takes into account a changing climate as a result of the change in ice sheet geometry by including the elevation-temperature feedback and an approximation to the albedo feedback. We outline (1) the GCM simulations performed, (2) the ice sheet model spin-up procedure and (3) details of the coupling method used between climate and ice sheet model.

LIG GCM simulations
GCM simulations representing 130, 125 and 120 ka are forced with insolation anomalies resulting from changes in the Earth's orbital parameters for the early to mid-part of the LIG. These time-slices are chosen because they cover the interval of peak LIG warmth as well as the maximum sea level highstand (Petit et al., 1999;Lisiecki and Raymo, 2005;Kopp et al., 2009). Compared with the pre-industrial period, larger eccentricity, obliquity and Northern Hemisphere summer (as opposed to winter) occurring at perihelion (see Table 2) result in greater seasonality, leading to pronounced high northern latitude summer insolation, consistent with warming observed in the geological record (North-GRIP, 2004;Kaspar et al., 2005;CAPE members, 2006) (see Fig. 2). This seasonal variation in insolation is important because ice sheet surface mass balance is particularly sensitive to summer warming.
The three LIG snapshot time-slices are run for 100 model years (70-yr spin-up and 30 yr for averaging) with the following Greenland boundary conditions: 1. modern-day GrIS present, 2. partial GrIS present derived from a tuned ice sheet model experiment forced with a 560-ppmv climate (Stone et al., 2010), 3. no GrIS present with bedrock in isostatic equilibrium.
This procedure gives a range of climate states between which the "expected" climate over a partially melted GrIS during the LIG might lie. One caveat of these climate simulations  (Berger and Loutre, 1991) and the maximum 65 • N June insolation anomaly. Also shown for comparison are the parameters for the pre-industrial period. concerns the use of isostatic equilibrium for the orography in the ice-free state. Obviously, if there was a substantial ice sheet present before the start of the LIG, as inferred from the eustatic sea level curve (Siddall et al., 2007), there would likely have been insufficient time for all the ice to melt, the bedrock to rebound fully and soil to develop on the bare rock surface. However, this provides the most contrasting climate scenario to a fully glaciated Greenland being present throughout the LIG (which is also unlikely). Another limitation of this approach (which was required for computational efficiency) is that the climate state is constrained to linearly interpolate between these states (see Sect. 2.3.3). For the LIG, the changed forcings from present day are the following: the modified trace gas concentrations and the seasonal and latitudinal insolation changes at the top of the atmosphere associated with Milankovitch orbital forcing (Milankovitch, 1941). For this modelling framework the perturbed astronomical values are consistent with those used in the standard Paleoclimate Modelling Intercomparison Project (PMIP) LIG simulations (see http://pmip.lsce. ipsl.fr/). Figure 2a shows the variation in insolation from 140 to 110 ka for the spring and summer months at three latitudes over Greenland: 65 • N, 74 • N and 80 • N. Insolation anomalies over Greenland relative to present day ( Fig. 2b) are at a maximum at ∼ 130 ka for May and June and decrease toward 120 ka. Smaller anomalies, for July and August, peak from ∼ 125 to 120 ka. Orbital parameters are taken from Berger and Loutre (1991) for the three snapshots at 130, 125 and 120 ka. Table 2 shows the obliquity, eccentricity and perihelion for these three scenarios. A further HadCM3 experiment at 136 ka is also included in order to spin up the ice sheet model sufficiently but differs slightly by including a MOSES 1 land surface scheme (Cox et al., 1999). This simulation is run for 500 model years with an averaging time of 30 yr.
Also shown in Fig. 2 is the atmospheric CO 2 concentration, reconstructed from ice cores, from 140 to 110 ka based on Lüthi et al. (2008). All CO 2 values are on the EDC3 gas age scale (Loulergue et al., 2007). There is a sharp rise in CO 2 concentration between 140 ka and 130 ka from ∼ 200 to 260 ppmv. Thereafter, this trace gas concentration stabilises between 260 and 290 ppmv. Since the greenhouse gases do Fig. 1. Distribution of 500 Glimmer parameter experiments produced by Latin hypercube sampling and projected onto two-dimensional slices through the five-dimensional space. The parameters are as follows: the geothermal heat flux (G), the positive degree day (PDD) factor for snow (α s ), the PDD factor for ice (α i ), the atmospheric vertical lapse rate (L G ) and the enhancement flow factor (f ). The experiments highlighted in red are those which are valid for the LIG. not markedly vary from the pre-industrial period during the LIG (Lüthi et al., 2008) and it has been shown that climate perturbations were predominantly orbitally driven at this time (Slowey et al., 1996;Loutre et al., 2007;Yin and Berger, 2012), gas concentrations are held constant and unchanged from the values used in the pre-industrial simulations. In this way any changes in LIG climate from the the pre-industrial period are due to changes in the orbital parameters of the Earth. CO 2 is, therefore, held constant at 280 ppmv for all experiments performed using HadCM3 between 130 and 120 ka. All other trace gases are equivalent to pre-industrial values. The exception is for the simula-tion at 136 ka where CO 2 , methane (CH 4 ) and nitrous oxide (N 2 O) are lower compared with the pre-industrial period at 200 ppmv, 413 ppbv and 229 ppbv, respectively. This is because differences in the trace gases compared with the preindustrial period are the driving mechanism for this earlier perturbed climate rather than changes in the orbital parameters compared with the pre-industrial period (see Fig. 2b where summer high-latitude insolation anomalies are small at 136 ka).
An additional simulation, the pre-industrial control, includes trace gas concentrations (280 ppmv for CO 2 , 760 ppbv for CH 4 and 270 ppbv for N 2 O) and orbital parameters  Petit et al. (1999) and Pépin et al. (2001) for the LIG (they are on the EDC3 gas age scale (Loulergue et al., 2007)). The colours correspond to the following months: May (light blue), June (blue), July (orange) and August (green). Line styles refer to different latitudes over Greenland.
(obliquity 23.45 • , perihelion occurs on day 2.6 of the year and eccentricity 0.01724) appropriate for 1850 AD. Outside of Greenland, global vegetation coverage is prescribed at present-day distributions. The simulations where the GrIS is removed/partially melted are prescribed bare soil coverage in place of Greenland ice, while the simulations with a full GrIS included use the present-day ice sheet mask with bare soil in ice-free regions. The land-sea mask remains unchanged from modern day since there were no significant tectonic changes to the continents between 130 ka and present and the estimated sea level change would result in negligible land-sea mask changes.
The 130-, 125-and 120-ka GCM simulations are continued from pre-industrial simulations of 100 model years with the appropriate bedrock and ice coverage. The spin-up time in large-scale atmosphere-ocean models is governed by the slow processes in the deep ocean and is usually on the order of several thousand years. However, due to computational expense this is not easily achievable. As such the ocean component of HadCM3 does not fully represent changes in Greenland temperature time series for the three LIG snapshots with a GrIS, partial GrIS and without a GrIS included. The first 100 yr represent pre-industrial greenhouse gases and orbital conditions. The last 100 yr are the temperature response to changed orbital parameters. The black line is the annual mean and the red line is the 10-yr mean. The thick black horizontal line shows the 30-yr annual Greenland temperature average required to force the ice sheet model. ocean circulation, but it does fully interact thermodynamically with the atmosphere in our modelling framework. Figure 3 shows the average temperature evolution over Greenland (one of the inputs into the Glimmer ice sheet model) including this pre-industrial spin-up. A 10-yr mean trend (red) is shown and indicates sufficient spin-up of the model nearsurface temperature in response to the changed orbits. This trend shows that compared with inter-annual variability the multi-year average temperature response of the simulations is close to equilibrium.

Obtaining a 136 ka GrIS
The size of the GrIS is not known at 130 ka (or at any other point during the LIG), although sea level was similar to present day (Siddall et al., 2007;Kopp et al., 2009), implying a substantial amount of ice must have been present at high northern and southern latitudes. The correct method of reconstructing the initial state of the GrIS under past climate forcings is unclear, but two main methods have been adopted in previous studies: (1) steady state simulations driven by present-day or past climatic conditions (e.g. Ritz et al., 1997;Stone et al., 2010) and (2) transient simulations driven by palaeoclimatic reconstructions (e.g. Applegate et al., 2012). Each method has its own caveats which have been investigated recently by Rogozhina et al. (2011). For example, they show that initialising from an ice-free state under glacial forcings is not a good choice for simulations that start under colder-than-modern conditions. Because it is not practically possible to spin up an ensemble of coupled HadCM3 ice sheet model configurations for several glacial-interglacial cycles, we use an approach that assumes the ice sheet is in equilibrium at the start of the transient ice sheet model simulations. We adopt a similar methodology to Rogozhina et al. (2011) by initialising from a modern GrIS geometry spun-up with a constant glacial climate forcing derived from HadCM3 then apply a time-dependent forcing into the interglacial period. We do not use palaeoclimatic reconstructions to obtain an initial state for the GrIS because, prior to the onset of the LIG, processes occurring in deeper parts of ice cores make them somewhat unreliable and extending beyond the LIG is, therefore, unrealistic (Grootes et al., 1993;Johnsen et al., 1997).
In order that changes in the ice sheet response to climate at 130 ka are not a result of inadequate spin-up of the ice sheet model, simulations begin at 136 ka when the climate was substantially colder. As a result, the ice sheet model is initiated with an ice sheet in equilibrium with the 136ka climate. The ice sheet model is spun up for 50 000 yr in anomaly mode using the 136 ka climatology. This method requires GCM monthly mean changes in precipitation and near-surface temperature (defined relative to a pre-industrial climate) to be superimposed onto a present-day reference climatology (see Sect. 2.2.2) used by the surface mass balance model in Glimmer. Anomaly coupling is used to reduce climate model bias both for precipitation and temperature which affects the ice sheet model output, as in previous studies (Lunt et al., 2008(Lunt et al., , 2009).

Coupling the climate and ice sheet models
We model a total of 16 000 yr, representing the time period from 136 to 120 ka. Figure 4 shows a diagram of the coupling process, which is outlined in detail below. The monthly average variables of temperature and precipitation, here denoted as C(t), are linearly interpolated along the time-axis from 136 to 130 ka where the notation C state of Greenland is used (i.e. state of Greenland in HadCM3 is either ice covered: ice; partial ice: pice; or ice-free: 0): The interpolation is between the 136 ka climate, C ice (136), and the 130 ka climate, C ice (130), where t 1 is 6000 model years. Glimmer is initiated with the equilibrated ice sheet geometry described above. At 130 ka the climate is allowed to evolve each year between the three climate scenarios (with a GrIS, a partial GrIS and without a GrIS) according to a weighting function defined by the ratio of the ice volume (V (t)) at time t and the ice volume predicted at 130 ka (V (130) the solid blue, orange and red arrows, respectively, in Fig. 4) similar to Eq. (1): and where C ice (125) is the 125-ka climate with the GrIS present; C pice (125) and C pice (130) are the 125-and 130-ka climates, respectively, with a partial GrIS; C 0 (125) and C 0 (130) are the 125-and 130-ka climates, respectively, with the GrIS removed; and t 2 is 5000 yr. Likewise, similar linear interpolations are also performed from 125 to 120 ka. If the ice volume, V (t), is greater than the partial ice volume (defined as V pice = 0.46V ice (130)), then the climate, C(t), at each year is now also weighted either towards the climate with a partial GrIS, C pice (t), or the GrIS climate, C ice (t), according to Overlain is the maximum observed LIG summer temperature anomalies from palaeotemperature proxies (terrestrial (circles) and marine (triangles)) (CAPE members, 2006;Kaspar et al., 2005). White regions are not statistically significant (at the 95 % confidence interval).
towards the climate with no GrIS, C 0 (t), or the partial GrIS climate, C pice (t) , according to 3 Results

The modelled climate of the LIG
The GCM-simulated annual average global temperature anomaly at 130 ka (with a modern-day fixed GrIS included) is only 0.13 • C relative to the pre-industrial period, consistent with the small mean annual forcing associated with the orbital configuration for the LIG. However, the seasonal temperature anomaly is −1.6 • C and 2.0 • C in the Northern Hemisphere for winter and summer, respectively. Figure 5 shows a comparison of the LIG simulated Northern Hemisphere maximum summer warming with reconstructed terrestrial temperature anomalies derived from ice cores, pollen and macrofossils and reconstructed marine peak summer sea surface temperature anomalies derived from foraminifera CAPE members, 2006). Overall, the agreement at high Northern Hemisphere latitudes is good with 65 % of the data points coinciding (within the uncertainty) with the 1 : 1 line in Fig. S1 (see also Table 3). However, during the summer months the maximum LIG average temperature anomaly over Greenland is 3.5 • C, cooler than values inferred (4 to 5 • C) from the temperature reconstructions over this region (CAPE members, 2006). This implies that the GrIS during the LIG was likely smaller than today and represents a minimum temperature anomaly estimate. Simulated LIG warmth in Greenland is sustained under a 130-and 125-ka climate but with significant cooling by 120 ka consistent with the change in summer insolation distribution (see Fig. 2). These changes are amplified by sea ice feedbacks discussed below. However, comparisons with proxy-derived estimates of temperature at the location of the NorthGRIP ice core show a simulated summer temperature of 4.2 • C ± 1.3 • C, and an annual precipitation-weighted temperature of 3.3 • C, lower than the 5 • C estimate obtained from the ice core oxygen isotope record (NorthGRIP, 2004). Over much of the Greenland region, predicted annual precipitation rate changes throughout the LIG are small. Since the ice sheet climate coupling requires a set of GCM simulations where the GrIS is removed and replaced with bare soil we can assess the climate of the extreme scenario of an ice-free Greenland under LIG climate conditions. At the location of the NGRIP ice core, simulated maximum annual precipitation-weighted temperature anomalies relative to the pre-industrial period are in excess of 20 • C (see Fig. S4c), and the average maximum summer Greenland anomaly ranges from 14 to 16 • C for the time period 125 to 130 ka. These are clearly greater than the annual proxy palaeodata estimate of 5 • C (CAPE members, 2006), which supports the ice core evidence that the GrIS did not completely disappear during the LIG (NorthGRIP, 2004).
The increased insolation, relative to the pre-industrial period, during the early part of the LIG results in spring/summer melting of Arctic sea ice with reduced concentrations, compared with the pre-industrial period, throughout the summer months. At 130 ka, sea ice concentration is reduced by up to 40 %, compared with the preindustrial period, in the central part of the Arctic Ocean, similar to results from Otto-Bliesner et al. (2006). This reduction of summer sea ice around the margins of Greenland results in a positive sea ice albedo feedback and contributes to the observed warming in this region, particularly in the Labrador Sea. At 125 ka there is still a reduction in sea ice in the Arctic, compared with the pre-industrial period, but only up to 20 % over the majority of the region. By 120 ka the summer sea ice concentration is similar if not greater than the pre-industrial period with over 50 % sea ice present again in the vicinity of the Labrador Sea. This increase in sea ice is attributed to the cooler climate as a result of reduced summer insolation forcings toward the termination of the LIG. Although this reduction in average sea ice over the Arctic Ocean implies a significant temperature difference relative to the pre-industrial period, the inter-annual variability over the averaging period of the simulations ranges from ∼ 0 to +1 • C and, therefore, results in the regional temperature differences being statistically insignificant (see Fig. 5). Table 3. Comparison of LIG temperature anomalies (in • C) derived from palaeoproxy reconstructions (CAPE members, 2006) with the simulated maximum LIG summer temperature anomalies from HadCM3. All locations described are shown in Fig. 5. The values in parentheses for comparison with ice core data on Greenland (NorthGRIP and Renland) refer to the warmest annual precipitation-weighted temperatures. Values derived from HadCM3 include ± 2σ .

Location
Observed T Modelled T

Ice sheet model ensemble results
In order to estimate the contribution of the GrIS to LIG sea level change, we drive 500 realisations of the ice sheet model with the GCM evolving climate from 136 to 120 ka. Consequently, ice sheet geometry is predicted throughout the LIG and compared with reconstructed ice surface extent data as implied from various ice cores drilled in Greenland. The impact of ice sheet model parametric uncertainty (Stone et al., 2010) on the evolution of the GrIS through the LIG is used to derive a probability density function of the Greenland contribution to LIG sea level rise, contingent on our modelling choices. This also takes into account the mismatch between present-day observed and modelled ice sheets, most likely due to missing higher-order physical ice dynamics and the inclusion of a parameterised surface mass balance scheme.
All 500 ice sheet model simulations show contraction of the ice sheet in response to peak LIG warming. Figure 6a shows the evolution of absolute ice volume throughout the 16 000yr ice sheet simulations. Also shown is the spin-up for the modern-day GrIS for each ensemble member and subsequent spin-up, using the 136-ka climatology to give an approximation of the initial GrIS state at 136 ka. It is possible to reject a number of the GrIS LHS experiments using palaeodata from the LIG. We use the criteria of the presence of ice persisting at the Summit (Raynaud et al., 1997) ice cores on Greenland throughout the LIG and ice present at NorthGRIP up to 123 ka (NorthGRIP, 2004). The Dye-3, Camp Century and Renland ice cores are not, however, used to reject/accept simulations, as the evidence for the presence of ice there is more equivocal. In addition, simulations which make a negative contribution to sea level change are also rejected. As a result a subset of 73 simulations are selected according to this evidence from the ice core data; that is, simulations where ice is absent at the NorthGRIP and Summit ice cores are rejected. The selected simulations are shown in Fig. 6b, including a representation of their ability to reproduce the modern-day GrIS according to a skill-score (for a given set of input parameters θ ) given by  (c) is the skill-score for the simulated modernday GrIS (see Eq. 7) with the simulated modern-day ice volume also shown on the right-hand axis of (b). The star represents the modern-day observed GrIS volume . The solid black line represents the simulation with the highest skill-score for the modern-day GrIS. The dashed black line represents the average for all accepted simulations. where n is the number of grid-points, x i is the observational ice thickness at each grid-point i, f i (θ ) is the experimental ice thickness at each grid-point for each ensemble member, σ is the ice thickness Root Mean Squared Error (RMSE) of the median parameter set experiment in terms of the LHS shown in Fig. 1 and τ 2 is the observational error variance at each grid-point. The observational error is assumed to be constant across all grid-points. This skill-score for modern ice thickness measures the spatial fit over the model domain, assuming the differences between model and observation at each grid-point location are independent and normally distributed. We calculate the differences with respect to the digital elevation model derived by , interpolated to a 20-km resolution.
The ice sheet retreats in all selected cases, compared with the pre-industrial period, in response to the orbitally induced warming, with minimum ice sheet volume reached between 125 ka and 120.5 ka. All simulations show recovery towards the end of the LIG in response to the reduction in summer insolation. This is also shown by the average surface temperature anomaly (as seen in the ice sheet model) over the Greenland region which peaks at around 2 to 5 • C for the selected members of the ensemble (see Fig. 7). Maximum GrIS contribution to LIG sea level rise ranges between 0.4 and 3.8 m (Fig. 6c). None of the accepted simulations show an absence of ice in the vicinity of the Dye-3 ice core in accordance with some evidence suggesting that ice persisted through the LIG at this location (NorthGRIP, 2004;Willerslev et al., 2007). However, there is large uncertainty in the dating of basal ice at this location (Willerslev et al., 2007), which is why it is not appropriate to use this data as a direct constraint on GrIS extent. Figure 8a-c shows the GrIS geometries for parameter sets resulting in the maximum, most likely (according to the skill-score) and minimum contribution to LIG sea level change. Also shown is the respective ensemble member modern-day GrIS geometry (Fig. 8e-g). The associated precipitation and temperature forcings for the simulation with the highest skill-score, derived from HadCM3 according to the coupling methodology, are shown in Fig. S2 in Supplement. The cases where minimum ice volume and maximum temperature anomaly are reached are given and illustrate the latitudinal gradient in temperature from the enhanced insolation forcing and the change in topographic height in response to the warming. The most likely extent of the GrIS shows retreat from the northern margins, but ice is still present over central and southern Greenland (Fig. 8b). This contrasts with several previous studies (Cuffey and Marshall, 2000;Tarasov and Peltier, 2003;Lhomme et al., 2005;Otto-Bliesner et al., 2006) where ice sheet retreat is sensitive in the south but not the north. However, this sensitivity of the northern margin agrees with other more recent GrIS simulations (Fyke et al., 2011;Greve et al., 2011;Born and Nisancioglu, 2012;Quiquet et al., 2013). An isolated cap remains in the vicinity of the Camp Century and Renland ice core locations for all simulations where ice also persists in the Summit region, in agreement with evidence suggesting ice also persisted here (Johnsen et al., 2001). The drawdown of the ice surface at the Summit core location in Fig. 8a and b is ∼ 450 m and ∼ 60 m, respectively, consistent with ice core data (Raynaud et al., 1997). In contrast, Fig. 8c shows little change from the modern-day ice sheet extent with an increase of ∼ 50 m at the location of Summit.

Probabilistic assessment of GrIS contribution to the LIG highstand
It is possible to derive a probabilistic assessment of GrIS contribution to LIG sea level rise by considering the LIG palaeo-evidence of the GrIS geometry, uncertainty in ice sheet model parameterisation and the ability of the ice sheet model to reproduce the modern-day ice sheet. In this section we outline our probabilistic method followed by an assessment of the likely contribution of the GrIS to LIG sea level rise, including a sensitivity analysis to the method used.

Probabilistic method
From Bayes' theorem for a continuous distribution, the skill-score given in Eq. (7): where A is a normalising constant such that the P [Y |θ ] = 1 and the logistic function, l(θ ), accounts for the uncertainty as to where the simulated ice sheet margin lies relative to the ice core locations at the resolution of the ice sheet model domain: Y (θ ) is the maximum sea level change for each member of the ensemble, Y max is the maximum contribution to LIG sea level rise from the accepted simulations (in this case 3.8 m) and l w is the logistic width. The prior probability distribution, P [θ], weights each ensemble member according to its parameter set probability. According to Stone et al. (2010) the parameter sets can reasonably be weighted as Gaussian 2-sigma ranges such that the extreme parameter choices are penalised. Hence, we model the prior probability distribution as a multivariate Gaussian distribution: where θ j is the value of each parameter j , σ j is the standard deviation for each parameter and µ j is the mean for each parameter range (see Table 1). A comparison of the derived probability density function between Gaussian and uniform prior probability distributions indicates the choice of prior probability distribution does not have a notable affect on the outcome of the overall probability density function. Subsequently, the posterior probability distribution of the ensemble and the associated maximum LIG sea level contribution are used to construct a probability density function using a kernel density estimator (Wand and Jones, 1995;Bowman and Azzalini, 1997). A probability density function is a function that describes the relative likelihood of a variable (in this case maximum sea level change) to take on a particular given value. The probability for the variable to fall within a particular region is given by the integral of this variable's density over the region. A kernel estimator is a nonparametric way of estimating the probability density function of a particular variable and is closely related to a histogram. Unlike a histogram, a smooth kernel function rather than a discrete box is used and each of these is centred directly over each model output in order to remove the dependence of end points of bins which occurs using a histogram method (Wand and Jones, 1995). In this way the kernel estimator smoothes out the contribution of each observed data point over a local www.clim-past.net/9/621/2013/  Fig. 10a (+1.5 m). Red spots show Greenland ice core locations. neighbourhood to that data point. The kernel density estimator at any point Y ,ĝ (Y ), is of the form where n is the number of ensemble members, K is a function satisfying K (Y ) dY = 1, the kernel, whose variance is controlled by the parameter, h (usually known as the window width or smoothing parameter). K is chosen to be a unimodal probability density function that is symmetric about zero. In this case we implement a normal density function K (Y ) = 1 √ 2π e − 1 2 Y 2 . The choice of h is important since structure in the data can be lost by over-smoothing. Scott (1992) shows that the reference rule bandwidth with a normal kernel is whereσ is the sample standard deviation, in this case for maximum LIG sea level, and n is the sample number. Alternatively, we can choose a kernel width based on the modern

Sea level equivalent height (m)
Probability density x h=0 h such that x=σ h such that x=1.5σ h such that x=2σ Fig. 9. Probability density functions constructed from the 500 member ensemble of modern-day GrIS sea level equivalent height. The red line denotes the observation from Bamber et al. (2001). The distance x represents the difference between the mean of the ensemble and the observation. The grey line shows the probability density function with no smoothing. The black lines show the cases where the smoothing parameter, h, results in a probability density function where x = σ (solid), x = 1.5σ (dotted) and x = 2σ (dashed).
ice sheet volume ensemble distribution. Figure 9 shows kernel widths that result in the measured ice volume lying 1, 1.5 and 2 standard deviations away from the mean of the ensemble. In this way the smoothing parameter accounts for the additional uncertainty in the ice sheet model resulting from overestimation of the modern-day GrIS volume (see Fig. 6a).

Probabilistic results and sensitivities
From the ensemble of 500 simulations, we have derived a probabilistic assessment of the likely contribution from the GrIS to LIG sea level change (Fig. 10) with the uncertainty in the ice model parameter distributions, modern-day GrIS observations and the location of the palaeodata constraints taken into account. Although the maximum contribution from the selected ensemble is 3.8 m, Fig. 10a shows that the most likely maximum GrIS contribution to LIG sea level change is 1.5 m with a 90 % probability that the maximum contribution falls between 0.3 and 3.6 m. Figure 8d shows the Glimmer simulation ice extent which most closely matches the 1.5-m sea level contribution for the LIG derived from this probability density function. This simulation shows a similar pattern of retreat to the ensemble member with the highest skill-score: predominantly from the north and south-west (Fig. 8c). We further show that the maximum contribution range varies from a maximum of 0.2-4.7 m to a minimum of 0.5-2.4 m depending on the parameters chosen in the formulation of the density function which takes into account ice sheet model uncertainty. There is a 90 % probability of the GrIS contribution exceeding 0.6 m during the LIG and a 67 % probability of exceeding 1.3 m. However, it is unlikely (< 33 % probability) the contribution exceeded 2.2 m and very unlikely (< 10 %) that it exceeded 3.2 m (Fig. 10b). Compared with estimates of the LIG sea level highstand (Rostami et al., 2000;Muhs et al., 2002;Kopp et al., 2009;Dutton and Lambeck, 2012) exceeding 4 m, we find that sources other than the GrIS are required to account for this high sea level, such as the West Antarctic Ice Sheet (Scherer et al., 1998;Huybrechts, 2002) and/or the Canadian ice fields (Otto-Bliesner et al., 2006).
In order to assess the sensitivity of our probability density function to various uncertainties in its construction, we first examined the effect of varying the kernel width. Figure 11a shows the case where the kernel width is applied to the LIG for the optimal width (0.40 m according to Eq. 13) and the modern-day observation lying one (h = 1.50 m) and two (h = 0.75 m) standard deviations away from the modern modelled ensemble mean. Although the peak of the probability density function does not change, the upper tail is sensitive to the kernel width with a very unlikely sea level contribution exceedance ranging between 3.1 and 4.1 m. The case with the optimal kernel width assumes the anomaly in ice volume between the LIG and present day being biased in a consistent way. The alternative extreme scenario is the case where the uncertainty in the anomaly is equivalent to the model error such that the modern-day ensemble lies only one standard deviation away from the observation (h = 1.50 m). We choose a kernel width of half this width, 0.75 m, as our most plausible case, described above and shown in Fig. 10.
In addition to assessing the effect of varying the kernel width, we also varied σ (Fig. 11b), the observational error on modern-day ice thickness (τ ) (Fig. 11c) (both given as input in Eq. 7) and the width of the logistic function (Fig. 11d). Figure 11b shows that when σ is equal to zero the peak of the probability density function coincides closely with the simulation with the highest skill-score. The spread shown is a result of the kernel smoothing method used. When all simulations have equal skill (equal weighting) the probability density function shows a similar response to when σ is equal to the RMSE of the median experiment. With regard to observational error, the vertical accuracy of observational ice thickness is between 10 and 100 m Layberry and Bamber, 2001), while Bogorodskiy (1985) reports that a typical radar-sounding survey has an inherent uncertainty of about 15 m for ice depth measurements. Figure 11c shows that the observation error between 10 and 100 m makes no noticeable change to the overall probability density function. Therefore, we use a value of 15 m. Finally, in terms of the variability in the logistic function, Figure 11d indicates that the choice of the logistic width parameter does show some sensitivity for the upper tail of the probability density function. For the most plausible case a value of 0.2 m is selected.
We further tested the robustness of our skill-score on the resultant probability density function by modifying Eq. (7) such that n = 1 and using only the average ice thickness as our metric. Figure 11e shows that this makes very little difference to the probability density function.
We also assess the effect of including an additional constraint on minimum GrIS contribution to LIG sea level derived from palaeo-sea-level data. According to the global sea level estimate for the LIG derived by Kopp et al. (2009) the distribution suggests a 95 % probability that the GrIS reached a minimum at which it was at least 2.5 m of equivalent sea level smaller than today. By including this constraint we show a shift in the probability density function with a peak contribution estimate of 3.2 m closer to the estimate of many recent studies (Fig. 11f). However, this estimate is based on a probabilistic analysis of local sea level indicators where (e) Comparison of two different methods of calculating the skill-score such that s 1 (θ) takes into account the error in ice thickness at each grid-point while s 2 (θ) uses the average ice thickness of the whole ice sheet. The parameters highlighted in bold are those used for the most plausible case shown in Fig. 10a. (f) The resultant probability density function (using the default parameters) when the constraint that the GrIS contributed at least 2.5 m to global sea level (estimated from Kopp et al., 2009) is taken into account. A logistic function (with a width of 0.2 m) on the lower tail of the probability density function is applied to constrain the minimum GrIS to sea level rise at 2.5 m. The most plausible scenario is also shown by the solid line (note the different y-axis scale). the prior and posterior probability distributions are assumed to be Gaussian. The authors urge caution when interpreting Northern Hemisphere ice volume from the posterior distribution as it can only be approximated using a Gaussian because it has a definite upper bound of 7 m constrained by the amount of ice in the Northern Hemisphere available to melt today. As such, we do not use this constraint in the majority of our analysis, but we use it here as an example of what this constraint would do to our probability distribution, if future analysis did confirm this initial minimum estimate from Kopp et al. (2009). Furthermore, new results from the NEEM ice core project (NEEM community members, 2013) suggest that no more than 2-m sea level rise could come from the GrIS, lower than this minimum suggested by Kopp et al. (2009).
In order to test the robustness of the coupling methodology we performed an ensemble of simulations where only two modelled climates (with and without the GrIS) were used in the coupling method illustrated in Fig. 4. We found that although this increased the number of accepted simulations it did not result in a notable difference in the overall structure of the probabilistic distribution of GrIS contribution to LIG sea level (see Fig. S3).
Finally, the very recent NEEM ice core drilling project (NEEM community members, 2013) suggests that ice persisted throughout the LIG at this location, allowing the GrIS contribution to LIG sea level rise to be constrained further (61 accepted simulations compared with 73 when NEEM is not included) with values very likely (> 90 % probability) greater than 0.5 m but very unlikely (< 10 % probability) greater than 2.8 m (see Fig. 10c-d).

Discussion and conclusions
There are several caveats that should be discussed in the context of this study. Firstly, the uncertainty in dating basal ice limits to an extent the usefulness of this binary criterion. With the advent of improved ice cores (such as NEEM) in the future it may be possible to preferentially weight the skill toward these ice cores. In the future other aspects of the new ice cores could also be used for model evaluation, e.g. downcore temperature profiles. However, uncertainties associated with these observations are currently quite large.
Secondly, these results, of course, are somewhat limited by the absence of climate model uncertainty. We use only one model where we linearly interpolate between three possible LIG climate states. It is difficult to estimate the uncertainty in the LIG climate since there is only limited palaeodata for this time. Future work could assess the impact of structural climate model error on LIG sea level change as part of the next phase of PMIP.
Thirdly, the PDD scheme used in calculating the surface mass balance, although efficient as it only needs temperature as an input and does not require the use of regional climate models, has been shown by van de Berg et al. (2011) to significantly underestimate melt for simulations which include LIG insolation forcing compared with an approach which takes insolation and albedo explicitly into account (Robinson et al., 2010). Van de Berg et al. (2011) show that surface melt is affected not only by higher ambient temperatures but also directly through stronger summertime insolation and associated non-linear feedbacks (melting snow absorbs twice as much solar radiation as dry snow). Temperature-melt relationships assume a fixed relation between near-surface air temperature and melt rate, but this relation is also dependent on insolation and, therefore, changes in orbital forcing parameters and the latitude. In essence, the PDD scheme fails to capture north-south melt gradients driven by insolation gradients. As a result, inclusion of this process could melt the GrIS further back during the LIG. Future improvements to the PDD scheme could be to use PDD factors which are a function of insolation change.
Fourthly, our climate model simulations did not include interactive vegetation. Inclusion of this feedback could partially explain the mismatch between data and model in terms of Arctic temperature response to enhanced solar insolation. Indeed, previous work with HadCM3 has shown that vegetation feedbacks can have a significant impact on the evolution of the GrIS (Stone and Lunt, 2013). In addition, other previous modelling studies have highlighted the positive feedback from vegetation changes in response to increased solar insolation during the Holocene and LIG (e.g. Foley et al., 1994;Harrison et al., 1995).
Fifthly, and perhaps most critically, the majority of the ensemble have an associated modern ice sheet which is too large (Fig. 6a-b), a feature of many ice sheet models Ridley et al., 2005;Robinson et al., 2011). This is partly due to additional ice at the margins not captured in the ice surface extent observation , which includes only the contiguous ice sheet. In common with many other studies (Lhomme et al., 2005;Robinson et al., 2011), we assume that the predicted LIG volume anomaly with respect to the predicted modern one is more robust. This is because the overestimation of volume, which is thought to result from the lack of higher-order terms in the ice flow equations, is likely to affect both modern and LIG ice sheets in a consistent manner. The omission of basal sliding may also result in simulations being biased towards higher values for modern ice volume. However, in order to account for this potential bias, we choose a plausible probability density function that takes into account this uncertainty. The skillscore used to generate the probability density function (Eq. 7) does also ensure that the simulations which have the best representation of the modern ice sheet contribute most to the probability density function.
We observe substantial retreat of the GrIS in the north, while the ice sheet remains relatively stable in the south in contrast with many previous studies using a different forcing methodology (e.g. Cuffey and Marshall, 2000;Tarasov and E. J. Stone et al.: Quantification of the Greenland ice sheet contribution to LIG sea level rise Peltier, 2003;Lhomme et al., 2005). One fundamental difference between LIG ice sheets derived using climate forcings reconstructed from ice core records (e.g. Letréguilly et al., 1991;Cuffey and Marshall, 2000;Lhomme et al., 2005) compared with a GCM is that the forcing fails to capture changes in atmospheric circulation patterns, precipitation changes and the heterogeneity of climate trends over Greenland. Figure S4 illustrates the temperature anomaly patterns over Greenland due to insolation changes and precipitation anomalies as a result of changing orography in the GCM, clearly demonstrating the spatial heterogeneity of climate trends over Greenland. The failure to capture these processes using the ice-core-derived climate forcings is because the method uses the present-day temperature pattern which is perturbed by a spatially homogenous anomaly of temperature derived from proxy data reconstructions (e.g. the GRIP ice core record). Precipitation anomalies are simply calculated using a standard relationship where precipitation is a function of temperature. Our method is similar to Born and Nisancioglu (2012), who partly explain the preferential LIG warming and melting of northern Greenland in their results (which we also observe), but absent from most previous studies, due to the impact of larger insolation changes in the north of Greenland not adequately captured using the proxy reconstruction forcing methods. Further differences between our study and previous work include the bedrock topography used (e.g. Cuffey and Marshall, 2000;Otto-Bliesner et al., 2006), which has been previously shown to considerably affect simulated present-day ice volume (Stone et al., 2010), and the use of the PDD scheme compared with a method which takes the impact of insolation on melt into account, such as that used by Robinson et al. (2011) (see discussion above).
Our climate model, when forced with LIG insolation anomalies, shows good agreement with maximum summer warmth from LIG proxy temperature estimates in the Arctic region. We show that the GrIS contribution to LIG sea level change, consistent with ice core data, is between 0.4 m and 3.8 m. However, it is very likely that the GrIS contributed between 0.3 and 3.6 m to LIG sea level rise, lower than the range of many recent estimates of 2.5 to 4.5 m (Cuffey and Marshall, 2000;Tarasov and Peltier, 2003;Kopp et al., 2009;Robinson et al., 2011) but similar to the lower bound of Robinson et al. (2011) and the estimate of 1.6 m from Colville et al. (2011). We also show that ice persists throughout the LIG at the Dye-3 ice core for all accepted simulations, consistent with the suggestion that ice at the base of Dye-3 may predate the beginning of the LIG (Willerslev et al., 2007;Colville et al., 2011) although dating of basal ice at this location is equivocal (Willerslev et al., 2007).
In conclusion, this study emphasises the importance of including ice sheet model parametric uncertainty and palaeodata as well as modern observations, in the context of a probabilistic assessment when evaluating the impact of climate change on ice sheets. Furthermore, we show that, in order for a full probabilistic analysis to effectively take into account robust skill-scores based on simulating the modern-day GrIS, efforts should be directed at improving the existing ice sheet model physics and representation of fast-flowing processes in models used by the palaeoclimate community whilst still minimising computational costs.