The Cryosphere, 13, 1911–1923, 2019 https://doi.org/10.5194/tc-13-1911-2019 © Author(s) 2019. This work is distributed under the Creative Commons Attribution 4.0 License. Submarine melt as a potential trigger of the North East Greenland Ice Stream margin retreat during Marine Isotope Stage 3 Ilaria Tabone1,2, Alexander Robinson1,2, Jorge Alvarez-Solas1,2, and Marisa Montoya1,2 1Universidad Complutense de Madrid, 28040 Madrid, Spain 2Instituto de Geociencias, Consejo Superior de Investigaciones Cientificas-Universidad Complutense de Madrid, 28040 Madrid, Spain Correspondence: Ilaria Tabone (itabone@ucm.es) Received: 22 October 2018 – Discussion started: 29 November 2018 Revised: 26 April 2019 – Accepted: 20 June 2019 – Published: 15 July 2019 Abstract. 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 tempera- tures in the subpolar North Atlantic, which enhance subma- rine basal melting and mass discharge. This demonstrates the high sensitivity of this region to oceanic changes. In addi- tion, a recent study suggested that the NEGIS grounding line was 20–40 km behind its present-day location for 15 ka dur- ing 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 expan- sion. To explain this anomalous retreat a combination of at- mospheric 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, po- tentially 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 dis- charging 12 % of the whole ice sheet through three outlet glaciers (Rignot and Mouginot, 2012): 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 re- cently 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 var- ious 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 calv- ing 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 Heimbach, 2013, 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 Published by Copernicus Publications on behalf of the European Geosciences Union. 1912 I. Tabone et al.: The NEGIS glacial retreat triggered by the ocean 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 conser- vative (Larsen et al., 2018). Reconstructions suggest that during the Last Glacial Max- imum (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 ex- posure and radiocarbon dating has recently facilitated the re- construction 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 assem- bled in the last 20 years (Arndt et al., 2015, 2017; Ben- nike and Weidick, 2001; Evans et al., 2009; Weidick et al., 1996; Winkelmann et al., 2010), suggest that its ice mar- gin considerably fluctuated in magnitude throughout this pe- riod. 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 to- ward the shelf break in the LGM and retreated again during the last deglaciation, at ca. 70 km behind its present-day po- sition, 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 forc- ings. 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 plausi- ble 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, assess- ing 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 peri- ods 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 con- siderable 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 ap- proximation (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 pro- portional 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 τ b =−β ub, (1) where β = cfNeff. (2) The term cf depends on the characteristics of the bedrock to- pography (e.g. presence of sediments); Neff is calculated as Neff = ρgH −pw, where ρ is the ice density, g the gravita- tional 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 adjust- ment of the bedrock due to variations in the ice load is repro- duced through the elastic lithosphere–relaxing asthenosphere model (Greve and Blatter, 2009). Unlike some recent hybrid models, the grounding-line position is defined through a pure flotation criterion involving ice thickness at the marine mar- gin 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): Tatm(t)= Tclim,atm+ (1−α(t))(TLGM,atm− TPD,atm). (3) The α(t) index is derived from the Greenland temperature reconstruction for the Holocene (Vinther et al., 2009), the The Cryosphere, 13, 1911–1923, 2019 www.the-cryosphere.net/13/1911/2019/ I. Tabone et al.: The NEGIS glacial retreat triggered by the ocean 1913 North Greenland Ice Core Project (NGRIP) reconstruction for the last glacial period (Kindler et al., 2014) and the North Greenland Eemian Ice Drilling (NEEM) reconstruc- tion for the LIG (NEEM, 2013), as in Tabone et al. (2018). The composed signal is then smoothed so that the spectral components below orbital frequencies are removed (i.e. pe- riods 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 Levermann, 2008). The precipitation field is obtained following a similar approach based on the ratio of LGM and present-day precipitation, scaled by α(t), as Pann(t)= Pclim,ann · ( α(t)+ (1−α(t)) · PLGM,ann PPD,ann ) , (4) 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 precipi- tation 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 Peltier, 2002; Marshall and Koutnik, 2006; Philippon et al., 2006; Zweck and Huybrechts, 2005). Surface ablation is calculated by the simple positive degree day (PDD) scheme (Reeh, 1989). 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): Bm(t)= Bref+ κ1Tocn(t), (5) 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 (1Tocn(t)) evolve as 1Tocn(t)= (1−α(t))(TLGM,ocn− TPD,ocn), (6) 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 val- ues to Bref, κ , TLGM,ocn and TPD,ocn. However, some simpli- fications 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 de- fined as Bref = κ(Tclim,ocn− Tf). (7) Tclim,ocn is the climatological mean of the oceanic tempera- ture considered at the grounding-line depth (K) and Tf is the freezing point temperature at the grounding line (K). The for- mer is depth-dependent; the latter also depends on the distri- bution of salinity in the water column. Introducing Bref in the equation is a simplification made to avoid the choice of val- ues to be assigned to these two variables, that might be chal- lenging and unconstrained (Beckmann and Goosse, 2003). 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 Hargreaves, 2013; MARGO, 2009). 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 Jacobs, 2002). Moreover, the resulting Bref val- ues are in the range of the submarine melt observed at the grounding line of PD Greenland glaciers that have floating ice shelves (Wilson and Heimbach, 2017). 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 Steffen, 2008; Wilson and Heimbach, 2017). However, this decrease is not parame- terised here as a function of the distance from the grounding line. Instead, submarine melt is assumed to have a binary be- haviour: 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 simplifica- tion of reality. Melting and refreezing may vary strongly at local scales, as we know from present observations in Antarc- tica (e.g. Rignot et al., 2013) and Greenland (e.g. Wilson and Heimbach, 2017). 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 spec- trum of resulting submarine melt rates leads to 11 different www.the-cryosphere.net/13/1911/2019/ The Cryosphere, 13, 1911–1923, 2019 1914 I. Tabone et al.: The NEGIS glacial retreat triggered by the ocean Figure 1. Evolution of the climatic index α and the resulting past oceanic temperature anomaly 1Tocn (K) (a). Potential submarine melt-rate evolution during the last glacial period for increasing Bref and κ values considered in the experiments (b). configurations, for which an increase in the oceanic sensi- tivity entails an increase of the melting rate during MIS-3 (Fig. 1). These configurations allow the role of the subma- rine 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 tran- sient 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. 3 Results We calculate the grounding-line distance from the PD po- sition 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 incep- tion (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 dis- Figure 2. Map 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 evo- lution of the grounding-line position are shown in purple. tance of 40 km from the shelf break, reducing the area of the floating ice shelf in the region (Fig. 4a–e). In all other simulations, the ocean forcing is switched on (κ,Bref > 0) and intensifies for increasing κ (Fig. 1). The lo- cation 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 simu- lations is the response of the grounding-line position right af- ter 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 loca- tion. 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 exper- iment. 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 further- The Cryosphere, 13, 1911–1923, 2019 www.the-cryosphere.net/13/1911/2019/ I. Tabone et al.: The NEGIS glacial retreat triggered by the ocean 1915 Figure 3. Simulated evolution of the NEGIS grounding line rela- tive 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). more 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 ad- vance towards the continental shelf after glacial inception (ca. 116 ka) and trigger a slight inland grounding-line re- treat by 80 km more, which culminates at around 45 ka. A higher oceanic sensitivity (κ = 5 m a−1 K−1) leads to a fur- ther 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 dur- ing 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 ad- vance. 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 mar- gin 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 buttress- ing, 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 maxi- mum 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. 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 sub- marine melt at the NEGIS marine margins enhances the re- treat and triggers it slightly earlier. This feature, then, satu- rates for Bref > 3 ma−1 K−1, as further inland retreat is con- strained 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 re- sponse 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 per- forming a precise one-to-one comparison with model data. Specifically, while the strong retreat during the Holocene www.the-cryosphere.net/13/1911/2019/ The Cryosphere, 13, 1911–1923, 2019 1916 I. Tabone et al.: The NEGIS glacial retreat triggered by the ocean Figure 4. Snapshots of U (ma−1) in total absence of submarine melting (a–e) and in the presence of active orbitally driven oceanic forcing (κ = 8 ma−1 K−1; Bref = 8 ma−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. Figure 5. Simulated ice flux at the NEGIS sector for different oceanic forcings. Colours refer to the colour scale of Fig. 3. 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 be- tween 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 pre- vious 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 ra- diation and increasing summer air temperatures operating to- gether. Since we have not investigated the sensitivity to these forcings separately, and our experiments do not show this ex- tended retreat, we can neither confirm nor discard their hy- pothesis. 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 (Huybrechts, 2002; 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 discrep- The Cryosphere, 13, 1911–1923, 2019 www.the-cryosphere.net/13/1911/2019/ I. Tabone et al.: The NEGIS glacial retreat triggered by the ocean 1917 Figure 6. Evolution 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 (κ = 8 ma−1 K−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 κ = 8 ma−1 K−1 experiment. ancy between our two studies, both performed using the same ice-sheet–shelf model, is likely due to differences in the dy- namics. 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 con- sequently 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 ad- vance 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 wa- ter 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 frac- tion 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 re- quire 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 re- treat for the last 45 kyr of the last glacial period do not pre- cisely correspond to those proposed by Larsen et al. (2018). In the experiments that show a significant retreat during MIS- 3 (κ > 4 ma−1 K−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 subma- rine melting signal representing oceanic temperature anoma- lies, 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 www.the-cryosphere.net/13/1911/2019/ The Cryosphere, 13, 1911–1923, 2019 1918 I. Tabone et al.: The NEGIS glacial retreat triggered by the ocean Figure 7. Simulated GrIS extent during the LGM for differ- ent 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 repre- sents 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). built through the same index α, any uncertainty in their evo- lution at orbital timescales affects the ice-sheet retreat during the Holocene, which is supposed to be a combination of at- mospheric and oceanic temperatures (Larsen et al., 2018). Although the Holocene maximum is quite well repro- duced 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 con- tinental borders, and the inaccurate position simulated at the PD is a direct consequence of this simplification. Gen- erally, the grounding-line retreat at the PD is proportional to the magnitude of the submarine melt rate imposed dur- ing 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 κ > 3 ma−1 K−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 Weidick, 2001; Larsen et al., 2018), its persistence until the present day is clearly unrealistic. Moreover, it is likely that imposing PD subma- rine melting peaks of 50–75 ma−1 as estimated at the 79N grounding line (Anhaus et al., 2019; Wilson and Heimbach, 2017), 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, en- hancing 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 wa- ter temperatures strongly fluctuated during the Holocene due to the variability of the oceanic currents. Surface (Sztybor and Rasmussen, 2017) and subsurface (Werner et al., 2013; Werner and Polyak, 2016) 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 Polyak, 2016) 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 or- bital timescales is available, and whether the best representa- tion of reality would be through oceanic temperatures vary- ing in phase or antiphase with the atmosphere is unclear. However, proxy-based temperature reconstructions indicate glacial–interglacial surface temperature anomalies to be be- tween 0 and −3 K (Annan and Hargreaves, 2013; MARGO, 2009), and the value chosen for TLGM,ocn− TPD,ocn is within this range (−1 K). The rapid occurrence of warm oceanic pulses on millen- nial 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 The Cryosphere, 13, 1911–1923, 2019 www.the-cryosphere.net/13/1911/2019/ I. Tabone et al.: The NEGIS glacial retreat triggered by the ocean 1919 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 be- fore cooling (Müller and Stein, 2014), while the opposite is found during stadials, when warmer Atlantic subsurface wa- ters 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 Arc- tic Ocean and the Fram Strait for the last 50 kyr, suggesting that several peaks of warmings occurred during MIS-3 reach- ing 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 pale- ooceanographic records of the north-eastern part of Green- land 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 comprehen- sively 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 sur- face 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 rep- resented 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 ma−1 (Joughin et al., 2018). However, the simulated ice flow of outlet glaciers and ice streams shows more discrepan- cies. 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, in- stead 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 prob- lem of ice-sheet models that investigate the evolution of the GrIS at large spatial scales. Most of these models underesti- mate the stream velocity and do not properly capture its out- line (Aschwanden et al., 2016; Calov et al., 2018; Golledge and Edwards, 2019; Greve and Herzfeld, 2013; Seddik et al., 2012). Greve and Otsu (2007) succeed in reproducing a cor- rect 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 geome- try. A good agreement between model and data is found in Price et al. (2011) and Peano et al. (2017), who use a spa- tially variable basal friction coefficient derived from an iter- ative 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 sig- nificant degree of freedom in ice-sheet models. A better rep- resentation of basal hydrology and sliding could help to im- prove the simulation of the ice stream. In parallel, new stud- ies 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 perfor- mance 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 thresh- old 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. How- ever, 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 kilo- metres. 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 re- sults showed that the overall effect of this parameter is to modulate the amplitude of the response to the oceanic per- turbations, 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 strik- ing 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 www.the-cryosphere.net/13/1911/2019/ The Cryosphere, 13, 1911–1923, 2019 1920 I. Tabone et al.: The NEGIS glacial retreat triggered by the ocean 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 ac- cumulation. 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 in- vestigated 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 rela- tive 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 per- form simulations of the GrIS for which basal melt follows a ice-core-proxy-derived curve assumed to represent the evo- lution of both atmospheric and oceanic temperatures at or- bital 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. Specif- ically, 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 pro- longed 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. Supplement. The supplement related to this article is available online at: https://doi.org/10.5194/tc-13-1911-2019-supplement. 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. Acknowledgements. 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 Span- ish 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 Promo- tion 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 re- viewed by Kerim Nisancioglu and two anonymous referees. References 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, https://doi.org/10.5194/cp-15-957-2019, 2019. Anhaus, P., Smedsrud, L. H., Årthun, M., and Straneo, F.: Sensitiv- ity of submarine melting on North East Greenland towards ocean forcing, The Cryosphere Discuss., https://doi.org/10.5194/tc- 2019-35, in review, 2019. Annan, J. D. and Hargreaves, J. C.: A new global reconstruction of temperature changes at the Last Glacial Maximum, Clim. Past, 9, 367–376, https://doi.org/10.5194/cp-9-367-2013, 2013. Arndt, J. E., Jokat, W., Dorschel, B., Myklebust, R., Dowdeswell, J. A., and Evans, J.: A new bathymetry of the Northeast Green- land continental shelf: Constraints on glacial and other processes, Geochem. Geophys. Geosyst., 16, 3733–3753, 2015. Arndt, J. E., Jokat, W., and Dorschel, B.: The last glaciation and deglaciation of the Northeast Greenland continental shelf re- vealed by hydro-acoustic data, Quaternary Sci. Rev., 160, 45–56, 2017. Aschwanden, A., Fahnestock, M. A., and Truffer, M.: Complex Greenland outlet glacier flow captured, Nat. Commun., 7, 10524, https://doi.org/10.1038/ncomms10524, 2016. The Cryosphere, 13, 1911–1923, 2019 www.the-cryosphere.net/13/1911/2019/ https://doi.org/10.5194/tc-13-1911-2019-supplement https://doi.org/10.5194/cp-15-957-2019 https://doi.org/10.5194/tc-2019-35 https://doi.org/10.5194/tc-2019-35 https://doi.org/10.5194/cp-9-367-2013 https://doi.org/10.1038/ncomms10524 I. Tabone et al.: The NEGIS glacial retreat triggered by the ocean 1921 Banderas, R., Alvarez-Solas, J., Robinson, A., and Montoya, M.: A new approach for simulating the paleo-evolution of the North- ern Hemisphere ice sheets, Geosci. Model Dev., 11, 2299–2314, https://doi.org/10.5194/gmd-11-2299-2018, 2018. Beckmann, A. and Goosse, H.: A parameterization of ice shelf- ocean interaction for climate models, Ocean Modell., 5, 157– 170, 2003. Bennike, O. and Weidick, A.: Late Quaternary history around Nioghalvfjerdsfjorden and Jøkelbugten, North-East Greenland, Boreas, 30, 205–227, 2001. 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 parameter- ization and relative sea level forcing in an ice-sheet-ice-shelf model, Clim. Past, 14, 619–635, https://doi.org/10.5194/cp-14- 619-2018, 2018. Buizert, C., Keisling, B., Box, J., He, F., Carlson, A., Sinclair, G., and DeConto, R.: Greenland-Wide Seasonal Temperatures Dur- ing the Last Deglaciation, Geophys. Res. Lett., 45, 1905–1914, 2018. Calov, R., Beyer, S., Greve, R., Beckmann, J., Willeit, M., Kleiner, T., Rückamp, M., Humbert, A., and Ganopolski, A.: Simula- tion of the future sea level contribution of Greenland with a new glacial system model, The Cryosphere, 12, 3097–3121, https://doi.org/10.5194/tc-12-3097-2018, 2018. Charbit, S., Ritz, C., and Ramstein, G.: Simulations of North- ern Hemisphere ice-sheet retreat:: sensitivity to physical mech- anisms involved during the Last Deglaciation, Quaternary Sci. Rev., 21, 243–265, 2002. 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, https://doi.org/10.5194/cp-3-15-2007, 2007. 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 cen- tury, Geophys. Res. Lett., 44, 11071–11079, 2017. 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. 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 incep- tions, Clim. Past, 10, 269–291, https://doi.org/10.5194/cp-10- 269-2014, 2014. Consolaro, C., Rasmussen, T. L., and Panieri, G.: Palaeoceano- graphic 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. 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. Evans, J., Ó Cofaigh, C., Dowdeswell, J. A., and Wadhams, P.: Ma- rine 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. Falardeau, J., de Vernal, A., and Spielhagen, R. F.: Paleoceanogra- phy of northeastern Fram Strait since the last glacial maximum: Palynological evidence of large amplitude changes, Quaternary Sci. Rev., 195, 133–152, 2018. 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 fu- ture sea level rise using the regional atmospheric climate model MAR, The Cryosphere, 7, 469–489, https://doi.org/10.5194/tc- 7-469-2013, 2013. 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. Golledge, N. R., Keller, E. D., Gomez, N., Naughten, K. A., Bernales, J., Trusel, L. D., and Edwards, T. L.: Global environ- mental consequences of twenty-first-century ice-sheet melt, Na- ture, 566, 65–72, 2019. Greve, R. and Blatter, H.: Dynamics of ice sheets and glaciers, Springer Science & Business Media, 2009. Greve, R. and Herzfeld, U. C.: Resolution of ice streams and out- let glaciers in large-scale simulations of the Greenland ice sheet, Ann. Glaciol., 54, 209–220, 2013. 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, https://doi.org/10.5194/tcd-1-41-2007, 2007. Howat, I. M., Negrete, A., and Smith, B. E.: The Green- land Ice Mapping Project (GIMP) land classification and surface elevation data sets, The Cryosphere, 8, 1509–1518, https://doi.org/10.5194/tc-8-1509-2014, 2014. 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. Joughin, I., Fahnestock, M., MacAyeal, D., Bamber, J. L., and Gogi- neni, P.: Observation and analysis of ice flow in the largest Green- land ice stream, J. Geophys. Res.-Atmos., 106, 34021–34034, 2001. Joughin, I., Smith, B. E., and Howat, I.: Greenland Ice Map- ping Project: ice flow velocity variation at sub-monthly to decadal timescales, The Cryosphere, 12, 2211–2227, https://doi.org/10.5194/tc-12-2211-2018, 2018. Keisling, B. A., Christianson, K., Alley, R. B., Peters, L. E., Chris- tian, 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. Khan, S. A., Kjær, K. H., Bevis, M., Bamber, J. L., Wahr, J., Kjeld- sen, 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.: Sus- tained mass loss of the northeast Greenland ice sheet triggered by regional warming, Nat. Clim. Change, 4, 292–299, 2014. 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, https://doi.org/10.5194/cp-10-887-2014, 2014. 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 North- east Greenland Ice Stream over the last 45,000 years, Nat. Com- mun., 9, 1872, 1–8, 2018. www.the-cryosphere.net/13/1911/2019/ The Cryosphere, 13, 1911–1923, 2019 https://doi.org/10.5194/gmd-11-2299-2018 https://doi.org/10.5194/cp-14-619-2018 https://doi.org/10.5194/cp-14-619-2018 https://doi.org/10.5194/tc-12-3097-2018 https://doi.org/10.5194/cp-3-15-2007 https://doi.org/10.5194/cp-10-269-2014 https://doi.org/10.5194/cp-10-269-2014 https://doi.org/10.5194/tc-7-469-2013 https://doi.org/10.5194/tc-7-469-2013 https://doi.org/10.5194/tcd-1-41-2007 https://doi.org/10.5194/tc-8-1509-2014 https://doi.org/10.5194/tc-12-2211-2018 https://doi.org/10.5194/cp-10-887-2014 1922 I. Tabone et al.: The NEGIS glacial retreat triggered by the ocean Lecavalier, B. S., Milne, G. A., Simpson, M. J., Wake, L., Huy- brechts, 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. Margo, P. M.: Constraints on the magnitude and patterns of ocean cooling at the Last Glacial Maximum, Nat. Geosci., 2, 127–132, 2009. Marshall, S. J. and Koutnik, M.: Ice sheet action versus reac- tion: Distinguishing between Heinrich events and Dansgaard- Oeschger cycles in the North Atlantic, Paleoceanography, 21, PA2021, https://doi.org/10.1029/2005PA001247, 2006. Marshall, S. J., Tarasov, L., Clarke, G., and Peltier, W. R.: Glacio- logical reconstruction of the Laurentide Ice Sheet: physical pro- cesses and modelling challenges, Can. J. Earth Sci., 37, 769–793, 2000. 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. 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. Mikkelsen, T. B., Grinsted, A., and Ditlevsen, P.: Influence of temperature fluctuations on equilibrium ice sheet volume, The Cryosphere, 12, 39–47, https://doi.org/10.5194/tc-12-39-2018, 2018. Montoya, M. and Levermann, A.: Surface wind-stress threshold for glacial Atlantic overturning, Geophys. Res. Lett., 35, L03608, https://doi.org/10.1029/2007GL032560, 2008. 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., Jakob- sson, M., Jordan, T. M., Kjeldsen, K. K., Millan, R., Mayer, L., Mouginot, J., Noël, B. P. Y., O’Cofaigh, C., Palmer, S., Rys- gaard, 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 sound- ing combined with mass conservation, Geophys. Res. Lett., 44, 11–051, 2017. 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. 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. Münchow, A., Padman, L., and Fricker, H. A.: Interannual changes of the floating ice shelf of Petermann Gletscher, North Green- land, from 2000 to 2012, J. Glaciol., 60, 489–499, 2014. NEEM: Eemian interglacial reconstructed from a Greenland folded ice core, Nature, 493, 489–494, 2013. Peano, D., Colleoni, F., Quiquet, A., and Masina, S.: Ice flux evolu- tion in fast flowing areas of the Greenland ice sheet over the 20th and 21st centuries, J. Glaciol., 63, 499–513, 2017. 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, https://doi.org/10.5194/cp-3-375-2007, 2007. 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. Poirier, R., Cronin, T., Briggs, W., and Lockwood, R.: Central Arc- tic paleoceanography for the last 50 kyr based on ostracode fau- nal assemblages, Mar. Micropaleontol., 88, 65–76, 2012. Price, S. F., Payne, A. J., Howat, I. M., and Smith, B. E.: Commit- ted 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. Rasmussen, T. L., Thomsen, E., and Nielsen, T.: Water mass ex- change between the Nordic seas and the Arctic Ocean on mil- lennial timescale during MIS 4-MIS 2, Geochem. Geophys. Geosyst., 15, 530–544, 2014. Rathmann, N., Hvidberg, C., Solgaard, A., Grinsted, A., Gud- mundsson, 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. Reeh, N.: Parameterization of melt rate and surface temperature on the Greenland ice sheet, Polarforschung, 59, 113–128, 1989. Rignot, E. and Jacobs, S. S.: Rapid bottom melting widespread near Antarctic Ice Sheet grounding lines, Science, 296, 2020–2023, 2002. Rignot, E. and Mouginot, J.: Ice flow in Greenland for the interna- tional polar year 2008–2009, Geophys. Res. Lett., 39, L11501, https://doi.org/10.1029/2012GL051634, 2012. Rignot, E. and Steffen, K.: Channelized bottom melting and sta- bility of floating ice shelves, Geophys. Res. Lett., 35, L02503, https://doi.org/10.1029/2007GL031765, 2008. Rignot, E., Jacobs, S., Mouginot, J., and Scheuchl, B.: Ice-shelf melting around Antarctica, Science, 341, 266–270, 2013. 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. 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. 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. 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. 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. The Cryosphere, 13, 1911–1923, 2019 www.the-cryosphere.net/13/1911/2019/ https://doi.org/10.1029/2005PA001247 https://doi.org/10.5194/tc-12-39-2018 https://doi.org/10.1029/2007GL032560 https://doi.org/10.5194/cp-3-375-2007 https://doi.org/10.1029/2012GL051634 https://doi.org/10.1029/2007GL031765 I. Tabone et al.: The NEGIS glacial retreat triggered by the ocean 1923 Seddik, H., Greve, R., Zwinger, T., Gillet-Chaulet, F., and Gagliar- dini, 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. Simpson, M. J., Milne, G. A., Huybrechts, P., and Long, A. J.: Cali- brating 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. Straneo, F. and Heimbach, P.: North Atlantic warming and the retreat of Greenland’s outlet glaciers, Nature, 504, 36–43, https://doi.org/10.1038/nature12854, 2013. Sztybor, K. and Rasmussen, T.: Late glacial and deglacial palaeo- ceanographic changes at Vestnesa Ridge, Fram Strait: Methane seep versus non-seep environments, Palaeogeogr. Palaeoclima- tol. Palaeoecol., 476, 77–89, 2017. 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, https://doi.org/10.5194/cp-14-455-2018, 2018. 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., Kipfs- tuhl, 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, https://doi.org/10.5194/tc-8-1275-2014, 2014. 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. Weidick, A., Andreasen, C., Oerter, H., and Reeh, N.: Neoglacial glacier changes around Storstrommen, North-East Greenland, Polarforschung, 64, 95–108, 1996. Werner, K., Spielhagen, R. F., Bauch, D., Hass, H. C., and Kandi- ano, 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. 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. Wilson, N., Straneo, F., and Heimbach, P.: Satellite-derived subma- rine melt rates and mass balance (2011–2015) for Greenland’s largest remaining ice tongues, The Cryosphere, 11, 2773–2782, https://doi.org/10.5194/tc-11-2773-2017, 2017. Winkelmann, D., Jokat, W., Jensen, L., and Schenke, H.-W.: Sub- marine end moraines on the continental shelf off NE Greenland– Implications for Lateglacial dynamics, Quaternary Sci. Rev., 29, 1069–1077, 2010. 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, https://doi.org/10.5194/tc-5-715-2011, 2011. Zweck, C. and Huybrechts, P.: Modeling of the northern hemi- sphere ice sheets during the last glacial cycle and glacio- logical sensitivity, J. Geophys. Res.-Atmos., 110, D07103, https://doi.org/10.1029/2004JD005489, 2005. www.the-cryosphere.net/13/1911/2019/ The Cryosphere, 13, 1911–1923, 2019 https://doi.org/10.1038/nature12854 https://doi.org/10.5194/cp-14-455-2018 https://doi.org/10.5194/tc-8-1275-2014 https://doi.org/10.5194/tc-11-2773-2017 https://doi.org/10.5194/tc-5-715-2011 https://doi.org/10.1029/2004JD005489 Abstract Introduction Methods Results Discussion Comparison between modelled and data-derived reconstructions Oceanic forcing at orbital timescales Model performance on the whole GrIS Future perspectives Conclusions Data availability Supplement Author contributions Competing interests Acknowledgements Financial support Review statement References