Interactive comment on “ Simulation of the Greenland Ice sheet over two glacial cycles : Investigating a sub-ice shelf melt parameterisation and relative sea level forcing in an ice sheet-ice shelf model

Observational evidence, including offshore moraines and sediment cores confirm that at the 10 Last Glacial maximum (LGM) the Greenland ice sheet (GrIS) expanded to a significantly larger spatial extent than seen at present, grounding into Baffin Bay and to the continental shelf break. Given this larger spatial extent and it is close proximity to the neighboring Laurentide (LIS) and Innuitian Ice sheet (IIS), it is likely these ice sheets will have had a strong non-local influence on the spatial and temporal behaviour of the GrIS. Most previous paleo ice sheet modelling simulations recreated an ice sheet that 15 either did not extend out onto the continental shelf; or utilized a simplified marine ice parameterisation and therefore did not fully include the effect of ice shelves, and or the sensitivity of the GrIS to this nonlocal bedrock signal from the surrounding ice sheets. In this paper, we investigated the evolution of the GrIS over the two most recent glacial-interglacial cycles (240 kyr BP to present day), using the ice sheet-ice shelf model, IMAU-ICE and investigated the 20 solid earth influence of the LIS and IIS via an offline relative sea level (RSL) forcing generated by a GIA model. This RSL forcing controlled via changes in the water depth below the ice shelves, the spatial and temporal pattern of sub-ice shelf melting, which was parameterized in as a function of changes in the water depth. In the ensemble of simulations, the GrIS at the glacial maximums coalesced with the IIS to the north, 25 expanded to the continental shelf break to the south west but remained too restricted to the north east. In terms of an ice-volume equivalent sea level contribution, at the Last Interglacial (LIG) and LGM the ice sheet added 1.46m and -2.59m to the budget respectively. This estimated lowering of the sea level by the GrIS is considerably higher (~1.26 m) than most previous studies indicated whereas the contribution to the LIG highstand is lower (~0.7 m). The spatial and temporal behaviour of the northern margin was 30 highly variable in all simulations, controlled by melting below the grounded ice which was dictated by the RSL forcing and the glacial history of the IIS and LIS. In contrast, the southwestern part of the ice sheet was insensitive to these forcings, with a uniform response in an all simulations controlled by the surface air temperature (SAT) forcing, derived from ice cores.

There have been many ice sheet modelling studies of the glacial-interglacial evolution of the Northern hemisphere ice sheets (NHIS) (including the Greenland Ice sheet (GrIS) and/or Laurentide Ice Sheet (LIS) (Charbit et al., 2007;Greve et al., 1999;Helsen et al., 2013;Ritz et al., 1996;Quiquet et al., 2013), in which there was no expanse of the ice sheet beyond the present day (PD) coastline during glacial periods.The ice sheet model in these studies modelled solely the evolution of grounded ice, where the edge of the grounded ice margin was determined based upon by the flotation criterion.However, the wealth of new observational data infers that at glacial maximums the GrIS extended beyond the PD coastline, grounding out onto the continental shelf (Vasskog et al., 2015 and references therein Sect. 2).
This shows there is a mismatch between the observed extent and the modelled extent in the abovementioned models.
A review publication by Dutton et al., 2015 stated that the exact magnitude and contribution of the various global ice sheets to ice-volume equivalent sea level (ESL) during the Last Interglacial (LIG,130-115 kyr BP) is still largely unresolved.From analysis of far-field sea level records, it is estimated to have reached a peak between 6-9 m above PD.However, the precise contribution from the GrIS is poorly constrained and ranges based on ice sheet models between 0.6 and 3.5m (Cuffey and Marshall, 2000;Dutton et al., 2015;Helsen et al., 2013;Robinson et al., 2011;Stone et al., 2013), within in addition a highly variable spatial extent (Vasskog et al., 2015).Also, Clark and Tarasov, 2014 highlight that closing the Last Glacial Maximum (LGM) ESL budget is now also becoming increasing problematic.This is mostly due to the reduction in the estimated contribution from the Antarctic Ice sheet (AIS), derived from both modelling and/or observational studies.In addition, Glacial Isostatic Adjustment (GIA) modelling studies have estimated the contribution of the GrIS to the LGM ESL budget to be ~ 5m (Lecavalier et al., 2014); whereas ice sheet modelling based studies indicate significantly less, typically < 2.5 m (average < 1m) (Fyke et al., 2011;Letreguilly et al., 1991;Ritz et al., 1996).The lower values are possibly caused by restricting the glacial maximum extent to the PD coastline.Consequently, ice sheet modelling based studies which simulate a sufficiently large GrIS during glacial periods, both in terms of maximum spatial extent and total contribution to the ESL sea level budget do not yet exist and as a consequence resolving the GrIS ESL contribution over the last two glacial remains problematic.
There have been two ice sheet modelling based approaches to address the expansion of the grounded ice sheet beyond the PD margin.In the first approach, often referred to as a marine parameterisation, ice is permitted to flow and ground beyond the PD coastline to a specified 'critical water depth', regardless of the ice thickness.This critical water depth is either a function of changes in ESL or constrained by a series of masks reconstructed from observational data sets (Zweck and Huybrechts, 2005).This approach has been adopted in many ice sheet only modelling studies, solving only for grounded ice and reconstructed an extended GrIS during glacial periods (i.e (Huybrechts, 2002;Lecavalier et al., 2014;Simpson et al., 2009).However, rather than the ice sheet evolving freely, it is preconditioned to match with the observational data and does not use any physically based principles.
The second approach includes ice shelves dynamics in combination with a calculation of subice shelf melting (SSM).The sub-ice shelf melt is calculated by a parameterisation which is typically based on changes in water depth, estimated using an ESL forcing.This heuristic approach allows the ice sheet to expand onto the continental shelf but not into the open ocean.There have been a number of publications which applied the second approach using, for instance the GRISLI ice sheet model (Ritz et al., 2001).For example, Colleoni et al., 2014 parameterised the SSM as a uniform value in relation to changes in water depth to examine the growth of the NHIS during MIS7 and MIS5.During glacial periods, the reconstructed GrIS grounded across the Nares Strait, the Smith Sound and out onto the continental shelf to the NE and SW (see Fig. 7, Colleoni et al., 2014).Although in this reconstruction the ice sheet retreated from the latter two offshore regions (NE and SW) by the LIG minimum (~ 115 kyr BP), it remained grounded across the Nares Strait, which is contrary to the observational data which are reviewed in Sect. 2. Implicit in both these approaches is that the changes in paleo water depth surrounding the ice sheet are driven by the ESL forcing, generally derived from a benthic ∂ 18 O record (Lisiecki and Raymo, 2005) or by an inverse forward modelling approach (Bintanja and van de Wal, 2008).However, sea level variations are in fact not simply ESL (i.e with no spatial variation), but vary spatially due to numerous processes that dominate over different time scales (Kopp et al., 2015;Rovere et al., 2016), with GIA the dominant processes on the time scales of this study.This study will advance the second approach, using an ice sheet-ice shelf model 'IMAU-ICE' (Sect.3.1).The GrIS will be simulated over two glacial cycles (240 kyr BP to PD), focusing on the parameterisations for SSM (Sect.3.2).Secondly, to investigate the influence of the spatial and temporal variability in sea level (or water depth) on the GrIS evolution, an offline forcing derived from a GIA model (Sect.3.3) will be included.The first goal is to investigate if a larger than present day GrIS can be simulated for glacial maximum conditions, which is coherent with the observational data (Sect.2), and thereby address the current mismatch between the ice sheet model and GIA based GrIS reconstructions.Secondly, we aim to evaluate the spatial and temporal sensitivity of the GrIS to changes in the SSM and the sea level forcing.Finally, we will address the question of the GrIS contribution to ESL over the last two glacial cycles.

Observational data
There have been numerous recent publications which have reviewed the wealth of (new) observational data that can be used to constrain the spatial and temporal history of the GrIS reconstructed from ice sheet model simulations (e.g Funder et al., 2011;Vasskog et al., 2015;Cofaigh et al., 2016).It is not the aim of this study to replicate this information, rather a selection of studies are outlined below which are particular useful to constrain the ice sheet model simulations (summarised in Table 1, Fig. 1).
There are currently six Greenland ice core records (Fig. 1, white circles) that contain evidence for LIG age ice, which can be used to constrain the minimum extent that the ice sheet reached during LIG (Fig. 1).Only simulations where these six sites remained glaciated at the LIG were considered valid.
From the NEEM record (Dahl-Jensen et al., 2013) it is inferred that at 122 kyr BP, the surface elevation thinned by 130 ± 300 metres.The other five ice core sites remained ice covered, including Dye-3 (Yau et al., 2016).Additionally, analysis of Sr-Nd-Pb isotope ratios in offshore material collected from Erik Drift (Colville et al., 2011) supports a smaller than present day retreat of the southern margin and indicated that the ice sheet did not undergo complete deglaciation during the LIG.
Constraining the offshore extent at the Penultimate Glacial maximum (PGM) or earlier glaciations is complicated as the older geomorphological evidence (i.e moraines) are overridden by the subsequent readvances.As a consequence, the preservation of offshore sediments is limited.Therefore, we assumed that ice extent during the PGM equals the extent during the LGM.The aim with all simulations within the study was to reproduce a spatially expanded grounded ice sheet which reached the constraints given below during these two glaciations.
Offshore geomorphological evidence collected from numerous geophysical surveys indicate that the ice sheet grounded out onto the continental shelf (Table 1, Fig. 1), specifically to the shelf break along the SW and north and central east at the LGM.This evidence includes moraines, grounding zone wedges (Hogan et al., 2016), large scale glacial lineations (Cofaigh et al., 2004) and glacio-marine sediments dated and analysed from offshore cores (Jennings et al., 2006).Table 1 provides an overview of the asynchronous nature of the timing of retreat from the glacial maximum towards the PD margin.
Through the expanse of Petermann and Humboldt Glacier at the NW margin out into the Kane and Hall Basin and the Nares Strait (Fig. 1), the ice sheet coalesced with the Innuitian Ice sheet (IIS) at the glacial maximums (LGM or PGM).The grounded ice margin reached south into the north of Baffin Bay and out along the Arctic coastline to the north.Dating from one of the few sediment cores from the offshore NW margin (Table 1), Jennings et al., 2011 documented that the retreat of the grounded ice from the Kane and Hall Basin initiated after 10.3 kyr BP.The margin retreated in an 'unzipping manner', first from west (Kane Basin) and later to the east (Hall Basin), driven in part by the retreat of IIS back onto Ellesmere Island.The final opening of the connection between the Arctic Ocean and Baffin Bay, via the Nares Strait did not occur until after 9 kyr BP, implying that this region was one of the last regions to deglaciate.The retreat of the grounded ice sheet across the Nares Strait and back to the PD margin was a key feature which was used to constrain the simulations, and if ice remained grounded across this margin at PD, the model simulation was rejected.
Along the NE margin, Evans et al., 2009 concluded that the ice sheet advanced out onto the middle to outer shelf, covering the Westwind trough (open blue square, Fig. 1).It grounded close to (as indicated by ice-rafted debris (IRD)) but did not extend as far as the Fram Strait, limited by the continental shelf break (see Fig. 1).No dated material was recovered so the timing is unresolved.
Progressing further south, the lateral extent and timing are better constrained (Table 1) due to the greater availability of data.Retreat from the central east outer shelf initiated by ~ 15 kyr BP (blue star, Fig. 1) stabilising on the inner shelf at 10 kyr BP (green star, Fig. 1) and reaching the PD margin by 7.4 kyr BP (red star, Fig. 1) (Cofaigh et al., 2004;Evans et al., 2002).Along the SE margin, the retreat from the outer to inner shelf is highly asynchronous, with the retreat from the outer Kangerdlugssuaq trough around 17.8 kyr BP (dark blue triangle, Fig. 1) (earlier than from the central east), reaching the inner fjord by 11.8 kyr BP at Kangerdlugssuaq (open blue triangle, Fig. 1), and by 10.8 kyr BP at Sermilik (red cross, Fig. 1).It is suggested that the timing of retreat across this region is strongly influenced by the warm incursion of Irminger current (Dyke et al., 2014).
The SW region of Greenland, around Disko Bugt and Uumannaq trough is one of the more extensively studied region of the ice sheets, with a range of observational data confirming the grounding of the ice sheet to the continental shelf break (Cofaigh et al., 2013;Jennings et al., 2014;Sheldon et al., 2016;Winsor et al., 2015).The retreat from the outer shelf (cluster of red triangles, Fig. 1) between 19.3-18.6 kyr BP is inferred to have been driven by either a change in sea level and/or an ongoing gradual rise in boreal summer insolation rather than changes in ocean temperatures.The margin stabilised at the middle shelf, Hellefiske moraine (open red circle, Fig. 1) retreating at 12.24 kyr BP, reaching the inner shelf by 10.9 kyr BP.The question of whether a change in sea level could initiate such a retreat is just one aspect that the use of a RSL forcing in this study will examine.The retreat from the outer shelf edge in the vicinity of the Uumannaq fjord (cluster of green triangles, Fig. 1) was later, after 14.9 kyr BP, reaching the outer Uumannaq trough by 10.8 kyr BP (Sheldon et al., 2016).Against this background of geological evidence, we evaluated our model results as presented in Sect. 4 and later.

IMAU-ICE: ice sheet-ice shelf model.
As the aim of this study is to simulate the expansion and retreat of the GrIS onto the continental shelf, it is essential to utilise an ice sheet model which includes the possibility for ice shelves to ground and thereby the ice sheet expanding beyond the PD margin.To achieve this, we used a 3D thermomechanical ice sheet-ice shelf model IMAU-ICE (previously known as ANICE) (de Boer et al., 2014).For regions of grounded ice, IMAU-ICE uses the commonly adopted SIA approximation (Hutter, 1983) to simulate ice velocities in combination with a 3D thermodynamical approach.For regions of floating ice, the ice shelf velocities are approximated using the SSA solution (Macayeal, 1989).The model does not accurately solve for grounding line dynamics, rather the grounding line is defined as the transition between ice sheet (grounded) and ice shelf (floating) points using the flotation criterion.The complex marginal topography of Greenland, with narrow troughs with steep gradients, can lead to complications when adopting the usual SSA approach.To address this problem, a 2D one-sided surface height gradient discretisation scheme was included for ice shelf points neighbouring the grounded line.
At regions within the ice sheet where the basal temperature reaches pressure melting point, the ice sheet is allowed to slide using a Weertman-type sliding law, which relates the sliding velocity (  ), to the basal shear stress (   ) such that Where   is defined as the sliding coefficient which can be taken as inversely proportional to the bed roughness, z is the reduced normal load and p and q are spatially uniform constants over the ice sheet domain.As the roughness at the base of ice sheet is a relatively unknown quantity, a range of sliding coefficients were investigated, between 0.04 10 -10 and 1.8 10 -10 m 8 N -3 yr -1 .
First, IMAU-ICE was ran for 100 kyr, under a constant present day climate (using the input climate fields taken from the RACMO2 dataset), to reach an equilibrium state with the aim of replicating the observed present day configuration (Bamber et al., 2013).The sensitivity to the flow enhancement factor (menh, varied between 1-5) was investigated for a range of sliding coefficients (As) values.The resultant ice volume varied over the range of menh and As by ~ 0.1210 15 m 3 , with all simulations producing an underprediction of ice thickness across the centre of the ice sheet and overprediction at the narrow outlet glaciers.Based on this preliminary evaluation a value of menh=3.5 was used in all simulations.The output of these simulations is used as the initial conditions for the simulation over the two glacial-interglacial cycles.
Secondly, each simulation was ran for 240 kyr using a spatially uniform SAT forcing taken from Helsen et al., 2013 (Fig. 2a) combined with a SSM parametrisation (Sect.3.2) and sea level forcing (derived from a GIA model, Sect.3.3), to simulate the GrIS over the two glacial-interglacial cycles.As there is no GrIS SAT record that extends beyond 128 kyr BP, this SAT forcing record was produced by combining the Vostok ice core (Petit et al., 1999) with the GRIP ice core record (Johnsen et al., 2001) using the glacial-index method (Greve, 2005).We note that using a SAT forcing record derived from ice cores will not account for any spatial variability in the SAT during these two glacial-interglacial cycle.
The SMB-gradient method (Helsen et al., 2012) was applied at each time step to calculate a new SMB field resulting from this SAT forcing.In this approach, first this uniform temperature forcing (Fig. 2a) is converted into a spatially variable climate-driven surface elevation change using an atmosphere lapse rate of -7.4 K km -1 .Second, the SMB gradient fields are calculated based on a linear regression between this new surface elevation field and the mean SMB in an area with a radius of 150km.With this approach, the spatially uniform temperature forcing (Fig. 2a) can be translated in the spatially varying SMB field and ensures that the local mass balance height feedback is captured.
The resultant suite of simulations was evaluated using the observational data defined in Sect.2.

Parameterisation of sub-ice shelf melt (SSM)
As full physical based models including sub ice shelf melting are still under development, we investigated a sub-ice shelf melt (SSM) parameterisations (Fig. 2b) primarily based around the assumption (Sect.1) that for an increase in paleo water depth (or sea level) there will be a corresponding increase in the amount of SSM.Hence the SSM does not depend on temperature: temperature changes only affect the surface mass balance changes.In this method, the SSM increases with water depth (WD) by a power law relation with a constant a, and exponent m.
In order to conveniently fit this power law through two points (SSM1, WD1) and (SSM2, WD2), we solve: The range of parameter values for SSM1, SSM2 and WD1 (water depth1) and WD2 (water depth2) adopted are listed in Table 2 (Fig. 2b).

Relative Sea Level or Water Depth forcing
In this study, the output from a GIA model is incorporated into IMAU-ICE to examine the influence of spatial and temporal variability in the relative sea level (RSL) forcing via the SSM parameterisation on the expansion and retreat of the GrIS.
Sea Level (or water depth), WD(, , ) can be defined as the vertical distance between the equilibrium ocean surface, the geoid G(, , ) and the solid earth surface R(, , ) (bed topography) (Mitrovica and Milne, 2003).A change in the water depth ∆  (, , ) can result from any vertical deformation in these two surfaces, and is defined as: Where   (, , ) and   (, , ) are the vertical perturbations in the geoid and solid earth surface, at  co-latitude,  east-longitude and time t.
In most ice sheet modelling studies of the GrIS, a spatially uniform, time varying ESL forcing is used to represent the perturbation in the geoid/ocean surface (G(, , )) and the deformation of the solid earth (R(, , )) calculated using the elastic lithosphere-relaxing asthenosphere (ELRA) method (Le Meur and Huybrechts, 1996).This method only includes the local changes in the solid earth surface resulting from the deglaciation of the GrIS but in reality, the water depth/sea level signal surrounding the GrIS is both highly spatially and temporally variable because of the neighbouring LIS across North America.On the time scales of this study the main processes driving this spatial and temporal variability is GIA (Rovere et al., 2016).The spatial variability in the signal results from the interplay between the GrIS driven local changes, as is typical for near field regions and the non-local changes driven by the LIS (Lecavalier et al., 2014).This is because Greenland sits on the resulting forebulge of the LIS.Ideally, the most complete method of incorporating this complex spatial relative sea level (water depth) signal would be with a coupled ice sheet-GIA model, as in de Boer et al., (2014).Instead, a simpler alternative method was adopted in this study, by incorporating an offline coupling with the output of a GIA model.
To incorporate the output from the GIA model, first Eq.( 4) was decomposed into (a) a local (subscript L) signal, driven by changes in the GrIS (∆  , ∆  ) and (b) a non-local signal (subscript NL) driven by the influence of all other ice sheets, primarily the LIS (∆  , ∆  ).
Hence the change in water depth is: ∆  (, , ) = (∆  (, , ) + ∆  (, , )) − (∆  (, , ) + ∆  (, , )) In order to solve this relationship for water depth a GIA model was used to calculate the non-local contributions (∆  (Fig. 3d) and ∆  (Fig. 3a)).This model has three user defined input components: a reconstruction of Late Pleistocene ice sheet history (Peltier, 2004), an Earth model that simulates the solid earth deformation due to changes in the surface mass redistribution between the oceans and ice sheets (Peltier, 1974) and a model of sea-level change (Farrell and Clark, 1976).The sea-level model adopted included perturbations to the rotation vector (Milne and Mitrovica, 1998;Mitrovica et al., 2001;Mitrovica et al., 2005), time-dependent shoreline migration and an accurate treatment of sea-level change in areas characterised by ablating marine based ice (Kendall et al., 2005;Mitrovica and Milne, 2003).
To run the GIA model over the two glacial cycles (240 kyr to PD) to determine the non-local signals, an input global ice reconstruction is required which reproduces the spatial and temporal history of all global ice sheets, apart from the GrIS during this interval.As a basis for this reconstruction, the ICE5G global ice model (Peltier, 2004) was adopted, which extends from 122 kyr BP to PD; one glacial cycle.As the history for two glacial cycles was required, the ice history over the 122 kyr was duplicated to represent the previous glacial-interglacial cycles (240 kyr to 122 kyr BP), resulting in an ice sheet reconstruction from 240 kyr to PD.The GrIS component was removed from ICE5G to produce the final 'non-local' input ice history.The adopted earth model is characterised by a 96 km lithosphere, and upper and lower mantle viscosity values of 5 ×10 20 Pa s and 1×10 22 Pa s, respectively.These viscosity parameters fall approximately within the middle of the range of values commonly inferred.Using this input ice history and earth model the GIA model was ran offline to produce the non-local geoid and deformation fields (∆  (Fig. 3d) and ∆  (Fig. 3a)) As the GIA model is used to produce the non-local components, the 'local' fields (∆  , ∆  ) driven by the GrIS are still required to solve Eq. ( 5).The local driven changes in the solid earth surface, ∆  were calculated internally within IMAU-ICE, using the elastic lithosphere-relaxing asthenosphere (ELRA) method (Le Meur and Huybrechts, 1996) (Fig. 3b).This local field is combined with ∆  (from the GIA model) to calculate the total deformation of the solid earth ∆  (Eq.4,Fig. 3c) which along with the non-local geoid signal,   is used to drive the ice sheet model at each time step.
Referring back to Eq. ( 5), using this method result in the following revised equation for ∆  .
∆  (, , ) = (∆  (, , )) − (∆  (, , ) + ∆  (, , )) (6) However, comparing Eq. ( 6) to Eq. ( 5), it is evident that the local geoid, GL is not calculated using this approach and can be defined as a missing signal.To calculate this local geoid signal GL, would require solving the sea level equation (as within the GIA model) resulting from these local GrIS driven changes within IMAU-ICE, which is not possible within the adopted approach.To estimate the magnitude of this missing 'local geoid' signal the difference between the total signal (Fig. 3f, derived from Eq.5) and the signal as obtained from Eq.6 (Fig. 3e) was calculated.This difference is small (contoured on Fig. 3f), but is a shortcoming of the modelling that is accepted given the simplicity of the other components of the model.It is noted that this approach neglects for example feedbacks between changes in the sea level and the marine-based ice and the stabilizing influence this may have on the evolution of the ice shelves (de Boer et al., 2014;Gomez et al., 2010).
As the Fig. 3a illustrates, at the PGM there is a significant non-local deflection in the solid earth surface.Across Ellesmere Island and the NW Greenland, the LIS and IIS produces significant subsidence, up to 200 m with central Greenland and Baffin Bay elevated by up to 100 m and 30 m respectively, due to the influence of the forebulge.In contrast, the non-local geoid signal is much smaller, with a range of ~ 40 m (Fig. 3d).Comparing these two signals, it is apparent that it is the deflection of the solid earth surface, which will be the main contributor to driving changes in RSL (Fig. 3e and 3f).

Results of simulations
There were only nine combinations of SSM1, SSM2, WD1 and As from the ensemble of simulations that resulted in glacialinterglacial retreat over the two glacial cycles (Table 3) and fulfilled the conditions defined in Sect.2, with the spatial extent at selected time periods illustrated for one simulation on Fig. 4.
Two additional simulations are included (LowAs_lowSSM1-0.25_deep and HighAs_highSS1) which only resulted in a glacial-interglacial expanse between the LGM and PD.
At the glacial maximums (PGM and LGM) the simulated ice sheet reached the inferred observational limits along the northern and eastern margin (Fig. 4a and 4c), however at the SW margin (see red and green triangles, Fig. 4a and 4c) the ice sheet remains too restricted, possibly related to a too strong mass balance height feedback in this region.The average LGM ESL contribution is -2.59 m, which is at the upper limit relative to previous ice sheet modelling based studies (Fyke et al., 2011), but remains still ~50% smaller than estimates from GIA modelling based studies (i.e.Lecavalier et al., 2014).
Therefore, closing the LGM ESL budget remains problematic.
The simulated LIG minimum extent in all 11 simulations complied with the spatial limits inferred from the LIG ice core data; of a thinning at NEEM ( ~ 250m) and a moderate inland retreat of the SW margin, but with DYE3 remaining covered with grounded ice (Fig. 4b).The average LIG ESL contribution was 1.46 m (Table 3), which lies within the most recent range for the GrIS contribution of between 0.6-3.5 m (Dutton et al., 2015).At the PD, the SW margin has retreated too far inland (Fig. 4d and e) and there is a pronounced overthickening (up to 500 m) along most of the coastline (Fig. 4e).
Preliminary simulations concluded that increasing the resolution to 10x10km reduced this misfit but a more detailed modelling of outlet glacier at scales down to kilometres is likely to be needed to fully resolve this misfit.

5: Forcing mechanisms controlling the spatial and temporal variability during deglaciation.
There is an evident correlation between the temporal variability of SAT forcing and the total ice volume in all simulations (Fig. 5); the periods of maximum ice volume (PGM, LGM) corresponding with the minimum in SAT and vice versa (Fig. 2a).This would imply that the timings of the glacial-interglacial variations are strongly dependent on the adopted SAT forcing.However, there is a spatially variable response between the NW and SW margins which alludes that the two regions are responding regionally to a different forcing mechanism or at least a different timing of the same mechanism.Therefore, the interplay between the SAT and RSL forcing and the spatial and temporal variability in ice sheet margin will be examined in greater detail for the last 20 kyr BP.
It is evident from Fig. 5 that there is minimal variation in total ice volume and spatial extent between the nine simulations from the LGM (~ 19 kyr BP) to 14.6 kyr BP (Fig. 6a).This corresponds to a period of relatively stable SAT ~ -15°C and minimal variations in the non-local RSL forcing (either the predicted bedrock depth or sea surface height (Fig. 3)) due to only minor changes in the glacial history of the LIS (Peltier, 2004).Following this, there are three time periods (highlighted on Fig. 5) where changes in the ice volume and SAT correlate with a significant retreat/readvance along the SW, SE and to a lesser extent NE (Fig. 6), but a static NW margin.Between 14.6 kyr BP (Fig. 5) and ~13.9 kyr BP there is a rapid retreat in the grounded SW margin (Fig. 5, Fig. 6a-b) and a fall in ice volume of ~ 1.0 × 10 15 m 3 (~ 0.24m ESL).This coincides with a warming (~ 10°C (Fig. 5)) and a strong non-local RSL signal due to a significant retreat of the LIS.As the LIS deglaciated, this will induce non-local subsidence of the bedrock (Fig. 3a) across this margin, increasing the water depth and in turn the SSM.Following this retreat, there is a ~ 1 kyr stillstand in the grounded ice extent (Fig. 6c), during which there is a slow gradual cooling (Fig. 5).From ~12.9 till 11.5 kyr BP (Fig. 5, Fig. 6d) during a period of pronounced cooling (~ 15°C, Fig. 5) the ice sheet readvances along the SW, producing a small increase in total ice volume (largest in simulations with high As), with the main period of retreat commencing at 11.5 kyr BP at the onset of the sharp rise in SAT (~ 12C).This readvance (12.9-11.5 kyr BP) coincides with the ongoing large non-local RSL signal (subsidence) which in turn results in an increase in SSM.This interplay implies that changes in the SSM (driven by the RSL signal) have only a secondary influence of the dynamics of the SW margin.This is emphasised by the minimal variations in the behaviour of the SW margins between the nine final simulations.In Sect.2, from analysis of observational data, it was inferred that the retreat from this margin may, in part be driven by the changes in RSL forcing.The simulation carried out in this study suggest that this is not the case, with the retreat driven primarily responding to the SAT forcing.
The spatial and temporal behaviour of the NW margin (blue box, Fig. 1) in all 11 simulations (Table 3) is highly variable, correlating with changes in the SSM, driven by the non-local and local RSL forcing.There is minimal correlation with the timings of the SAT forcing.In all simulations, the timing of final deglaciation of this NW margin was too late (~ 10-9 kyr BP from observations), but the spatial pattern as inferred by Jennings et al., (2011), of an inwards retreat initiated at the western margin first and later to the east is replicated.This is due to the faster ice velocity within the narrow outlet fjords to the west, i.e.Humboldt glacier, which feed into the grounded ice sheet across the Kane Basin (relative to the eastern grounded margin cross the Hall Basin).The initiation of this retreat (which is at the earliest 8.9 kyr BP, Table 3) was controlled in part by the timing of the final deglaciation of the LIS within the ICE5G model (Peltier, 2004) but also by the influence of the IIS which develops within IMAU-ICE.In ICE5G the LIS retreats across Hudson Bay at 10 kyr BP with complete deglaciation of the high Arctic by ~ 8 kyr BP.This drives the onset of the non-local subsidence of the bedrock (∆  ) across this region (Fig. 3a and 3c), as the LIS forebulge collapses.It is noted that changes in the choice of earth model and/or the spatial and temporal deglaciation history of the LIS during this final deglaciation interval will of course directly impact on the timing of the GrIS retreat.
The non-local influence of the IIS (which develops across Ellesmere Island) also strongly governed the timing of the retreat of the NW margin, which can be seen by comparing the results from the two simulations: HighA s _lowSSM1-0.25 with HighA s _lowSSM1 (see Table 3).It could be assumed given the lower SSM1 (0.25 m/yr c.f 0.5 m/yr) in HighAs_lowSSM1-0.25 which results in a lower SSM close to the edge of the grounded ice margin that the onset of the retreat would be later.However, the retreat is in fact 1 kyr earlier.In this simulation (HighAs_lowSSM1-0.25) the IIS is considerably thicker (> 1500m), increasing the depression of the bedrock (due to the increased ice loading), the water depth and in turn producing a higher SSM which drives the earlier deglaciation.This highlights the influence of the IIS on controlling the deglaciation of the GrIS across this region.
The amount of basal sliding (via the choice of As) also influences the timing of the onset of this retreat: with a lower amount of basal sliding generally promoting an earlier retreat (comparing the average As (labelled AvAs) simulations to the High As simulations (labelled HighAs) on Table 3).This is examined in detail for two simulations, AvAs_lowSSM1-0.5_shallowand HighAs_lowSSM1-0.25 on Fig. 7a and b.The retreat is initiated 5 kyr earlier in the simulation with a lower As value, AvAs_lowSSM1-0.5_shallow.In part the earlier onset of the retreat with a lower As is due to more restricted and thinner grounded ice sheet across the NW margin, so there is a smaller volume of ice to retreat (compare the red contours in Fig. 7a and Fig. 7b) but also the different SSM parameters.In AvAs_lowSSM1-0.5_shallow the combination of a higher SSM1 (0.5 m/yr c.f 0.25 m/yr) and a shallower WD1 (300m c.f 475m) results in SSM that is higher at all water depths; > 0.5 m/yr and > 10 m/yr at water depths up to 300m and > 600m respectively.It is this higher SSM combined with the lower As drives the earlier onset of retreat and more restricted glacial maximum extent (Fig. 7a).
The SSM at deeper water depths (> WD1), controlled by SSM2, also strongly influences the behaviour of the NW margin via the impact on the PGM to LIG glacial history.Fig. 7c-d compares the difference in the simulated water depth between two simulations (AvAs+AvSSM1 and AvAs+AvSSM1_redSSM2) where the SSM2 is reduced by 25 m/yr (from 100 m/yr to 75 m/yr).It could be assumed given the reduction in SSM at deeper water depth, that the retreat would be later.However, the onset of retreat is 2 kyr earlier (8.9 kyr BP c.f 6.9 kyr BP).This is due to the influence of the PGM to LIG glacial history (first glacial-interglacial cycle) on the dynamics of the LGM to PD retreat.In the AvAs+AvSSM1_redSSM2 simulation, during the first advance of the ice sheet, the lower SSM at water depths > 400m results in a thicker ice sheet across the Nares Strait and eastern Ellesmere Island.This increases the bedrock subsidence and the water depth (Fig. 7c) resulting in a higher SSM surrounding the retreated ice margin during the subsequent glacial-interglacial cycle (after the LIG minimum).This higher SSM restricts the maximum spatial extent that the grounded ice margin reaches during the subsequent LGM-PD cycle (compare Fig. 7d to Fig. 7a and 7b).Therefore, with a smaller ice extent, surrounded by a region of higher SSM, this induces an earlier onset of retreat.

Conclusions
In this study using the ice sheet-ice shelf model IMAU-ICE, the evolution of the GrIS over the two most recent glacial-interglacial cycles (240 kyr BP to PD) was investigated.The sensitivity of the spatial and temporal behaviour of the ice sheet to an offline RSL forcing, generated by a GIA model was incorporated.Through this, the influence of the glacial history of the LIS and IIS was explored.This RSL forcing governed via changes in the water depth below the GrIS developing ice shelves, the spatial and temporal pattern of SSM.The SSM was parameterisation in relation to the water depth, where for an increase in water depth, the SSM increased.We note that we do not investigate the influence of these two ice sheets (LIS and IIS) on the atmospheric circulation; there was no climate model used within our study.
At the LIG minimum, all of the LIG ice cores remain ice covered, with a ~ 250 m thinning at NEEM and an inland retreat of the SW margin, contributing 1.46 m to the LIG highstand; a reduction of ~0.7m relative to previous studies.At the glacial maximums, the ice sheet has expanded offshore to coalesces with the IIS reaching the Smith Sound at the north of Baffin Bay and out onto the continental shelf along the SW but is still too restricted to the NE.A LGM ESL contribution of -2.59m is considerably higher that most previous studies (~ 1.26m), but closing the LGM ESL budget remains problematic.
The temporal response of the SW margin was primarily controlled by the adopted SAT forcing (taken from ice core records).The RSL forcing and the choice of SSM parameterisation were of secondary influence.However, the inclusion of the RSL forcing improved the reconstructed present day GrIS by reducing a underprediction along the SW margin (relative to observations).Conversely the NW margin, where the ice sheet coalesced with the IIS was relatively insensitive to the imposed SAT forcing; instead the spatial and temporal response were controlled by variation in resulted SSM patterns driven by the spatial and temporal variability in the RSL forcing and the glacial history of the LIS and IIS.RSL and temperature change combined to generate a highly variable temporal response, where optimum parameters were found to be a sliding coefficient As in the range of 1.0×10 -10 to 1.2×10 -10 m 8 N -3 yr -1 ; a relatively low SSM close to the grounded ice margin for an expanse to occur, and a higher SSM at deeper water depths to promote retreat.
There is minimal change between the extent at 5.9 and present day.Observed data constraining the timing of retreat are summarised on Table .1. Small floating ice shelves formed at the edge of the grounded ice sheet, but these are not shown. 710

Figure 1 :
Figure 1: Summary of the main place names and regions referred to in the main text and location of

Figure 4 :
Figure 4: Simulated extent of the grounded ice sheet and relative water depth using the parameters highlighted on Table 3 (AvAs+AvSSM1) at four time periods (a) Penultimate Glacial Maximum (PGM), 135

Figure 5 :
Figure 5: Comparison of the Grounded Ice Volume (10 15 m 3 ) from the nine optimum simulations on Table 3 (see Table for colours) and the Surface Air Temperature forcing (°C) (black line with crosses, Fig. 2a) from 15 to 10 kyr BP.Highlighted are the timings of the retreat-readvance-retreat of the SW margin (red box, Fig.1), the spatial pattern of which is illustrated for AvAs+AvSSM1 (light red line) on Fig. 6. Results for the NW are illustrated on Fig. 7.

Figure 6 :
Figure 6: Retreat of the grounded ice sheet along the SW margin (region bounded by the red box on Fig.1) for AvAs+AvSSM1 (kyr BP).The red contour marks the edge of the grounded ice sheet at an earlier time step

Figure 7 :Fig
Figure 7: Examples of the influence of the choice of SSM parameter and RSL forcing on the spatial and temporal retreat pattern across the NW margin (region bounded by the blue box on Fig.1).Panels (a) and (b) compare the grounded ice sheet extent in two simulations (AvAs_lowSSM1-0.5_shallow)and