Articles | Volume 13, issue 7
Research article
15 Jul 2019
Research article |  | 15 Jul 2019

Submarine melt as a potential trigger of the North East Greenland Ice Stream margin retreat during Marine Isotope Stage 3

Ilaria Tabone, Alexander Robinson, Jorge Alvarez-Solas, and Marisa Montoya

The Northeast Greenland Ice Stream (NEGIS) has been suffering a significant ice mass loss during the last decades. This is partly due to increasing oceanic temperatures in the subpolar North Atlantic, which enhance submarine basal melting and mass discharge. This demonstrates the high sensitivity of this region to oceanic changes. In addition, a recent study suggested that the NEGIS grounding line was 20–40 km behind its present-day location for 15 ka during Marine Isotope Stage (MIS) 3. This is in contrast with Greenland temperature records indicating cold atmospheric conditions at that time, expected to favour ice-sheet expansion. To explain this anomalous retreat a combination of atmospheric and external forcings has been invoked. Yet, as the ocean is found to be a primary driver of the ongoing retreat of the NEGIS glaciers, the effect of past oceanic changes in their paleo evolution cannot be ruled out and should be explored in detail. Here we investigate the sensitivity of the NEGIS to the oceanic forcing during the last glacial period using a three-dimensional hybrid ice-sheet–shelf model. We find that a sufficiently high oceanic forcing could account for a NEGIS ice-margin retreat of several tens of kilometres, potentially explaining the recently proposed NEGIS grounding-line retreat during Marine Isotope Stage 3.

1 Introduction

The Northeast Greenland Ice Stream (NEGIS) is the largest ice stream in the Greenland Ice Sheet (GrIS), extending more than 600 km inland (Joughin et al.2001) and discharging 12 % of the whole ice sheet through three outlet glaciers (Rignot and Mouginot2012): Nioghalvfjerdsfjord Gletscher (79N), Zachariae Isstrøm (ZI), and Storstrømmen Gletscher (SG), which today is a surging glacier. These marine-terminating glaciers have suffered huge changes in the last decades. In less than 15 years the ZI floating tongue has lost 95 % of its size as a result of an enhanced mass loss (Mouginot et al.2015). Concurrently, since 1999 the 79N ice shelf has lost 30 % of its thickness at the grounding line (Mouginot et al.2015), contributing to its inland retreat by 2 km (Mayer et al.2018). However, since 79N is retreating over an upward-sloping bed (Mouginot et al.2015), it may be less prone than ZI to an unstable retreat. This has been recently examined through an ice-flow model pointing out that its floating tongue has to lose several tens of kilometres of ice before the glacier becomes unstable (Rathmann et al.2017). Enhanced stability of 79N has been recently tested under various future warming scenarios by another modelling study (Choi et al.2017), suggesting that it may be related to the presence of pinning points (such as ice rises) near the calving front. Ice loss from these two marine-terminating glaciers is thought to be partly related to the increasing temperature of North Atlantic waters (Khan et al.2014; Mouginot et al.2015), which increases the oceanic heat flux and accelerates the submarine melting (Mayer et al.2018). This hypothesis is supported by the three-decade-long observed warming in the subpolar North Atlantic (Straneo and Heimbach2013, and references therein). Moreover, warmer oceanic waters in Fram Strait could directly reach 79N, further increasing its basal melting and potentially causing the loss of its floating ice tongue (Schaffer et al.2017). Although 79N has been suggested to be more resistant to increasing basal and frontal melt than ZI (Choi et al.2017), new evidence revealing that both glaciers retreated beyond their PD margins during the Holocene indicates that this conclusion may be too conservative (Larsen et al.2018).

Reconstructions suggest that during the Last Glacial Maximum (LGM), ca. 21 ka, the north-eastern region of the GrIS considerably advanced 250–300 km from the present-day coastline, likely reaching the continental shelf break (Arndt et al.2015, 2017; Evans et al.2009; Winkelmann et al.2010). Although the age of these LGM reconstructions is still poorly constrained, the combination of cosmogenic exposure and radiocarbon dating has recently facilitated the reconstruction of the position of the NEGIS over the last 45 kyr (Larsen et al.2018). The paleo records emerging from this study, combined with a collection of geological data assembled in the last 20 years (Arndt et al.2015, 2017; Bennike and Weidick2001; Evans et al.2009; Weidick et al.1996; Winkelmann et al.2010), suggest that its ice margin considerably fluctuated in magnitude throughout this period. Around 41–26 ka during Marine Isotope Stage 3 (MIS-3, ca. 60–25 ka) the NEGIS front was ca. 20–40 km farther inland than today, then advanced by more than 250 km toward the shelf break in the LGM and retreated again during the last deglaciation, at ca. 70 km behind its present-day position, where it stopped most of the mid-Holocene and late Holocene (7.8–1.2 ka). The Holocene retreat was likely due to an increase in both atmospheric and oceanic temperatures, whilst the retreat during MIS-3 was attributed by Larsen et al. (2018) to a combination of atmospheric and external forcings. However, the potential role of oceanic forcing in this retreat has not been explicitly investigated from a modelling perspective. In the light of the ongoing changes in the GrIS attributed to ice–ocean interactions, this appears as a plausible mechanism that needs to be investigated. Moreover, since it is expected that warmer Atlantic waters entering the fjords will strongly affect the NEGIS margin in the future, assessing its response to similar past warm oceanic conditions will provide new insights into the future stability of its glaciers' front.

Here we use an ice-sheet–shelf model to investigate the sensitivity of the NEGIS grounding line to changing oceanic conditions during the last glacial period. The submarine melting at the grounding line is parameterised in such a way that basal melt is allowed during relatively warm time periods such as the present, the last interglacial (LIG; ca. 128–116 ka) or MIS-3, whereas it reaches zero at the onset of the LGM. We study the NEGIS marine margin response to increasing basal melting rates during MIS-3 to show that a sufficiently high oceanic sensitivity could have driven a considerable NEGIS grounding-line retreat during MIS-3 from its former glacial position.

2 Methods

To simulate the NEGIS response to past oceanic forcing, we use the three-dimensional, hybrid ice-sheet–shelf model GRISLI-UCM (Alvarez-Solas et al.2019; Tabone et al.2018), adapted from the extensively used GRISLI model (Ritz et al.2001). Grounded, slow-moving ice-sheet regions and floating shelves are treated through the shallow-ice approximation (SIA) and shallow-shelf approximation (SSA), respectively. In the transition between these two regimes (i.e. fast-moving, grounded ice), the dynamics is solved by the simple addition of the SIA and SSA velocity solutions (Winkelmann et al.2011). The SSA boundary condition is provided by basal sliding below the ice streams following a linear friction law, in which the basal shear stress τb is proportional to the basal velocity ub and to a friction coefficient β dependent on the effective pressure of the water at the base of the ice sheet Neff, as

(1) τ b = - β u b ,


(2) β = c f N eff .

The term cf depends on the characteristics of the bedrock topography (e.g. presence of sediments); Neff is calculated as Neff=ρgH-pw, where ρ is the ice density, g the gravitational acceleration and H the ice thickness. The subglacial water pressure pw comes from a simple basal hydrological model based on a Darcy-type law, for which water flows at the base of temperate ice as driven by a gradient of hydraulic pressure. Despite the simplicity of this hydrology scheme, it provides a fair description of the outflow systems at the base of the ice sheet (Peyaud et al.2007). Glacial isostatic adjustment of the bedrock due to variations in the ice load is reproduced through the elastic lithosphere–relaxing asthenosphere model (Greve and Blatter2009). Unlike some recent hybrid models, the grounding-line position is defined through a pure flotation criterion involving ice thickness at the marine margin and prescribed sea level. Calving occurs whenever a two-constraint thickness rule is satisfied at the ice–ocean interface (Colleoni et al.2014): first, the ice-front thickness must be lower than a fixed threshold (H=200 m here); second, the upstream ice advection does not succeed in preserving the ice-front thickness above that threshold.

The atmospheric temperature forcing applied to the model follows an anomaly method according to which the present-day climatological temperature Tclim,atm is perturbed by past anomalies obtained from a spatially uniform proxy-derived index α(t):

(3) T atm ( t ) = T clim , atm + ( 1 - α ( t ) ) ( T LGM , atm - T PD , atm ) .

The α(t) index is derived from the Greenland temperature reconstruction for the Holocene (Vinther et al.2009), the North Greenland Ice Core Project (NGRIP) reconstruction for the last glacial period (Kindler et al.2014) and the North Greenland Eemian Ice Drilling (NEEM) reconstruction for the LIG (NEEM2013), as in Tabone et al. (2018). The composed signal is then smoothed so that the spectral components below orbital frequencies are removed (i.e. periods less than 16 kyr). By construction, α=1 in the present day (PD) and α=0 in the LGM. Tclim,atm is taken from the regional climate model MAR forced by ERA-Interim (Fettweis et al.2013), averaged over the years 1981–2010. TLGM,atm-TPD,atm is the glacial minus present-day (meaning preindustrial) atmospheric temperature anomaly simulated by the climate model of intermediate complexity CLIMBER-3α (Montoya and Levermann2008). The precipitation field is obtained following a similar approach based on the ratio of LGM and present-day precipitation, scaled by α(t), as

(4) P ann ( t ) = P clim , ann α ( t ) + ( 1 - α ( t ) ) P LGM , ann P PD , ann ,

where PLGM,ann and PPD,ann are the LGM and PD annual precipitation provided by the same climate simulations as TLGM,atm and TPD,atm. This approach has been adopted by many ice-sheet models to represent transient past precipitation when they are not coupled to a climate model (e.g. Banderas et al.2018; Charbit et al.2002, 2007; Colleoni et al.2014; Marshall et al.2000; Marshall and Peltier2002; Marshall and Koutnik2006; Philippon et al.2006; Zweck and Huybrechts2005). Surface ablation is calculated by the simple positive degree day (PDD) scheme (Reeh1989). Although this method does not account for past insolation changes, since here we primarily investigate the sensitivity of the NEGIS to the oceanic forcing during glacial times, the choice of this melt scheme should only have second-order effects on the overall results of this work.

The oceanic forcing is prescribed at the grounding line through a parameterisation of the submarine melt rate based on an anomaly method for which the present-day melt rate is perturbed by its past changes associated with variations in the oceanic temperature (Tabone et al.2018):

(5) B m ( t ) = B ref + κ Δ T ocn ( t ) ,

where Bm(t) is the melt rate at the grounding line at a given time (m a−1), Bref is the present-day basal melting rate at the grounding line (m a−1) and κ is a coefficient representing the heat flux exchanged between water and ice at the ice–ocean front (m a−1 K−1). Past oceanic temperatures below the ice (ΔTocn(t)) evolve as

(6) Δ T ocn ( t ) = ( 1 - α ( t ) ) ( T LGM , ocn - T PD , ocn ) ,

where the α(t) index is that of Eq. (3), and TLGM,ocn-TPD,ocn is the glacial minus interglacial oceanic temperature anomaly (K). The system of Eqs. (5)–(6) can be solved assigning values to Bref, κ, TLGM,ocn and TPD,ocn. However, some simplifications can be considered during the parameter assignment. Following Tabone et al. (2018), the reference melting rate Bref is proportional to the oceanic sensitivity κ, as it is defined as

(7) B ref = κ ( T clim , ocn - T f ) .

Tclim,ocn is the climatological mean of the oceanic temperature considered at the grounding-line depth (K) and Tf is the freezing point temperature at the grounding line (K). The former is depth-dependent; the latter also depends on the distribution of salinity in the water column. Introducing Bref in the equation is a simplification made to avoid the choice of values to be assigned to these two variables, that might be challenging and unconstrained (Beckmann and Goosse2003). For the sake of simplicity, Tclim,ocn-Tf can be considered to be spatially (horizontally and vertically) constant, in the way that Bref is defined to scale directly with κ. Here, we prescribe Tclim,ocn-Tf=1 K; thus Bref=κ1 K. Also, the glacial–interglacial temperature anomaly TLGM,ocn-TPD,ocn is considered here to be spatially uniform and is set to a value of −1 K (Annan and Hargreaves2013; MARGO2009). With these simplifications, the system of Eq. (5)–(6) is thus reduced to a problem of 1 degree of freedom (κ). Investigated values of κ range from 0 to 10 m a−1 K−1; thus Bref ranges from 0 to 10 m a−1. These κ values are consistent with the inference from the Antarctic Ice Sheet that a change of 1 K in the oceanic temperature varies the melt rate by 10 m a−1 (Rignot and Jacobs2002). Moreover, the resulting Bref values are in the range of the submarine melt observed at the grounding line of PD Greenland glaciers that have floating ice shelves (Wilson and Heimbach2017). Melting at the base of the ice shelves is defined to be 10 % of that calculated at the grounding line which reflects the decrease of melting rate observed towards the ice shelves (Anhaus et al.2019; Münchow et al.2014; Rignot and Steffen2008; Wilson and Heimbach2017). However, this decrease is not parameterised here as a function of the distance from the grounding line. Instead, submarine melt is assumed to have a binary behaviour: it is equal to Bm at the grounding line and to 10 % of Bm at all floating grid cells. Since the submarine melting rate at the grounding line is calculated to be spatially constant along the whole domain, the resulting value of the sub-shelf melt rate is also spatially uniform and is shared by all the ice-shelf grid cells of the domain. Note that refreezing below the grounding line is not allowed, and it is cut off to zero; thus there is neither melting nor refreezing during the LGM for the whole set of experiments, which is probably a simplification of reality. Melting and refreezing may vary strongly at local scales, as we know from present observations in Antarctica (e.g. Rignot et al.2013) and Greenland (e.g. Wilson and Heimbach2017). However due to the lack of data for basal melt along the NEGIS margins for the last glacial and the coarse resolution of our model (10 km), this assumption is considered to be the most reasonable approach. The spectrum of resulting submarine melt rates leads to 11 different configurations, for which an increase in the oceanic sensitivity entails an increase of the melting rate during MIS-3 (Fig. 1). These configurations allow the role of the submarine melting rate in the NEGIS margin position during the last glacial period to be investigated. Model simulations of the whole GrIS are initialised at 250 ka using the PD GrIS topography from Schaffer et al. (2016) and run under transient climatic conditions for two full glacial cycles. The first glacial cycle has been considered a spin-up. The analysis of the results focuses on the NEGIS sector.

Figure 1Evolution of the climatic index α and the resulting past oceanic temperature anomaly ΔTocn (K) (a). Potential submarine melt-rate evolution during the last glacial period for increasing Bref and κ values considered in the experiments (b).


3 Results

We calculate the grounding-line distance from the PD position on 48 transects intersecting the ZI and the continental shelf break (Fig. 2). Then we average the results to create one transient evolution for the grounding line for each oceanic forcing. The experiment with submarine melt prescribed to zero (κ=0, Bref=0), which is hereafter referred to as the unperturbed experiment, shows the NEGIS margin rapidly advancing towards the continental shelf during glacial inception (Fig. 3). In less than 20 kyr after the peak of the LIG, the grounding line advances through the inner sector of the continental shelf, extending offshore to a distance of about 250 km from the PD NEGIS margin at around 65 ka. During MIS-3, the ice-margin position gradually advances towards its maximum glacial extent, which is reached at about 20 ka (LGM), when the ice sheet becomes grounded at a mean distance of 40 km from the shelf break, reducing the area of the floating ice shelf in the region (Fig. 4a–e).

Figure 2Map of the NEGIS sector showing the location of its three outlet glaciers (79N, ZI and SG), the observed present-day grounding-line position (solid black line), the observed present-day surface velocities (from Joughin et al.2018), the offshore bathymetry and the onshore ice elevation (both from Schaffer et al.2016) and the maximum (dotted black line) and minimum (dashed black line) grounding-line positions reconstructed for the LGM (Funder et al.2011). The 48 transects used to calculate the evolution of the grounding-line position are shown in purple.


Figure 3Simulated evolution of the NEGIS grounding line relative to its observed present-day position for the set of experiments (coloured lines). The grounding-line distance has been calculated along 48 transects which approximately follow the flow direction of the NEGIS ZI glacier towards the shelf break (Fig. 2). The dashed black line shows the reconstruction by Larsen et al. (2018). Shaded regions represent the time periods corresponding to the LIG, MIS-3 and Holocene. The three dotted lines show the PD NEGIS grounding-line position (0 km) and the maximum (300 km ± 50 km) expected advance of the north-east GrIS to the continental shelf break during the LGM according to Funder et al. (2011).


Figure 4Snapshots of U (m a−1) in total absence of submarine melting (a–e) and in the presence of active orbitally driven oceanic forcing (κ = 8 ma-1K-1; Bref=8m a−1) (f–j) at different times through MIS-3 and the LGM. The sector shown spans an area of about 600 km by 600 km. The black line represents the position of the simulated grounding line. The grey thin solid line represents the observed PD grounding-line position (Schaffer et al.2016). Maximum (dotted black line) and minimum (dashed black line) grounding-line positions reconstructed for the LGM (Funder et al.2011) are also shown.


In all other simulations, the ocean forcing is switched on (κ,Bref>0) and intensifies for increasing κ (Fig. 1). The location of the grounding line at the LIG is the same in all simulations and thus insensitive to κ and set mainly by the atmospheric forcing. Another common feature of these simulations is the response of the grounding-line position right after the peak of the LIG (Fig. 3): the inclusion of positive melt rates before 70 ka somewhat constrains the NEGIS margins 300 km upstream of the grounding-line position obtained for the unperturbed experiment, remaining close to its LIG location. At about 70 ka, the ice margin starts to move towards the continental shelf break, stopping ca. 170 km from its PD position after 10 kyr of rapid advance.

The strongest reaction of the NEGIS grounding line to the applied submarine melting rate is found during MIS-3. By including a basal melt rate of 0–0.5 m a−1 during MIS-3 (κ=1 m a−1 K−1), the location of the NEGIS margin moves 100 km further inland with respect to the unperturbed experiment. Additionally, increasing the oceanic forcing not only helps to preclude the grounding-line advance (as compared to the unperturbed case with no oceanic forcing) but furthermore triggers its retreat. Submarine melt rates of 0–1.2 m a−1 during MIS-3 (κ=4 m a−1 K−1) constrain the NEGIS advance towards the continental shelf after glacial inception (ca. 116 ka) and trigger a slight inland grounding-line retreat by 80 km more, which culminates at around 45 ka. A higher oceanic sensitivity (κ=5 m a−1 K−1) leads to a further and earlier retreat during MIS-3. The minimum extent of grounded ice during MIS-3 is reached at around 50 ka, when the grounding line retreats by more than 100 km inland from its position simulated at 60 ka. The ice margin then remains steady until the end of MIS-3 (Fig. 3). This value of κ and the resulting basal melt configuration (MIS-3 values above 1.6 m a−1) act as a threshold above which the submarine melt rate forces the grounding line to retreat by several kilometres inland during MIS-3, stopping only 40 km far from the PD position. Grounding-line advance and retreat are often very rapid, especially during the first advance after the LIG or during the MIS-3. This is primarily due to the oceanic forcing applied, since the large advance and retreat follow submarine melt-rate evolution well. Part of this stepped nature may be due to the bathymetry too. The area connecting ZI and 79N to the inland shows a bedrock depth of 100–200 m (Morlighem et al.2017). In periods of relatively high sea level, such as the first thousand years (kyr) after the last interglacial, this deep bathymetry may be crucial in driving the grounding line evolution (in our model through the flotation criterion), since it hampers the floating ice to ground and thus the ice sheet advance. This is in line with recent work suggesting that deep bathymetry combined with warmer waters entering the fjord may have important consequences in the destabilisation of 79N (Schaffer et al.2017).

The effect of submarine melt applied to the NEGIS marine margin during MIS-3 is also perceived far inland. The basal melt imposed at the ice–ocean interface causes the ice margin to retreat inland, strongly enhancing ice discharge (Figs. 5 and 6). The reduction of buttressing previously ensured by the presence of ice on the continental shelf increases margin velocities, which propagate inland (Fig. 4f), causing a decrease of ice thickness in the ice-sheet interior (Fig. 6). An initial strong peak in ice discharge is observed, following the initial increase of submarine melting and loss of buttressing, but the effect persists with further ice discharge until the end of MIS-3. At this moment, the absence of melt imposed through the LGM allows the grounding line to advance again towards the continental shelf break (Fig. 4g–j). The maximum distance reached at the peak of the LGM and the time of the onset of the advance are inversely proportional to the melt rate suffered in the previous millennia (Fig. 1). A strong melt rate imposed during MIS-3 leads to a delayed triggering and spatially constrained grounding-line advance, and vice versa.

Figure 5Simulated ice flux at the NEGIS sector for different oceanic forcings. Colours refer to the colour scale of Fig. 3.


Figure 6Evolution of the ice thickness (a), basal melt (b) and ice velocity (c) averaged within the regions A (red lines) and B (blue lines), simulated in the presence of submarine melt during MIS-3 (κ=8ma-1K-1 experiment). Grey lines in (b) show the contribution to surface mass balance (accumulation minus ablation) simulated by the model and averaged over the regions A and B. Dashed lines in the same panel show the potential contributions that would be observed if regions A and B were ice-covered. The black line on the side map represents the LGM maximum extent for the κ=8ma-1K-1 experiment.


Figure 7Simulated GrIS extent during the LGM for different oceanic forcings compared to other glacial reconstructions. Coloured lines follow the colour scale of Fig. 3. The solid black line refers to the maximum glacial extent simulated by Lecavalier et al. (2014), calibrated to match the minimum LGM configuration (Funder et al.2011) in the north-east. The dashed black line represents the expected maximum glacial extent in the north-east sector as inferred from various geological data (Arndt et al.2015, 2017; Evans et al.2009; Winkelmann et al.2010).


By construction, submarine melt occurs again after 20 ka, when both atmospheric and oceanic temperatures increase, contributing to the grounding line being pushed back towards the ice-sheet interior (Fig. 3). This retreat is also simulated in the unperturbed experiment, which demonstrates that the Holocene ice loss is driven by both increasing atmospheric and oceanic temperatures. Nevertheless, the presence of submarine melt at the NEGIS marine margins enhances the retreat and triggers it slightly earlier. This feature, then, saturates for Bref>3ma-1K-1, as further inland retreat is constrained by the bathymetry. However, to specifically assess the relative role between atmospheric and oceanic forcings in the evolution of the NEGIS margin, an equal sensitivity test on the atmospheric forcing, and/or further experiments with another melt scheme, should be carried out.

4 Discussion

4.1 Comparison between modelled and data-derived reconstructions

The NEGIS grounding-line fluctuations simulated in response to a high oceanic forcing in this set of experiments are similar to those suggested by Larsen et al. (2018) for the last 45 kyr. This reconstruction is a result of averaging the evolution of three NEGIS outlet glacier fronts (79N, ZI and SG) inferred from the various geological records with respect to their position at 2014 (Howat et al.2014). Although it is a valuable tool providing a rough idea of the margin fluctuation during the last 45 kyr, caution should be taken before performing a precise one-to-one comparison with model data. Specifically, while the strong retreat during the Holocene is documented for all these glaciers, records showing their margin position during MIS-3 are available only for ZI and SG, which were behind their present location by ca. 20 and 40 km, respectively. However, since they all shared the same behaviour during the Holocene, it is likely that the 79N front was as far inland as the others during MIS-3 (Larsen et al.2018). On this premise, there are some major differences between our results and theirs that deserve further attention.

First, we do not simulate the MIS-3 retreat farther inland than the PD position (20–40 km), although our simulations do show a retreat of more than 100 km with respect to the previous millennia. The observed MIS-3 retreat behind the PD NEGIS margin position has been attributed by Larsen et al. (2018) to lower accumulation rates, high incoming solar radiation and increasing summer air temperatures operating together. Since we have not investigated the sensitivity to these forcings separately, and our experiments do not show this extended retreat, we can neither confirm nor discard their hypothesis. However, our work has demonstrated that orbitally driven oceanic warming during MIS-3 is enough to cause a substantial retreat of the NEGIS margin during part of the last glacial period.

Second, our simulated grounding-line advance during the LGM is smaller than the maximum extension suggested by reconstructions from geological records (Fig. 7). This bias furthermore increases with increasing oceanic forcing. Even in the unperturbed experiment, which allows the largest ice-sheet expansion due to the absence of melting at the marine margins, the grounding line does not reach the continental shelf break. Nevertheless, our simulated extent is still one of the best reconstructions of north-east Greenland during the LGM obtained with an ice-sheet model (Bradley et al.2018; Lecavalier et al.2014; Simpson et al.2009; Tabone et al.2018). This discrepancy in the LGM extents is reflected in the transient GrIS sea-level contribution from the LGM to the present (Fig. S4 in the Supplement), that is underestimated as compared to other recent modelling work (Lecavalier et al.2014; Tabone et al.2018). Nevertheless, our estimation is not far from others (Huybrechts2002; Simpson et al.2009) and well within the range proposed by Buizert et al. (2018). Note that although the LGM extent simulated by Lecavalier et al. (2014) is smaller than ours in the north-east, their ice volume contribution during the glacial maximum is about 1 m sea-level equivalent (SLE) higher. This could be partly due to their larger grounding-line advance in the north-west, but it might also be related to a more active dynamics in our simulations. This is plausible since also the volume discrepancy between our two studies, both performed using the same ice-sheet–shelf model, is likely due to differences in the dynamics. The main reason seems to be related to the fact that SIA and SSA velocities are simply summed up here instead of being mixed through a weighting function as in Tabone et al. (2018). This increases the velocities in the transition zones, promoting discharge of ice from the interior and consequently limiting the ice volume accretion. Second, Tabone et al. (2018) accounted for refreezing processes at the base of the ice shelves, which allowed the grounding line to advance easily, leading to a glacial state in which almost all the GrIS margins were able to reach the shelf break. It is clear that this larger extent could account for a substantial part of the ice volume discrepancy. Another possible reason could be that here we increased the basal drag at the base of grounded temperate ice (by increasing its coefficient cf; see Eq. 2). More friction at the base may foster the production of water at the ice–bed interface through heat release, lubricating the bed and causing the ice flow to accelerate. However, we expect that this process is responsible for only a small fraction of the ice volume discrepancy, since it is counteracted by the increase in basal friction itself. Increasing the total ice volume during the glacial (and its extent) would probably require a substantial tuning effort that is beyond the scope of this study. Our goal is not to provide a perfect match with the LGM but to illustrate a plausible mechanism behind the retreated ice margin at MIS-3 and its subsequent advance.

Third, the timings of the grounding-line advance and retreat for the last 45 kyr of the last glacial period do not precisely correspond to those proposed by Larsen et al. (2018). In the experiments that show a significant retreat during MIS-3 (κ>4ma-1K-1), we simulate both the grounding-line advance at the end of this stage and the retreat at the onset of the Holocene earlier than expected. This is due to the submarine melting signal representing oceanic temperature anomalies, which saturates at zero at around 35 ka and is switched on again at 20 ka, assuming that the LGM starts and ends at these times. Since both atmospheric and oceanic forcings are built through the same index α, any uncertainty in their evolution at orbital timescales affects the ice-sheet retreat during the Holocene, which is supposed to be a combination of atmospheric and oceanic temperatures (Larsen et al.2018).

Although the Holocene maximum is quite well reproduced in our submarine melting configurations (Fig. 1) and in the atmospheric temperature evolution, the slight basal melt decrease applied in the late Holocene is not sufficient to make the grounding line advance back towards the continental borders, and the inaccurate position simulated at the PD is a direct consequence of this simplification. Generally, the grounding-line retreat at the PD is proportional to the magnitude of the submarine melt rate imposed during the mid-Holocene and late Holocene, which is related to the value of Bref used. However, this correspondence is very weak, and the retreat quickly saturates about 70 km away from the PD position along the glacier flow direction for κ>3ma-1K-1, where it is stopped by the presence of bedrock above sea level (Schaffer et al.2016; Morlighem et al.2017). Although such a retreat is supported by proxies for the mid-Holocene (Bennike and Weidick2001; Larsen et al.2018), its persistence until the present day is clearly unrealistic. Moreover, it is likely that imposing PD submarine melting peaks of 50–75 m a−1 as estimated at the 79N grounding line (Anhaus et al.2019; Wilson and Heimbach2017), which are even higher than the Bref values considered in this work, could cause a further inland retreat. This bias could be related to (1) the low spatial resolution of the model (10 km), which does not allow for a precise treatment of the grounding-line zone and may trigger non-linear effects, enhancing grounding-line retreat farther inland than expected; (2) the design of the submarine melt signal itself during the Holocene, which shows a constant increase from 0 to the set Bref value through the last 20 kyr. It is unlikely that such a monotonic increase could have persisted for most of the Holocene, since several records inferred from sediment cores in the Arctic Ocean and in the Fram Strait indicate that water temperatures strongly fluctuated during the Holocene due to the variability of the oceanic currents. Surface (Sztybor and Rasmussen2017) and subsurface (Werner et al.2013; Werner and Polyak2016) oceanic temperatures increased by 3–4 K from the beginning of the Holocene due to the inflow of Atlantic warming waters. After 9–8 kyr, however, these records report a drop in temperatures, gradually (Falardeau et al.2018; Werner and Polyak2016) or interrupted by some peaks of warming (Consolaro et al.2018). These different oceanic conditions between the early and mid-Holocene and late Holocene suggest that such a continuously high melting rate during the whole Holocene is likely overestimated.

4.2 Oceanic forcing at orbital timescales

The oceanic forcing is defined here to be in phase with the atmospheric forcing, as they are both set to evolve in time through the same NGRIP-derived index α. To the best of our knowledge, little evidence on oceanic changes at orbital timescales is available, and whether the best representation of reality would be through oceanic temperatures varying in phase or antiphase with the atmosphere is unclear. However, proxy-based temperature reconstructions indicate glacial–interglacial surface temperature anomalies to be between 0 and −3 K (Annan and Hargreaves2013; MARGO2009), and the value chosen for TLGM,ocn-TPD,ocn is within this range (−1 K).

The rapid occurrence of warm oceanic pulses on millennial timescales is an important characteristic of MIS-3, which is not taken into account here. Available records based on sediment cores of the Arctic Ocean suggest rapid temperature fluctuations as a result of large changes in water masses at different depths. Given the non-linear response of subglacial melting to temperature variations (e.g. Mikkelsen et al.2018), this effect could potentially modulate the orbitally driven response on shorter timescales. Generally, strong oceanic variations are found between glacial–interglacial but also between larger stadial–interstadial transitions (Poirier et al.2012). High sea surface temperature (SST) and low intermediate water temperatures are typical of interstadials, with warmer surface waters generally lasting for 3–4 kyr before cooling (Müller and Stein2014), while the opposite is found during stadials, when warmer Atlantic subsurface waters enter the Nordic Seas up to the Arctic Ocean (Rasmussen et al.2014). Such a strong oceanic temperature variability is also documented by a stack of sediment cores of the Arctic Ocean and the Fram Strait for the last 50 kyr, suggesting that several peaks of warmings occurred during MIS-3 reaching temperatures 1–3 K higher than those recorded for the Holocene (Cronin et al.2012). Nevertheless, the evolution of this temperature record at long (orbital) timescales agrees qualitatively well with that of the melting rate signal used in this work: high melting during MIS-3, prolonged cooling during the LGM and resumed melting during the Holocene.

Thus, even though we remove some degree of realism by not considering the millennial-scale variability in the ocean, our experimental design could represent the evolution of northern Greenland oceanic conditions over long timescales fairly well. A complete treatment of the problem from this perspective is difficult, however, because of the lack of paleooceanographic records of the north-eastern part of Greenland that provide information on the oceanic state during most of the last glacial period at high temporal resolution.

4.3 Model performance on the whole GrIS

The comparison of our results with observations is a good strategy to assess the model performance and to comprehensively evaluate the robustness of our results. At large spatial scales our simulations represent the present state of the GrIS reasonably well (Fig. S1). The maximum differences in surface elevations are found in the south-west and in the east due to a mismatch in ice cover. There, the ice sheet ends in many steep and narrow fjords, which are not adequately represented by the 10 km resolution model. Also, the NEGIS front is located farther inland than observed. The velocity field shows a pretty good agreement in the interior of the ice sheet, where ice speeds are expected to be lower than 50 m a−1 (Joughin et al.2018). However, the simulated ice flow of outlet glaciers and ice streams shows more discrepancies. The speed of the inland flow is generally overestimated, whilst the velocities of streams as they extend far inland is underestimated. By zooming into our domain of interest we see that this pattern is also shared by the NEGIS (Fig. S2 and left panel of Fig. S3). The stream geometry is not adequately resolved, although the spatial distribution of the velocities is somewhat consistent with observations (faster flow at the margins and reduced speed in the interior, as seen in Fig. S2). However, the fast-flow features of the tributary that feed 79N are not reproduced; the SG is faster than expected, and, instead of the long penetrating tongue of ice that characterises the NEGIS, the model simulates a stream draining a wider area. Properly modelling the NEGIS is a well-known problem of ice-sheet models that investigate the evolution of the GrIS at large spatial scales. Most of these models underestimate the stream velocity and do not properly capture its outline (Aschwanden et al.2016; Calov et al.2018; Golledge and Edwards2019; Greve and Herzfeld2013; Seddik et al.2012). Greve and Otsu (2007) succeed in reproducing a correct magnitude of its speed by increasing the basal sliding under the NEGIS by 3 orders of magnitude relative to the rest of the ice sheet, but they fail in reproducing its geometry. A good agreement between model and data is found in Price et al. (2011) and Peano et al. (2017), who use a spatially variable basal friction coefficient derived from an iterative inverse method to match the observed velocities. Our imperfect reproduction of the NEGIS is probably related to a combination of low spatial resolution (10 km) and problems in capturing the dynamics at the base of the ice sheet. Our basal friction coefficient β is a function of the effective water pressure at the base of the ice sheet (Eq. 2), which is a significant degree of freedom in ice-sheet models. A better representation of basal hydrology and sliding could help to improve the simulation of the ice stream. In parallel, new studies on the origin of the stream (following Rogozhina et al.2016), its basal characteristics (following e.g. Keisling et al.2014; Christianson et al.2014; Riverman et al.2019) and new data from the EGRIP ice core (following e.g. Vallelonga et al.2014) will bring new insights in this direction.

Also, the model behaves satisfactorily in simulating the advance and retreat of the GrIS margins throughout the last 120 kyr (see also Tabone et al.2018). Part of this performance is related to the two-condition calving law, that is a function of the critical ice thickness Hf below which the ice edge is calved. Thus, depending on its value, this law may be more or less conservative. Here, with an imposed threshold of 200 m, both 79N and ZI floating tongues are lost at present. It could be that by increasing the value for Hf the model would show slightly more resistance to calve. However, the impact of the calving law is limited to the grid cells at the ice front, while the retreat caused by the submarine melt involves fluctuations in the margin of hundreds of kilometres. Alvarez-Solas et al. (2019) assessed this issue in a different context (the sensitivity of the Eurasian ice sheet to the oceanic forcing during the last glacial period). Their results showed that the overall effect of this parameter is to modulate the amplitude of the response to the oceanic perturbations, but its value did not qualitatively affect their main results. Thus, we expect that changes in the calving critical thickness may cause only second-order effects on the retreat.

4.4 Future perspectives

This work represents the first attempt to simulate the striking margin retreat reconstructed for the NEGIS during MIS-3 (Larsen et al.2018) only accounting for changes in the oceanic forcing. However, such an ocean-driven retreat of the ice margin may have triggered feedbacks on the local climate that are not taken into account in this work. For example, it is possible that this large ice retreat would have caused changes in the albedo, affecting surface air temperatures and snow accumulation. Other feedbacks related to the freshwater flux into the ocean could have led to variations in sea ice and local oceanic circulation. All these processes, not included here, could have additionally contributed to variations in the ice thickness and grounding-line position and should be investigated in the future for a complete understanding of the conundrum. Further experiments accounting for changes in the atmospheric temperatures and precipitation or variations in the external forcing (i.e. insolation) should be carried out for a full understanding of the mechanisms involved in this retreat, here explained by considering the sole impact of the ocean. Particularly, a sensitivity study on climatic variations performed with a prescribed ocean could help to constrain the effect of the atmosphere to eventually evaluate the relative role between the forcings in driving the NEGIS margin.

5 Conclusions

We have studied the sensitivity of the NEGIS ice margin to oceanic forcing during the last glacial period. To this end, we used a three-dimensional, hybrid ice-sheet–shelf model in which the submarine melt rate is parameterised to perform simulations of the GrIS for which basal melt follows a ice-core-proxy-derived curve assumed to represent the evolution of both atmospheric and oceanic temperatures at orbital scales. The increase in basal melt during MIS-3 reflects a relatively warm oceanic state, whereas the lack of basal melt during the LGM corresponds to the associated expected minimum in oceanic temperatures. We showed that in the absence of submarine melting during the entire last glacial period, the grounding line advances towards the continental shelf just after the LIG. On the other hand, switching on the oceanic forcing helps to limit the ice margin advance. Specifically, sufficiently high submarine melt rates during MIS-3 eventually trigger its retreat by more than 100 km from its former position. The lack of basal melt during the LGM then resumes the grounding-line advance by 200 km towards the continental shelf break. Our results robustly show that a prolonged presence of submarine melt at the NEGIS ice margin is enough to substantially contribute to grounding-line retreat there, which helps to explain the recently suggested NEGIS ice margin retreat during MIS-3.

Data availability

The GRISLI-UCM model outputs analysed herein are available upon request.


The supplement related to this article is available online at:

Author contributions

IT performed the experiment, analysed the results and wrote the manuscript. All the authors of this work helped to conceptualise the experiment and write the paper.

Competing interests

The authors declare that they have no conflict of interest.


The model simulations were performed in the HPC of Climate Change of the International Campus of Excellence of Moncloa (EOLO), supported by MECD and MICINN. We are grateful to two anonymous referees, Kerim Nisancioglu and the handling editor Andreas Vieli for their valuable help in improving this work. Also, we are thankful to Catherine Ritz for providing the original GRISLI model.

Financial support

This research has been supported by the Spanish Ministry of Science and Innovation (grant nos. MOCCA (CGL2014-59384-R), RIMA (CGL2017-85975-R)). Ilaria Tabone was funded by the Spanish National Programme for the Promotion of Talent and Its Employability (grant no. BES-2015-074097). Alexander Robinson was funded by the Ramón y Cajal Programme of the Spanish Ministry for Science, Innovation and Universities (grant no. RYC-2016-20587).

Review statement

This paper was edited by Andreas Vieli and reviewed by Kerim Nisancioglu and two anonymous referees.


Alvarez-Solas, J., Banderas, R., Robinson, A., and Montoya, M.: Ocean-driven millennial-scale variability of the Eurasian ice sheet during the last glacial period simulated with a hybrid ice-sheet-shelf model, Clim. Past, 15, 957–979,, 2019. a, b

Anhaus, P., Smedsrud, L. H., Årthun, M., and Straneo, F.: Sensitivity of submarine melting on North East Greenland towards ocean forcing, The Cryosphere Discuss.,, in review, 2019. a, b

Annan, J. D. and Hargreaves, J. C.: A new global reconstruction of temperature changes at the Last Glacial Maximum, Clim. Past, 9, 367–376,, 2013. a, b

Arndt, J. E., Jokat, W., Dorschel, B., Myklebust, R., Dowdeswell, J. A., and Evans, J.: A new bathymetry of the Northeast Greenland continental shelf: Constraints on glacial and other processes, Geochem. Geophys. Geosyst., 16, 3733–3753, 2015. a, b, c

Arndt, J. E., Jokat, W., and Dorschel, B.: The last glaciation and deglaciation of the Northeast Greenland continental shelf revealed by hydro-acoustic data, Quaternary Sci. Rev., 160, 45–56, 2017. a, b, c

Aschwanden, A., Fahnestock, M. A., and Truffer, M.: Complex Greenland outlet glacier flow captured, Nat. Commun., 7, 10524,, 2016. a

Banderas, R., Alvarez-Solas, J., Robinson, A., and Montoya, M.: A new approach for simulating the paleo-evolution of the Northern Hemisphere ice sheets, Geosci. Model Dev., 11, 2299–2314,, 2018. a

Beckmann, A. and Goosse, H.: A parameterization of ice shelf-ocean interaction for climate models, Ocean Modell., 5, 157–170, 2003. a

Bennike, O. and Weidick, A.: Late Quaternary history around Nioghalvfjerdsfjorden and Jøkelbugten, North-East Greenland, Boreas, 30, 205–227, 2001. a, b

Bradley, S. L., Reerink, T. J., van de Wal, R. S. W., and Helsen, M. M.: 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, Clim. Past, 14, 619–635,, 2018. a

Buizert, C., Keisling, B., Box, J., He, F., Carlson, A., Sinclair, G., and DeConto, R.: Greenland-Wide Seasonal Temperatures During the Last Deglaciation, Geophys. Res. Lett., 45, 1905–1914, 2018. a

Calov, R., Beyer, S., Greve, R., Beckmann, J., Willeit, M., Kleiner, T., Rückamp, M., Humbert, A., and Ganopolski, A.: Simulation of the future sea level contribution of Greenland with a new glacial system model, The Cryosphere, 12, 3097–3121,, 2018. a

Charbit, S., Ritz, C., and Ramstein, G.: Simulations of Northern Hemisphere ice-sheet retreat:: sensitivity to physical mechanisms involved during the Last Deglaciation, Quaternary Sci. Rev., 21, 243–265, 2002. a

Charbit, S., Ritz, C., Philippon, G., Peyaud, V., and Kageyama, M.: Numerical reconstructions of the Northern Hemisphere ice sheets through the last glacial-interglacial cycle, Clim. Past, 3, 15–37,, 2007. a

Choi, Y., Morlighem, M., Rignot, E., Mouginot, J., and Wood, M.: Modeling the response of Nioghalvfjerdsfjorden and Zachariae Isstrøm Glaciers, Greenland, to ocean forcing over the next century, Geophys. Res. Lett., 44, 11071–11079, 2017. a, b

Christianson, K., Peters, L. E., Alley, R. B., Anandakrishnan, S., Jacobel, R. W., Riverman, K. L., Muto, A., and Keisling, B. A.: Dilatant till facilitates ice-stream flow in northeast Greenland, Earth Planet. Sci. Lett., 401, 57–69, 2014. a

Colleoni, F., Masina, S., Cherchi, A., Navarra, A., Ritz, C., Peyaud, V., and Otto-Bliesner, B.: Modeling Northern Hemisphere ice-sheet distribution during MIS 5 and MIS 7 glacial inceptions, Clim. Past, 10, 269–291,, 2014. a, b

Consolaro, C., Rasmussen, T. L., and Panieri, G.: Palaeoceanographic and environmental changes in the eastern Fram Strait during the last 14,000 years based on benthic and planktonic foraminifera, Mar. Micropaleontol., 139, 84–101, 2018. a

Cronin, T. M., Dwyer, G. S., Farmer, J., Bauch, H. A., Spielhagen, R. F., Jakobsson, M., Nilsson, J., Briggs Jr., W., and Stepanova, A.: Deep Arctic Ocean warming during the last glacial cycle, Nat. Geosci., 5, 631–634,, 2012. a

Evans, J., Ó Cofaigh, C., Dowdeswell, J. A., and Wadhams, P.: Marine geophysical evidence for former expansion and flow of the Greenland Ice Sheet across the north-east Greenland continental shelf, J. Quatern. Sci. Published for the Quatern. Res. Assoc., 24, 279–293, 2009. a, b, c

Falardeau, J., de Vernal, A., and Spielhagen, R. F.: Paleoceanography of northeastern Fram Strait since the last glacial maximum: Palynological evidence of large amplitude changes, Quaternary Sci. Rev., 195, 133–152, 2018. a

Fettweis, X., Franco, B., Tedesco, M., van Angelen, J. H., Lenaerts, J. T. M., van den Broeke, M. R., and Gallée, H.: Estimating the Greenland ice sheet surface mass balance contribution to future sea level rise using the regional atmospheric climate model MAR, The Cryosphere, 7, 469–489,, 2013. a

Funder, S., Kjeldsen, K. K., Kjær, K. H., and Cofaigh, C. Ó.: The Greenland Ice Sheet during the past 300,000 years: A review, in: Developments in Quaternary Sciences, 15, 699–713, 2011. a, b, c, d

Golledge, N. R., Keller, E. D., Gomez, N., Naughten, K. A., Bernales, J., Trusel, L. D., and Edwards, T. L.: Global environmental consequences of twenty-first-century ice-sheet melt, Nature, 566, 65–72, 2019. a

Greve, R. and Blatter, H.: Dynamics of ice sheets and glaciers, Springer Science & Business Media, 2009. a

Greve, R. and Herzfeld, U. C.: Resolution of ice streams and outlet glaciers in large-scale simulations of the Greenland ice sheet, Ann. Glaciol., 54, 209–220, 2013. a

Greve, R. and Otsu, S.: The effect of the north-east ice stream on the Greenland ice sheet in changing climates, The Cryosphere Discuss., 1, 41–76,, 2007. a

Howat, I. M., Negrete, A., and Smith, B. E.: The Greenland Ice Mapping Project (GIMP) land classification and surface elevation data sets, The Cryosphere, 8, 1509–1518,, 2014. a

Huybrechts, P.: Sea-level changes at the LGM from ice-dynamic reconstructions of the Greenland and Antarctic ice sheets during the glacial cycles, Quaternary Sci. Rev., 21, 203–231, 2002. a

Joughin, I., Fahnestock, M., MacAyeal, D., Bamber, J. L., and Gogineni, P.: Observation and analysis of ice flow in the largest Greenland ice stream, J. Geophys. Res.-Atmos., 106, 34021–34034, 2001. a

Joughin, I., Smith, B. E., and Howat, I.: Greenland Ice Mapping Project: ice flow velocity variation at sub-monthly to decadal timescales, The Cryosphere, 12, 2211–2227,, 2018. a, b

Keisling, B. A., Christianson, K., Alley, R. B., Peters, L. E., Christian, J. E., Anandakrishnan, S., Riverman, K. L., Muto, A., and Jacobel, R. W.: Basal conditions and ice dynamics inferred from radar-derived internal stratigraphy of the northeast Greenland ice stream, Ann. Glaciol., 55, 127–137, 2014. a

Khan, S. A., Kjær, K. H., Bevis, M., Bamber, J. L., Wahr, J., Kjeldsen, K. K., Bjørk, A. A., Korsgaard, N. J., Stearns, L. A., Van Den Broeke, M. R., Liu, L., Larsen, N. K., and Muresan, I.: Sustained mass loss of the northeast Greenland ice sheet triggered by regional warming, Nat. Clim. Change, 4, 292–299, 2014. a

Kindler, P., Guillevic, M., Baumgartner, M., Schwander, J., Landais, A., and Leuenberger, M.: Temperature reconstruction from 10 to 120 kyr b2k from the NGRIP ice core, Clim. Past, 10, 887–902,, 2014. a

Larsen, N. K., Levy, L. B., Carlson, A. E., Buizert, C., Olsen, J., Strunk, A., Bjørk, A. A., and Skov, D. S.: Instability of the Northeast Greenland Ice Stream over the last 45,000 years, Nat. Commun., 9, 1872, 1–8, 2018. a, b, c, d, e, f, g, h, i, j, k

Lecavalier, B. S., Milne, G. A., Simpson, M. J., Wake, L., Huybrechts, P., Tarasov, L., Kjeldsen, K. K., Funder, S., Long, A. J., Woodroffe, S., Dyke, A. S., and Larsenk, N. K.: A model of Greenland ice sheet deglaciation constrained by observations of relative sea level and ice extent, Quaternary Sci. Rev., 102, 54–84, 2014. a, b, c, d

Margo, P. M.: Constraints on the magnitude and patterns of ocean cooling at the Last Glacial Maximum, Nat. Geosci., 2, 127–132, 2009. a, b

Marshall, S. J. and Koutnik, M.: Ice sheet action versus reaction: Distinguishing between Heinrich events and Dansgaard-Oeschger cycles in the North Atlantic, Paleoceanography, 21, PA2021,, 2006. a

Marshall, S. J., Tarasov, L., Clarke, G., and Peltier, W. R.: Glaciological reconstruction of the Laurentide Ice Sheet: physical processes and modelling challenges, Can. J. Earth Sci., 37, 769–793, 2000. a

Marshall, S. J., T. L. C. G. K. and Peltier, W. R.: North American ice sheet reconstructions at the Last Glacial Maximum, Quaternary Sci. Rev., 21, 175–192, 2002. a

Mayer, C., Schaffer, J., Hattermann, T., Floricioiu, D., Krieger, L., Dodd, P. A., Kanzow, T., Licciulli, C., and Schannwell, C.: Large ice loss variability at Nioghalvfjerdsfjorden Glacier, Northeast-Greenland, Nat. Commun., 9, 2768, 1–11, 2018. a, b

Mikkelsen, T. B., Grinsted, A., and Ditlevsen, P.: Influence of temperature fluctuations on equilibrium ice sheet volume, The Cryosphere, 12, 39–47,, 2018. a

Montoya, M. and Levermann, A.: Surface wind-stress threshold for glacial Atlantic overturning, Geophys. Res. Lett., 35, L03608,, 2008. a

Morlighem, M., Williams, C. N., Rignot, E., An, L., Arndt, J. E., Bamber, J. L., Catania, G., Chauché, N., Dowdeswell, J. A., Dorschel, B., Fenty, I., Hogan, K., Howat, I., Hubbard, A., Jakobsson, M., Jordan, T. M., Kjeldsen, K. K., Millan, R., Mayer, L., Mouginot, J., Noël, B. P. Y., O’Cofaigh, C., Palmer, S., Rysgaard, S., Seroussi, H., Siegert, M. J., Slabon, P., Straneo, F., van den Broeke, M. R., Weinrebe, W., Wood, M., and Zinglersen, K. B.: BedMachine v3: Complete bed topography and ocean bathymetry mapping of Greenland from multibeam echo sounding combined with mass conservation, Geophys. Res. Lett., 44, 11–051, 2017. a, b

Mouginot, J., Rignot, E., Scheuchl, B., Fenty, I., Khazendar, A., Morlighem, M., Buzzi, A., and Paden, J.: Fast retreat of Zachariæ Isstrøm, northeast Greenland, Science, 350, 1357–1361, 2015. a, b, c, d

Müller, J. and Stein, R.: High-resolution record of late glacial and deglacial sea ice changes in Fram Strait corroborates ice–ocean interactions during abrupt climate shifts, Earth Planet. Sci. Lett., 403, 446–455, 2014. a

Münchow, A., Padman, L., and Fricker, H. A.: Interannual changes of the floating ice shelf of Petermann Gletscher, North Greenland, from 2000 to 2012, J. Glaciol., 60, 489–499, 2014. a

NEEM: Eemian interglacial reconstructed from a Greenland folded ice core, Nature, 493, 489–494, 2013. a

Peano, D., Colleoni, F., Quiquet, A., and Masina, S.: Ice flux evolution in fast flowing areas of the Greenland ice sheet over the 20th and 21st centuries, J. Glaciol., 63, 499–513, 2017. a

Peyaud, V., Ritz, C., and Krinner, G.: Modelling the Early Weichselian Eurasian Ice Sheets: role of ice shelves and influence of ice-dammed lakes, Clim. Past, 3, 375–386,, 2007. a

Philippon, G., Ramstein, G., Charbit, S., Kageyama, M., Ritz, C., and Dumas, C.: Evolution of the Antarctic ice sheet throughout the last deglaciation: A study with a new coupled climate-north and south hemisphere ice sheet model, Earth Planet. Sci. Lett., 248, 750–758, 2006. a

Poirier, R., Cronin, T., Briggs, W., and Lockwood, R.: Central Arctic paleoceanography for the last 50 kyr based on ostracode faunal assemblages, Mar. Micropaleontol., 88, 65–76, 2012. a

Price, S. F., Payne, A. J., Howat, I. M., and Smith, B. E.: Committed sea-level rise for the next century from Greenland ice sheet dynamics during the past decade, P. Natl. Acad. Sci. USA, 108, 8978–8983, 2011. a

Rasmussen, T. L., Thomsen, E., and Nielsen, T.: Water mass exchange between the Nordic seas and the Arctic Ocean on millennial timescale during MIS 4-MIS 2, Geochem. Geophys. Geosyst., 15, 530–544, 2014. a

Rathmann, N., Hvidberg, C., Solgaard, A., Grinsted, A., Gudmundsson, G. H., Langen, P. L., Nielsen, K., and Kusk, A.: Highly temporally resolved response to seasonal surface melt of the Zachariae and 79N outlet glaciers in northeast Greenland, Geophys. Res. Lett., 44, 9805–9814, 2017. a

Reeh, N.: Parameterization of melt rate and surface temperature on the Greenland ice sheet, Polarforschung, 59, 113–128, 1989. a

Rignot, E. and Jacobs, S. S.: Rapid bottom melting widespread near Antarctic Ice Sheet grounding lines, Science, 296, 2020–2023, 2002. a

Rignot, E. and Mouginot, J.: Ice flow in Greenland for the international polar year 2008–2009, Geophys. Res. Lett., 39, L11501,, 2012. a

Rignot, E. and Steffen, K.: Channelized bottom melting and stability of floating ice shelves, Geophys. Res. Lett., 35, L02503,, 2008. a

Rignot, E., Jacobs, S., Mouginot, J., and Scheuchl, B.: Ice-shelf melting around Antarctica, Science, 341, 266–270, 2013. a

Ritz, C., Rommelaere, V., and Dumas, C.: Modeling the evolution of Antarctic Ice Sheet over the last 420,000 years. Implications for altitude changes in the Vostok region, J. Geophys. Res., 106, 31943–31964, 2001. a

Riverman, K., Alley, R. B., Anandakrishnan, S., Christianson, K., Holschuh, N., Medley, B., Muto, A., and Peters, L.: Enhanced Firn Densification in High-Accumulation Shear Margins of the NE Greenland Ice Stream, J. Geophys. Res.-Earth Surf., 124, 365–382, 2019. a

Rogozhina, I., Petrunin, A. G., Vaughan, A. P., Steinberger, B., Johnson, J. V., Kaban, M. K., Calov, R., Rickers, F., Thomas, M., and Koulakov, I.: Melting at the base of the Greenland ice sheet explained by Iceland hotspot history, Nat. Geosci., 9, 366–269, 2016. a

Schaffer, J., Timmermann, R., Arndt, J. E., Kristensen, S. S., Mayer, C., Morlighem, M., and Steinhage, D.: A global, high-resolution data set of ice sheet topography, cavity geometry, and ocean bathymetry, Earth Syst. Sci. Data, 8, 543–557, 2016. a, b, c, d

Schaffer, J., von Appen, W.-J., Dodd, P. A., Hofstede, C., Mayer, C., de Steur, L., and Kanzow, T.: Warm water pathways toward Nioghalvfjerdsfjorden Glacier, northeast Greenland, J. Geophys. Res.-Oceans, 122, 4004–4020, 2017. a, b

Seddik, H., Greve, R., Zwinger, T., Gillet-Chaulet, F., and Gagliardini, O.: Simulations of the Greenland ice sheet 100 years into the future with the full Stokes model Elmer/Ice, J. Glaciol., 58, 427–440, 2012. a

Simpson, M. J., Milne, G. A., Huybrechts, P., and Long, A. J.: Calibrating a glaciological model of the Greenland ice sheet from the Last Glacial Maximum to present-day using field observations of relative sea level and ice extent, Quaternary Sci. Rev., 28, 1631–1657, 2009. a, b

Straneo, F. and Heimbach, P.: North Atlantic warming and the retreat of Greenland's outlet glaciers, Nature, 504, 36–43,, 2013. a

Sztybor, K. and Rasmussen, T.: Late glacial and deglacial palaeoceanographic changes at Vestnesa Ridge, Fram Strait: Methane seep versus non-seep environments, Palaeogeogr. Palaeoclimatol. Palaeoecol., 476, 77–89, 2017. a

Tabone, I., Blasco, J., Robinson, A., Alvarez-Solas, J., and Montoya, M.: The sensitivity of the Greenland Ice Sheet to glacial-interglacial oceanic forcing, Clim. Past, 14, 455–472,, 2018. a, b, c, d, e, f, g, h, i

Vallelonga, P., Christianson, K., Alley, R. B., Anandakrishnan, S., Christian, J. E. M., Dahl-Jensen, D., Gkinis, V., Holme, C., Jacobel, R. W., Karlsson, N. B., Keisling, B. A., Kipfstuhl, S., Kjær, H. A., Kristensen, M. E. L., Muto, A., Peters, L. E., Popp, T., Riverman, K. L., Svensson, A. M., Tibuleac, C., Vinther, B. M., Weng, Y., and Winstrup, M.: Initial results from geophysical surveys and shallow coring of the Northeast Greenland Ice Stream (NEGIS), The Cryosphere, 8, 1275–1287,, 2014. a

Vinther, B. M., Buchardt, S. L., Clausen, H. B., Dahl- Jensen, D., Johnsen, S. J., Fisher, D., Koerner, R., Raynaud, D., Lipenkov, V., Andersen, K., Blunier, T., Rasmussen, S. O., Steffensen, J. P., and Svensson, A. M.:: Holocene thinning of the Greenland ice sheet, Nature, 461, 385–388, 2009.  a

Weidick, A., Andreasen, C., Oerter, H., and Reeh, N.: Neoglacial glacier changes around Storstrommen, North-East Greenland, Polarforschung, 64, 95–108, 1996. a

Werner, K., Spielhagen, R. F., Bauch, D., Hass, H. C., and Kandiano, E.: Atlantic Water advection versus sea-ice advances in the eastern Fram Strait during the last 9 ka: Multiproxy evidence for a two-phase Holocene, Paleoceanography, 28, 283–295, 2013. a

Werner, K., Müllerb, J., Husum, K., Spielhagend, R. F., Kandiano, E. S., and Polyak, L.: Holocene sea subsurface and surface water masses in the Fram Strait – Comparisons of temperature and sea-ice reconstructions, Quaternary Sci. Rev., 147, 194–209, 2016. a, b

Wilson, N., Straneo, F., and Heimbach, P.: Satellite-derived submarine melt rates and mass balance (2011–2015) for Greenland's largest remaining ice tongues, The Cryosphere, 11, 2773–2782,, 2017. a, b, c, d

Winkelmann, D., Jokat, W., Jensen, L., and Schenke, H.-W.: Submarine end moraines on the continental shelf off NE Greenland–Implications for Lateglacial dynamics, Quaternary Sci. Rev., 29, 1069–1077, 2010. a, b, c

Winkelmann, R., Martin, M. A., Haseloff, M., Albrecht, T., Bueler, E., Khroulev, C., and Levermann, A.: The Potsdam Parallel Ice Sheet Model (PISM-PIK) – Part 1: Model description, The Cryosphere, 5, 715–726,, 2011. a

Zweck, C. and Huybrechts, P.: Modeling of the northern hemisphere ice sheets during the last glacial cycle and glaciological sensitivity, J. Geophys. Res.-Atmos., 110, D07103,, 2005. a

Short summary
Recent reconstructions show that the North East Greenland Ice Stream (NEGIS) retreated away from its present-day position by 20–40 km during MIS-3. Atmospheric and external forcings were proposed as potential causes of this retreat, but the role of the ocean was not considered. Here, using a 3-D ice-sheet model, we suggest that oceanic warming is sufficient to induce a retreat of the NEGIS margin of many tens of kilometres during MIS-3, helping to explain this conundrum.