Simulation of the Greenland Ice Sheet over two glacial–interglacial cycles: investigating a sub-ice- shelf melt parameterization and relative sea level forcing in an ice-sheet–ice-shelf model

Observational evidence, including offshore moraines and sediment cores, confirm that at the 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 out onto the continental shelf break. Given this larger spatial extent and its close proximity to the neighbouring Laurentide Ice Sheet (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 either did not extend out onto the continental shelf or utilized a simplified marine ice parameterization which did not fully include the effect of ice shelves or neglected the sensitivity of the GrIS to this non-local 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 ka BP to the present day) using the ice-sheet–ice-shelf model IMAU-ICE. We investigated the solid earth influence of the LIS and IIS via an offline relative sea level (RSL) forcing generated by a glacial isostatic adjustment (GIA) model. The RSL forcing governed the spatial and temporal pattern of sub-ice-shelf melting via changes in the water depth below the ice shelves.

In the ensemble of simulations, at the glacial maximums, the GrIS coalesced with the IIS to the north and expanded to the continental shelf break to the southwest but remained too restricted to the northeast. In terms of the global mean sea level contribution, at the Last Interglacial (LIG) and LGM the ice sheet added 1.46 and −2.59 m, respectively. This LGM contribution by the GrIS is considerably higher (∼ 1.26 m) than most previous studies whereas the contribution to the LIG highstand is lower (∼ 0.7 m). The spatial and temporal behaviour of the northern margin was highly variable in all simulations, controlled by the sub-ice-shelf melting 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 all simulations controlled by the surface air temperature, derived from ice cores.

1 Introduction

There have been many ice-sheet modelling studies of the glacial–interglacial evolution of the Northern Hemisphere ice sheets (NHISs) (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 expansion of the ice sheet beyond the present-day (PD) coastline during glacial periods. The ice-sheet model in these studies solely modelled the evolution of grounded ice, where the edge of the grounded ice margin was determined 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 and the modelled extents.

A review publication by Dutton et al. (2015) stated that the exact magnitude and contribution of the various global ice sheets to global mean sea level (GMSL) during the Last Interglacial (LIG; 130–115 ka BP) is still largely unresolved. From the analysis of far-field sea level records, it is estimated to have reached a peak between 6 and 9 m above PD levels. However, the contribution from the GrIS is poorly constrained and its reconstructed spatial extent highly variable (Vasskog et al., 2015). Estimates from ice-sheet modelling-based studies of the contribution to the LIG highstand range between 0.6 and 3.5 m (Cuffey and Marshall, 2000; Dutton et al., 2015; Helsen et al., 2013; Robinson et al., 2011; Stone et al., 2013). Also, Clark and Tarasov (2014) highlight that closing the Last Glacial Maximum (LGM) GMSL budget is becoming increasingly problematic. This is mostly due to the reduction in the estimated contribution from the Antarctic Ice Sheet (AIS), derived from both modelling and observational studies. In addition, glacial isostatic adjustment (GIA) modelling studies have estimated the contribution of the GrIS to the LGM GMSL budget to be ∼ 5 m (Lecavalier et al., 2014), whereas most ice-sheet modelling-based studies indicate significantly less (typically < 2.5 m; average < 1 m) (Fyke et al., 2011; Letreguilly et al., 1991; Ritz et al., 1996). These lower estimates are possibly caused by restricting the glacial maximum extent to the PD coastline. Consequently, the number of 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 GMSL budget, are limited, and as a consequence resolving the GrIS GMSL contribution over the last two glacial–interglacial cycles remains problematic. We note, however, the recent Tabone et al. (2018) study, which does address this.

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 parameterization, 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 GMSL or constrained by a series of masks reconstructed from observational datasets (Zweck and Huybrechts, 2005). This approach has been adopted in many ice-sheet 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-shelf dynamics in combination with a calculation of sub-ice-shelf melting (SSM). The sub-ice-shelf melt is calculated by a parameterization which is typically based on changes in water depth, estimated using a GMSL 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) parameterized 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. 8 in 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 ka 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 GMSL forcing, generally derived either from a benthic δ¹⁸O record (Lisiecki and Raymo, 2005) or by inverse forward modelling (Bintanja and van de Wal, 2008). However, sea level variations are in fact not simply GMSL (i.e. with no spatial variations), but vary spatially due to numerous processes that dominate over different timescales, with GIA the dominant process on the timescales of this study (Kopp et al., 2015; Rovere et al., 2016).

Table 1 Summary of a selection of the observational data which were used to constraint the timing and spatial extent of the grounded ice sheet (Fig. 1).

Figure 1 Summary of the place names and regions referred to in the main text and locations of observational data. All information describing the symbols and references for the observational data is listed in Table 1. The red and blue boxes highlight the regions shown in Figs. 7 and 8, respectively.

This study advances the second approach, using the ice-sheet–ice-shelf model IMAU-ICE (Sect. 3.1). The GrIS will be simulated over two glacial–interglacial cycles (240 ka BP to the PD), focusing on the parameterizations adopted 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 GrIS larger than that of the PD can be simulated for glacial maximum conditions, which is coherent with the observational data (Sect. 2) and thereby addresses the current mismatch between ice-sheet model and GIA-based GrIS reconstructions. Thirdly, 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 GMSL over the last two glacial–interglacial cycles.

2 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 simulated by ice-sheet models (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 were useful to constrain the ice-sheet model simulations (summarized in Table 1 and Fig. 1).

There are currently six Greenland ice core records (Fig. 1, white circles) that contain evidence of LIG age ice, and so were used to constrain the minimum extent that the ice sheet reached during the 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 ka BP, the surface elevation thinned by 130 ± 300 m. 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) infers that the southern margin retreated to a smaller extent than in the PD but 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) is overridden by the subsequent readvances. As a consequence, the preservation of offshore sediments is limited. Therefore, we assumed that the ice extent during the PGM and the LGM was equal. 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 glacial maximums.

Offshore geomorphological evidence collected from numerous geophysical surveys indicates that the ice sheet grounded out onto the continental shelf (Table 1 and Fig. 1), specifically to the shelf break along the SW, 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 this expanded glacial maximum towards the PD margin.

Through the expansion of the Petermann and Humboldt glaciers at the NW margin into the Kane and Hall basins 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 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 basins initiated after 10.3 ka BP. The margin retreated in an “unzipping manner”, first from the 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 ka 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 in the 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 ka BP (blue star, Fig. 1), stabilizing on the inner shelf at 10 ka BP (green star, Fig. 1), and reaching the PD margin by 7.4 ka 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, retreating from the outer Kangerdlugssuaq Trough at ∼ 17.8 ka BP (dark blue triangle, Fig. 1) (earlier than from the central east), reaching the inner fjord by 11.8 ka BP at Kangerdlugssuaq (open blue triangle, Fig. 1) and by 10.8 ka 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 the Irminger current (Dyke et al., 2014).

The SW region of Greenland, around Disko Bugt and the Uumannaq Trough, is one of the more extensively studied regions of the ice sheet, with a range of observational data confirming that the ice sheet grounded out onto 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 and 18.6 ka BP is inferred to have been driven by either a change in sea level and/or the ongoing gradual rise in the boreal summer insolation rather than changes in ocean temperatures. The margin stabilized at the middle shelf near the Hellefiske moraine (open red circle, Fig. 1), retreating at 12.24 ka BP and reaching the inner shelf by 10.9 ka BP. The question of whether a change in sea level could initiate such a retreat is just one aspect that the inclusion of a relative sea level (RSL) forcing in this study will address. 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 ka BP, reaching the outer Uumannaq Trough by 10.8 ka BP (Sheldon et al., 2016). Against this background of geological evidence, we evaluated our model results as presented in Sect. 4.

3 Method

3.1 IMAU-ICE: ice-sheet–ice-shelf model

As the aim of this study is to simulate the expansion onto and retreat from the continental shelf of the GrIS, it is essential to utilize an ice-sheet model which includes the possibility for ice shelves to ground and thereby for the ice sheet to expand beyond the PD margin. To achieve this, we used a 3-D 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 shallow-ice approximation (SIA) (Hutter, 1983) to simulate ice velocities in combination with a 3-D thermodynamical approach. For regions of floating ice, the ice-shelf velocities are approximated using the shallow-shelf approximation (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 2-D one-sided surface height gradient discretization scheme was included for ice-shelf points neighbouring the grounded line.

In 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 (ν_b) to the basal shear stress (τ_b) such that

$$\begin{array}{}\text{(1)}& {\mathit{\nu }}_{\mathrm{b}}=A_{\mathrm{s}}{\displaystyle \frac{{\mathit{\tau }}_{\mathrm{b}}^p}{Z^q}},\end{array}$$

where A_s 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 (A_s) were investigated between 0.04 × 10⁻¹⁰ and 1.8 × 10⁻¹⁰ m⁸ N⁻³ yr⁻¹.

Figure 2 (a) Surface air temperature forcing (SAT, ^∘C) taken from Helsen et al. (2013), with the solid red line a 2 kyr mean. (b) Grounded ice volume (10¹⁵ m³) from the ensemble of simulations (grey lines) and the nine optimum simulations (see Table 3 for colours). The solid orange line marks the present-day ice volume (Bamber et al., 2013).

Present-day input fields of the ice thickness, the surface elevation, and the bed topography are taken from Bamber et al. (2013) with input climate fields (surface mass balance (SMB), refreezing, surface air temperature (SAT)) adapted from the RACMO2 dataset (van Angelen et al., 2014). All these external datasets are interpolated and projected onto the 20 × 20 km ice model grid using the mapping software OBLIMAP2.0 (Reerink et al., 2010, 2016). The adopted OBLIMAP grid projection parameters were λ=371.5, φ=71.8, and α=7.15.

First, IMAU-ICE was run for 100 kyr, under a constant PD climate (using the input climate fields taken from the RACMO2 dataset) to reach an equilibrium state with the aim of replicating the observed PD configuration (Bamber et al., 2013). The sensitivity to the flow enhancement factor (m_enh, varied between 1 and 5) was investigated for the range of sliding coefficients (A_s). The resultant ice volume varied over the range of m_enh and A_s by ∼ 0.12 × 10¹⁵ m³, with all simulations producing an underprediction of the ice thickness across the centre of the ice sheet and overprediction at the narrow outlet glaciers. Based on this preliminary evaluation, a value of m_enh=3.5 was used in all simulations. Also, as simulations with an A_s=1.8 × 10⁻¹⁰ m⁸ N⁻³ yr⁻¹ resulted in a significant retreat of the SW margin, the range of sliding coefficients was revised to 0.04 × 10⁻¹⁰ and 1.2 × 10⁻¹⁰ m⁸ N⁻³ yr⁻¹ in subsequent simulations. The output of these simulations was used as the initial conditions for the subsequent simulations of the two glacial–interglacial cycles.

Secondly, each simulation was run for 240 kyr using a spatially uniform SAT forcing taken from Helsen et al. (2013) (Fig. 2a), combined with a SSM parameterization (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 ka 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 cycles. The SMB gradient method (Helsen et al., 2012, 2013) 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⁻¹. 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 150 km. With this approach, the spatially uniform temperature forcing (Fig. 2a) can be translated into 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.

Table 2 Range of parameters adopted in the sub-ice-shelf melting (SSM) parameterizations, illustrated in Fig. 3.

3.2 Parameterization of SSM

As full physically based models including SSM are still under development, we investigated a SSM parameterization (Fig. 3) primarily based around the assumption 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. In this method, the SSM increases with water depth (W_D) by a power law relation with a constant a and exponent m.

$$\begin{array}{}\text{(2)}& \mathrm{SSM}=a{W}_{\mathrm{D}}^m\end{array}$$

In order to conveniently fit this power law through two points (SSM1, W_D1) and (SSM2, W_D2), we solve

$$\begin{array}{}\text{(3)}& m={\displaystyle \frac{\ln \left(\frac{\mathrm{SSM}\mathrm{2}}{\mathrm{SSM}\mathrm{1}}\right)}{\ln \left(\frac{W_{\mathrm{D}}\mathrm{2}}{W_{\mathrm{D}}\mathrm{1}}\right)}};a={\displaystyle \frac{\mathrm{SSM}\mathrm{1}}{W_{\mathrm{D}}{\mathrm{2}}^m}}.\end{array}$$

The range of parameter values for SSM1, SSM2, W_D1 (water depth1), and W_D2 (water depth2) are listed in Table 2, with three examples illustrated in Fig. 3.

Figure 3 Schematic illustration of three examples of the SSM parameterization.

3.3 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 RSL forcing via the SSM parameterization on the expansion and retreat of the GrIS.

Sea level (or water depth), $W_{\mathrm{D}}\left(\mathit{\theta },\mathit{\psi },t\right)$, can be defined as the vertical distance between the equilibrium ocean surface, the geoid G(θψt), and the solid earth surface R(θψt) (bed topography) (Mitrovica and Milne, 2003). A change in the water depth $\mathrm{\Delta }{W}_{\mathrm{D}}\left(\mathit{\theta },\mathit{\psi },t\right)$ can result from any vertical deformation in these two surfaces and is defined as

$$\begin{array}{}\text{(4)}& \mathrm{\Delta }{W}_{\mathrm{D}}\left(\mathit{\theta },\mathit{\psi },t\right)=\mathrm{\Delta }{G}_{\mathrm{T}}\left(\mathit{\theta },\mathit{\psi },t\right)-\mathrm{\Delta }{R}_{\mathrm{T}}\left(\mathit{\theta },\mathit{\psi },t\right),\end{array}$$

where $G_{\mathrm{T}}\left(\mathit{\theta },\mathit{\psi },t\right)$ and $R_{\mathrm{T}}\left(\mathit{\theta },\mathit{\psi },t\right)$ are the vertical perturbations in the geoid and solid earth surface at θ co-latitude, ψ eastern longitude, and time t.

In most ice-sheet modelling studies of the GrIS, a spatially uniform, time varying GMSL forcing is used to represent the perturbation in the geoid/ocean surface (G(θψt)), and the deformation of the solid earth (R(θψt)) is 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. In reality, the water depth/sea level signal surrounding the GrIS is highly spatially and temporally variable due to the influence of the neighbouring LIS and IIS. On the timescales of this study, the main processes driving this spatial and temporal variability is GIA (Rovere et al., 2016). The variability 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 and the IIS (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 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 coupling offline the output from a GIA model into IMAU-ICE.

Figure 4 Illustration of the various components used within the calculation of the offline RSL forcing at 135 ka BP. Panels (a–c) are the predicted (solid earth deformation) bedrock depth, ΔR: (a) non-local component, ΔR_NL (No GrIS); (b) local component, ΔR_L (GrIS only); and (c) total signal. ΔR_T (d) Predicted non-local geoid signal, ΔG_NL (No GrIS). The predicted water depth (W_D) signal is illustrated for (e) ΔG_NL−ΔR_T: combination of total predicted bedrock depth and non-local geoid. It is this signal which is used within all simulations. Contoured is the local signal, ΔG_L, which is not included. (f) Total water depth signal, ΔW_D. See text for extra details.

Table 3 Set of optimum parameters which resulted in a growth beyond the PD margin during glacial maximums (PGM and LGM) and a retreat by the present day (PD). Note that W_D2 is constant in all simulations (1000 m). The simulation highlighted in grey is shown in Fig. 5. The timing of the retreat across the Nares Strait for two interglacials is given in ka BP: PGM–LIG and LGM–PD, along with the total global mean sea level (GMSL) contribution from the GrIS only for the Last Interglacial (LIG) and the Last Glacial Maximum (LGM). The GMSL was calculated using the simulated ice volume.

Figure 5 Simulated extent of the grounded ice sheet and relative water depth using the parameters highlighted in Table 3 (AvA_s+ AvSSM1) at four time periods: (a) Penultimate Glacial Maximum (PGM), 135 ka BP; (b) Last Interglacial (LIG) minimum, 123 ka BP; (c) Last Glacial Maximum (LGM), 19 ka BP; and (d) the present day (PD). Panel (e) is the difference in the present-day surface elevation (relative to the observed surface elevation; Bamber et al., 2013), where positive value indicates an overprediction. The black circles mark the location of the GrIS ice core sites. Observed data constraining the timing of retreat is summarized in Table 1. Note that the colour bar extended beyond (±) 150. Small floating ice shelves formed at the edge of the grounded ice sheet, but these are not shown.

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 (ΔG_L,ΔR_L) and (b) a non-local signal (subscript NL) driven by the influence of all other ice sheets, primarily the LIS (ΔG_NL,ΔR_NL).

Hence the relationship for the change in water depth is written as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathrm{\Delta }{W}_{\mathrm{D}}\left(\mathit{\theta },\mathit{\psi },t\right)=\left(\mathrm{\Delta }{G}_{\mathrm{L}}\left(\mathit{\theta },\mathit{\psi },t\right)+\mathrm{\Delta }{G}_{\mathrm{NL}}\left(\mathit{\theta },\mathit{\psi },t\right)\right)\\ \text{(5)}& {\displaystyle }& {\displaystyle }-(\mathrm{\Delta }{R}_{\mathrm{L}}\left(\mathit{\theta },\mathit{\psi },t\right)+\mathrm{\Delta }{R}_{\mathrm{NL}}\left(\mathit{\theta },\mathit{\psi },t\right)).\end{array}$$

In order to solve this relationship, a GIA model was used to calculate the non-local contributions (ΔG_NL (Fig. 4d) and ΔR_NL (Fig. 4a)). This model has three key 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 included perturbations to the rotation vector (Milne and Mitrovica, 1998; Mitrovica et al., 2001, 2005), time-dependent shoreline migration, and an accurate treatment of sea level change in areas characterized by ablating marine-based ice (Kendall et al., 2005; Mitrovica and Milne, 2003).

To run the GIA model over the two glacial–interglacial cycles (240 ka to the PD) to produce 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 ka BP to the PD: one glacial–interglacial cycle. As the history for two glacial–interglacial cycles was required, the ice history over the 122 kyr was duplicated to represent the previous glacial–interglacial cycles (240 ka to 122 ka BP), resulting in an ice-sheet reconstruction from 240 ka BP to the PD. The GrIS component was removed from the ICE5G global ice model to produce the final non-local input ice history. The adopted earth model is characterized by a 96 km lithosphere and an upper and lower mantle viscosity of 5 × 10²⁰ Pa s and 1 × 10²² 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 run offline to produce the non-local geoid and deformation fields (ΔG_NL (Fig. 4d) and ΔR_NL (Fig. 4a)).

Figure 6 Comparison of the grounded ice volume (10¹⁵ m³) from the nine optimum simulations (see Table 3 for colours) and the surface air temperature forcing (^∘C) (black line; Fig. 2a) from 15 to 10 ka 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 AvA_s+AvSSM1 (light red line) in Fig. 7. Results for the NW are illustrated in Fig. 8.

Figure 7 Retreat of the grounded ice sheet along the SW margin (region bounded by the red box in Fig. 1) for AvA_s+AvSSM1 (ka BP). The red contour marks the edge of the grounded ice sheet at an earlier time step (given in red): (a) 14.9 (20.9), (b) 13.9 (14.9), (c) 12.9 (13.9), (d) 11.9 (12.9), (e) 9.9 (10.9), (f) 8.9 (9.9), and (g) 5.9. There is minimal change between the extent at 5.9 and the present day. Observed data constraining the timing of retreat are summarized in Table 1. Small floating ice shelves formed at the edge of the grounded ice sheet, but these are not shown.

Figure 8 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 in Fig. 1). Panels (a) and (b) compare the grounded ice-sheet extent between two simulations AvA_s_lowSSM1-0.5_shallow and HighAs_lowSSM1-0.25 to highlight the impact of the choice of A_s in controlling the onset of the LGM-to-PD retreat. The red contour marks the edge of the grounded ice sheet at an earlier time step. Panels (c) and (d) illustrate the difference in the simulated water depth and maximum grounded ice-sheet extent (shaded grey region) between simulations AvA_s+ AvSSM1 and AvA_s+ AvSSM1_redSSM2. The main difference between these two simulations is the choice of SSM2 as 100 m yr⁻¹ (AvA_s+ AvSSM1) compared to 75 m yr⁻¹ (AvA_s+ AvSSM1_redSSM2).

As the GIA model is used to produce the non-local components, the local fields (ΔR_L, ΔG_L) driven by the GrIS are still required to solve Eq. (5). The locally driven changes in the solid earth surface, ΔR_L, were calculated internally within IMAU-ICE, using the ELRA method (Le Meur and Huybrechts, 1996) (Fig. 4b). This local field (ΔR_L) is combined with ΔR_NL(from the GIA model) to calculate the total deformation of the solid earth surface ΔR_T (Fig. 4c). This is combined with the non-local geoid signal, G_NL (Fig. 4d), to produce the final ΔW_D, which is used to force the ice-sheet model at each time step (Fig. 4e).

Referring back to Eq. (5), using this method results in the following revised equation for ΔW_D.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathrm{\Delta }{W}_{\mathrm{D}}\left(\mathit{\theta },\mathit{\psi },t\right)=\left(\mathrm{\Delta }{G}_{\mathrm{NL}}\left(\mathit{\theta },\mathit{\psi },t\right)\right)\\ \text{(6)}& {\displaystyle }& {\displaystyle }-(\mathrm{\Delta }{R}_{\mathrm{L}}\left(\mathit{\theta },\mathit{\psi },t\right)+\mathrm{\Delta }{R}_{\mathrm{NL}}\left(\mathit{\theta },\mathit{\psi },t\right))\end{array}$$

However, comparing Eqs. (6) to (5), it is evident that the local geoid, ΔG_L, is not calculated using this approach and can be defined as a missing signal. To calculate this local geoid signal, ΔG_L, would require solving the sea level equation (as within the GIA model; see de Boer et al., 2014) 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, which is calculated using the GIA model (Fig. 4f, derived from Eq. 5) and the signal as obtained from Eq. (6) (Fig. 4e) was calculated. This difference is small (contoured in Fig. 4e) 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, for example, neglects 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 Fig. 4a illustrates, at the PGM there is a significant non-local deflection in the solid earth surface. Across Ellesmere Island and NW Greenland, the LIS and IIS produce significant subsidence of up to 200 m. Central Greenland and Baffin Bay are elevated by up to 100 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. 4d). Comparing these two signals, it is apparent that the deflection of the solid earth surface will be the main contributor to driving changes in water depth/sea level in this study (Fig. 4e and f).

4 Results of simulations

There were only nine combinations of SSM1, SSM2, W_D1, and A_s from the ensemble of simulations that resulted in glacial–interglacial retreat over the two glacial–interglacial cycles (Table 3) and fulfilled the conditions defined in Sect. 2. Two additional simulations are included in Table 3 (LowA_s_lowSSM1-0.25_deep and HighA_s_highSS1) which only resulted in a glacial–interglacial expansion between the LGM and the PD, one glacial–interglacial cycle. The spatial extent at selected time periods is illustrated for one example simulation in Fig. 5.

At the glacial maximums (PGM and LGM) the simulated ice sheet reached the inferred observational limits along the northern and eastern margin (Fig. 5a and c); however, at the SW margin (see red and green triangles, Fig. 5a and c) the ice sheet remains too restricted, possibly related to too strong a mass balance height feedback in this region. The average LGM GMSL contribution is −2.59 m, which is still ∼ 50% smaller than estimates from GIA modelling-based studies (i.e. Lecavalier et al., 2014). Therefore, closing the LGM GMSL budget remains problematic.

The simulated LIG minimum extent in all nine simulations complied with the spatial limits inferred from the LIG ice core data, with a thinning at NEEM (∼ 250 m) and a moderate inland retreat of the SW margin, but with Dye-3 remaining covered with grounded ice (Fig. 5b). The average LIG GMSL contribution was 1.46 m (Table 3), which lies within the most recent estimated range of between 0.6 and 3.5 m (Dutton et al., 2015). At the PD, the SW margin has retreated too far inland (Fig. 5d and e) and there is a pronounced overthickening (up to 500 m) along most of the coastline (Fig. 5e). Preliminary simulations concluded that increasing the resolution to 10 × 10 km reduced this misfit, but a more detailed modelling of outlet glaciers at scales down to kilometres is likely 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, the periods of maximum ice volume (PGM, LGM) corresponding with the minimum in SAT and vice versa (Fig. 2). 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 implies that the two regions respond 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 these two margins is examined in greater detail for the last 20 kyr.

It is evident from Fig. 6 that there is minimal variation in total ice volume and spatial extent between the nine simulations from the LGM (∼ 19 ka BP) to 14.6 ka BP (Fig. 7a). 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 (similar to that illustrated in Fig. 4)) due to only minor changes in the glacial history of the LIS (Peltier, 2004). Following this, there are three time periods (highlighted in Fig. 6) where changes in the ice volume and SAT correlate with a significant retreat/readvance along the SW, SE, and to a lesser extent NE margins (Fig. 7), but with the NW margin remaining stable. Between 14.6 ka BP (Fig. 6) and ∼ 13.9 ka BP there is a rapid retreat in the grounded SW margin (Figs. 6, 7a–b) and a fall in ice volume of ∼ 1.0 × 10¹⁵ m³ (∼ 0.24 m GMSL). This coincides with a warming (∼ 10 ^∘C; Fig. 6) and a strong non-local RSL signal due to a significant retreat of the LIS. As the LIS deglaciated, it produced a non-local subsidence of the bedrock (Fig. 4a) 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. 7c), during which there is a slow gradual cooling (Fig. 6). From ∼ 12.9 to 11.5 ka BP (Figs. 6, 7d), during a period of pronounced cooling (∼ 15 ^∘C, Fig. 6), the ice sheet readvances along the SW margin, producing a small increase in total ice volume (largest in simulations with high A_s), with the main period of retreat commencing at 11.5 ka BP at the onset of the sharp rise in SAT (∼ 12 ^∘C). This readvance (12.9–11.5 ka 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 on the dynamics of the SW margin. This is emphasized 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 simulations carried out in this study suggest that this is not the case, with the retreat driven primarily by SAT forcing.

The spatial and temporal behaviour of the NW margin (blue box, Fig. 1) in all nine 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 the NW margin was too late compared to observations (∼ 10–9 ka BP), but the spatial pattern as inferred by Jennings et al. (2011) of a retreat initiated first at the western margin 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. in the Humboldt Glacier, which feed into the grounded ice sheet across the Kane Basin (relative to the eastern grounded margin across the Hall Basin). The initiation of this retreat (at 8.9 ka BP at the earliest, Table 3) was controlled in part by the timing of the final deglaciation of the LIS within the ICE5G global ice model (Peltier, 2004) but also by the influence of the IIS which was simulated within IMAU-ICE. In ICE5G, the LIS retreats across Hudson Bay at 10 ka BP with complete deglaciation of the high Arctic by ∼ 8 ka BP. This drives the onset of the non-local subsidence of the solid earth surface (bedrock) (ΔR_NL) across this region (Fig. 4a and c), 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 two simulations: HighA_s_lowSSM1-0.25 to HighA_s_lowSSM1 (see Table 3). It could be assumed – given the lower SSM1 (0.25 m yr⁻¹ compared to 0.5 m yr⁻¹) in HighA_s_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 (HighA_s_lowSSM1-0.25), the IIS is considerably thicker (> 1500 m), increasing the subsidence of the solid earth surface (bedrock) (due to the increased ice loading) and 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 A_s) also influences the timing of the onset of the NW margin retreat, with a lower amount of basal sliding generally promoting an earlier retreat (comparing the average A_s (labelled AvA_s) simulations to the High A_s simulations (labelled HighA_s) in Table 3). This is examined in detail for two simulations: AvA_s_lowSSM1-0.5_shallow and HighA_s_lowSSM1-0.25 (Fig. 8a and b). The retreat is initiated 5 kyr earlier in the simulation with a lower A_s value (AvA_s_lowSSM1-0.5_shallow). The earlier onset of the retreat with a lower A_s is due in part to the 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. 8a and b), and is also due to the different SSM parameters. In AvA_s_lowSSM1-0.5_shallow the combination of a higher SSM1 (0.5 m yr⁻¹ compared to 0.25 m yr⁻¹) and a shallower W_D1 (300 m compared to 475 m) results in SSM that is higher at all water depths. It is this combination of a higher SSM with the lower A_s which drives the earlier onset of retreat and more restricted glacial maximum extent (Fig. 8a).

The SSM at deeper water depths (> W_D1), controlled by SSM2, also strongly influences the behaviour of the NW margin via the impact on the PGM-to-LIG glacial history, i.e. the first glacial–integlacial cycle. Fig. 8c and d compare the difference in the simulated water depth between two simulations (AvA_s+AvSSM1 and AvA_s+AvSSM1_redSSM2) where the SSM2 is reduced by 25 m yr⁻¹ (from 100 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 ka BP compared to 6.9 ka 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 (second glacial–interglacial cycle). In the AvA_s+AvSSM1_redSSM2 simulation, during the first advance of the ice sheet, the lower SSM at water depths > 400 m results in a thicker ice sheet across the Nares Strait and eastern Ellesmere Island (part of the IIS). This increases the bedrock subsidence and the water depth (Fig. 8c) 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-to-PD glacial–interglacial cycle (compare Fig. 8d to a and b). Therefore, with a smaller ice extent, surrounded by a region of higher SSM, this induces an earlier onset of retreat.

6 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 ka BP to the 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 the spatial and temporal pattern of SSM via changes in the water depth below the ice shelves that developed around the ice sheet. The SSM was parameterized 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.7 m relative to previous studies. At the glacial maximums, the ice sheet expanded offshore to coalesce with the IIS, reaching the Smith Sound at the north of Baffin Bay as well as the continental shelf along the SW. However, it is still too restricted to the NE. A LGM GMSL contribution of −2.59 m is considerably higher than most previous studies (∼ 1.26 m), but closing the LGM GMSL 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 parameterization were of secondary influence. However, the inclusion of the RSL forcing improved the reconstructed PD GrIS by reducing an 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 was controlled by variations in the resultant SSM patterns that are driven by the variability in the RSL forcing and the glacial history of the LIS and IIS. The combined RSL and temperature changes generate a highly variable temporal response, where optimum parameters were found to be a sliding coefficient A_s in the range of 1.0 × 10⁻¹⁰ to 1.2 × 10⁻¹⁰ m⁸ N⁻³ yr⁻¹; a relatively low SSM close to the grounded ice margin to allow glacial expansion; and a higher SSM at deeper water depths to promote interglacial retreat.