Articles | Volume 15, issue 2
Research article
28 Mar 2019
Research article |  | 28 Mar 2019

Impact of millennial-scale oceanic variability on the Greenland ice-sheet evolution throughout the last glacial period

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

Temperature reconstructions from Greenland ice-sheet (GrIS) ice cores indicate the occurrence of more than 20 abrupt warmings during the last glacial period (LGP) known as Dansgaard-Oeschger (D-O) events. Although their ultimate cause is still debated, evidence from both proxy data and modelling studies robustly links these to reorganisations of the Atlantic Meridional Overturning Circulation (AMOC). During the LGP, the GrIS expanded as far as the continental shelf break and was thus more directly exposed to oceanic changes than in the present. Therefore oceanic temperature fluctuations on millennial timescales could have had a non-negligible impact on the GrIS. Here we assess the effect of millennial-scale oceanic variability on the GrIS evolution from the last interglacial to the present day. To do so, we use a three-dimensional hybrid ice-sheet–shelf model forced by subsurface oceanic temperature fluctuations, assumed to increase during D-O stadials and decrease during D-O interstadials. Since in our model the atmospheric forcing follows orbital variations only, the increase in total melting at millennial timescales is a direct result of an increase in basal melting. We show that the GrIS evolution during the LGP could have been strongly influenced by oceanic changes on millennial timescales, leading to oceanically induced ice-volume contributions above 1 m sea level equivalent (SLE). Also, our results suggest that the increased flux of GrIS icebergs as inferred from North Atlantic proxy records could have been triggered, or intensified, by peaks in melting at the base of the ice shelves resulting from increasing subsurface oceanic temperatures during D-O stadials. Several regions across the GrIS could thus have been responsible for ice mass discharge during D-O events, opening the possibility of a non-negligible role of the GrIS in oceanic reorganisations throughout the LGP.

1 Introduction

Ice cores from the Greenland ice sheet (GrIS) indicate that the climate of Greenland experienced multiple abrupt temperature increases during the last glacial period (LGP) known as Dansgaard-Oeschger (D-O) events (Dansgaard et al.1993). D-Os have waiting times between consecutive events of around 1500 years and, with decreasing probability, 3000 and 4500 years (Alley et al.2001) and are characterised by an initial abrupt warming of 10–16 C (interstadial state) followed by a slow cooling which finally ends with an abrupt return to glacial background conditions (stadial state) (Johnsen et al.1992; Kindler et al.2014). These millennial-scale climate fluctuations were first observed in Greenland ice cores (Anklin et al.1993; Grootes et al.1993), but strong evidence is also found in numerous marine sediments records, especially in the Nordic Seas and North Atlantic (Bond et al.1993; Rasmussen et al.1996; Voelker et al.1998), tropical and subtropical stalagmite proxy archives (Fleitmann et al.2009; Wang et al.2001), and in Antarctic ice-sheet (AIS) ice cores (Barbante et al.2006), suggesting a worldwide imprint of D-O events (Voelker2002).

The numerous and extensive evidence from paleo records (e.g. Shackleton et al.2000; Böhm et al.2015; Gottschalk et al.2015; Henry et al.2016) and from results from a variety of models, ranging from Earth system models of intermediate complexity (EMICs) to comprehensive general circulation models (GCMs) (e.g. Ganopolski and Rahmstorf2001; Shaffer et al.2004; Peltier and Vettoretti2014; Vettoretti and Peltier2016, 2018) widely attribute the cause of D-O events to reorganisations of the Atlantic Mean Overturning Circulation (AMOC), which are likely modulated by variations in the oceanic convection over the Nordic Seas and North Atlantic (Dokken and Jansen1999; Rasmussen and Thomsen2004; Dokken et al.2013; Barker et al.2015; Hoff et al.2016; Rasmussen et al.2016); a recent review on the available evidence was provided by Lynch-Stieglitz (2017). However, the final cause of these AMOC reorganisations is still unclear. Numerous modelling studies have explored the problem. Much work attributes AMOC instability to freshwater discharge from the Northern Hemisphere (NH) ice sheets (Ganopolski and Rahmstorf2001; Vellinga and Wood2002; Menviel et al.2014; Bagniewski et al.2017) directly connected to changes in the strength (Skinner and Elderfield2007) and in the location (Sévellec and Fedorov2015) of deep convection. Other possible mechanisms link the origin of D-O events to sea-ice cover variability (Li et al.2005, 2010; Sime et al.2019) or to linked sea-ice–ice-shelf fluctuations (Boers et al.2018; Petersen et al.2013). Still others connect AMOC reorganisations to climatic perturbations in the atmosphere associated with changes in ice-sheet dynamics (Wunsch2006; Zhang et al.2014a), to progressive CO2 atmospheric variations (Zhang et al.2017), to changes in atmospheric heat transport (Wang et al.2015), or to combined changes in wind and atmospheric CO2 concentrations driven by the Southern Ocean (Banderas et al.2015). Recently, D-O events have been explained as the result of a non-linear internal salt oscillator in the Atlantic (Peltier and Vettoretti2014; Vettoretti and Peltier2016, 2018) without the need to invoke any external forcing.

Despite the evidence for millennial-scale climate variability during the LGP in Greenland, the specific role of the GrIS in this framework has not been investigated in depth. Whether the GrIS was an active or a passive player in this oscillatory system is unknown and its evolution between the last interglacial (LIG) and the Last Glacial Maximum (LGM) is still debated (Vasskog et al.2015). Evidence of ice-rafted debris (IRD) deposition in marine sediment cores suggests that the GrIS could have been a non-negligible source of iceberg discharge during the LGP (Jonkers et al.2010). Cores drilled in the Irminger Sea (Bond and Lotti1995; Van Kreveld et al.2000), Denmark Strait (Voelker and Haflidason2015; Andrews et al.2017), close to Scoresby Sund (Stein et al.1996) and in the North Atlantic close to the Reykjanes Ridge (Barker et al.2015; Hodell et al.2010; Rasmussen et al.2016) show evidence of iceberg transport sourced in East Greenland and Northeast Greenland, among other origins; also, sediments in Fram Strait show imprints of iceberg discharge from the northern GrIS (Darby et al.2002). Still, other work links the IRD deposition found close to the east–southeast marine margin of the GrIS to local ice-sheet instability (Verplanck et al.2009) and associated sediments found in the Labrador Sea to iceberg discharge coming from Baffin Bay during the Younger Dryas (Andrews et al.2012). Also, other work based on paleo reconstructions explicitly affirms that ocean temperature increase played a fundamental role in the GrIS retreat in the last deglaciation (Jennings et al.2017). Recent proxy data collected by morain-derived marine shells suggest that the margin of the Northeast Greenland Ice Stream (NEGIS) region may have fluctuated throughout the LGP by more than 200 km (Larsen et al.2018). All these examples suggest that the GrIS may have experienced substantial variability during the LGP and beyond.

From a modelling perspective, only a few studies have tackled the response of the GrIS to millennial-scale climate variability throughout the LGP. Charbit et al. (2007) simulated the evolution of NH ice sheets by forcing their ice-sheet model with a glacial–interglacial climate anomaly scaled by a climatic index taken from the GRIP reconstruction, which retained millennial-scale temperature fluctuations. However, their model accounted only for grounded ice; ice shelves were not included and the effect of the ocean was thus not taken into account. Huybrechts (2002) also investigated the response of the GrIS to millennial-scale climate variability during the LGP with an ice-sheet model. However, the GrIS extent was controlled by orbital-only variations in past eustatic sea level, and millennial-scale fluctuations, both in sea level and in the ocean temperature, were omitted. Marshall and Koutnik (2006) assessed the response of the NH ice sheets to climate variability at millennial timescales. They simulated the GrIS evolution using an ice-sheet model that included a calving parameterisation and that was forced by millennial-scale temperature variations. Their results showed very little response of the GrIS to the imposed climate variability, with a weak increase in iceberg flux only during interstadials. However, this study did not explicitly investigate the effect of the millennial-scale variability in the ocean either.

A recent study has demonstrated the important role of oceanic conditions in the evolution of the GrIS throughout the past two glacial cycles (Tabone et al.2018). This study, however, focused on orbital oceanic variations. Thus, the effect of oceanic millennial-scale variations on the GrIS remains to be explored. Sea surface temperatures (SSTs) inferred from planktonic foraminifera in North Atlantic sediment cores typically vary from 7–8.5 C during interstadials to 3–3.5 C during stadials (Rasmussen et al.2016). Recent SST estimates from a stack of sediment cores of the North Atlantic suggest a broad range of interstadial–stadial temperatures fluctuating between a maximum of 15 C and a minimum of 0 C, but typical SST anomalies in each location are lower (2–6 C) (Jensen et al.2018). Moreover, oceanic temperatures reconstructed from sediment cores in Denmark Strait (Sessford et al.2018) and the Nordic Seas (Ezat et al.2014) suggest stadial–interstadial subsurface temperature anomalies of 3–4 and 2–5 C, respectively, during the central part of the LGP. Oceanic temperature variations associated with D-O events may thus have an appreciable impact on the GrIS.

Here we use a three-dimensional, thermomechanical, hybrid ice-sheet–shelf model to investigate this issue. Stadial–interstadial temperature variations in the ocean are expressed in the model through a linear parameterisation of the submarine melting rate at the grounding line and below the ice shelves. We characterise the millennial-scale variability in the ocean by performing two large ensembles of simulations: one forced by both orbital and millennial-scale components in the oceanic temperature signal and the other forced solely by orbital frequencies. In both ensembles, the atmospheric forcing only includes orbital changes; thus comparison between the two large ensembles returns the impact of millennial-driven oceanic fluctuations on the GrIS evolution. Moreover, we investigate the impact of the oceanic-induced perturbations at the marine margin on the GrIS interior. This gives further insight into the potential role of the ocean in dynamic reorganisations of the ice sheet and enhanced ice discharge.

The paper is structured as follows: the ice-sheet–shelf model, the forcing methods and the methodology used to study the sensitivity of the GrIS to millennial fluctuations in the submarine melting are described in Sect. 2; in Sect. 3, we first characterise the millennial-scale variability effect in the ocean on the GrIS evolution; then we describe the impact of this oceanic oscillation on the GrIS transient dynamics. We then discuss results and caveats of the model, and this is followed by the conclusions.

2 Tools and methods

2.1 GRISLI-UCM ice-sheet–shelf model

Here we use the three-dimensional, thermo-mechanical, hybrid ice-sheet–shelf model GRISLI-UCM (Alvarez-Solas et al.2018; Tabone et al.2018), developed from the GRISLI ice-sheet model (Ritz et al.2001), which has been substantially modified through diverse parameterisations of ice-sheet boundary conditions (e.g. surface mass balance, basal sliding and basal melt). The model solves deformational-flow-driven grounded areas via the shallow-ice approximation (SIA; Hutter1983) and floating shelves predominantly moving under plug flow via the shallow-shelf approximation (SSA; Morland et al.1987). Transitional areas where these two flow regimes coexist are solved by adding the velocities of the SIA and SSA solutions (Winkelmann et al.2011). The ice base of SIA-solved areas is supposed to be in a cold phase (below the pressure melting point) and does not admit sliding. By contrast, transitional zones conditionally allow for basal sliding as long as the ice base is temperate (at the pressure melting point) and the effective pressure at the bed is below a critical threshold. Once this condition is satisfied, the basal shear traction is calculated to be proportional to the ice-base velocity and to a coefficient depending on the basal water pressure and the bed roughness. Glacial isostatic adjustment (GIA) of the bedrock related to changes in the ice pressure is calculated through the Elastic Lithosphere Relaxing Asthenosphere model (Le Meur and Huybrechts1996). Grounding-line migration is expressed through a flotation condition between the temporally variable sea level prescribed in the model and the ice thickness calculated at the ice–ocean interface. In this study the position of the grounding line is diagnosed at sub-grid-scale precision by interpolating the ice thickness over the grid cell including the grounding line (adapted from Gladstone et al.2010). The calculation returns the grounded percentage of the grid cell and avoids the simplistic assumption that the grounding line lies on the last grounded coarse-grid point before floating ice.

The total ice mass balance is defined at any time as the difference between the mass balance at the ice surface (accumulation minus ablation), basal melting at the ice base and ice calving. Surface ablation is calculated by the positive degree day (PDD) scheme (Reeh1989). This scheme is known to be overly simplistic for both ice-sheet models (Robinson and Goelzer2014) and EMICs (Bauer and Ganopolski2017) in the paleo context, as it does not incorporate the effect of incoming solar radiation changes. Nevertheless, since here we focus on the sensitivity of the ice sheet to millennial-scale oceanic variations during the LGP, the choice of this scheme should be sufficient for our purposes. Surface precipitation is exponentially proportional to atmospheric temperatures, which vary through an index approach (Banderas et al.2018; Blasco et al.2019; Tabone et al.2018) (Sect. 2.2). Geothermal heat flux exchanged at the base of grounded ice is taken from Shapiro and Ritzwoller (2004). Basal melting at the grounding line and at the ice-shelf base is described in detail in Sect. 2.3. Calving is calculated following a two-rule criterion: the ice front must not exceed a critical thickness (H=200 m) and ice advected from upstream must fail to maintain the ice thickness above that imposed threshold (Peyaud et al.2007; Colleoni et al.2014).

The GRISLI-UCM model is applied here to the GrIS domain, with a spatial resolution of 20 km by 20 km. Since the goal of this work is to investigate the sensitivity of the GrIS to past millennial-scale variability in the ocean, we force the ice-sheet model through spatially variable surface atmospheric temperatures and precipitation that reflect orbital variations only. This allows us to study the direct effect of millennial-scale fluctuations on the GrIS evolution due to the ocean only, unperturbed by the millennial-scale atmospheric variability.

2.2 Atmospheric forcing

As in Tabone et al. (2018), the atmospheric temperature forcing applied to the model follows an index-anomaly method in which the present-day temperature Tatmclim is modified by fluctuations of past temperature anomalies ΔTatm(t):

(1) T atm ( t ) = T atm clim + Δ T atm ( t ) ,


(2) Δ T atm ( t ) = ( 1 - α ( t ) ) Δ T atm orb .

The present-day temperature climatology Tatmclim is obtained from the regional climate model MAR forced by the ERA-Interim reanalysis (Fettweis et al.2013) for the period 1981–2010. The resulting atmospheric temperature Tatm(t) is a 2-D field which varies in time according to α(t), which is a spatially uniform climatic index derived from a composite series of various proxy-derived temperature anomalies from the LIG to the present day (Tabone et al.2018) (Fig. 1a). This temperature anomaly signal is then smoothed to remove the spectral components below the orbital frequencies (1/f<18 ka) and is normalised between 0 and 1 to build the index α (α=1 at present day (PD) and α=0 at the LGM) (Fig. 1b). Thus, by construction, ΔTatmorb is the glacial–interglacial atmospheric temperature anomaly, which is here taken from glacial minus present-day anomaly snapshots of the EMIC CLIMBER-3α model (Montoya and Levermann2008). This orbital forcing method is likewise applied to the precipitation field, which is parameterised through the ratio of glacial and present-day precipitation anomalies δPannorb, as in Tabone et al. (2018):

(3) P ann ( t ) = P ann clim α ( t ) + ( 1 - α ( t ) ) δ P ann orb .

Figure 1North Greenland Ice Core (NGRIP)-temperature anomaly reconstruction (a) used to construct the orbital-index α(t) (b) and the millennial-index β(t) (c) for the last glacial period. α(t) is built from the NGRIP-derived temperature anomalies (Tabone et al.2018), filtered to remove the spectral components below the orbital period (1/f<18 ka) and normalised between 0 and 1. β(t) is derived by subtracting the orbital α(t) index from the NGRIP-derived temperature signal, normalised between 0 and 1 and then filtered below 1 ka to remove sub-millennial periods.


2.3 Oceanic forcing

The methodology used to force the ocean is similar to that of Tabone et al. (2018), except that here we include the millennial-scale variability in the ocean which was omitted in the previous work. This approach is analogous to that of Blasco et al. (2019) for the Antarctic domain. Here, the basal melting rate at the grounding line is parameterised as

(4) B ( t ) = κ ( T ocn clim - T f ) + κ Δ T ocn ( t ) ,


(5) Δ T ocn ( t ) = ( 1 - α ( t ) ) Δ T ocn orb + β Δ T ocn mil .

The formulation is adapted from Beckmann and Goosse (2003) and permits the translation of changes in ocean temperature into changes in basal melt through an ocean–ice heat-flux exchange scaling factor κ (ma-1K-1). The first term of Eq. (4) is the prescribed basal melting rate for the present day, where Tf is the temperature at the base of the ice shelf, assumed to be at the freezing point; Tocn is the temperature of the ocean at the present day; and Tocnclim is the mean climatological temperature of the ocean. Since the present-day basal melt rate for the whole Greenland domain is largely unconstrained, we assume it to be spatially constant (hereafter referred to as Bref=κ(Tocnclim-Tf) (m a−1)). The reference basal melt is then perturbed by a time-dependent anomaly, which includes both orbital and millennial variability. β (Fig. 1c) is the millennial index, derived from the same normalised temperature anomaly signal as α. Specifically, it is built by subtracting the orbital α index (Fig. 1b) from the unfiltered multi-proxy temperature signal (Fig. 1a) normalised between 0 and 1 and then filtering to remove periods below 1 ka to eliminate sub-millennial components of the signal. ΔTocnorb and ΔTocnmil are the glacial–interglacial and interstadial–stadial oceanic subsurface temperature anomalies (K), respectively. The term (1-α(t))ΔTocnorb reflects the subsurface temperature variations resulting from orbital changes. The term βΔTocnmil expresses the millennial-scale temperature changes at the subsurface assumed to be in antiphase with respect to the Greenland atmospheric temperature inferred from Greenland ice cores (e.g. Johnsen et al.2001 and Kindler et al.2014). Thus, we are assuming that subsurface water temperatures increase during stadials and decrease during interstadials. This is in agreement with the presence of warming subsurface waters during stadial periods as suggested by several records of the North Atlantic and Nordic Seas (Ezat et al.2014; Jonkers et al.2010; Rasmussen and Thomsen2004; Rasmussen et al.2016; Sessford et al.2018; Dokken et al.2013) and supported by modelling work (Vettoretti and Peltier2015; Brady and Otto-Bliesner2011; Mignot et al.2007; Knutti et al.2004; Marcott et al.2011). The result is a submarine melting signal that peaks during D-O stadials. This is consistent with the oceanic forcing of ice-sheet models used previously to inspect the effect of subsurface warming during the coldest stadials, i.e. Heinrich events (Alvarez-Solas et al.2010, 2013; Bassis et al.2017) or, as done recently, to investigate the origin of D-O events through a conceptual model (Boers et al.2018).

As in Tabone et al. (2018) and Blasco et al. (2019), basal melt under the ice shelves is 10 % of that near the grounding line, as supported by recent estimates in Greenland (Rignot and Steffen2008; Münchow et al.2014; Wilson et al.2017). To prevent unconstrained accretion below ice shelves and at the grounding line, negative basal melt rate (freezing) is cut off to 0 m a−1. Changes in global sea level at orbital timescales are prescribed in the model and are used to compute the grounding-line position. The applied sea-level variation time series is taken from Bintanja and Van de Wal (2008). The signal is inferred from a marine benthic oxygen isotope record reconstructed for the last 3 Myr through an ice-sheet model coupled to a simple marine temperature model.

2.4 Perturbed basal melting parameters and LE description

Following the discussion above, we can rewrite the basal melting equation (Eq. 4) as

(6) B ( t ) = B ref + κ ( 1 - α ( t ) ) Δ T ocn orb + β Δ T ocn mil .

Written in this form, the basal melting formulation depends on the choice of four parameter values: Bref, κ, ΔTocnorb and ΔTocnmil. These variables are here all considered as spatially uniform around Greenland for the sake of simplicity, leading to a spatially homogeneous basal melting rate. To assess the GrIS response to millennial-scale variability in the ocean, we could simply consider varying the value of κ, which is the sensitivity of the oceanic forcing (Tabone et al.2018). However, by construction of Eq. (6), increasing κ does not necessarily mean increasing the millennial-scale oceanic forcing alone, since this would enhance concurrently both the millennial and the orbital-scale components in the ocean. Therefore, investigating the oceanic millennial-scale variability effect on the past GrIS is not as straightforward as expected. Moreover, none of the four parameters of Eq. (6) is perfectly constrained in reality, and a sensitivity study on the influence of their chosen values on the GrIS evolution would be required.

For these reasons, it is first useful to characterise the impact of millennial-scale variability in the ocean on the GrIS evolution by testing a broad range of values of the key parameters in Eq. (6): ΔTocnmil, set to vary between 0 and 3 K; Bref, between 0 and 10 m a−1 (chosen as a reasonable climatic mean between Rignot et al.2010, Rignot et al.2016, Straneo et al.2012 and Wilson et al.2017 for the largest tidewater glaciers around the GrIS and Liu et al.2015 and Rignot et al.2013 for Antarctica); and κ, between 0 and 10 ma-1K-1 (following Rignot and Jacobs2002 for Antarctica). The reduction of the problem to 3 degrees of freedom is supported by estimates of subsurface temperatures at both short (millennial) and long (orbital) timescales that have the same amplitudes. Reconstructed LGM–present-day subsurface anomalies, which at orbital timescales are considered to follow those at the surface, are between 0 and 3 K around Greenland (Annan and Hargreaves2013; Margo2009). A similar range of values is found for the interstadial–stadial subsurface temperature anomalies (Alvarez-Solas et al.2018; Brady and Otto-Bliesner2011; Vettoretti and Peltier2015; Zhang et al.2014b). Also, in our experimental set-up the α and β indices share the same normalisation. Thus, interstadial–stadial and glacial–interglacial anomalies share the same amplitude of forcing. Also, it should be noted that Bref depends on κ (Eq. 4); thus any change in κ results in a change in Bref too. However, Bref also depends on the water and freezing-point temperatures at the grounding zone (Tocnclim and Tf, respectively), which are largely unconstrained, since they mostly depend on the characteristics of the considered grounding line (e.g. on its depth; Beckmann and Goosse2003). To take into account the uncertainty of these quantities, here, Tocnclim-Tf is treated as an additional degree of freedom that, although not explicitly perturbed in the equation, allows us to change Bref independently from κ. Of course, decoupling these two variables requires additional attention in considering a range of values for Bref and κ that allow for realistic Tocnclim-Tf. The simulations that do not fulfill this requirement will be discarded a posteriori in the analysis.

Table 1Range of parameter values sampled by the LHS technique for generating the TOT and ORB ensembles.

Download Print Version | Download XLSX

To test a wide range of combinations between the three key parameters, we perform a large ensemble (LE) of model simulations using the near-random Latin hypercube sampling (LHS) technique (McKay et al.1979), which allows us to efficiently explore the phase space of the key parameters minimising the LE computational cost with respect to the full-factorial sampling technique. The LHS method has already been used to constrain different ice-sheet model parameters and to assess their influence on the model's behavior (Stone et al.2010; Applegate et al.2012; Stone et al.2013; Robinson et al.2017). The parameter values are sampled from the specified ranges and, assuming, as discussed, that they are independent of each other, they are randomly combined to generate a total LE of 100 simulations, named TOT simulations. At the same time, we perform another set of identical simulations, except for the fact that the climatic index associated with the millennial-scale variability (β) is set to zero. These are named ORB simulations and are used for direct comparison with the TOT simulations, as discussed in Sect. 3 (see Table 1 for a full list of the phase space of parameters investigated within the two LEs). To initialise the model we use the present-day topography of Greenland from Schaffer et al. (2016). All the simulations of the LE cover the last two glacial cycles, with the first considered as a spin-up and therefore not analysed.

Figure 2(a) Ice-volume evolutions simulated by ORB (black curves) and TOT (red curves) throughout the LGP (m SLE); (b) TOT-ORB ice-volume anomaly, intended as the oceanic millennial-scale variability contribution to the GrIS ice-volume variation (m SLE; positive values indicate ice gained by TOT with respect to ORB). The colour scale indicates the root mean square deviation (ε) calculated in terms of ice-volume residuals between TOT and ORB simulations (m SLE). Black line in (b) refers to the simulation corresponding to the maximum GrIS response to millennial-scale variability in the ocean.


3 Results

3.1 Characterisation of oceanic millennial-scale variability effect

A first evaluation of the effect of millennial-scale oceanic variability on the whole GrIS can be made by looking at the evolution of its simulated ice volume throughout the LGP (Fig. 2a). Apart from low oceanic forcing configurations for which both ORB (black curves) and TOT (red curves) present small ice-volume deviations from their simulated present-day value, TOT simulations show higher ice-volume fluctuations than ORB ones, especially during Marine Isotope Stage 3 (MIS3, ca. 60–25 ka). A better representation of the impact of millennial-scale variability on the GrIS evolution can be obtained by considering deviations in ice volume from the orbital-only forced reference simulation (ORB) (Fig. 2b). Those anomalies can reach peaks of almost 2 m SLE at certain times during the LGP, with changing rates of more than 1 m SLE in less than 2 kyr. However, isolated peaks in volume anomalies do not necessarily imply a high millennial-scale ice-volume variation over the whole LGP since they may result from sparse melting peaks, and the contribution of the oceanic variability to ice-volume changes needs to be quantified for the entire time series. This can be done by calculating the root mean square deviation (here specified by the greek letter ε) of the ice-volume (V) residuals between TOT and ORB averaged in time for each simulation of the LE (m SLE) as

(7) ε V = 1 N t - 1 t = 1 N t V TOT ( t ) - V ORB ( t ) 2 .

Nt is the total number of time steps in each simulation and VTOT(t)−VORB(t) is the ice-volume residual calculated as the difference between the two simulations at each time step. This quantity tells us how much the volume in the TOT simulation, including millennial-scale variability, deviates from that in the background throughout the whole LGP. The larger its value, the larger the impact of oceanic millennial-scale variability on the GrIS.

A second and complementary approach is to calculate the deviation of the ice speed (U) simulated by TOT from its background state ORB. The methodology is similar to that used for the ice volume, except for the fact that now ε is first calculated for each grid point ij as εU,ij=1Nt-1t=1Nt(Uij,TOT(t)-Uij,ORB(t))2 and then averaged for the entire domain:

(8) ε U = i j ε i j 2 N ,

where N is the total number of grid points of the GrIS domain.

Ice-volume evolution is a good benchmark for understanding the overall effect of the oceanic variability throughout the time, but U is the direct expression of how the ice dynamics are affected by the forcing. Therefore, we jointly examine the ε calculated for ice volume and velocities. The higher the millennial-scale variability in the ocean is, the more the ice volume and ice velocity in TOT are expected to deviate from the ORB simulation. Specifically, we are interested in understanding which combination of perturbed parameters in the basal melting equation (Eq. 6) leads to the highest millennial variability in terms of both quantities. We find that the effect of a high oceanic millennial-scale variability is fairly well constrained in the parameter phase space, at least per each couple of analysed variables (Fig. 3). Very low oceanic-driven variability is generally associated with low values of κ and ΔTocn (weak oceanic forcing). These configurations lead to basal melting evolutions that do not deviate much from the reference-state basal melting throughout the LGP (blue curves of Fig. S1 in the Supplement). Low millennial-scale variability is also found for combinations of low Bref with high κ (or high ΔTocnmil) and for high κ with high ΔTocnmil. The reason is that these are associated with a basal-melt evolution that rapidly saturates at the cut-off value 0 m a−1 (no melting, since freezing is not allowed) after the LIG and maintains that state throughout the LGP (light blue curves in Fig. S1). It is between these two configurations that the resulting basal melting signal allows for a high millennial-scale variability in the ocean. This is usually ensured by medium–high Bref and high κ, but also by combinations of high Bref and relatively low κ (red curves showing a positive volume anomaly during MIS3 in Fig. 2b). In these cases, the high ε is a result of the large ice-volume variation between a GrIS configuration in ORB that is confined inland by a persistent high melting during the whole LGP and a less confined one in TOT as a result of a basal melting that saturates at zero at certain intervals. In any case, the highest millennial variability in the ocean is associated with a signal that contains sufficiently high melting peaks at millennial timescales and that sometimes saturates at 0 m a−1 (red curves in Fig. S1). A similar logic can be followed considering ΔTocnmil instead of κ. High κ must be associated with low ΔTocnmil and vice versa to produce high millennial-scale variability in the basal melt, as reflected in Fig. 3 as well.

Figure 3Evaluation of the effect of the oceanic millennial-scale variability for each pair of parameters perturbed in the basal melting equation (Eq. 6). The ranges of εV (m SLE, colour scale) from Eq. (7) and εU (m a−1, circle size) from Eq. (8) are shown. The simulation with the highest millennial-scale variability is represented by the largest red circle; the simulation with the lowest variability is represented by the smallest blue circle.


3.2 Oceanic millennial-scale variability impact on transient dynamics

The aim of this section is to investigate the effect of oceanic millennial-scale variability on the transient GrIS dynamics and ice-sheet evolution during the LGP. To this end, we choose the simulation from the LE which corresponds to the maximum millennial-scale variability response in terms of both ice volume and velocity and we compare it to its corresponding background simulation. Specifically, we associate the largest millennial-scale variability in the ocean to the maximum value given by the product of the two scores (εTOT=εVεU). As expected from the characterisation of the millennial variability discussed above, the chosen simulation has medium–high Bref (7.6 m a−1), high κ (8.3 ma-1K-1) and consequently medium ΔTocnmil (1.2 K) values. This experiment shows a millennial-scale contribution to ice-volume variation of more than 1 m SLE during the highest peaks of melting (black curve of Fig. 2b). From now on, this simulation and its corresponding orbital-only reference simulation are referred to as TOTmax and ORBmax, respectively.

We compare ice thickness, velocity and migration of the grounding-line position simulated by TOTmax with those of ORBmax for the whole LGP (movie in Supplement). Since the millennial-scale variability in the ocean is built to represent temperature variations at the subsurface (Fig. 4a), each submarine melting peak is associated with a stadial interval occurring between two near D-O events, also referred to as Greenland Stadial (GS) (Rasmussen et al.2014). Here we show the comparison specifically for GS 15 (ca. 55 ka), that is linked to one of the most pronounced basal melting peaks during MIS3 (Fig. 4). Note that when a comparison is made between two GrIS extensions that show a mismatch in ice cover, the resulting velocity difference refers to the ice velocity of the simulation in which the grounding line has not retreated yet. In TOTmax the high melting peak leads to a large ice-thickness decrease, associated with a substantial grounding-line retreat, especially in the Kangerlussuaq and Sermilik fjords (SE Greenland) and the NEGIS region (Fig. 5). The area near Baffin Bay (NW Greenland) is also affected by inland grounding-line migration, but the retreat is limited if compared to the other zones. The outlet glacier at Jakobshavn Isbrae (W Greenland) also strongly responds to the forcing, but the response, compared to the ORBmax simulation, decreases throughout the melting peak. This is related to the large ice thickness simulated by ORBmax in the 20 kyr previous to GS 15, which increases the amount of glacial water produced at the base of the ice stream, eventually enhancing basal sliding and ice discharge when a persistent submarine melting is applied. On the other hand, the inclusion of millennial-scale oceanic variability allows for frequent ice discharge, limiting the basal-water amount and the ice stream basal velocity, leading to discharge ice more promptly at the beginning of the peak. During the slow subsurface cooling after the GS melting peak (at ca. 52 ka), the grounding line is then free to readvance.

Figure 4Evolution of (a) millennial-scale subsurface oceanic temperature as defined in Eq. (6) to force the TOTmax simulation and (b) the resulting submarine melting used in both TOTmax (solid line) and ORBmax (dashed line) simulations for the last 100 kyr. Submarine melt peaks are associated with their corresponding Greenland stadials (GSs).


Figure 5TOTmax ORBmax ice thickness (m, a, b, c) and velocities (m a−1, d, e, f) during submarine melting rate peak corresponding to GS 15. TOTmax and ORBmax grounding-line positions are shown by solid and dashed lines, respectively.


Simulations show that increased velocities in the NEGIS and Kangerdlugssuaq regions associated with enhanced ice discharge are found also during GS 12 (ca. 43 ka) (movie in Supplement). Appreciable but less pronounced grounding-line variations are likewise observed during GS 13 (ca. 47 ka) and GS 16 (ca. 56 ka), especially in SE Greenland. In GS 16 the weak response to the submarine melting peak is related to a small grounding-line migration in the south. Around GS 13, the pattern is even reversed, since more ice is lost in ORBmax than in TOTmax. This is due to the presence of the persistent melting previous to GS 13 that has already led the SE Greenland to be more retreated and to discharge more ice.

Generally, ORBmax is much more static than TOTmax, showing smaller grounding-line and ice-discharge changes. Evidence from marine records suggest that stadial–interstadial oceanic temperature variations of about 2–5 K at intermediate depths (Ezat et al.2014; Sessford et al.2018) likely occurred around Greenland during D-O events. Even larger anomalies are found at the surface (e.g. Jensen et al.2018; Rasmussen et al.2016). It is reasonable that such consistent temperature fluctuations were associated with large variations in melting at the base of ice shelves present in Greenland during the last glacial, when grounded ice advanced over the continental shelf. Therefore, a more variable submarine melt configuration such as that of TOTmax is expected to better reflect the oceanic temperature variation found during D-O events. However, there is no certainty as to whether more ice variability associated with this melting fluctuation is actually a better representation of reality than the ORBmax case.

Figure 6Transient dynamics specified for three locations along Baffin Bay. Ice thickness (a, b, c) and velocity (d, e, f) LGP evolution is shown for the glacial marine margin (point A), present-day marine margin (point B) and far in the interior of the ice sheet (point C). Black curves indicate the dynamics of the ORBmax simulation, red curves represent the dynamics of the TOTmax simulation. The green line on the map on the left-hand side shows the maximum GrIS glacial extent simulated in TOTmax.


The inclusion of millennial-scale variability in the forcing results in ice-thickness variations of hundreds of meters at the marine margin when compared to the corresponding orbital-only case. These ocean-driven fluctuations in ice mass are first exhibited at the ice–ocean front but penetrate through the ice-sheet interior by several tens of kilometers. This is related to the propagation of velocities along fast-flowing ice streams, which, especially around the west GrIS, Sermilik, Kangerlussuaq and NEGIS regions, can penetrate far inland promoting significant ice discharge into the ocean. This effect can be seen by analysing the LGP ice-thickness and velocity evolution for a single grid point of the domain, at three strategic locations: close to the glacial marine margin (A), close to the PD marine margin (B, glacial ice-sheet interior) and far inland in the ice sheet (C). We show the results of this analysis explicitly for the Baffin Bay region (Fig. 6), but another three areas around the Greenland domain have been investigated (Fig. S2 for Jakobshavn Isbrae, Fig. S3 for the NEGIS and Fig. S4 for Kangerlussuaq fjord). Around Baffin Bay the highest differences in transient dynamics between ORBmax and TOTmax are found at the beginning of the LGP. For almost 30 kyr, the margin of ORBmax is constrained close to the continental border, whereas intermittent periods of submarine melt in TOTmax allow the ice edge to rapidly advance and retreat, showing high-velocity fluctuations (600–800 m a−1) (point A). A similar response is found in the proximity of the PD marine margin, about 160 km far from the glacial ice–ocean border (point B), where the effect of the millennial-scale variability in the ocean is found in ice-thickness variations of 500 m. This effect is diffused inland, still reaching differences of 200 m 300 km away from the ocean border (point C), despite minor or no observable difference in velocities. These estimates further corroborate the results found for the Baffin Bay region in the 2-D plots (movie in Supplement). The region including the Jakobshavn Isbrae (Fig. S2) seems to be slightly less affected by the millennial fluctuation in basal melting than the region of Baffin Bay. Still, appreciable ice-thickness variations are found close to the marine margin, especially at the beginning of the glacial period and around 75–50 ka, when the ice shelf repeatedly disappears in the presence of rapid melting peaks. Little evidence of this ice fluctuation is found in the ice interior, however, although slightly increased velocities with respect to ORBmax around 60 ka limit ice growth, culminating in ice ∼700 m lower than that of the ORBmax (point B). The oceanically induced variability rapidly decays further inland (point C). However, it must be noted that this latter location is more than 600 km away from the marine margin; thus this attenuation is expected. Millennial-scale fluctuations in submarine melting also impact the NEGIS glacial evolution (Fig. S3). Intermittent periods of no melting at the beginning of the glacial period allow the ice margin to strongly fluctuate back and forth over the continental shelf (point B). By contrast, persistent basal melting in ORBmax associated with very high velocity (higher than 3000 m a−1) constrains ice expanding offshore (point C). The absence of ice at location A for almost the whole LGP as simulated by TOTmax precludes the buttressing effect that limits ice discharge from the ice-sheet interior, as seen in the ORBmax simulation. By contrast, the presence of millennial-scale fluctuations in the ocean helps to maintain ice advection from the interior toward the ice–ocean margin, favouring ice discharge and limiting the ice increase by 300 m (point B) and 200 m inland (point C). Finally, oceanic temperature variability at millennial timescales strongly affects the region around Kangerlussuaq fjord (Fig. S4), where ice growth and loss clearly induce ice-thickness differences of more than 1000 m up to 300 km away from the ice-shelf edge.

Figure 7Root mean square deviation ε for ice thickness H (m) and velocity U (m a−1) calculated between TOTmax and ORBmax, following Eq. (8).


The response of the GrIS to millennial-scale variability in the ocean throughout the whole LGP is summarised in Fig. 7 at a large spatial scale. Ice thickness and velocity ε both calculated following Eq. (8) for the whole time period show that regions that exhibit a strong response to the ocean in terms of ice-thickness variations are also subjected to strong changes in velocity. This is highlighted in regions such as Baffin Bay, Jakobshavn Isbrae, Kangerlussuaq and Sermilik fjords, and NEGIS, confirming the results discussed in this section for specific times and locations.

4 Discussion

4.1 Comparison of model results to proxy data

Our simulations show that small ocean temperature variations can result in significant ice-volume fluctuations from the background glacial GrIS configuration. Also, the GrIS may have contributed to IRD discharge during the LGP, as suggested by proxy records (e.g. Andrews et al.2017; Barker et al.2015; Bond and Lotti1995; Hodell et al.2010; Jonkers et al.2010; Rasmussen et al.2016; Van Kreveld et al.2000; Voelker and Haflidason2015), and it is therefore interesting to outline its possible imprint associated with millennial-scale variability in the ocean.

Figure 8Simulated submarine melt (a) and simulated ice flux averaged over the southeast GrIS region (b) are compared to ice-discharge proxy data taken from the SO82-2 (c) and the JPC-13 (d) sediment cores in the North Atlantic.


We compare our model results to proxy data taken from sediment cores drilled in locations that may have been partly affected by recurrent ice discharge from the GrIS throughout the LGP. The analysed sites are located in the central North Atlantic, one along the eastern side of the Reykjanes Ridge (core JPC-13; Hodell et al.2010), close to the so-called IRD belt (Ruddiman1977), and the other on its western flank, further north in the ocean (core SO82-2; Rasmussen et al.2016). IRD from both sediment cores are thought to have a possible east GrIS source (Hodell et al.2010; Van Kreveld et al.2000). Therefore, we compare those IRD data to the modelled ice flux averaged over the area between Kangerlussuaq and Sermilik Fjords (Fig. 8). A very good agreement is found between the ice flux modelled by TOTmax and the proxy data, especially during MIS 3. It is suprising how well the ice-dicharge events corresponding to GS 16 and GS 13 (and GS 12 for SO82-2) are reproduced by the model when oceanic millennial-scale variability is considered. Some other peaks in ice discharge, such as those linked to GS 6 and GS 7, are also recognizable. These results suggest that high peaks in submarine melting, resulting from strong stadial–interstadial oceanic warming, may have triggered massive iceberg discharge from the east GrIS ice shelf to both locations. Those GrIS-sourced icebergs eventually reached the central North Atlantic, at least during certain stadial phases of D-O cycles. Although these results agree with the evidence that most of the IRDs found close to Reykjanes Ridge have eastern GrIS and Icelandic origins (Van Kreveld et al.2000), the hypothesis that part of the IRDs, especially those from sediment core JPC-13 (Hodell et al.2010), have a Laurentide Ice Sheet source cannot be discarded. Also, although evidence supports the presence of an ice shelf in Denmark Strait during certain periods of the LGP (Barker et al.2015; Van Kreveld et al.2000; Voelker et al.1998), this assumption is still debated (Andrews et al.2017) and a more profound analysis of sediment cores and models is required to confirm this hypothesis.

Another two locations situated in the Labrador Sea (MD95-2024; Weber et al.2001) and in the northwestern margin of Iceland (PS2644-5; Voelker and Haflidason2015) have also been investigated (Figs. S5 and S6). Gamma ray density (g m−3) reconstructed from sediment core MD95-2024 (Labrador Sea, Fig. S5) is compared to the simulated ice flux averaged over the Baffin Bay region. All main Heinrich events (HEs) and some of the GS of the LGP can be recognised in the gamma ray density signal (Weber et al.2001), which follows the IRD data extracted from the same sediment core. A very weak model–data agreement is found. The most important source of IRD found in the Labrador Sea is thought to be the Hudson Bay, which is likely responsible for the majority of iceberg discharge during HEs; thus this could explain the very low correspondence between model and data for this location. Finally, the total lithic fragment (IRD and thepra grains) extracted from sediment core PS2644-5 (NW Iceland, Fig. S6) is compared to the simulated ice flux averaged over the ice-covered Denmark Strait area close to Scoresby Sund. Some correspondence between modelled ice flux from the east GrIS and IRD deposition is found, especially during MIS3. Thus, we cannot rule out that millennial-scale variability in the ocean is partly responsible for enhanced iceberg production from the east–southeast GrIS bringing IRD to the PS2644-5 site by the East Greenland Current (Voelker and Haflidason2015).

Despite individual disagreements and with certain nuances, all of these comparisons show that a more responsive GrIS, as obtained when the model is forced with millennial-scale oceanic variability at the subsurface (TOTmax), is in better agreement with proxies compared to the orbital-only simulation (ORBmax). Rapid (millennial-scale) basal melting fluctuations allow for ice-flux increases that match the timing of some of the peaks of the proxy data, meaning that the GrIS could have contributed to the ice discharge at those times. With this analysis, we do not aim to precisely reconstruct the timing and spatial distribution of the iceberg discharge during D-O events. Instead, we aim to explore the implications of considering the oceanic millennial-scale variability on the GrIS ice fluctuation during the LGP. It is remarkable how much the sole presence of millennial-scale variability in the ocean can influence the GrIS evolution and ice discharge during the LGP interstadials, notwithstanding that millennial-scale variability is neglected in the atmosphere. However, tracing the origin of IRD deposition is a complex and open problem that needs a deeper understanding of changes in the ocean column during D-O events, as well as a better knowledge of the oceanic and glaciological processes that ultimately determine the deposition of IRD on the ocean floor due to iceberg discharge.

4.2 Model limitations and caveats in basal melting parameterisation

Some of the GRISLI-UCM model limitations already discussed by Tabone et al. (2018), such as the coarse model resolution (20 km by 20 km), which impedes the correct solution of small and steep fjords, the relatively simple GIA scheme, which takes into account only local changes in the ice load, and the PDD ablation scheme, which does not consider changes in past insolation, are also present in this work.

The simplicity of the submarine melt parameterisation adopted here allows for a direct investigation of the effect of past oceanic temperature fluctuations on the GrIS evolution. However it may disregard complex processes occurring at the ice–ocean interface. Moreover, here we have introduced a cut-off value to prevent freezing at the ice-shelf base and thus limit ice accretion during the glacial period. This assumption is reasonable for the low spatial resolution of the model; however, it ignores an existing process, since freezing is observed at the base of many marine glaciers, even for the present day (at least in Antarctica; Rignot et al.2013). A parameterisation accounting for variations in salinity (Holland et al.2008) or the geometry of the ice-shelf cavity (Gladstone et al.2012) or the slope of the base of the ice shelf (Little et al.2012) or more complex schemes describing meltwater buoyancy and water convection in the ice-shelf cavity (Jenkins2011) could help to describe in detail the processes below the ice shelves and along the grounding line. Recently, a fairly simple parameterisation applied for the AIS resolving ice-shelf cavities and using a boundary-layer theory for the melt has been capable of capturing the effect of ocean circulation below the ice shelves (Reese et al.2018). With this sensitivity study we do not aim to perfectly reproduce reality but rather to analyse the possible effects of an ocean driven by millennial variability on the ice-sheet evolution in Greenland. The parameterisation used here was explicitly developed for paleoclimatic simulations and fits this purpose well.

The basal melting evolution reproduced in the TOTmax simulation (solid line in Fig. 4b) is strongly related to a specific combination of the perturbed parameters in the basal melting equation (Eq. 6). This combination is reasonably chosen from the range of possible values as it produces the maximum oceanic millennial variability in the LE, but the same analysis could be performed for other basal melting configurations. A different magnitude of the melting peaks during MIS3 (Fig. S1) or a different distribution of the melting peaks (resulting from the use of a different temperature reconstruction to create the α and β indices) may vary the timing of grounding-line advance and retreat or ice growth and loss. This is also true for the orbitally driven basal melting forcing applied to the ORBmax simulation (dashed line in Fig. 4), which shows submarine melting rates higher than 0 m a−1 for only 10 kyr during MIS3 (between 55 and 45 ka). A different melt configuration could vary the timing and spatial distribution of grounding-line migration and ice cover.

Many authors associate D-O events to changes in sea-ice cover over the Nordic Seas (Dokken et al.2013; Hoff et al.2016; Jensen et al.2018; Li et al.2010; Sime et al.2019), which were disregarded in our set-up. The use of spatially variable fields could add a further degree of realism to the simulations. Perturbing the model with spatially variable subsurface oceanic temperatures, which may reflect variations in sea-ice cover, instead of a simple temperature time series, would likely affect the simulation of local ice mass variations throughout the last glacial cycle. Both stadial–interstadial and glacial–interglacial temperature anomalies could be taken from existing transient model outputs, for instance. However, a complete map of the observed PD basal melting rates for the whole Greenland domain does not exist yet, in contrast to Antarctica (Rignot et al.2013), thus limiting the effectiveness of including additional complexity at this time. The fact that the few available model reconstructions of stadial–interstadial oceanic temperatures suggest that subsurface warming increases almost homogeneously across the ocean during the stadials, at least in the Nordic Seas (Zhang et al.2014b; Alvarez-Solas et al.2018), supports the choice of ignoring spatial oceanic variations as a first approach. This is also in agreement with large-scale subsurface warming suggested by Rasmussen and Thomsen (2004) and Marcott et al. (2011). Thus, considering a spatial variation in oceanic temperatures (and sea ice) around the GrIS is not considered fundamental for this study, which primarily aims to look at the response of the whole GrIS to D-O events.

Finally, atmospheric precipitation and temperature used to perturb the ice-sheet model vary only at long (orbital) timescales, while shorter (millennial-scale) variations related to the D-O events are not considered here. This simplified experimental design is chosen following previous results demonstrating the important role of the ocean in the GrIS evolution over the last glacial periods at orbital timescales (Tabone et al.2018). The experiments here allow us to directly investigate how millennial-scale variability in the ocean can impact the evolution of the GrIS. However, we have omitted an important component of the climate system, and to comprehensively understand the effect of millennial-scale fluctuations on Greenland, further experiments should be done by including this variability in air temperature and precipitation.

5 Conclusions

We have assessed the effect of the millennial-scale oceanic variability on the evolution of the GrIS during the LGP. To do so, we used an ice-sheet–shelf model, in which the millennial-scale variability in the ocean was imposed as a fluctuation in the basal melting rate at the grounding line and below the ice shelves, resulting from oceanic temperature anomalies at the subsurface, where warmer (cooler) waters are associated with D-O stadials (interstadials). We first characterised the millennial variability through a sensitivity test for a broad range of values of the perturbed parameters in the submarine melting equation. We showed that the millennial-scale contribution to ice-volume variations during the LGP could have reached peaks of more than 1 m SLE. The southeastern area around Kangerlussuaq fjord, Baffin Bay and the NEGIS regions were found to be very sensitive to millennial-scale variability in the ocean. Ice thicknesses simulated at the marine margin differed by 500–1000 m from those simulated by orbital-only driven oceanic variations. Moreover, imprints of these differences are still found for several tens (hundreds – in certain regions) of kilometers far from the ice–ocean interface due to the velocity-driven upstream propagation of the ice-flow perturbation. Although the aim of this work was far from assessing the true timing and spatial distribution of any GrIS ice discharge that occurred during the D-O events, we showed that considering the millennial-scale variability in the ocean is necessary to reproduce some of the IRD peaks observed in proxy data. Specifically, the basal melting increase during D-O stadials, associated with warmer oceanic waters at the subsurface, could have been responsible for the enhanced release of Greenland icebergs as inferred from North Atlantic sediment cores. Our work thus suggests that millennial-scale-induced changes in ocean circulation and temperature may be important drivers of the GrIS evolution during the LGP, advancing the hypothesis of a potential role of the GrIS in oceanic reorganisations at millennial timescales.

Code availability

The GRISLI-UCM model code is available from the authors upon request.

Data availability

This work used published available data. Data from SO82-2 core were digitised using (Rohatgi, 2018) from Fig. 4 of Rasmussen et al. (2016).


The supplement related to this article is available online at:

Author contributions

IT performed the simulations, analysed the results and wrote the paper. All the authors of this work contributed to conceiving the experiment and writing the paper.

Competing interests

The authors declare that they have no conflict of interest.


The work was supported by the Spanish Ministry of Science and Innovation in the framework of the project MOCCA (Modelling Abrupt Climate Change, grant no. CGL2014-59384-R). Ilaria Tabone is funded by the Spanish National Programme for the Promotion of Talent and Its Employability through grant no. BES-2015-074097. Alexander Robinson is funded by the Ramón y Cajal Programme of the Spanish Ministry for Science, Innovation and Universities. The model simulations were carried out in the HPC of Climate Change of the International Campus of Excellence of Moncloa (EOLO), supported by MECD and MICINN. Finally, we are thankful to Catherine Ritz for providing the original model GRISLI.

Review statement

This paper was edited by Marit-Solveig Seidenkrantz and reviewed by Anders Svensson and one anonymous referee.


Alley, R., Anandakrishnan, S., and Jung, P.: Stochastic resonance in the North Atlantic, Paleoceanogr. Paleocl., 16, 190–198, 2001. a

Alvarez-Solas, J., Charbit, S., Ritz, C., Paillard, D., Ramstein, G., and Dumas, C.: Links between ocean temperature and iceberg discharge during Heinrich events, Nat. Geosci., 3, 122–126, 2010. a

Alvarez-Solas, J., Robinson, A., Montoya, M., and Ritz, C.: Iceberg discharges of the last glacial period driven by oceanic circulation changes, P. Natl. Acad. Sci. USA, 110, 41, 16350–16354, 2013. a

Alvarez-Solas, J., Banderas, R., Robinson, A., and Montoya, M.: Oceanic forcing of the Eurasian Ice Sheet on millennial time scales during the Last Glacial Period, Clim. Past Discuss.,, in review, 2018. a, b, c

Andrews, J., Dunhill, G., Vogt, C., and Voelker, A.: Denmark Strait during the Late Glacial Maximum and Marine Isotope Stage 3: Sediment sources and transport processes, Mar. Geol., 390, 181–198, 2017. a, b, c

Andrews, J. T., Barber, D., Jennings, A., Eberl, D., Maclean, B., Kirby, M., and Stoner, J.: Varying sediment sources (Hudson Strait, Cumberland Sound, Baffin Bay) to the NW Labrador Sea slope between and during Heinrich events 0 to 4, J. Quaternary Sci., 27, 475–484, 2012. a

Anklin, M., Barnola, J. M., Beer, J., Blunier, T., Chappellaz, J., Clausen, H. B., Dahl-Jensen, D., Dansgaard, W., De Angelis, M., Delmas, R. J., Duval, P., Fratta, M., Fuchs, A., Fuhrer, K., Gundestrup, N., Hammer, C., Iversen P., Johnsen, S., Jouzel, J., Kipfstuhl, J., Legrand, M., Lorius, C., Maggi, V., Miller, H., Moore, J. C., Oeschger, H., Orombelli, G., Peel, D. A., Raisbeck, G., Raynaud, D., Schott-Hvidberg, C., Schwander, J., Shoji, H., Souchez, R., Stauffer, B., Steffensen, J. P., Stievenard, M., Sveinbjornsdottir, A., Thorsteinsson, T., and Wolff, E. W.: Climate instability during the last interglacial period recorded in the GRIP ice core, Nature, 364, 203–207, 1993. a

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

Applegate, P. J., Kirchner, N., Stone, E. J., Keller, K., and Greve, R.: An assessment of key model parametric uncertainties in projections of Greenland Ice Sheet behavior, The Cryosphere, 6, 589–606,, 2012. a

Bagniewski, W., Meissner, K. J., and Menviel, L.: Exploring the oxygen isotope fingerprint of Dansgaard-Oeschger variability and Heinrich events, Quaternary Sci. Rev., 159, 1–14, 2017. a

Banderas, R., Alvarez-Solas, J., Robinson, A., and Montoya, M.: An interhemispheric mechanism for glacial abrupt climate change, Clim. Dynam., 44, 2897–2908, 2015. 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

Barbante, C., Barnola, J.-M., Becagli, S., Beer, J., Bigler, M., Boutron, C., Blunier, T., Castellano, E., Cattani, O., Chappellaz, J., Dahl-Jensen, D., Debret, M., Delmonte, B., Dick, D., Falourd, S., Faria, S., Federer, U., Fischer, H., Freitag, J., Frenzel, A., Fritzsche, D., Fundel, F., Gabrielli, P., Gaspari, V., Gersonde, R., Graf, W., Grigoriev, D., Hamann, I., Hansson, M., Hoffmann, G., Hutterli, M. A., Huybrechts, P., Isaksson, E., Johnsen, S., Jouzel, J., Kaczmarska, M., Karlin, T., Kaufmann, P., Kipfstuhl, S., Kohno, M., Lambert, F., Lambrecht, A., Landais, A., Lawer, G., Leuenberger, M., Littot, G., Loulergue, L., Lüthi, D., Maggi, V., Marino, F., Masson-Delmotte, V., Meyer, H., Miller, H., Mulvaney, R., Narcisi, B., Oerlemans, J., Oerter, H., Parrenin, F., Petit, J.-R., Raisbeck, G., Raynaud, D., Röthlisberger, R., Ruth, U., Rybak, O., Severi, M., Schmitt, J., Schwander, J., Siegenthaler, U., Siggaard-Andersen, M.-L., Spahni, R., Steffensen, J. P., Stenni, B., Stocker, T. F., Tison, J.-L., Traversi, R., Udisti, R., Valero-Delgado, F., van den Broeke, M. R., van de Wal, R. S. W., Wagenbach, D., Wegner, A., Weiler, K., Wilhelms, F., Winther, J.-G., and Wolff, E.: One-to-one coupling of glacial climate variability in Greenland and Antarctica, Nature, 444, 195–198, 2006. a

Barker, S., Chen, J., Gong, X., Jonkers, L., Knorr, G., and Thornalley, D.: Icebergs not the trigger for North Atlantic cold events, Nature, 520, 333–336, 2015. a, b, c, d

Bassis, J. N., Petersen, S. V., and Mac Cathles, L.: Heinrich events triggered by ocean forcing and modulated by isostatic adjustment, Nature, 542, 332–334, 2017. a

Bauer, E. and Ganopolski, A.: Comparison of surface mass balance of ice sheets simulated by positive-degree-day method and energy balance approach, Clim. Past, 13, 819–832,, 2017. a

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

Bintanja, R. and Van de Wal, R.: North American ice-sheet dynamics and the onset of 100,000-year glacial cycles, Nature, 454, 869–872, 2008. a

Blasco, J., Tabone, I., Alvarez-Solas, J., Robinson, A., and Montoya, M.: The Antarctic Ice Sheet response to glacial millennial-scale variability, Clim. Past, 15, 121–133,, 2019. a, b, c

Boers, N., Ghil, M., and Rousseau, D.-D.: Ocean circulation, ice shelf, and sea ice interactions explain Dansgaard-Oeschger cycles, P. Natl. Acad. Sci. USA, 115, E11005–E11014,, 2018. a, b

Böhm, E., Lippold, J., Gutjahr, M., Frank, M., Blaser, P., Antz, B., Fohlmeister, J., Frank, N., Andersen, M., and Deininger, M.: Strong and deep Atlantic meridional overturning circulation during the last glacial cycle, Nature, 517, 73–76, 2015. a

Bond, G., Broecker, W., Johnsen, S., McManus, J., Labeyrie, L., Jouzel, J., and Bonani, G.: Correlations between climate records from North Atlantic sediments and Greenland ice, Nature, 365, 143–147, 1993. a

Bond, G. C. and Lotti, R.: Iceberg discharges into the North Atlantic on millennial time scales during the last glaciation, Science, 267, 1005–1010, 1995. a, b

Brady, E. C. and Otto-Bliesner, B. L.: The role of meltwater-induced subsurface ocean warming in regulating the Atlantic meridional overturning in glacial climate simulations, Clim. Dynam., 37, 1517–1532, 2011. a, b

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

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

Dansgaard, W., Johnsen, S., Clausen, H., Dahl-Jensen, D., Gundestrup, N., Hammer, C., Hvidberg, C., Steffensen, J., Sveinbjörnsdottir, A., Jouzel, J., and Bond, G.: Evidence for general instability of past climate from a 250-kyr ice-core record, Nature, 364, 218–220, 1993. a

Darby, D. A., Bischof, J. F., Spielhagen, R. F., Marshall, S. A., and Herman, S. W.: Arctic ice export events and their potential impact on global climate during the late Pleistocene, Paleoceanography, 17, 15-1–15-17, 2002. a

Dokken, T. M. and Jansen, E.: Rapid changes in the mechanism of ocean convection during the last glacial period, Nature, 401, 458–461, 1999. a

Dokken, T. M., Nisancioglu, K. H., Li, C., Battisti, D. S., and Kissel, C.: Dansgaard-Oeschger cycles: Interactions between ocean and sea ice intrinsic to the Nordic seas, Paleoceanogr. Paleocl., 28, 491–502, 2013. a, b, c

Ezat, M. M., Rasmussen, T. L., and Groeneveld, J.: Persistent intermediate water warming during cold stadials in the southeastern Nordic seas during the past 65 ky, Geology, 42, 663–666, 2014. a, b, c

Fettweis, X., Franco, B., Tedesco, M., Van Angelen, J., Lenaerts, J., Van den Broeke, M., and Gallée, H.: Estimating 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

Fleitmann, D., Cheng, H., Badertscher, S., Edwards, R., Mudelsee, M., Göktürk, O. M., Fankhauser, A., Pickering, R., Raible, C., Matter, A., Kramers, J., and Tüysüz, O.: Timing and climatic impact of Greenland interstadials recorded in stalagmites from northern Turkey, Geophys. Res. Lett., 36, L19707,, 2009. a

Ganopolski, A. and Rahmstorf, S.: Rapid changes of glacial climate simulated in a coupled climate model, Nature, 409, 153–158, 2001. a, b

Gladstone, R. M., Payne, A. J., and Cornford, S. L.: Parameterising the grounding line in flow-line ice sheet models, The Cryosphere, 4, 605–619,, 2010. a

Gladstone, R. M., Lee, V., Rougier, J., Payne, A. J., Hellmer, H., Le Brocq, A., Shepherd, A., Edwards, T. L., Gregory, J., and Cornford, S. L.: Calibrated prediction of Pine Island Glacier retreat during the 21st and 22nd centuries with a coupled flowline model, Earth Planet. Sc. Lett., 333, 191–199, 2012. a

Gottschalk, J., Skinner, L. C., Misra, S., Waelbroeck, C., Menviel, L., and Timmermann, A.: Abrupt changes in the southern extent of North Atlantic Deep Water during Dansgaard–Oeschger events, Nat. Geosci., 8, 950–954, 2015. a

Grootes, P. M., Stuiver, M., White, J., Johnsen, S., and Jouzel, J.: Comparison of oxygen isotope records from the GISP2 and GRIP Greenland ice cores, Nature, 366, 552–554, 1993. a

Henry, L., McManus, J. F., Curry, W. B., Roberts, N. L., Piotrowski, A. M., and Keigwin, L. D.: North Atlantic ocean circulation and abrupt climate change during the last glaciation, Science, 353, 470–474, 2016. a

Hodell, D. A., Evans, H. F., Channell, J. E., and Curtis, J. H.: Phase relationships of North Atlantic ice-rafted debris and surface-deep climate proxies during the last glacial period, Quaternary Sci. Rev., 29, 3875–3886, 2010. a, b, c, d, e

Hoff, U., Rasmussen, T. L., Stein, R., Ezat, M. M., and Fahl, K.: Sea ice and millennial-scale climate variability in the Nordic seas 90 kyr ago to present, Nat. Commun., 7, 12247,, 2016. a, b

Holland, P. R., Jenkins, A., and Holland, D. M.: The response of ice shelf basal melting to variations in ocean temperature, J. Clim., 21, 2558–2572, 2008. a

Hutter, K.: Theoretical Glaciology: Material Science of Ice and the Mechanics of Glaciers and Ice Sheets (Mathematical Approaches to Geophysics), edited by: Reidel, D., Dordrecht, the Netherlands, 1983. 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

Jenkins, A.: Convection-driven melting near the grounding line of ice shelves and tidewater glaciers, J. Phys. Oceanogr., 41, 2279–2294,, 2011. a

Jennings, A. E., Andrews, J. T., Cofaigh, C. Ó., Onge, G. S., Sheldon, C., Belt, S. T., Cabedo-Sanz, P., and Hillaire-Marcel, C.: Ocean forcing of Ice Sheet retreat in central west Greenland from LGM to the early Holocene, Earth Planet. Sc. Lett., 472, 1–13, 2017. a

Jensen, M. F., Nummelin, A., Nielsen, S. B., Sadatzki, H., Sessford, E., Risebrobakken, B., Andersson, C., Voelker, A., Roberts, W. H. G., Pedro, J., and Born, A.: A spatiotemporal reconstruction of sea-surface temperatures in the North Atlantic during Dansgaard–Oeschger events 5–8, Clim. Past, 14, 901–922,, 2018. a, b, c

Johnsen, S. J., Clausen, H. B., Dansgaard, W., Fuhrer, K., Gundestrup, N.and Hammer, C. U. I. P., Jouzel, J., Stauffer, B., and Steffensen, J. P.: Irregular glacial interstadials recorded in a new Greenland ice core, Nature, 359, 311–313, 1992. a

Johnsen, S. J., Dahl-Jensen, D., Gundestrup, N., Steffensen, J. P., Clausen, H. B., Miller, H., Masson-Delmotte, V., Sveinbjörnsdottir, A. E., and White, J.: Oxygen isotope and palaeotemperature records from six Greenland ice-core stations: Camp Century, Dye-3, GRIP, GISP2, Renland and NorthGRIP, J. Quaternary Sci., 16, 299–307, 2001. a

Jonkers, L., Moros, M., Prins, M. A., Dokken, T., Dahl, C. A., Dijkstra, N., Perner, K., and Brummer, G.-J. A.: A reconstruction of sea surface warming in the northern North Atlantic during MIS 3 ice-rafting events, Quaternary Sci. Rev., 29, 1791–1800, 2010. a, b, c

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, b

Knutti, R., Flückiger, J., Stocker, T., and Timmermann, A.: Strong hemispheric coupling of glacial climate through freshwater discharge and ocean circulation, Nature, 430, 851–856, 2004. 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,, 2018. a

Le Meur, E. and Huybrechts, P.: A comparison of different ways of dealing with isostasy: examples from modelling the Antarctic ice sheet during the last glacial cycle, Ann. Glaciol., 23, 309–317, 1996. a

Li, C., Battisti, D. S., Schrag, D. P., and Tziperman, E.: Abrupt climate shifts in Greenland due to displacements of the sea ice edge, Geophys. Res. Lett., 32, L19702,, 2005. a

Li, C., Battisti, D. S., and Bitz, C. M.: Can North Atlantic sea ice anomalies account for Dansgaard–Oeschger climate signals?, J. Clim., 23, 5457–5475, 2010. a, b

Little, C. M., Goldberg, D., Gnanadesikan, A., and Oppenheimer, M.: On the coupled response to ice-shelf basal melting, J. Glaciol., 58, 203–215, 2012. a

Liu, Y., Moore, J. C., Cheng, X., Gladstone, R. M., Bassis, J. N., Liu, H., Wen, J., and Hui, F.: Ocean-driven thinning enhances iceberg calving and retreat of Antarctic ice shelves, P. Natl. Acad. Sci. USA, 112, 3263–3268, 2015. a

Lynch-Stieglitz, J.: The Atlantic meridional overturning circulation and abrupt climate change, Ann. Rev. Mar. Sci., 9, 83–104, 2017. a

Marcott, S. A., Clark, P. U., Padman, L., Klinkhammer, G. P., Springer, S. R., Liu, Z., Otto-Bliesner, B. L., Carlson, A. E., Ungerer, A., Padman, J., Feng , H., Cheng, J., and Schmittner A.: Ice-shelf collapse from subsurface warming as a trigger for Heinrich events, P. Natl. Acad. Sci. USA, 108, 13415–13419, 2011. a, b

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

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

McKay, M. D., Beckman, R. J., and Conover, W. J.: Comparison of three methods for selecting values of input variables in the analysis of output from a computer code, Technometrics, 21, 239–245, 1979. a

Menviel, L., Timmermann, A., Friedrich, T., and England, M. H.: Hindcasting the continuum of Dansgaard–Oeschger variability: mechanisms, patterns and timing, Clim. Past, 10, 63–77,, 2014. a

Mignot, J., Ganopolski, A., and Levermann, A.: Atlantic subsurface temperatures: Response to a shutdown of the overturning circulation and consequences for its recovery, J. Clim., 20, 4884–4898, 2007. a

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

Morland, L., Van der Veen, C., and Oerlemans, J.: Dynamics of the West Antarctic ice sheet., Proc. Workshop on the Dynamics of the West Antarctic Ice Sheet, edited by: Reidel, D., Dordrecht, the Netherlands, 99–116, 1987. 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

Peltier, W. R. and Vettoretti, G.: Dansgaard-Oeschger oscillations predicted in a comprehensive model of glacial climate: A “kicked” salt oscillator in the Atlantic, Geophys. Res. Lett., 41, 7306–7313, 2014. a, b

Petersen, S. V., Schrag, D. P., and Clark, P. U.: A new mechanism for Dansgaard-Oeschger cycles, Paleoceanography, 28, 24–30, 2013. 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

Rasmussen, S. O., Bigler, M., Blockley, S. P., Blunier, T., Buchardt, S. L., Clausen, H. B., Cvijanovic, I., Dahl-Jensen, D., Johnsen, S. J., Fischer, H., Gkinis, V., Guillevic , M., Hoek, W., Lowe, J., Pedro, J., Popp, T., Seierstad, I., Steffensen, J. P., Svensson A. M., Vallelonga, P., Vinther, B., Walker, M., Wheatley, J. J., and Winstrup, M.: A stratigraphic framework for abrupt climatic changes during the Last Glacial period based on three synchronized Greenland ice-core records: refining and extending the INTIMATE event stratigraphy, Quaternary Sci. Rev., 106, 14–28, 2014. a

Rasmussen, T. L. and Thomsen, E.: The role of the North Atlantic Drift in the millennial timescale glacial climate fluctuations, Palaeogeogr. Palaeocl., 210, 101–116, 2004. a, b, c

Rasmussen, T. L., Thomsen, E., Labeyrie, L., and van Weering, T. C.: Circulation changes in the Faeroe-Shetland Channel correlating with cold events during the last glacial period (58–10 ka), Geology, 24, 937–940, 1996. a

Rasmussen, T. L., Thomsen, E., and Moros, M.: North Atlantic warming during Dansgaard-Oeschger events synchronous with Antarctic warming and out-of-phase with Greenland climate, Sci. Rep., 6, 20535,, 2016. a, b, c, d, e, f, g, h

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

Reese, R., Albrecht, T., Mengel, M., Asay-Davis, X., and Winkelmann, R.: Antarctic sub-shelf melt rates via PICO, The Cryosphere, 12, 1969–1985,, 2018. 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 Steffen, K.: Channelized bottom melting and stability of floating ice shelves, Geophys. Res. Lett., 35, L02503,, 2008. a

Rignot, E., Koppes, M., and Velicogna, I.: Rapid submarine melting of the calving faces of West Greenland glaciers, Nat. Geosci., 3, 187–191,, 2010. a

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

Rignot, E., Xu, Y., Menemenlis, D., Mouginot, J., Scheuchl, B., Li, X., Morlighem, M., Seroussi, H., Broeke, M. v., Fenty, I., Cai, C., An, L., and de Fleurian, B.: Modeling of ocean-induced ice melt rates of five west Greenland glaciers over the past two decades, Geophys. Res. Lett., 43, 6374–6382,, 2016. 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

Robinson, A. and Goelzer, H.: The importance of insolation changes for paleo ice sheet modeling, Geophys. Res. Lett., 8, 1419–1428, 2014. a

Robinson, A., Alvarez-Solas, J., Calov, R., Ganopolski, A., and Montoya, M.: MIS-11 duration key to disappearance of the Greenland ice sheet, Nat. Commun., 8, 16008,, 2017. a

Rohatgi, A.: WebPlotDigitizer Version: 4.1, available at:, last access: January, 2018. 

Ruddiman, W. F.: Late Quaternary deposition of ice-rafted sand in the subpolar North Atlantic (lat 40 to 65 N), Geol. Soc. Am. Bull., 88, 1813–1827, 1977. 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

Sessford, E., Tisserand, A., Risebrobakken, B., Andersson, C., Dokken, T., and Jansen, E.: High-Resolution Benthic Mg/Ca Temperature Record of the Intermediate Water in the Denmark Strait Across D-O Stadial-Interstadial Cycles, Paleoceanogr. Paleocl., 33, 1169–1185, 2018. a, b, c

Sévellec, F. and Fedorov, A. V.: Unstable AMOC during glacial intervals and millennial variability: The role of mean sea ice extent, Earth Planet. Sc. Lett., 429, 60–68, 2015. a

Shackleton, N. J., Hall, M. A., and Vincent, E.: Phase relationships between millennial-scale events 64,000–24,000 years ago, Paleoceanography, 15, 565–569, 2000. a

Shaffer, G., Olsen, S. M., and Bjerrum, C. J.: Ocean subsurface warming as a mechanism for coupling Dansgaard-Oeschger climate cycles and ice-rafting events, Geophys. Res. Lett., 31, L24202,, 2004. a

Shapiro, N. M. and Ritzwoller, M. H.: Inferring surface heat flux distributions guided by a global seismic model: particular application to Antarctica, Earth Planet. Sc. Lett., 223, 213–224,, 2004. a

Sime, L., Hopcroft, P., and Rhodes, R.: Impact of abrupt sea ice loss on Greenland water isotopes during the last glacial period, P. Natl. Acad. Sci. USA, 116, 4099–4104,, 2019. a, b

Skinner, L. and Elderfield, H.: Rapid fluctuations in the deep North Atlantic heat budget during the last glacial period, Paleoceanography, 22, PA1205,, 2007. a

Stein, R., Nam, S.-i., Grobe, H., and Hubberten, H.: Late Quaternary glacial history and short-term ice-rafted debris fluctuations along the East Greenland continental margin, Geol. Soc. Lond. Spec. Publ., 111, 135–151, 1996. a

Stone, E., Lunt, D., Annan, J., and Hargreaves, J.: Quantification of the Greenland ice sheet contribution to Last Interglacial sea level rise, Clim. Past, 9, 621–639,, 2013. a

Stone, E. J., Lunt, D. J., Rutt, I. C., and Hanna, E.: Investigating the sensitivity of numerical model simulations of the modern state of the Greenland ice-sheet and its future response to climate change, The Cryosphere, 4, 397–417,, 2010. a

Straneo, F., Sutherland, D. A., Holland, D., Gladish, C., Hamilton, G. S., Johnson, H. L., Rignot, E., Xu, Y., and Koppes, M.: Characteristics of ocean waters reaching Greenland's glaciers, Ann. Glaciol., 53, 202–210,, 2012. 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, j, k, l

Van Kreveld, S., Sarnthein, M., Erlenkeuser, H., Grootes, P., Jung, S., Nadeau, M., Pflaumann, U., and Voelker, A.: Potential links between surging ice sheets, circulation changes, and the Dansgaard-Oeschger cycles in the Irminger Sea, 60–18 kyr, Paleoceanography, 15, 425–442, 2000. a, b, c, d, e

Vasskog, K., Langebroek, P. M., Andrews, J. T., Nilsen, J. E., and Nesje, A.: The Greenland Ice Sheet during the last glacial cycle: Current ice loss and contribution to sea-level rise from a palaeoclimatic perspective, Geophys. Res. Lett., 43, 5336–5344, 2015. a

Vellinga, M. and Wood, R. A.: Global climatic impacts of a collapse of the Atlantic thermohaline circulation, Climatic Change, 54, 251–267, 2002. a

Verplanck, E. P., Farmer, G. L., Andrews, J., Dunhill, G., and Millo, C.: Provenance of Quaternary glacial and glacimarine sediments along the southeast Greenland margin, Earth Planet. Sc. Lett., 286, 52–62, 2009. a

Vettoretti, G. and Peltier, W. R.: Interhemispheric air temperature phase relationships in the nonlinear Dansgaard-Oeschger oscillation, Geophys. Res. Lett., 42, 1180–1189, 2015. a, b

Vettoretti, G. and Peltier, W. R.: Thermohaline instability and the formation of glacial North Atlantic super polynyas at the onset of Dansgaard-Oeschger warming events, Geophys. Res. Lett., 43, 5336–5344, 2016. a, b

Vettoretti, G. and Peltier, W. R.: Fast Physics and Slow Physics in the Nonlinear Dansgaard–Oeschger Relaxation Oscillation, J. Clim., 31, 3423–3449, 2018. a, b

Voelker, A. H.: Instability of the Atlantic overturning circulation during Marine Isotope Stage 3, Quaternary Sci. Rev., 21, 1185–1212, 2002. a

Voelker, A. H. and Haflidason, H.: Refining the Icelandic tephrachronology of the last glacial period–the deep-sea core PS2644 record from the southern Greenland Sea, Glob. Planet. Change, 131, 35–62, 2015. a, b, c, d

Voelker, A. H., Sarnthein, M., Grootes, P. M., Erlenkeuser, H., Laj, C., Mazaud, A., Nadeau, M.-J., and Schleicher, M.: Correlation of marine 14C ages from the Nordic Seas with the GISP2 isotope record: Implications for 14C calibration beyond 25 ka BP, Radiocarbon, 40, 517–534, 1998. a, b

Wang, Y., Wu, J., Wu, J., Mu, X., Xu, H., and Chen, J.: Correlation between high-resolution climate records from a Nanjing stalagmite and GRIP ice core during the last glaciation, Sci. China Ser. D, 44, 14–23, 2001.  a

Wang, Z., Zhang, X., Guan, Z., Sun, B., Yang, X., and Liu, C.: An atmospheric origin of the multi-decadal bipolar seesaw, Sci. Rep., 5, 8909,, 2015. a

Weber, M., Mayer, L. A., Hillaire-Marcel, C., Bilodeau, G., Rack, F., Hiscott, R., and Aksu, A.: Derivation of δ18O from sediment core log data: Implications for millennial-scale climate change in the Labrador Sea, Paleoceanogr. Paleocl., 16, 503–514, 2001. 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

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

Wunsch, C.: Abrupt climate change: An alternative view, Quaternary Res., 65, 191–203, 2006. a

Zhang, X., Lohmann, G., and Knorr, G.and Purcell, C.: Abrupt glacial climate shifts controlled by ice sheet changes, Nature, 3512, 290–294, 2014a. a

Zhang, X., Prange, M., Merkel, U., and Schulz, M.: Instability of the Atlantic overturning circulation during Marine Isotope Stage 3, Geophys. Res. Lett., 41, 4285–4293,, 2014a. a, b

Zhang, X., Knorr, G., Lohmann, G., and Barker, S.: Abrupt North Atlantic circulation changes in response to gradual CO2 forcing in a glacial climate state, Nat. Geosci., 10, 518–523, 2017. a

Short summary
By using a 3-D hybrid ice-sheet–shelf model, we investigate the impact of millennial-scale oceanic variability on the Greenland Ice Sheet (GrIS) evolution during the last glacial period (LGP). We show that the GrIS may have strongly reacted to oceanic temperature fluctuations associated with Dansgaard–Oeschger cycles, contributing to sea-level variations of more than 1 m. Our results open the chance for a non-negligible role of the GrIS in millennial-scale oceanic reorganisations during the LGP.