Oceanic forcing of the Eurasian Ice Sheet on millennial time scales during the Last Glacial Period

Abstract. The last glacial period (LGP; ca.110–10 ka BP) was marked by the existence of two types of abrupt climatic changes, Dansgaard-Oeschger (DO) and Heinrich (H) events. Although the mechanisms behind these are not fully understood, it is generally accepted that the presence of ice sheets played an important role in their occurrence. While an important effort has been made to investigate the dynamics and evolution of the Laurentide Ice Sheet (LIS) during this period, the Eurasian Ice Sheet (EIS) has not received much attention, in particular from a modeling perspective. However, meltwater discharge from this and other ice sheets surrounding the Nordic Seas is often implied as a potential cause of ocean instabilities that lead to glacial abrupt climate changes. Thus, a better understanding of its variations during the LGP is important to understand its role in glacial abrupt climate changes. Here we investigate the response of the EIS to millennial-scale climate variability during the LGP. We use a hybrid, three-dimensional, thermomechanical ice-sheet model that includes ice shelves and ice streams. The model is forced offline through a novel perturbative approach that includes the effect of both atmospheric and oceanic variations and provides a more realistic treatment of millennial-scale climatic variability than conventional methods. Our results show that the EIS responds with enhanced ice discharge in phase with interstadial warming in the North Atlantic when forced with surface ocean temperatures. Conversely, when subsurface ocean temperatures are used, enhanced ice discharge occurs both during stadials and at the beginning of the interstadials. Separating the atmospheric and oceanic effects demonstrates the major role of the ocean in controlling the dynamics of the EIS on millennial time scales. While the atmospheric forcing alone is only able to produce modest iceberg discharges, warming of the ocean leads to higher rates of iceberg discharges as a result of relatively strong basal melting at the margins of the ice sheet. Together with previous work, our results provide a consistent explanation for the response of the LIS and the EIS to glacial abrupt climate changes, and highlight the need for stronger constraints on the local North Atlantic behavior in order to improve our understanding of the ice sheet's glacial dynamics.


to relatively warm (interstadial) conditions within decades (Dansgaard et al., 1993) followed by a gradual cooling interval lasting from centuries to millennia and an ultimate phase of rapid cooling back to stadial conditions (Steffensen et al., 2008).
Superimposed on the millennial-scale variability associated with DO-events, an additional lower-frequency climatic cycle is identified. So-called Bond cycles are flanked by prolonged stadials ending with prominent DO-events within about 7-10 kyr (Bond et al., 1993). Preceding these, and concomitant with the culmination of the prolonged stadials, H-events are registered 5 in North Atlantic marine sediments as layers of remarkably high concentrations of ice-rafted debris (IRD) (Heinrich, 1988) as a result of massive iceberg discharges from the Laurentide ice-sheet (LIS) (Hemming, 2004).
While significant effort has been invested in understanding the role of the LIS in glacial abrupt climate changes, the dynamics of the Eurasian Ice Sheet (EIS) during the LGP has received comparatively less attention from a modeling perspective.
However, improving our understanding of its evolution and response to past climate changes is important for a number of 10 reasons. First, constraining freshwater inputs into the North Atlantic Ocean is crucial for a better understanding of the driving mechanisms of glacial abrupt climate changes (Rasmussen and Thomsen, 2013), since meltwater discharge from the ice sheets surrounding the Nordic Seas is often implied as a cause of ocean instabilities. Precursor events could possibly have originated from the European and Icelandic ice sheets (Grousset et al., 2000;Scourse et al., 2000). Meltwater peaks in the Norwegian Sea as well as in the southern border of the EIS during Marine Isotopic Stage 3 (MIS 3) have been associated with H events and 15 millennial-scale climate variability (Lekens et al., 2006;Toucanne et al., 2015). From a broader perspective, the EIS, consisting of the Fennoscandian, the British Isles and the Barents-Kara ice sheets (FIS, BIIS and BKSIS, respectively) contained a large marine-based sector at its maximum extension (Hughes et al., 2016) that was exposed to oceanic variations, and the BKSIS is often considered as an analog for the current West Antarctic ice sheet (WAIS). At the LGM both had a similar size, but while the WAIS endured the deglaciation, the BKSIS completely disappeared (Andreassen and Winsborrow, 2009). Understanding 20 the underlying mechanisms would provide important insights into the future evolution of the WAIS (Gudlaugsson et al., 2013(Gudlaugsson et al., , 2017. Reconstructing the EIS response to past glacial abrupt climate changes prior to the LGM has been difficult, in part because, in reaching its maximum extent, the ice sheet eroded and removed nearly all older deposits. Nevertheless, the available paleodata indicate that during MIS 3 the EIS was highly dynamic, with its advance and retreat closely linked to stadials and interstadials 25 (Toucanne et al., 2015). In this line, records from Norway (Mangerud et al., 2003(Mangerud et al., , 2010Olsen et al., 2002), Finland (Helmens and Engels, 2010) and Sweden (Wohlfarth, 2010) indicate rapid and rythmic ice-sheet variations in western Scandinavia, with advances and retreats during stadials and interstadials, respectively. Recent records also indicate enhanced meltwater discharges during interstadials from the Svalbard-Barents Sea ice sheet and probably also from the Scandinavian ice sheet (Rasmussen and Thomsen, 2013). The resolution and quality of geophysical data across marine sectors has improved considerably in the 30 past decade (Hughes et al. (2016) and references therein). The results confirm substantial variations of the EIS volume, with the largest uncertainties in marine sectors of the ice sheets. Strong variations in the deposition of IRD suggest high co-variability of the BIIS with changes in ocean sea surface temperature (Hall et al., 2011;Scourse et al., 2009) and variations in EIS ice streams (Becker et al., 2017). North Atlantic marine sediment records register widespread variations of IRD input throughout the LGP indicating variations of iceberg rafting from virtually all surrounding ice sheets. Sources and timing differ among different 35 2 Clim. Past Discuss., https://doi.org/10.5194/cp-2018-89 Manuscript under review for journal Clim. Past Discussion started: 9 October 2018 c Author(s) 2018. CC BY 4.0 License. sites. A dominant periodicity equal to that of DO-events was identified in the Irminger Sea, with the largest IRD peaks at the end of stadials originating in the Iceland and Greenland ice sheets (von Kreveld et al., 2000). Strong millennial-scale iceberg rafting variability of the BIIS has been documented as well in the North Sea (Hall et al., 2011;Peck et al., 2007;Scourse et al., 2009), but enhanced IRD seems to occur both during interstadials and stadials. For the FIS, IRD records in the Norwegian Sea show the characteristic DO periodicity, with IRD discharge occurring just before interstadial transitions (Lekens et al., 2006). 5 More recently, however, an increase in IRDs from Fennoscandia during interstadials has been reported (Dokken et al., 2013;Becker et al., 2017). Correlating IRD occurrence with temperature changes registered in Greenland remains difficult, however, because it requires an extremely well dated chronology to assess the phasing between ocean sediments and ice cores.
Progress has been achieved also in the past decade using ice-sheet models. Siegert and Dowdeswell (2004) used inverse modelling to simulate the EIS evolution during the second part of the LGP, matching the geological evidence presented by 10 optimizing the fit with data. Forsström and Greve (2004) used subsequent versions of a three-dimensional, polythermal icesheet model to simulate the EIS evolution throughout the LGP. Important variations in the EIS ice volume in response to temperature and precipitation variations were simulated. Clason et al. (2014) additionally included a parameterisation of surface meltwater enhanced sliding. In both cases too much ice was simulated in the northeastern EIS. Gudlaugsson et al. (2017) used the same model but introducing a simple representation of the subglacial hydrological system, focusing on its role in the 15 temporal evolution of the EIS. Recently, an ice-sheet model constrained by data has been used to simulate the EIS evolution throughout part of LGP (Patton et al., 2016). The model targets the most probable EIS distribution at different time slices and reproduces substantial ice-volume variations. However, all of these models suffer from limitations, such as the use of the shallow-ice approximation (SIA) and its associated lack of an explicit treatment of the oceanic forcing. Marshall and Koutnik (2006) investigated the production of icebergs from all the North American ice sheets with a parameterized calving model.

20
They found different behaviors on millennial time-scales depending on the local glaciological and climatic characteristic, with increased iceberg production both during stadials (e.g. from Iceland) or during interstadials (e.g. from Barents Sea). Nonetheless, sub-marine melting at the grounding line has not been explicitely considered until now and its impacts on millennial-scale variability have not been investigated up to now from a modelling perspective.
Here, we investigate the response of the EIS to millennial-scale climate variability during MIS 3 using a three-dimensional 25 ice-sheet model. To this end, a novel offline approach is used that provides a better representation of millennial-scale climate variability (Banderas et al., 2018). In addition, for the first time, both the atmospheric and oceanic effects of millennial scale climate variability associated with glacial abrupt climate changes are considered. This facilitates the quantification of the relative contribution of surface (ablation) and dynamic processes related to ice-ocean interactions.
The paper is organized as follows: in Section 2 the ice-sheet model, the forcing method and the experimental setup are 30 described. In Section 3 the response of the EIS to the imposed forcing is shown, the focus being the evolution of its ice volume, its impact on sea level and the mechanisms behind meltwater and ice discharge. Finally, the main conclusions are summarised in Section 5. The model used in this study is the ice-sheet model GRISLI-UCM, an extension of the original model GRISLI developed by Ritz et al. (2001). GRISLI-UCM is a hybrid three-dimensional thermomechanical ice-sheet model. Inland ice flows through deformation under the Shallow Ice Approximation (SIA, Hutter, 1983). Ice shelves and ice streams are described following 5 the Shallow Shelf Approximation (SSA, MacAyeal, 1989). Ice streams (areas of fast flow, typically faster than 10 2 m a −1 ) are considered as dragging ice shelves, allowing for basal movement of the ice (Bueler and Brown, 2009). Basal stress under ice streams is proportional to ice velocity and to the effective pressure of ice. The effects of varying this proportionality factor on the simulated ice streams are discussed in Alvarez-Solas et al. (2011). The locations of the ice streams are determined by the presence of basal water within areas where the sediment layer is saturated.

10
The criterion to activate SSA inland relies on the presence of water above 1 meter in places of soft sediments (Laske, 1997) and above 400 meters in absence of such sediments. The grounding line position dynamically evolves following the flotation criterion after the mass conservation equation is solved. Calving takes place at the-shelf front following a double criterium. First, its thickness must first fall below a threshold. (H calv = 150 m, in the standard setup). This is a semiempirical parameter reflecting the fact that this is the typical thickness of observed ice-shelf fronts. Second, the upstream advection 15 must fail to maintain the ice thickness above this threshold following a semi-Lagrangian approach (Peyaud et al., 2007) to account for the fact that ice-flux divergence fosters the formation of crevasses (Levermann et al., 2012). GRISLI-UCM thus explicitly calculates grounding line migration, ice-stream and ice-shelf velocities. This allows the model to properly represent both grounded and floating ice. GRISLI-UCM uses finite differences on a staggered Cartesian grid at a 40 km resolution, corresponding to 224×208 grid points for the Northern Hemisphere domain, including the EIS, with 21 vertical levels. By 20 default, initial topographic conditions are provided by surface and bedrock elevations built from the ETOPO1 dataset (Amante and Eakins, 2009) and ice thickness (Bamber et al., 2001). The surface mass balance is given by the sum of accumulation and ablation, both of which are calculated from monthly surface air temperatures (SATs) and monthly total precipitation. Accumulation is calculated by assuming that the fraction of solid precipitation is proportional to the fraction of the year with mean daily temperature below 2 • C. The daily temperature is computed from monthly SATs assuming that the annual 25 temperature cycle follows a cosine function. Ablation is calculated using the positive-degree-day (PDD) method (Reeh, 1989).
Its main parameters are the standard deviation of daily temperature, σ, and the conversion factors from PDDs to melt for snow and ice, f PDDsnow and f PDDice . Here, σ = 5 K, f PDDsnow = 0.003 mwe PDD −1 and f PDDice = 0.008 mwe PDD −1 . Refreezing is considered, with a value of C si = 60% (see Section 2 in the Supplementary Information). GRISLI-UCM accounts for changes in elevation at each time step considering a linear atmospheric vertical profile for temperature with different lapse rates in 30 summer and in the annual mean (0.0065 and 0.0080 K m −1 , respectively) to account for the smaller summer atmospheric vertical stability.
Basal melting inland depends on pressure and water content at the base of the ice sheet (Ritz et al., 2001) as well as on the geothermal heat flux, which is prescribed from the reconstruction by Shapiro and Ritzwoller (2004). Basal melting for floating ice is computed using a linear temperature anomaly with respect to the freezing point. The details of the implementation of the boundary conditions (SMB and oceanic basal melting) in this particular study are given below (Section 2.2).

Offline forcing method
SMB and oceanic basal melting are obtained through a time-varying synthetic climatology built through a novel method that is found to provide a more realistic offline forcing for ice-sheet models than classical offline methods (Banderas et al., 2018). 5 The method follows a perturbative approach in the sense that the forcing combines the present-day climatology, obtained from observational data, together with simulated anomalies. But in contrast to usual offline forcing methods, orbital and millennial scale variabilities are not lumped in a sole anomaly pattern but differentiated. The method thus combines present-day observations, simulated Last Glacial Maximum (LGM) anomalies relative to present, scaled by an orbital-timescale index, and simulated stadial-interstadial anomalies, scaled by a millennial-timescale index: Here, T atm (t) and P (t) are the SAT and precipitation fields at time t. T atm 0 and P 0 are the ERA-INTERIM present-day SAT and precipitation climatologies (Dee et al., 2011). ∆T atm orb = T atm lgm −T atm pd and δP orb = P lgm /P pd are the orbital temperature anomaly and precipitation ratio relative to the present day (not shown, see Banderas et al. (2018)), respectively, obtained from 15 previous equilibrium simulations for the preindustrial and LGM climates performed with the CLIMBER-3α model (Montoya and Levermann, 2008). ∆T atm mil = T atm is − T atm st and δP mil = P is /P st are the millennial temperature anomaly and precipitation ratio, respectively, for the interstadial relative to the stadial state (Section 2.3). The key differences between these climate modes are that in the stadial, North Atlantic Deep Water (NADW) formation is relatively weak and takes place south of Iceland. Accordingly the sea-ice front in the North Atlantic reaches 40 • N. In the interstadial state there is a northward shift and 20 intensification of NADW formation. Northward oceanic heat transport increases, and the North Atlantic and surrounding areas warm relative to the stadial state, in particular the Nordic Seas. The simulated interstadial state is thus characterised by a more vigorous NADW formation and AMOC together with reduced sea ice in the Nordic Seas, and a temperature increase of up to 10 K in the North Atlantic relative to the stadial state, with a maximum anomaly in the Nordic Seas. Note bold symbols indicate two-dimensional spatial fields. The stadial mode in our study is represented by a climate simulation of the LGM with 25 CLIMBER-3α (Montoya and Levermann, 2008). The interstadial mode is taken from a recent glacial transient simulation performed with the same model under glacial climatic conditions, but with intensified NADW formation (Banderas et al., 2015).
α and β are two indices that separately modulate the contribution of the orbital and millennial anomalies. Both were built based on two recent complementary temperature reconstructions over Greenland, one from the NGRIP ice-core record for the LGP (Kindler et al., 2014), and the other one from several ice-core records for the Holocene (Vinther et al., 2009). Their 30 combination (hereafter, the KV reconstruction) results in a continuous temperature reconstruction for Greenland for the past 120 ka (Banderas et al., 2018). α is obtained after applying a low-pass frequency filter (f c = 1/18 ka −1 ) to the original KV reconstruction based on a spectral decomposition; β is obtained following a similar procedure but retaining the high frequency signal. Both indices are tuned in such a way that the resulting synthetic temperature time series at the NGRIP site exactly matches the KV reconstruction (this distinguishes α and β from the raw α and β indices previous to this tuning ;Banderas et al. (2018)).
The net basal melting rate for floating parts B is assumed to follow a linear relation: where T ocn is the oceanic temperature close to the grounding line, T f is the temperature at the ice base, assumed to be at the freezing point, and κ is the heat flux exchange coefficient between ocean water and ice at the ice-ocean interface. Several marine-shelf basal melting parameterizations can be found in the literature. The submarine melt rate is thought to be directly influenced by the oceanic temperature variations below the ice shelves. Accordingly, most basal melting parameterizations are built as a function of the difference between the oceanic temperature at the ice-ocean boundary layer and the temperature at 10 the ice-shelf base, generally assumed to be at the freezing point. The dependence on this temperature difference can be linear (Beckmann and Goosse, 2003) or quadratic (Holland et al., 2008;Pollard and DeConto, 2012;DeConto and Pollard, 2016;Pattyn, 2017). The linear marine-shelf basal melting parameterization used in this study is the simplest case that allows testing of the ice-sheet sensitivity to past oceanic temperature changes. Nevertheless, it accounts separately for sub-ice-shelf areas near the grounding line and for purely floating ice (ice shelves). The basal melting rate for purely floating ice shelves (B sh ) is 15 given by the grounding-line basal melt B gl scaled by a constant factor In this study, γ is set to 0.1. Thus, we consider that the submarine melting rate for ice shelves is 10 times lower than that close to the grounding zone, which is qualitatively in agreement with observations in some Greenland glaciers (Münchow et al., 2014;Rignot and Jacobs, 2002;Wilson et al., 2017). The melt rate in the open ocean, that is considered as being beyond the 20 continental shelf break, is prescribed to a high value (20 m a −1 ) to avoid unrealistic ice growth beyond 750 m of ocean depth, following Peyaud et al. (2007).
Following the approach described above, T ocn (t) is assumed to be given by an expression analogous to Eq. 1. Thus Eq. 3 can be rewritten as: The specific details of the experimental setup used are described below.

Experimental setup
We herein investigate the response of the EIS to millennial-scale climate variability during MIS 3. The starting point of our 30 experiments is a control-run ice-sheet simulation with constant bounday conditions for MIS 3 that provides a representative configuration of the EIS for that time period (Figure 1). To this end, α was set to its value at 40 ka BP, that is, α = α 40K = −0.1, and β = 0 to preclude millennial-scale variations. Note however these values are to a certain extent arbitrary; they are intended to provide a stable mean background state similar but not necccessarily identical to background MIS 3 conditions. Thus: Note that although Eq. 8 is formally correct and consistent with the scheme used, in contrast to the present-day SAT or precipitation the present-day rate of oceanic basal melting cannot be determined. Thus, in practice we replace this equation by directly tuning the value of B 40K to obtain a reasonable ice-sheet configuration at 40 ka BP given the atmospheric forcing 10 fields expressed by equations 6-7. To this end, a constant basal melting rate of 0.1 m a −1 is assumed. The ice sheet was forced with the resulting climatologies for 100 kyr previous to the starting of the perturbations described below. This allows the vertical temperature profile within the ice sheet to be equilibrated with the climate. This procedure was found to facilitate the growth of European ice-sheets to an extent that satisfactorily agrees with previous reconstructions (Svendsen et al., 2004;Kleman et al., 2013). 15 Our forcing method allows to investigate the response of the EIS solely to millennial-scale climate variability at MIS 3 by keeping constant the orbital component of the forcing (α = α 40K ) and letting β vary throughout the LGP (eqs. 1, 2 and 5). In order to assess the relative roles of the atmosphere and the ocean, three independent experiments have been carried out. First, an atmospheric-only forced simulation (ATM) in which the time evolution of SAT and precipitation on millennial time scales is considered, while the oceanic forcing is kept constant to MIS 3 (i.e., 40 ka BP) background climatic conditions. Thus: Second, an oceanic-only forced simulation OCN in which the atmospheric forcing is kept constant while the oceanic basal melting is allowed to vary at millennial timescales around its background MIS 3 value: The magnitude and sign of oceanic temperature anomalies ∆T ocn depends on the depth at which T ocn is considered. In our simulations, a large part of the NE sector of the EIS is marine based with shallow bedrock depths between 500 m and less than  Finally, a simulation ALL combining both the atmospheric and the oceanic forcings: In all experiments β (t) dictates the millennial-scale variability of the forcings (Figure 3). Because our simulated stadial-to-10 interstadial transition results from an intensification of the AMOC, positive β values imply an increase in T atm relative to its background MIS 3 value (e.g., Eq. 15 and Figures 2 and 3). As a consequence, the atmosphere warms at interstadials relative to stadial periods, as reflected by the ∆T atm mil millennial-scale anomaly field (Figure 2). Note that refreezing is not allowed to occur in our model in the current setup. If κ β ∆T ocn mil < −B 40K (which would imply B(t) < 0) we simply impose the value B(t) = 0.

Results
Substantial differences are found in the response of the EIS to the forcing scenarios. Under constant forcing, the CTRL run 20 shows negligible millennial-scale sea-level equivalent (SLE) variations, although a lower frequency SLE fluctuation is found as a result of internal ice-sheet variability (Figure 3). When the model is forced only by changes in sea level (SL run), a slight response is observed on millennial-scales. These changes appear not be sufficient to cause a substantial migration of the grounding line, thus not affecting ice velocities (not shown). In ATM, the atmospheric forcing alone causes a sequence of enhanced ablation episodes resulting in modest ice volume variations (up to 1.5 m SLE) during the most prominent stadial- 2d). Thus, the out-of-phase relationship found in the dynamic response of the EIS between these two oceanic experiments results from the opposed sign of their spatial forcing patterns (Figure 2). When considering the forcing at the subsurface of the ocean together with the atmosphere (ALLsub), slight reductions of the EIS volume (less than 1 m of s.l.e) during interstadials are superimposed onto the previous behavior (Figure 3).
The magnitude of these changes in terms of sea-level rise rate and discharge, specifically for the MIS 3 period, is illustrated 5 in Figure 4. The simulations forced with the surface of the ocean (OCNsrf and ALLsrf) show the largest amplitudes, with peaks of sea-level rise above 4 mm yr −1 during DO-events and sustained contributions well above 1 mm yr −1 during entire interstadial periods. In ATM, a decline of the EIS during stadial-to-interstadial transitions is still observed but presents a smaller amplitude of 1-2 mm yr −1 . The simulations in which the ice sheet is forced with the subsurface of the ocean (OCNsub and ALLsub) present a decline of their volume during stadial periods and regrowth during interstadials as a consequence of 10 the inverted spatial pattern of temperature anomalies with respect to the surface. In OCNsub (and ALLsub) the amplitude of these changes is smaller than in OCNsrf (and ALLsrf), on the order of 0.5-1 mm yr −1 , reaching more than 1 mm yr −1 during pronounced stadials (as ca. at 44 ka BP). The ALLsrf and ALLsub simulations show a similar or slightly larger volume loss during interstadials, as a consequence of the additional atmospheric forcing, that is superimposed onto the OCNsrf and OCNsub behaviour. 15 The response of the EIS has been analyzed in terms of its mass balance decomposition for the all-forcing runs ( Figure 5).
In ALLsrf the surface ocean temperature varies in phase with the atmosphere (Figure 2). Thus, during stadial-to-interstadial transitions the high negative values of dV/dt can be explained by the conjunction of an initial sharp increase in ablation together with pronounced increases in basal melting and calving, which allow for a large grounding line retreat in the Bjørnøyrenna basin ( Figure 5 mid panel). The rate of ice loss by basal melting is similar to that resulting from the increase in ablation (as 20 reflected in the surface mass balance, SMB) during the peak of a stadial-to-interstadial period. However, basal melting is much more efficient than surface mass balance in decreasing volume along the whole duration of an interstadial. This is due to the fact that ablation is restricted to the southern borders of the EIS. Thus, when the ice sheet has retreated to areas of no ablation, in spite of a slight further loss provided by the elevation feedback it rapidly equilibrates and a negative surface mass balance cannot propagate further inland. In contrast, when enhanced basal melting from higher oceanic temperatures is applied, the 25 associated retreat can propagate further inland occupying a large proportion of the Bjørnøyrenna basin and facilitating high rates of volume loss (although similar in amplitude with respect to SMB) during the whole interstadial period (see the animation corresponding to ALLsrf in the Supplementary Information). During stadial periods, both the enhanced positive mass balance and the absence of basal melting (favored by the negative oceanic anomalies) favor the regrowth of the EIS. Subsurface ocean temperatures evolve also in phase with the atmosphere in the SW part of the EIS but in anti-phase in its NE part. In other words, 30 when forcing with the subsurface of the ocean, a slight warming (cooling) is observed around the Britain ice sheet while cooling (warming) of the Bjørnøyrenna basin is simulated during interstadial (stadial) periods (see Figure 2). Therefore, the ALLsub simulation presents volume declines during stadial-to-interstadial transitions due to an increase in ablation and basal melting in the SW part. Subsequently, reduced basal melting in the NE part of the EIS favors regrowth of the Bjørnøyrenna basin during interstadial periods. Finally, shifting to pronounced stadial periods (as in ca. 44 ka BP) favors the penetration of warm subsurface waters that increase basal melting enough to produce an ice-sheet retreat in the NE part in spite of the enhanced positive surface mass balance (Figures 4 and 5). When considering the atmosphere and the subsurface ocean forcing together in ALLsub, these competing processes translate into a smaller amplitude of millennial-scale EIS changes as compared to the case with surface ocean forcing (ALLsrf). Furthermore, declines of the EIS can be observed both during the beginning of interstadial periods and during pronounced stadial periods in ALLsub (Figures 4 and 5). 5 Focusing on the OCN and ATM simulations separately facilitates isolating the effects of the ocean on this complex pattern. To this end, the simulated ice-sheet distribution and velocities of OCNsrf, OCNsub and ATM are shown in Figure 6 for the period around DO-event 12, at ca. 47 ka BP. As expected, OCNsrf shows a widespread retreat both in the NE and the SW of the EIS from the stadial to the interstadial period ( Figure 6, bottom left). This is accompanied by an acceleration of the Bjørnøyrenna basin due to its grounding line thinning and retreat ( Figure 6, left panels). OCNsub presents a collapsed Bjørnøyrenna basin 10 during the stadial period previous to DO-event 12 due to enhanced basal melting from warmer subsurface waters. The transition to the interstadial period favors a slight regrowth of this NE part of the EIS due to decreased basal melting, while its SW section slightly retreats (Figure 6 upper mid panel). Concerning ATM, only in the southwestern (SW) part of the EIS is the atmospheric forcing capable of generating an important reduction in the EIS volume in response to the stadial-interstadial transition ( Figure   6 right bottom panel). This is a result of the spatial pattern of the forcing, with the largest SAT anomalies located around the 15 Nordic seas (Figure 2). Therefore, the ice volume reduction of the EIS in ATM during interstadials is due to the positive SAT anomaly, which leads to enhanced ablation in the SW part of the EIS. In turn, reduced SATs during stadials allow the regrowth of the ice sheet up to the continental margin of the Nordic seas. The more active dynamic response of the EIS in the OCN simulations can be attributed to the increase in oceanic temperatures by 2-4 • C (Figure 2) within the margins of the ice sheet during interstadial (in the case of OCNsrf) and stadial (OCNsub case) periods, which translates into enhanced basal melting 20 at the margins of the EIS. The SW sector of the EIS also responds to the warmer SSTs, actually with a larger reduction of ice volume than in ATM ( Figure 6).
The spatial patterns shown in Figure 6 are representative of the all other stadial-to-interstadial transitions. In OCNsrf, the EIS reacts to every abrupt surface warming with a substantial ice-flow acceleration, especially in the Bjørnøyrenna basin ( Figure   7). Ice shelves that are present during stadial periods suddenly retreat during DO-events and together with enhanced basal 25 melting favor thinning and retreat of the grounding line that translate in large iceberg discharges up to ca. 0.06 Sv. In OCNsub, ice velocities in the Bjørnøyrenna basin increase during stadials, when enhanced basal melting erodes the grounding line and favors its retreat. Peaks in calving are recorded accordingly during pronounced stadial periods. These peaks are however of smaller amplitude than in OCNsrf. This can be explained by the fact that transitions to stadials are usually more gradual than transitions to interstadials, thus the incursion of warmer (subsurface) waters happens in this case in a smoother manner. High 30 velocities reach their maxima at the end of the stadial and beginning of the interstadials. The latter are however not accompanied by an increase in calving due to the fact that ice shelves are expanding and thickening during this period thanks to reduced basal melting (Figure 8). In general, the extension of ice shelves is greatly reduced during periods of enhanced basal melting (Figures 7, 8), with no large unconfined ice shelves surviving during these episodes. Some thinner ice shelves remain, in spite of the enhanced basal melting, thanks to an increase in advection from the Bjørnøyrenna ice stream triggered by a grounding line retreat (Figure 6).
Note that changes in the position of the calving front are usually accompanied by a grounding line displacement (not shown).
For some minor ice-shelf breakups this close relationship can be broken, but with almost no effects upstream inland. Thus we consider that the grounding line position is the best indicator for characterizing the dynamic behavior of the marine part of 5 the EIS. Inspection of the temporal evolution of the grounding line position in OCN simulations confirms that ice dynamics control the majority of ice-volume variations in the EIS as opposed to the SMB processes involved in ATM (Figure 9). The migration of the grounding line through time has been characterized by means of an index (µ) that weighs the proportion of non-grounded points in the region of the Bjørnøyrenna basin: where N g (t) represents the evolution of the number of points of grounded ice within a fixed area of N points in the Barents Sea region. An increase (decrease) in µ thus indicates a retreat (advance) of the grounding line. While in ATM µ barely changes ( Figure 9), OCN runs show a large dynamic behavior of the basin. In OCNsrf, µ reflects a synchronous evolution of the grounding line position and the oceanic forcing, with major retreats coinciding with interstadial states (Figure 9). Conversely, the Bjørnøyrenna basin is generally much closer to a full retreat in OCNsub during stadials due to a larger penetration of warm 15 subsurface waters (Figure 2; OCNsub) compared to the surface waters (Figure 2;OCNsrf). However, the grounding line is able to advance and reach Svalbard during episodes of reduced basal melting at the interstadials.
The direct coupling between the oceanic forcing and the response of the Bjørnøyrenna ice stream is also evident from the relatively high negative correlation (r −0.9) found between µ and ice thickness in this area (Figure 9). In essence, in response to the grounding-line retreat (

25
In order to investigate the sensitivity of the results to the model parameters, eight additional OCN simulations, both for the surface and the subsurface, have been carried out with different κ parameters between 1-10 m a −1 K −1 , i.e., bracketing our standard case of κ = 5 m a −1 K −1 . This choice reflects the inferences based on measurements made on Antarctic ice shelves that a variation of 1 K in the effective oceanic temperature changes the melt rate by ca. 10 m a −1 (Rignot and Jacobs, 2002;Shepherd et al., 2004). A robust response of the EIS is found, with a more reactive EIS response for increasing κ values Our results suggest a highly dynamic Eurasian ice sheet at millennial time-scales largely responding to changes in the ocean temperatures. Some authors present the marine based Kara-Barents complex as an analogue for present-day West Antarctic ice sheet for which bedrock topography is a major control for stability. We have shown, in this sense, that the Bjørnøyrenna basin is highly susceptible to changes in the oceanic temperatures. The timing of this response with respect to changes registered in 5 Greenland depends, however, on whether the surface or the subsurface of the ocean is considered as the relevant forcing of the ice sheet.
Recently, IRD peaks of Fennoscandian origin reported from a high-resolution marine sediment core from the Norwegian Sea indicate the presence of more frequent IRD deposition and thus calving during interstadials than during stadials (Dokken et al., 2013). This result has been corroborated in a compilation of new and previously published data (Becker et al., 2017) clearly 10 showing that within MIS 3, the IRD deposition increases within interstadials. The coeval deposition of carbonate-rich, sorted fine sands and near-surface warming suggests the presence of Atlantic water along the margin, and is interpreted by the authors as the effects of winnowing due to an intensified AMOC during interstadials. This interpretation results in concordance with our results when considering the surface waters as the oceanic forcing. Thus, this agreement would play in favor of considering that the EIS was primarily responding to changes in the surface of the ocean.

15
Our results also provide a mechanism to explain the pervasive presence of IRD in the North Atlantic during MIS 3, both during stadials and interstadials, and originating both in the LIS and the EIS. During stadials, the simultaneous appearance of IRD across the wider North Atlantic Ocean can be explained through the build-up of subsurface heat in the high-latitude North Atlantic leading to increased iceberg calving in the presence of large, thick ice shelves, together with lower surface temperatures allowing for wider dispersal of icebergs (Barker et al., 2015). According to our results interstadials could lead 20 to enhanced calving of the EIS through oceanic surface subglacial melting as a result of the warmer surface conditions and relatively shallow grounding line of this ice sheet.
The identification of IRD layers with increased calving through ice-sheet instabilities must be taken with caution, since it is based on several untested assumptions (Clark and Pisias, 2000): (i) delivery of IRD to a specific site is caused solely by iceberg calving, versus transport by sea-ice; (ii) an increase in IRD represents an increase in the iceberg flux, versus a greater amount 25 of debris incorporated at the base of the ice sheet that delivers the icebergs, or a greater distance of iceberg transport; (iii) the amount of IRD carried by all the icebergs is similar, assuming therefore a direct relationship between IRD concentration and iceberg flux. However, the former assumptions have not been confirmed and, thus, the calving-IRD relationship might not be so direct. In addition, ocean temperatures affect melting of icebergs and thus their release of IRD. Variations in ocean temperatures can alter the IRD released by an iceberg at a certain site, causing variations in IRD deposition even for a constant 30 amount of icebergs produced at the source.
Our experimental setup is not intended to match the paleorecord, but to provide insight into the response of the EIS to millennial-scale variability. The EIS variations simulated here represent the upper-end amplitude of potential responses during the whole glacial cycle, due to its large size. Extending the study to cover the whole LGP would require the consideration of orbital variability as part of the forcing. In this case, the EIS would likely be smaller during the mildest phase of MIS 3, thus limiting its contact with the ocean and the production of iceberg discharges.
Also, our results depend somewhat on the particular SAT and oceanic temperature anomaly patterns simulated by our climate model, the magnitudes of the resulting forcing, and the initial size of the simulated EIS. The use of different atmospheric realisations is subject to the availability of climate simulations with different models for the three climate states needed: glacial 5 (stadial), present, and interstadial. The latter is only available for a reduced number of models. This makes the assessment of this issue difficult in the present study. Assessing the sensitivity to these features should be in the scope of future work, and illustrates the need for carrying out new simulations of both the interstadial and the stadial states with more sophisticated climate models. Nonetheless, our results indicate that the ocean is the major driver of the EIS ice-volume changes during MIS-3. Note the temporal index used is the same for the atmosphere and the ocean and the amplitude is given by an OGCM 10 simulation of two different oceanic states mimicking stadial and interstadial periods. We then translate those fields into ablation (through PDD, whose uncertainty has been extensively explored) and into basal melting (through a linear equation). The values of the oceanic sensitivity parameter (κ) we used here are in the range (or even below in most cases) of those suggested by data in Antarctica (Rignot et al. 2002). Note, in particular, that even from low-mid values of κ of 2 m a K −1 the response to the ocean begins to be of greater amplitude than that to the atmosphere, making our main conclusions robust. 15 Finally, our study lacks bi-directional coupling between the ice sheet, the atmosphere and the ocean. Eventually the goal is to investigate this matter with fully coupled climate-ice sheet models.
Meltwater discharge from the EIS and other ice sheets surrounding the Nordic Seas is often implied as a cause of ocean instabilities. The same would be the case for iceberg discharges. This issue is beyond the scope of this study; its assessment would require investigating the impact of these freshwater perturbations in deep water formation and the AMOC. Again, proper 20 assessment requires the use of a coupled climate-ice sheet model.

Conclusions
We have investigated the response of the EIS to millennial-scale climate variability associated with DO-events through a series of simulations with a three-dimensional, hybrid ice-sheet model that represents inland ice flow under the SIA and floating ice shelves and ice streams through the SSA. The model makes use of an offline forcing method that separately accounts for 25 orbital and millennial-scale climate variability during the LGP, improving the representation of the latter (Banderas et al., 2018).
Atmospheric and ocean forcings associated with millennial-scale variability were considered both separately and together.
Separating the effects of atmospheric and oceanic forcing during the glacial period has allowed us to quantify the contribution of each to EIS variability. Atmospheric forcing during stadial-interstadial transitions has a modest effect on the ice sheet, which is a consequence of the largest SMB changes being confined to SW sector of the EIS, where the forcing is strongest. 30 In contrast, the oceanic forcing has a larger effect, through changes in the ice dynamics in the Bjørnøyrenna basin of the EIS.
Ocean warming is able to induce a retreat of grounded ice in this part of the EIS through dynamic processes. As a consequence,  after the spinup was completed (right). This ice sheet represents the initial state previous to the applied perturbations. Bjørnøyrenna basin, as referenced in the text, is shown by the black rectangle.