The sensitivity of the Greenland Ice Sheet to glacial–interglacial oceanic forcing

Observations suggest that during the last decades the Greenland Ice Sheet (GrIS) has experienced a gradually accelerating mass loss, in part due to the observed speed-up of several of Greenland's marine-terminating glaciers. Recent studies directly attribute this to warming North Atlantic temperatures, which have triggered melting of the outlet glaciers of the GrIS, grounding-line retreat and enhanced ice discharge into the ocean, contributing to an acceleration of sea-level rise. Reconstructions suggest that the influence of the ocean has been of primary importance in the past as well. This was the case not only in interglacial periods, when warmer climates led to a rapid retreat of the GrIS to land above sea level, but also in glacial periods, when the GrIS expanded as far as the continental shelf break and was thus more directly exposed to oceanic changes. However, the GrIS response to palaeo-oceanic variations has yet to be investigated in detail from a mechanistic modelling perspective. In this work, the evolution of the GrIS over the past two glacial cycles is studied using a three-dimensional hybrid ice-sheet–shelf model. We assess the effect of the variation of oceanic temperatures on the GrIS evolution on glacial–interglacial timescales through changes in submarine melting. The results show a very high sensitivity of the GrIS to changing oceanic conditions. Oceanic forcing is found to be a primary driver of GrIS expansion in glacial times and of retreat in interglacial periods. If switched off, palaeo-atmospheric variations alone are not able to yield a reliable glacial configuration of the GrIS. This work therefore suggests that considering the ocean as an active forcing should become standard practice in palaeo-ice-sheet modelling.

1 Introduction

Recent observations show that the Greenland Ice Sheet (GrIS) has lost mass at an accelerated rate over the past decades (Rignot et al., 2011; Zwally et al., 2011; Sasgen et al., 2012; Shepherd et al., 2012; van den Broeke et al., 2016). On average, the GrIS contributed to 0.47 ± 0.23 mm a⁻¹ of sea-level rise from 1991 to 2015 (van den Broeke et al., 2016), with an accelerated rate of 0.89 ± 0.09 mm a⁻¹ from 2010 to 2014 (Yi et al., 2015). In the future, the GrIS is expected to continue losing mass, contributing to a sea-level rise relative to the 20th century between 90 and 280 mm by 2100 in the worst-case scenario (RCP8.5) (Bindschadler et al., 2013; IPCC, 2013; Clark et al., 2015; Fürst et al., 2015). This accelerated ice loss is due to a combination of increased surface melting and enhanced ice discharge from marine-terminating glaciers to the ocean (van den Broeke et al., 2016). High surface melting has been attributed to rising atmospheric Greenland temperatures (Box et al., 2009; Hall et al., 2008; Tedesco et al., 2016), which may also increase crevassing and calving at the ice front. Conversely, the recently enhanced discharge of ice into the ocean is thought to be directly connected to warmer Atlantic waters entering Greenland's fjords (Holland et al., 2008a; Rignot et al., 2010; Straneo et al., 2010; Straneo and Heimbach, 2013). Higher oceanic temperatures increase the submarine melting at the calving front of tidewater glaciers, contributing to their acceleration, ice mass discharge into the ocean and potentially grounding-line retreat. This acceleration-and-retreat mechanism has been found in several Greenland glaciers that terminate in the ocean (Rignot and Kanagaratnam, 2006). Jakobshavn Isbræ, West Greenland's fastest glacier, experienced a high rate of basal melting (Motyka et al., 2011) initially induced by the intrusion of warmer waters from the Irminger Sea (Holland et al., 2008a), more than doubling its speed in the last 25 years (Joughin et al., 2012) and suffering a rapid retreat of its terminus. Following enhanced subglacial melting observed in the early 2000s, the Helheim Glacier (southeast Greenland) also doubled its speed (Howat et al., 2005; Sutherland and Straneo, 2012) and suffered peak thinning rates of 90 m a⁻¹ (Stearns and Hamilton, 2007), with its terminus retreating by about 7 km over just 3 years (Howat et al., 2007; Straneo et al., 2016).

The complex mechanisms that lead to ice-shelf thinning, loss of buttressing and potential grounding-line instability have been studied largely for the Antarctic Ice Sheet (AIS) (DeConto and Pollard, 2016; Favier et al., 2014; Hanna et al., 2013; Joughin et al., 2014a; Pritchard et al., 2012; Rignot et al., 2004; Shepherd et al., 2004; Wouters et al., 2015). The thinning of the Larsen C ice shelf (Holland et al., 2015) and its recent calving event (Hogg and Gudmundsson, 2017; Jansen et al., 2015), the collapse of Larsen B and the melting of the Antarctic Peninsula glaciers (Cook et al., 2016), the widespread retreat of Pine Island and other glaciers in West Antarctica (Alley et al., 2015; Joughin et al., 2014b; Rignot et al., 2014) and the thinning of some East Antarctica ice shelves (Rignot et al., 2013) are notable examples of the direct connection between changes in oceanic forcing and glacier-termini adjustment (Alley et al., 2015). Only in the last several years has the scientific community also focused its attention on the ice–ocean interaction in Greenland, motivated by the observed acceleration and retreat of major GrIS outlet glaciers. Although marine-terminating glaciers cover only a small fraction of the entire GrIS, modifications at the ice–ocean boundaries due to oceanic changes may considerably affect the inland ice geometry. The effects induced by outlet-glacier acceleration are transferred onshore by ice-flow dynamics, causing adjustments to the entire inland ice-mass configuration (Nick et al., 2009; Fürst et al., 2013; Golledge et al., 2012). For this reason, a full understanding of the interaction between ice and ocean is crucial to assess the response of the GrIS to past and future climate changes.

Various numerical models have been used to simulate current submarine melt rates (Jenkins, 2011; Rignot et al., 2016; Sciascia et al., 2013; Xu et al., 2012, 2013) and dynamic retreat (Morlighem et al., 2016; Vieli and Nick, 2011) of the GrIS marine-terminating glaciers, as well as ice-dynamic future projections of the whole GrIS (Fürst et al., 2015; Nowicki et al., 2013), due to changes in the oceanic temperatures. However, how this thermal forcing affected the past GrIS configuration has not been explored from a modelling perspective so far. Recently, Bradley et al. (2017) simulated the GrIS evolution for the two last glacial cycles by considering a sub-shelf melt parameterization which is a function of the water depth below the ice shelves. Under this assumption, the submarine melt rate increases when the past sea-level rises. However, their approach does not take into account ocean temperature changes. Other studies have reconstructed the GrIS past evolution as driven essentially by atmospheric forcing (Langebroek and Nisancioglu, 2016; Quiquet et al., 2012, 2013; Robinson et al., 2011; Stone et al., 2013), while the dynamic evolution of the entire GrIS including the influence of the past oceanic forcing too has only been investigated in a simplified manner. To this end, Huybrechts (2002) used a three-dimensional ice-sheet model in which marine extent is controlled by changes in water depth based on past eustatic sea-level variations, while Tarasov and Peltier (2002), Simpson et al. (2009) and Lecavalier et al. (2014) performed a palaeo-reconstruction of the entire GrIS constraining their ice-sheet models with past relative sea-level (RSL) reconstructions. However, submarine melting was not taken into account as an active forcing in these studies.

The main purpose of this work is to assess the impact of ice–ocean interaction on the evolution of the whole GrIS throughout the last two glacial cycles. By implementing a submarine melting rate parameterization suitable for palaeo-climatic studies into a three-dimensional hybrid ice-sheet–shelf model, we evaluate the sensitivity of the GrIS to past climatic variations, including changes in oceanic temperatures (in terms of heat-flux variations), and investigate their capability to trigger grounding-line advance and retreat through time. First, we describe the ice-sheet–shelf model used to simulate the GrIS evolution, focusing on the implementation of the submarine melt rate parameterization, and the sensitivity tests performed for this study (Sect. 2). In Sect. 3, we show the results obtained in each experiment and we compare them with data for the last interglacial (LIG), the Last Glacial Maximum (LGM) and the present day (PD) found in the literature. After discussing the main model uncertainties and caveats (Sect. 4), we summarize the main conclusions of this work (Sect. 5).

2 Model description and experimental design

2.1 Model

To investigate the oceanic sensitivity of the GrIS throughout the last two glacial cycles, we use the three-dimensional, hybrid, ice-sheet–shelf model GRISLI-UCM. The model is an extension of the GRISLI ice-sheet model (Ritz et al., 2001), which has already been successfully used to simulate the evolution of the past Greenland (Quiquet et al., 2012, 2013) and Antarctic ice sheets (Ritz et al., 2001; Philippon et al., 2006; Alvarez-Solas et al., 2010), as well as the Laurentide Ice Sheet (Alvarez-Solas et al., 2011, 2013). GRISLI-UCM combines the Shallow Ice Approximation (SIA) for slow inland deformational flow and the Shallow Shelf Approximation (SSA) over fast-flowing areas, that is, ice streams and ice shelves, where plug flow is dominant. Since we assume deformational ice-sheet regions to be frozen at the bedrock, no basal sliding is considered for SIA-dominated areas. Basal sliding for ice streams is determined through a basal drag term (τ_b), defined as a function of the effective pressure (N_eff) between ice and water pressure, and the basal horizontal velocity (u_b), considered as follows:

$$\begin{array}{}\text{(1)}& {\mathit{\tau }}_{\text{b}}=-\mathit{\beta }{\mathit{u}}_{\text{b}},\end{array}$$

where

$$\begin{array}{}\text{(2)}& \mathit{\beta }\propto N_{\text{eff}}.\end{array}$$

Dragging at the floating ice-shelf base is considered to be zero. The position of the grounding line is evaluated following a flotation rule dependent on the current sea level and the ice thickness. Calving at the ice front is based on a two-condition thickness criterion (Peyaud et al., 2007; Colleoni et al., 2014). First, an ice-shelf front must have a thickness lower than 200 m to potentially contribute to ice discharge. This threshold is in agreement with the thickness of many observed shelves at their ice–ocean interface. Second, if the ice advected from each upstream point is not sufficient to maintain the ice-front thickness higher than that threshold, the grid point at the front calves. The entire GrIS ice dynamics is solved on a computational grid of 20 km × 20 km horizontal resolution and 21 vertical layers. The glacial isostatic adjustment (GIA) is described by the elastic lithosphere – relaxed asthenosphere method (Le Meur and Huybrechts, 1996), for which the viscous asthenosphere responds to the ice load with a tunable characteristic relaxation time (see Sect. 2.4).

The hybrid scheme uses a weighting function to combine the non-sliding horizontal SIA velocities (u_SIA) with the SSA horizontal velocities (u_SSA) and is defined as (Bueler and Brown, 2009)

$$\begin{array}{}\text{(3)}& \mathit{u}=(\mathrm{1}-f(u_{\text{SSA}}\left)\right){\mathit{u}}_{\text{SIA}}+f\left(u_{\text{SSA}}\right){\mathit{u}}_{\text{SSA}},\end{array}$$

where the weighting function f(u_SSA) depends on the module of the SSA component (u_SSA) through

$$\begin{array}{}\text{(4)}& f\left(u_{\text{SSA}}\right)={\displaystyle \frac{\mathrm{2}}{\mathit{\pi }}}\text{arctan}{\left({\displaystyle \frac{u_{\text{SSA}}}{u_{\text{ref}}}}\right)}^{\mathrm{2}},\end{array}$$

ranging between 0 and 1. In this work, the reference velocity u_ref is set to 100 m a⁻¹. For small values of u_SSA, f(u_SSA)≈0 and the horizontal velocities are calculated within the SIA, while f(u_SSA)≈1 for u_SSA≫u_ref, for which the contribution of SSA dominates.

2.2 Atmospheric forcing

The surface mass balance (SMB) is calculated by the positive degree-day (PDD) scheme (Reeh, 1989) forced by surface atmospheric temperatures and precipitation. This melting scheme is admittedly too simple for palaeo-simulations as it omits the contribution of insolation-induced effects on surface melting, which are important in past warmer periods such as the Eemian (Robinson and Goelzer, 2014). However, since this study focuses on the melting effects induced by past ocean temperature variations, the PDD melt model is sufficient to give a first approximation of surface melt that allows the ice sheet to retreat during interglacial periods in a realistic way. The atmospheric temperature forcing is a spatially and temporally variable field. It is retrieved using an index-anomaly approach in which the present-day climatological field (T_{clim, atm}) is perturbed by past temperature anomalies derived through a spatially uniform climatic index α(t) (Fig. 1), as follows:

$$\begin{array}{}\text{(5)}& T_{\text{atm}}\left(t\right)=T_{\text{clim, atm}}+(\mathrm{1}-\mathit{\alpha }(t\left)\right)(T_{\text{LGM, atm}}-T_{\text{PD, atm}}).\end{array}$$

The index α(t) is built through a multiproxy approach. First, we combine the temperature reconstruction for Greenland by Vinther et al. (2009) from 11.7 ka BP to present, the North Greenland Ice Core Project (NGRIP) reconstruction (Kindler et al., 2014) for 115–11.7 ka BP and the North Greenland Eemian Ice Drilling (NEEM) reconstruction (NEEM, 2013) for 135–115 ka BP, and we generate a synthetic temperature anomaly time series for 250–135 ka BP based on Antarctic isotope records following Barker et al. (2011). Second, the composite signal undergoes a windowed low-pass frequency filter ($f_{\text{c}}=\mathrm{1}/\mathrm{16}$ ka⁻¹) in order to remove the spectral components associated with millennial timescales and below. Finally, the index α is obtained by normalizing the resulting signal to be in agreement with Eq. (5), i.e. α=0 for the LGM and α=1 for the present day. The present-day climatological field is taken from the regional climate model MAR forced by ERA-Interim (Fettweis et al., 2013). T_{LGM, atm}−T_{PD, atm} is the 2-D surface atmospheric temperature (SAT) difference between the LGM and the present, as simulated by the climatic model of intermediate complexity CLIMBER-3α (Montoya and Levermann, 2008). For the precipitation rate, the procedure is similar, but the annual present-day precipitation is scaled by the ratio of the past precipitation to its present value, as in Banderas et al. (2017). At the base of the grounded ice, the melt rate is calculated as a function of the geothermal heat flux, which is prescribed following Shapiro and Ritzwoller (2004), and the local pressure melting point. The submarine melt rate is described in detail in the next subsection.

Figure 1 The 250 ka Greenland annual temperature anomaly signal built through a multiproxy approach based on the reconstruction by Vinther et al. (2009) from 11.7 ka BP to present, the NGRIP reconstruction (Kindler et al., 2014) for 115–11.7 ka BP, the NEEM reconstruction (NEEM, 2013) for 135–115 ka BP and a synthetic temperature anomaly time series for 250–135 ka BP following Barker et al. (2011) (black line). The red line shows the filtered and normalized climatic index α used to correct the present-day climatological fields when forcing the model. The same signal can be interpreted as the palaeo-oceanic temperature anomaly of Eq. (8) (in blue).

2.3 Oceanic forcing

Several marine basal melting rate parameterizations can be found in the literature. Generally, the submarine melt rate is thought to be directly influenced by the oceanic temperature variations below the ice shelves. Accordingly, most basal melting parameterizations are built as function of the difference between the oceanic temperature at the ice–ocean boundary layer and the temperature at the ice-shelf base, generally assumed to be at the freezing point. The dependence on this temperature difference can be linear (Beckmann and Goosse, 2003) or quadratic (Holland et al., 2008b; Pollard and DeConto, 2012; DeConto and Pollard, 2016; Pattyn, 2017). Because of the increasing temperature anomaly approaching the onshore ice-shelf limit, both schemes ensure a higher basal melting rate close to the grounding line, as suggested by observations (Dutrieux et al., 2013; Rignot and Jacobs, 2002; Wilson et al., 2017).

The marine basal melting rate parameterization used in this work follows a linear approach that accounts separately for sub-ice-shelf areas near the grounding line and for purely floating ice (ice shelves). A linear scheme is the simplest case that allows testing of the GrIS sensitivity to past oceanic temperature changes. The formulation is derived from the net basal melt rate B_gl (m a⁻¹) for ice-shelf cavities close to the grounding line and terminating in shallow ocean zones, expressed by Beckmann and Goosse (2003) as

$$\begin{array}{}\text{(6)}& B_{\text{gl}}\left(t\right)=\mathit{\kappa }\left(T_{\text{ocn}}\right(t)-T_{\text{f}}),\end{array}$$

where T_ocn is the oceanic temperature close to the grounding line (K), T_f is the temperature at the ice base, assumed to be at the freezing point (K), and κ is the heat flux exchanged between ocean water and ice at the ice–ocean interface ($\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$). Since knowledge of past T_ocn and T_f is challenging for the complex heat-flux transfer between ice shelves and the surrounding water, we opted for substituting these quantities and rearranging the equations to make them more suitable for palaeo-studies. A representation of the transient oceanic temperature T_ocn can be given by the climatological oceanic temperature T_{clim, ocn} corrected by the LGM–present temperature anomaly (T_{LGM, ocn}−T_{PD, ocn}) scaled by the same climatic index α=α(t) used to correct the atmospheric climatological fields (Fig. 1). Under this assumption, Eq. (6) can be rewritten as

$$\begin{array}{}\text{(7)}& B_{\text{gl}}\left(t\right)=\mathit{\kappa }(T_{\text{clim, ocn}}+\mathrm{\Delta }{T}_{\text{ocn}}(t)-T_{\text{f}}),\end{array}$$

where

$$\begin{array}{}\text{(8)}& \mathrm{\Delta }{T}_{\text{ocn}}\left(t\right)=(\mathrm{1}-\mathit{\alpha }(t\left)\right)(T_{\text{LGM, ocn}}-T_{\text{PD, ocn}}).\end{array}$$

Combining and reorganizing these equations as

$$\begin{array}{}\text{(9)}& B_{\text{gl}}\left(t\right)=\mathit{\kappa }(T_{\text{clim, ocn}}-T_{\text{f}})+\mathit{\kappa }\mathrm{\Delta }{T}_{\text{ocn}}\left(t\right),\end{array}$$

we can finally retrieve the expression for the basal melting rate at the grounding line B_gl (m a⁻¹) as used in this work:

$$\begin{array}{}\text{(10)}& B_{\text{gl}}\left(t\right)=B_{\text{ref}}+\mathit{\kappa }\mathrm{\Delta }{T}_{\text{ocn}}\left(t\right).\end{array}$$

B_ref (m a⁻¹) $=\mathit{\kappa }(T_{\text{clim, ocn}}-T_{\text{f}})$ is assumed to represent the present-day basal melting rate around the ice sheet, κ represents the sensitivity of the basal melting rate to changes in the oceanic temperature, and ΔT_ocn (K) expresses the temporal evolution of the melting at the ice base. In this way, B_gl coincides with the present-day melt (B_ref) for α=1 and its LGM (21 ka BP) value for α=0. When B_gl is negative, the model allows for refreezing, and the grounding line can advance offshore.

In a more realistic setup, all parameters in Eq. (10) could be described by 2-D spatially variable fields. However, for the sake of simplicity, here we considered all the parameters to be spatially uniform around all the GrIS marine borders, as described in Sect. 2.4. The glacial–interglacial temperature anomaly T_{LGM, ocn}−T_{PD, ocn} (Eq. 8) is set constant to −3 K, which corresponds to the mean value of the reconstructed LGM sea surface temperature (SST) anomalies for the Atlantic Ocean between 60 and 80^∘ N latitude (MARGO Project Members, 2009). This value slightly differs from the LGM mean SST anomaly reconstructed by Annan and Hargreaves (2013) (between −1 and −2 K). However, a variation in κ or an identical change in ΔT_ocn equally affect the oceanic forcing applied to the model. Therefore, considering a different value for ΔT_ocn would not alter the magnitude of the oceanic sensitivity applied to the GrIS. These simplifications allow here for a spatially uniform, but time-dependent, B_gl.

Table 1 Summary of all parameter values used to perturb the basal melting rate equation (Eq. 10) in each sensitivity test.

The basal melting rate for purely floating ice shelves (B_sh) is given by the grounding-line basal melt B_gl scaled by a constant factor γ:

$$\begin{array}{}\text{(11)}& B_{\text{sh}}\left(t\right)=\mathit{\gamma }{B}_{\text{gl}}\left(t\right).\end{array}$$

In this study, γ is set to 0.1. Hence, we consider that the basal melting rate for ice shelves is 10 times lower than that close to the grounding zone, which is qualitatively in agreement with melt rates observed in some Greenland glaciers (Münchow et al., 2014; Rignot and Steffen, 2008; Wilson et al., 2017). Conversely, the melt rate in the open ocean, that is considered as being beyond the continental shelf break, is prescribed to a high value (50 m a⁻¹) to avoid unrealistic ice growth beyond 1500 m of ocean depth.

2.4 Experimental design

To study how oceanic changes impact the evolution of the GrIS over the last glacial cycles, we performed a set of sensitivity tests by perturbing the two key parameters of the basal melting rate equation (Eq. 10): the estimated present-day submarine melting B_ref and the heat-flux coefficient κ. For each experiment, we ran an ensemble of simulations over the GrIS domain throughout the last 250 ka. In this study, the model is initialised with the present-day Greenland topography (Bamber et al., 2013), the characteristic relaxation time for the lithosphere is set to 3000 years, and the model is forced by the past relative sea-level reconstruction of Grant et al. (2014). The first ∼ 100 ka of the simulation are considered as a spin-up and are not analysed. A summary of all the parameter values used in each sensitivity test is shown in Table 1.

First, we analyse the sensitivity of the model to different constant (in space and time) B_ref values applied at the base of the ice-sheet marine margins. Due to the scarcity of submarine melt observations along the GrIS coasts, and since the only available estimates have focused on few very rapid tidewater Greenland glaciers that cannot be representative of the basal melt rate for the entirety of GrIS marine areas (Rignot et al., 2010; Motyka et al., 2011; Straneo et al., 2012; Xu et al., 2013; Enderlin and Howat, 2013; Fried et al., 2015; Rignot et al., 2016; Wilson et al., 2017), we assume present-day basal melting rates for Greenland comparable to those from Antarctic ice shelves (Rignot et al., 2013). The range of values of B_ref is set between 0 and 40 m a⁻¹, while κ is set to zero to make the ocean contribution constant in time. The resulting basal melting rate is thus equal to the tested B_ref value and a condition of no oceanic basal melting around the GrIS is achieved only when both B_ref and κ are set to zero.

Second, we study the sensitivity of the GrIS to the basal melt rate sensitivity κ at the ice–ocean interface. The range of tested values for κ is between 0 (expressing a temporally constant basal melting rate) and 10 $\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$. The choice of this range reflects the inference made in Antarctica by Rignot and Jacobs (2002) that a variation of 1 K in the effective oceanic temperature changes the melt rate by 10 m a⁻¹ (Eq. 6). Due to the lack of data for Greenland, as a first approximation, we can assume such a value is also realistic there. This is surely a simplification of the problem, as the relation between ocean temperature and melt rate is not universal but depends on many factors, such as the water salinity, the depth, the conformation of the cavity, the water velocity below the ice shelf and subglacial discharge. The sensitivity test for κ is firstly done for B_ref=1 m a⁻¹ and then for other B_ref values to show that the GrIS response to the melting rate sensitivity κ depends on the chosen reference basal melting rate (see Table 1).

Figure 2 Time evolution of (a) grounded ice volume (million km³) and (b) ice-covered area (million km²) simulated for different values of B_ref (m a⁻¹) having set κ=0 (see Table 1). Dashed lines show the present-day estimated volume and area of the GrIS (Bamber et al., 2013).

Figure 3 Glacial maximum GrIS surface elevation (km) simulated at Termination II (a–c) and Termination I (d–f) for different values of the reference basal melting rate (B_ref=0, 5, 10 m a⁻¹) under constant oceanic conditions (κ=0). The timing at which the ice volume reaches its maximum value during a glacial cycle depends on the experiment and is stated in black for each snapshot. Blue lines indicate the GrIS extension at the following peak of deglaciation with its corresponding timing reported in blue. Red zones represent the ice shelves extending beyond the glacial maximum grounding line (black line). Black circles indicate the locations of the Camp Century, NEEM, NGRIP, GRIP and Dye3 ice cores (from north to south).

3 Results

In this section, we present the results of each sensitivity study aiming to assess the impact of the ocean on the evolution of the GrIS throughout the last two glacial cycles, especially focusing on the LIG, the LGM and the PD GrIS. The present work involved a total of 110 model simulations, although only the most representative cases for each sensitivity study are discussed.

3.1 Sensitivity to the reference submarine melting

In this experiment, the maximum ice volume reached in glacial times ranges between 3.4 and 4.3 million km³ (Fig. 2), 15–45 % higher than the observed current value (Bamber et al., 2013), suggesting that under constant oceanic forcing, the GrIS is limited to a configuration close to that of the present day (Fig. 3). The highest glacial ice volume is reached by imposing a null basal melting to the GrIS margins (B_ref=0), which corresponds to a simulation forced solely by palaeo-atmospheric variations. The varying SMB throughout the cycles still results in a changing GrIS ice volume over time. However, during glacials, most grounded ice remains on land above sea level, and only small ice shelves are able to grow (Fig. 3a, d).

For B_ref>0, a positive basal melt rate is applied to the marine margins of the whole GrIS throughout the two glacial cycles. The submarine melting not only inhibits the grounding-line advance during the glacials but contributes to thinning the few marine-terminating glaciers still present, constraining the grounding line further inland and resulting in a GrIS extent close to the observed present-day configuration (Fig. 3b, c, e, f). This mechanism can still be quite active during glacial times, such that the ice volume can be even lower than that simulated at the present (Fig. 2). Note that the ice volume is more sensitive to B_ref during the glacial periods, as during the interglacial periods the effect of the ocean is limited by the topography of the Greenland itself. Thus, the retreat is almost entirely driven by the surface ablation and the elevation–melt feedback. For low B_ref values, the ice lost in a deglaciation is to a large extent determined by the GrIS configuration in the preceding glacial. As high basal melting rates inhibit the ice growth during the cold phase, the higher the B_ref is applied to the marine margins, the lower ice loss is simulated in the following interglacial (Fig. 4). However, for B_ref=5 m a⁻¹, ice loss becomes insensitive to the melting applied, since the GrIS is also totally land based during glacial periods and any subsequent ice mass loss is therefore uniquely driven by ablation (compare Fig. 3a and b or Fig. 3d and e).

Figure 4 Distribution of the ice volume (a) and area (b) lost during the LIG (triangles) and the Holocene (circles) as a function of B_ref. Grey and yellow shades show the range between the maximum glacial and the minimum interglacial ice volumes (area) for LIG and Holocene, respectively. The loss is calculated between the time at which the ice volume reaches its maximum value simulated before the deglaciation (between 132 and 128 ka BP for TII and between 13 and 9 ka BP for TI) and its following ice minimum (between 122 and 121 ka BP for the Eemian and between 8 and 0 ka BP for the Holocene). The colours of the points follow the legend of Fig. 2 for clarity. Each ice volume loss has been converted to value of sea-level-equivalent anomaly (m s.l.e. with respect to its simulated present-day volume.

3.2 Sensitivity to the heat-flux coefficient

We next study the sensitivity to the ocean for a fixed B_ref value of 1 m a⁻¹ (Fig. 5). This value is within present-day submarine melting rates estimates, between those found in the largest remaining outlet glaciers in Greenland (Wilson et al., 2017) and those of smaller marine-terminating glaciers with presumably much lower ocean-induced melt. Under this assumption, the maximum ice volume simulated in both glacial periods for different κ values ranges between 4 and 5.4 million km³, greatly exceeding the range found for the case with constant oceanic forcing. Prescribing positive or zero uniform submarine melting to the marine boundaries limits the glacial expansion of the GrIS (Figs. 6a and 7a), as discussed in Sect. 3.1. Conversely, by intensifying the oceanic forcing applied to the margins (with increasing values of κ), the glacial ice volume increases. For κ=1 $\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$, the model simulates a GrIS glacial expansion to the continental shelf break in which the grounding line has already advanced from the present-day continental boundaries, and large ice shelves are generated in the eastern GrIS, especially in the northeast (Figs. 6b and 7b). The maximum expansion is simulated for the last glaciation, where the grounding line has almost reached the continental shelf break and large ice shelves in the east cover the remaining shallower zones of the bathymetry. For κ=10 $\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$, the GrIS extends all the way to the continental shelf break at its glacial maximum, while only a few small floating ice shelves are present (Figs. 6c and 7c).

Figure 5 (a) Time evolution of GrIS grounded ice volume (million km³) and (b) ice area (million km²) simulated for different values of the heat-flux coefficient κ, having set B_ref=1 m a⁻¹. Dashed lines shows the GrIS ice volume and area estimated for the present day (Bamber et al., 2013).

Figure 6 GrIS surface elevation (km) simulated at the penultimate glacial maximum (TII) (a–c) and at the LIG minimum (Eemian) (d–f) for three values of the melting rate sensitivity κ having set B_ref=1 m a⁻¹. The timing of these snapshots depends on the experiment and is stated in black for each snapshot. Corresponding ice volume (in s.l.e.) is shown in blue. Red zones represent the ice shelves extending beyond the glacial maximum grounding line (black line). Black circles indicate the locations of the Camp Century, NEEM, NGRIP, GRIP and Dye3 ice cores (from north to south).

Figure 7 GrIS surface elevation (km) simulated at the LGM (TI) (a–c) and present-day GrIS (d–f) for three values of the heat-flux coefficient κ having set B_ref=1 m a⁻¹. The timing of the snapshots depends on the experiment and is stated in black for each snapshot. Corresponding ice volume (in s.l.e.) is shown in blue. Red zones represent the ice shelves extending beyond the LGM grounding line (black line). Black circles indicate the locations of the Camp Century, NEEM, NGRIP, GRIP and Dye3 ice cores (from north to south).

A larger ice sheet loses more ice during a deglaciation, leading to an interglacial state that is almost independent of κ (Fig. 5). This response is related to the saturation of the oceanic forcing in warm peaks, when the GrIS is almost totally land-based and the ice loss is hence mostly due to the increase in atmospheric temperature and precipitation. Since glacial accretion affects the ice growth much more than basal melting during the retreat, the ice loss during a deglaciation monotonically increases with increasing κ (Fig. 8). Thus, for larger κ values, more ice grows during glacial periods and more ice is lost, and faster, during the subsequent deglaciation. Mass loss is mostly due to the large number of grounded-ice zones that are converted into ice-free areas during the deglaciation (Fig. 8b). The percentage of grounded points lost until the peak of an interglacial period saturates for κ above 3 $\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ in correspondence with preceding glacial GrIS configurations which present a grounding-line expansion to the continental shelf break. The slightly increasing ice loss still observed for higher oceanic sensitivities is mostly related to the ice lost in the GrIS interiors due to the positive elevation–melt feedback.

Figure 8 Distribution of ice volume (a) and area (b) lost during the LIG and the Holocene as a function of κ, for B_ref=1 m a⁻¹. Grey and yellow shades show the deviation between the maximum and the minimum ice volumes (area) for LIG and Holocene, respectively (see Fig. 5). The loss is calculated between the time at which the ice volume reaches its maximum value simulated before deglaciation (between 140 and 128 ka BP for TII and between 19 and 10 ka BP for TI) and the subsequent ice minimum (between 122 and 121 ka BP for the Eemian and between 8 and 0 ka BP for the Holocene). The colours of the points follow the legend of Fig. 5 for clarity.

Due to our melting parameterization (Eq. 10) and to the B_ref value chosen, water below the ice shelves is allowed to freeze for κ>0.5 $\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$, favouring ice growth and GrIS expansion (Fig. 5). Below this threshold, the model still allows for submarine melting rates across the margins in glacial times and the GrIS expansion is almost totally driven by surface accumulation. However, the sensitivity with respect to κ strictly depends on the value of B_ref, as it defines the positive threshold that the glacial GrIS has to overcome to start reacting to the oceanic forcing imposed at the margins (Fig. 9). For B_ref=10 m a⁻¹, the GrIS responds to the ocean only for κ>3 $\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$, while for B_ref=30 m a⁻¹ the GrIS starts to expand only for κ>8 $\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$. For high B_ref, since a constant high submarine melting is applied overall, the glacial GrIS is almost constrained to the PD configuration and exposure to the ocean is reduced. Only a sufficiently high κ to counteract this strong melting is able to make the GrIS expand and then retreat during the interglacial. Once the reaction has started, the sensitivity of the GrIS to κ increases with increasing B_ref; i.e. small variations in the magnitude of κ lead to a fast and large growth of ice during glacials and consequently to a fast and large loss of ice during the deglaciation. Similar results are found for the LIG (not shown).

3.3 Last interglacial

The amount of ice lost during the LIG period increases with the oceanic sensitivity κ. High κ values lead to higher glacial ice volumes and to larger ice losses during the consequent deglaciation (Fig. 8). The range of observed volume changes spans between 4.2 m s.l.e. (for κ=0) and 8 m s.l.e. (for κ=10 $\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$), above the present-day GrIS ice volume. Despite this large ice-loss range, all GrIS configurations simulated at the LIG ice minimum (Eemian) present a similar extension (Fig. 6d–f). In all experiments, a large retreat is observed in the north (especially in the northeast), where melting overcomes the low accumulation rates, and in the southwest, where the ice discharge from the interior is enhanced by the presence of fast ice streams and, in some areas, by the fact that the bedrock is below sea level. Although the position of the land-ice borders at the Eemian is not very sensitive to κ, the corresponding surface elevation fields show some differences depending on κ. For high values of κ, a lower ice elevation is simulated over the GrIS (compare Fig. 6d and f), a tendency that is reflected in a slightly lower ice volume too (Fig. 5).

It is interesting to note that even when imposing a very high κ, the complete disappearance of the GrIS is not simulated. The GrIS is only partly deglaciated and all ice-core sites are still covered by ice (including the discussed ice core locations of Dye3 and NEEM). Since the oceanic-driven retreat is limited by the land-based configuration observed in the interglacials, the retreat during the LIG is mainly controlled by the atmospheric temperatures and precipitations with which the model is forced.

Figure 9 Distribution of the ice volume lost in the Holocene as a function of the heat-flux coefficient κ, simulated for three selected reference basal melting rates (B_ref=1 m a⁻¹ in green, B_ref=10 m a⁻¹ in blue and B_ref=30 m a⁻¹ in red). The ice volume loss is calculated between the time at which the ice volume reaches its maximum value before the deglaciation and the present day. The green points are the same as the circles of Fig. 8a (for the Holocene).

The amount of ice lost during the Eemian relative to the present day (Fig. 10), which ranges between 2.9 and 3.2 m s.l.e., is within the uncertainty range of ice volumes suggested by some previous studies (e.g. 1.2–3.5 m s.l.e. for Helsen et al. (2013), 0.4–4.4 m s.l.e. for Robinson et al. (2011) and 0.4–3.8 m s.l.e. for Stone et al., 2013). Also, the timing at which the peak of deglaciation occurs, which spans between 122.3 and 121.6 ka BP in all the simulations, agrees with the timing proposed in many previous studies (Calov et al., 2015; Langebroek and Nisancioglu, 2016; Robinson et al., 2011; Stone et al., 2013; Yau et al., 2016). The time at which the peak of the Eemian occurs in our experiments depends partly on the timing of the atmospheric temperature peak and partly on the duration of the post-glacial rebound, which controls the intrusion of warm waters into the GrIS bays enhancing the ocean-driven retreat. However, the Eemian peak does not depend on the maximum insolation since the PDD scheme used does not account for past insolation changes.

Figure 10 GrIS ice volume evolution simulated for different values of the melting rate sensitivity κ during the last interglacial (see Fig. 5 for the line colour legend). The ice volumes have been converted to values of s.l.e. anomaly with respect to the present-day volumes estimated in each specific simulation. Grey shading represents the reference basal melting rates B_ref investigated for the case of constant-in-time oceanic forcing (κ=0 $\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$). The upper bound refers to B_ref=40 m a⁻¹ and the lower bound to B_ref=0 m a⁻¹. Black and white symbols indicate the LIG minimum ice volumes estimated by previous studies. The tight clustering of our estimates compared to previous work is due to the fact that the sole uncertainty is here related to the oceanic forcing through κ.

3.4 Last Glacial Maximum

Although many uncertainties about the GrIS configuration during the last glacial period still exist, several estimates of the sea-level contribution from the GrIS during the last deglaciation can be found in the literature: 2.6 m s.l.e. (Bradley et al., 2017), 2.7 m s.l.e. (Huybrechts, 2002), between 2 and 3 m s.l.e. (Clark and Mix, 2002), 3.1 m s.l.e. (Fleming and Lambeck, 2004), 4.1 m s.l.e. (Simpson et al., 2009), between 3.1 and 4.5 m s.l.e. (Buizert et al., 2018) and 4.7 m s.l.e. (Lecavalier et al., 2014). These estimates come from ice-sheet models of different complexity, with their own dynamics and boundary conditions. Particularly the ice-sheet model used by Simpson et al. (2009) and Lecavalier et al. (2014) is run in combination with a GIA and RSL model and then constrained by past surface elevations derived from ice-core data, observations of past changes in RSL and the present-day GrIS configuration. These models do not solve the dynamics of the ice shelves or the grounding-line migration, which is parameterized. However, their estimates of the GrIS spatial extent can be considered as the most realistic reconstructions of the recent past glacial GrIS so far.

Under constant oceanic conditions, the LGM-PD ice excess simulated by our model at 21 kyr BP spans between 0 and 1.4 m s.l.e. for B_ref ranging from 0 to 40 m a⁻¹, increasing with decreasing B_ref values (grey shaded region – Fig. 11). This range is well below previous LGM ice volume reconstructions found in the literature (grey points). However, slightly larger ice volumes (0.6–2 m s.l.e.) are found at the peak simulated further in time in the glaciation (∼ 13–10 ka BP). For the case with no submarine melting (B_ref=0 m a⁻¹), the maximum ice volume (lower bound of grey shadow, at ∼ 12 kyr BP) is close to those of Huybrechts (2002) and Bradley et al. (2017). In this simulation, the GrIS increases moderately as its extension surpasses its PD borders and the grounding line approaches the continental shelf (Fig. 12a). Nevertheless, the atmospheric forcing alone is not sufficient to make the GrIS expand as expected during the LGM. According to reconstructions, the GrIS extended as far as the continental shelf break in every direction, except in the northeast region where the grounding line remains closer to the coast (Lecavalier et al., 2014). In our simulations, the GrIS reaches a glacial expansion consistent with the literature only for κ≥1 $\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ (Fig. 12b). However, the ice volume reached for this oceanic sensitivity is still smaller than the LGM volumes of Simpson et al. (2009) and Lecavalier et al. (2014) (Fig. 11), since only with κ>3 $\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ does the model simulate a maximum ice volume comparable to those ranges. The discrepancy in volumes, despite the same extension, could be related to the different dynamics and boundary conditions applied in the two models. Nevertheless, our simulated ice volumes are in agreement with recent estimates corrected for seasonal surface air temperatures in Greenland during the LGM (Buizert et al., 2018).

Figure 11 GrIS ice volume evolution simulated for different values of the melting rate sensitivity κ during the last deglaciation. The ice volumes have been converted to values of s.l.e. anomaly with respect to the present-day volumes estimated in each specific simulation. As in Fig. 10, grey shading represents the simulations for the different reference basal melting rates (B_ref) investigated for the case of constant-in-time oceanic forcing (κ=0). The upper bound refers to B_ref=40 m a⁻¹ and the lower bound to B_ref=0 m a⁻¹. Grey dots and orange shading indicate estimates of the GrIS ice volume at the LGM (21 ka BP) and at the maximum ice volume reached before the last deglaciation (16.5 ka BP), as suggested by previous work.

The timing of the reconstructed deglaciation can also provide information for comparison. The maximum increases suggested by Simpson et al. (2009) (4.6 m s.l.e.) and Lecavalier et al. (2014) (5.1 m s.l.e.) occur at 16.5 ka BP, while our simulations suggest a timing dependent on κ ranging from 20 to 10 ka BP for very low κ values (Fig. 11). The magnitudes of the oceanic sensitivity that best approximate the evolution of the GrIS before the Holocene are thus between 5 $\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ (4.6 m s.l.e. at 17.4 ka BP) and 10 $\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ (5.3 m s.l.e. at 14 ka BP). However, some discrepancies between our GrIS glacial extension and that of Lecavalier et al. (2014) are still present (Fig. 12b).

Figure 12 GrIS total extent (ice shelves are included) simulated at the peak of the last glaciation for (a) no melting/freezing at the grounding line (cyan) and (b) κ=0, 1 and 10 $\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ (red, green and purple lines, respectively) for B_ref=1 m a⁻¹. The timing of the glacial maximums is (a) 12 kyr BP and (b) 10, 20 and 14 kyr BP for κ=0, 1 and 10 m a⁻¹, respectively. LGM (21 kyr BP) GrIS grounding-line position estimated by Lecavalier et al. (2014) is shown for comparison (black line).

3.5 Present-day GrIS

Given that the topography of the present-day GrIS is one of the trustworthy measures used to assess the reliability of an ice-sheet model, we compare our present-day GrIS ice thickness and extent simulated for κ=10 $\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ to those estimated by Bamber et al. (2013) (Fig. 13). The choice of this particular κ value is based on the discussion above (Sect. 3.4) for the LGM and is supported by the good agreement between the simulated present-day ice volume and observations (Bamber et al., 2013) (Fig. 7f). The simulated extent of the GrIS matches reasonably well the observations. However, notable discrepancies are observed in some sectors. The main differences are found in the northeast, where GRISLI-UCM predicts an ice margin somewhat too far inland, and in the southwest, where our model is not able to make the GrIS retreat as expected. The ice loss in the north is a known problem that appears in many studies when simulating the GrIS during an interglacial (Stone et al., 2010; Born and Nisancioglu, 2012). In the interior, the difference in ice thickness is relatively low. However, the GrIS simulated by our model generally shows thicker ice along the margins, a tendency that propagates inland. Other areas in which our simulated ice thickness is lower than that observed are located in the centre of the continent and in the very southeast corresponding to a mountainous region. However, the focus of our work is not to exactly reproduce the observed present-day GrIS ice volume at the end of the simulations but rather to demonstrate the impact of the ocean on the GrIS past evolution. From this perspective, the simulations arrive at a reasonable representation of the present day and within the range of other models.

Figure 13 Modelled minus observed surface elevation for the present day. Modelled data are taken from the GRISLI-UCM simulation which best estimates the presumed LGM extension (B_ref=1 m a⁻¹ and κ=10 $\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$) while the observed surface elevation is taken from Bamber et al. (2013). Purple and black lines represent simulated and observed GrIS extensions, respectively. Black circles indicate the locations of the Camp Century, NEEM, NGRIP, GRIP and Dye3 ice cores (from north to south).

4 Discussion

Our model simulates the advance and retreat of the GrIS during the last two glacial cycles. Transient simulations reflect the ice-sheet response to the specific oceanic forcing applied to the model. This reaction is different for glacial and interglacial periods (Fig. 5). Since during the interglacial periods the GrIS is almost totally land based and therefore less exposed to the ocean, the minimum ice volume reflects the oceanic imprint only mildly and is limited to a small range of possible values. On the other hand, the volume reached in glacial periods is much more sensitive to κ. Although the maximum ice volume loss is constrained by the imposed limited extension to the continental shelf break, the large ice loss observed for a high oceanic sensitivity is closely related to the GrIS configuration in the previous glacial, which is essentially marine-based at the margins and therefore more subjected to oceanic changes (Fig. 6c). As water temperatures rise at the beginning of the deglaciation, the basal melting rate increases too (Eq. 10), thinning ice shelves at the boundaries, enhancing outflow of ice and triggering grounding-line retreat. The effects of this ocean-driven retreat are not locally confined but are propagated inland through a dynamic response of the grounded ice sheet. The ice loss at the margins triggers ice advection from the interior which further increases the ice discharge into the ocean, and, as the thickness of the inland ice decreases, the elevation–melt feedback begins. At a given stage of the deglaciation, when the whole ice sheet starts to become land based, this atmosphere-driven retreat becomes the sole driver of ice mass loss. The simulated retreat during this phase is influenced by the choice of the surface melt scheme used in the model. At the peak of the Eemian, the melt determined by the PDD scheme can be 20–50 % lower than the melt calculated if past insolation changes are taken into account (Robinson and Goelzer, 2014). This inaccuracy therefore influences the GrIS contributions to sea-level rise for the last interglacial (Fig. 10), which could be underestimated.

As discussed in Sect. 3.4 and 3.5, the oceanic forcing that seems to best reconstruct the past (LGM) and the present GrIS is achieved for a heat-flux coefficient of 10 $\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$. However, the submarine melt scheme used and some simplifications made in its treatment may partly influence our results. Firstly, only a limited range of reference submarine melting rates has been investigated, since only two of the system model parameters have been explored (Table 1). Secondly, our melting parameterization is highly conditioned by the B_ref value assumed to represent the present-day submarine melting rate around the GrIS (Fig. 8), as it consequently determines the minimum κ value needed to allow the GrIS to respond to the ocean (Fig. 3). Using a single value for this term is a coarse approximation to reality, but since the detailed distribution of the present-day sub-shelf melt along the coasts does not yet exist for Greenland, the retrieval of a 2-D field would be complex and highly uncertain. However, we have considered the same order of magnitude of melt rates as is proposed in the literature for the AIS, which spans values from negative to above 40 m a⁻¹ in some very active regions (Rignot and Jacobs, 2002). In support of this evidence, similar basal melting rates have been found recently in some GrIS ice tongues (Wilson et al., 2017). Thirdly, the basal melting equation strongly depends on the oceanic temperature anomaly T_LGM,ocn−T_{PD, ocn}, which has been prescribed to a spatially constant value of −3 K. Since this term impacts the oceanic sensitivity through κ (Eq. 10), it is clear that the same results obtained in this work would have been reached by fixing one value of κ and instead examining the influence of different levels of the ΔT_ocn on the GrIS past evolution. Considering a spatially constant SST anomaly represents an idealized simplification of the oceanic forcing for two reasons: the temperature of the water is clearly not uniform along the GrIS coasts and the melt at the grounding line is presumably controlled by water temperature deeper in the ocean column (between 100 and 1000 m in Greenland; Rignot et al., 2016). These issues could be avoided, for example, by using spatially variable (horizontally and vertically) oceanic temperatures from available model outputs for Greenland. To see whether this simplification could influence our results, some tests using 2-D temperatures from CLIMBER-3α snapshots (Montoya and Levermann, 2008) have been run (not shown). Despite some differences in the ice distribution and the time of the retreat, the main results obtained in this work did not change. Finally, another simplification made here is the assignation of the same climatic index α to both atmosphere and ocean. In principle, forcing the ocean with an index derived from past ocean temperatures could be more appropriate. To this end, we ran additional simulations by applying the multiproxy index α for the atmosphere and another index for the ocean calculated from benthic-retrieved ocean temperatures (Waelbroeck et al., 2002). The results of the new simulations show very little differences from the ones reported here, while the same sensitivity to the ocean is preserved (not shown). Thus, such a distinction in forcing does not affect the main outcomes of this work.

Our results may well be model dependent, and some model limitations should be noted. As described in Sect. 2.1, our ice-sheet–shelf model is provided with an internal GIA scheme which accounts for bedrock deformation due to changes in the GrIS ice load. However, since the GrIS rests on the peripheral forebulge of the North American ice sheets (NAIS), such as the Laurentide Ice Sheet, variations in the NAIS ice load induce consequent vertical motions of the lithosphere beneath the GrIS (Lecavalier et al., 2014). The resulting GrIS isostatic adjustment is therefore the combination of these local and non-local responses which make the GIA treatment rather complex. In principle, these non-local effects should be taken into account as they contribute to the sea-level variability, becoming especially relevant at the beginning of deglaciations when the ice mass loss is significantly induced by sea-level rise (Lecavalier et al., 2017). However, for the sake of simplicity, the GrIS isostatic adjustment is assumed here to be only due to local ice mass variations, as other works have done in the past (Greve and Blatter, 2009; Helsen et al., 2013; Huybrechts, 2002; Langebroek and Nisancioglu, 2016; Stone et al., 2013).

The simulated ice volume at the present day is overestimated for all investigated values of κ (Fig. 8). This fact suggests that our model has a tendency to overestimate the ice thickness of the GrIS, especially in the marginal zones of the domain, a well-known phenomenon (Calov et al., 2015). These discrepancies are partly linked to the relatively low model resolution (20 km × 20 km), which limits the accuracy in estimating the margins especially along the fjords, and partly due to the boundary conditions applied to the ice-sheet model, such as the basal sliding. The coarse model resolution prevents the model from resolving fine-scale physical processes at the marine-terminating outlet glaciers that end in narrow fjords, although they are considered as the primary sources of ice discharge today due to oceanic changes. Such an inability of our model may be more relevant when modelling the GrIS retreat during the LIG and the Holocene. The lack of a sub-grid fjord treatment does not allow for a proper analysis of the ice front processes which become relevant when the retreat has reached the continental area above the sea level. Especially when, as in our case, the submarine melt goes abruptly to a high value at the grounding line, the implementation of a sub-grid-scale parameterization would allow the small processes at the fjords to be accurately resolved (Calov et al., 2015; Favier et al., 2016; Gladstone et al., 2017). However, these limitations lead to only second-order effects given the scope of our work.

The parameterization used for the submarine melting rate at the GrIS marine margins is a simplification compared to other temperature-dependent submarine melting schemes. We are aware that the melting rate depends on many regional factors such as the temperature and salinity of the ocean at the ice-shelf margin, the shape of the ice-shelf cavity and the depth of the grounding line, which our equations do not take into account. However, our simple construction allows us to test the sensitivity of the GrIS to the oceanic forcing in a straightforward manner and is found to be particularly suitable for palaeo-studies.

Our basal melting scheme is implemented in such a way that the melting at the grounding line (Eq. 10) is higher than the one set below the ice shelves (Eq. 11). This approach is supported by sub-shelf melting rate estimates (Dutrieux et al., 2013; Reese et al., 2017; Rignot and Jacobs, 2002; Wilson et al., 2017). Moreover, we assume that the ratio between the two is 1:10, which is valid for the present day, but could be inaccurate for glacial times. However, some experiments done with ratios of 1:5 and 1:15 differ very little from the results presented in this work (not shown). Therefore, our parameterization is much less sensitive to the melting rate below the ice shelves with respect to that at the grounding line. On the other hand, a recent study shows the need to make the basal melting decrease smoothly to zero when approaching the grounding line from the ice shelf to avoid resolution-dependent performances (Gladstone et al., 2017). This can be achieved, for example, by considering the submarine melt to be dependent on the water-column thickness beneath the ice shelf, as Bradley et al. (2017) suggested in their work. It is interesting to compare our results with theirs, as we address the same scientific problem, i.e. the impact of submarine melting on the evolution of the past GrIS, from two different points of view. Our submarine melt scheme is implicitly a linear function of the water depth, as, going down through the water column, the melt rate maintains the same value until it reaches a critical zone at which the sub-shelf melt is set to 50 m a⁻¹ to avoid improbable ice expansion (Sect. 2.3). Our work shows that without melting/freezing at the grounding line (for B_ref=0 and κ=0), the GrIS is not able to reach the continental shelf break (Fig. 12a). However, it is able to extend past the present-day coastline, similar to the simulations presented by Bradley et al. (2017). Moreover, experiments performed under the same oceanic conditions with increased basal sliding at the margins show that our model allows further expansion during the glacial periods (not shown). On the other hand, the model used by Bradley et al. (2017) has the capability of making the GrIS retreat during interglacial periods only if the submarine meltwater depth relation is exponential and if RSL variations due to both local and non-local effects are considered. On the contrary, a proper retreat during the deglaciations is always achieved in our simulations (Figs. 3 and 6–7), although the GIA does not account for global effects. These discrepancies are probably due not only to the different submarine melt schemes considered in each model but also to the features of the model dynamics, such as the sliding law and the grounding-line migration scheme. Following these assumptions, a sub-grid treatment of the small-scale processes taking place at the grounding line, such as basal sliding, sub-shelf melting, hydrology and migration, will be added in our model in the future. This will provide a more realistic description of grounding-line processes such as the enhanced submarine melting as well as the basal drag at the margin of fast grounded ice.

We should finally remark that the GrIS evolution during the last two glacial cycles has been assessed here only from an oceanic point of view, while the influence of different atmospheric forcings has not been investigated. This simplification may be especially important for the results shown for the LIG and the Holocene, in which the retreat is mostly induced by surface ablation. However, this point will be in the scope of future work.

5 Conclusions

Here, we assessed the impact of palaeo-oceanic temperature variations on the evolution of the GrIS on a glacial–interglacial timescale. By using a three-dimensional hybrid ice-sheet–shelf model including a parameterization of the basal melting rate at the GrIS marine margins, the model simulates the evolution of the whole ice sheet under temporally variable oceanic conditions. Firstly, the magnitude of the oceanic forcing applied at the ice–ocean interface triggers and drives the grounding-line advance (through water freezing) and retreat (through ice melting). Secondly, it induces a dynamic adjustment of the grounded ice sheet, determining the amount of ice grown (lost) during the cold (warm) stages. Although the GrIS evolution is a result of the atmospheric and oceanic forcings operating together, we have shown that the ocean is a primary driver of the GrIS glacial advance. Not only must the oceanic forcing be activated, but it must be strong enough to reproduce a reliable GrIS evolution throughout the glacial cycles. It is important to remark that other factors which could affect the GrIS evolution have not yet been explored in detail. Sensitivity tests to the atmospheric forcing, glacial isostatic adjustment effects and spatially non-uniform submarine melt rates should be taken into account in the future to analyse the scientific problem from a broad range of points of view. Nevertheless, we have shown that changing oceanic conditions is a fundamental contributor to the evolution of the whole GrIS, suggesting that the oceanic component should be included as an active forcing in palaeo-ice-sheet models.