The Antarctic Ice Sheet response to glacial millennial-scale variability

The Antarctic Ice Sheet (AIS) is the largest ice sheet on Earth and hence a major potential contributor to future global sea-level rise. A wealth of studies suggest that increasing oceanic temperatures could cause a collapse of its marine-based western sector, the West Antarctic Ice Sheet, through the mechanism of marine ice-sheet instability, leading to a sea-level increase of 3–5 m. Thus, it is crucial to constrain the sensitivity of the AIS to rapid climate changes. The last glacial period is an ideal benchmark period for this purpose as it was punctuated by abrupt Dansgaard–Oeschger events at millennial timescales. Because their center of action was in the North Atlantic, where their climate impacts were largest, modeling studies have mainly focused on the millennial-scale evolution of Northern Hemisphere (NH) paleo ice sheets. Sea-level reconstructions attribute the origin of millennial-scale sea-level variations mainly to NH paleo ice sheets, with a minor but not negligible role of the AIS. Here we investigate the AIS response to millennial-scale climate variability for the first time. To this end we use a three-dimensional, thermomechanical hybrid, ice sheet–shelf model. Different oceanic sensitivities are tested and the sea-level equivalent (SLE) contributions computed. We find that whereas atmospheric variability has no appreciable effect on the AIS, changes in submarine melting rates can have a strong impact on it. We show that in contrast to the widespread assumption that the AIS is a slow reactive and static ice sheet that responds at orbital timescales only, it can lead to ice discharges of around 6 m SLE, involving substantial grounding line migrations at millennial timescales.

1 Introduction

The Antarctic Ice Sheet (AIS) presently stores around 60 m of potential sea-level rise (Fretwell et al., 2013). It is divided into two parts, the East Antarctic Ice Sheet (EAIS) and the West Antarctic Ice Sheet (WAIS), including the Antarctic Peninsula (AP). Present-day observations show that the mass balance of the AIS is negative due to mass loss from the WAIS, whereas the EAIS maintains a positive mass balance (Martín-Español et al., 2016; Shepherd et al., 2018). Because ablation in the AIS is almost negligible except in the small region of the AP, the mechanisms that contribute to mass loss are submarine melting of floating ice shelves and calving processes at the ice front (Paolo et al., 2015; Rignot et al., 2013). The WAIS is a marine ice sheet, i.e., most of it is grounded below sea level, and it contains several large ice shelves that are thinning or calving more rapidly than the storage provided by surface accumulation. The positive mass balance of the EAIS can be explained by the fact that the amount of floating ice is considerably smaller than in the WAIS, and thus the mass loss via calving and basal melting does not surpass the accumulation.

Rising oceanic temperatures in the coming century in response to climate change can boost basal melt and reduce ice shelves. Although thinning of floating ice shelves does not directly contribute to sea-level rise, it can lead to a reduction of ice-shelf buttressing, enhancing inland ice flow as seen after the collapse of the Larsen B ice shelf (Fürst et al., 2016; Rignot et al., 2004) and Pine Island Glacier (Favier et al., 2014; Jacobs et al., 2011). In addition, most parts of the WAIS lie on a retrograde bed slope. Conceptual models suggest the existence of an inherent instability in such ice sheets, the marine ice-sheet instability (MISI; Schoof, 2007; Weertman, 1974), that could lead to a collapse of the marine grounded zones in the WAIS region. Mercer (1978) speculated about the fact that this instability could be triggered through a rise in oceanic temperatures. Collapse of the WAIS sector could cause a sea-level increase of 3–5 m (Bamber et al., 2009; Feldmann and Levermann, 2015; Sutter et al., 2016), with major implications for coastal zones (Nicholls and Cazenave, 2010). From a modeling perspective, projections differ considerably in future sea-level contributions depending on the model used and the process parameterizations therein (Bakker et al., 2017a, b; DeConto and Pollard, 2016; Golledge et al., 2015).

Improving our understanding of the AIS sensitivity is thus essential to constrain future projections (Bakker et al., 2017a). Some of the most remarkable abrupt climate changes of the near past are those of the last glacial period (LGP; 110–10 ka). Thus, one way to gain insight in this respect is to assess the response of the AIS to these past rapid climate changes. In addition, understanding the AIS behavior during these millennial-scale abrupt events will help in identifying the ultimate causes of the Dansgaard–Oeschger (DO) events. Ice-core records from the Greenland Ice Sheet (GrIS) during the LGP show the characteristic signal of DO events: a rapid warming of more than 10 K on decadal timescales followed by a slow cooling that can last from several centuries to thousands of years (e.g., Dansgaard et al., 1993). Modeling studies (e.g., Ganopolski and Rahmstorf, 2001; Rahmstorf, 2002; Shaffer et al., 2004) and reconstructions (Barker et al., 2015; Böhm et al., 2015; Henry et al., 2016; McManus et al., 2004) support the hypothesis that DO events were caused by reorganizations of the Atlantic Meridional Overturning Circulation (AMOC), with enhanced (reduced) North Atlantic Deep Water (NADW) formation during interstadials (stadials) transporting more (less) heat into high northern latitudes. In addition, marine sediment records across large areas of the North Atlantic show quasi-periodic deposition of ice-rafted detritus (IRD) (Hemming, 2004) known as Heinrich (H) events. H events are thought to have been caused by massive iceberg discharges from the paleo Laurentide Ice Sheet (LIS), possibly in response to reductions in NADW formation that, through positive feedbacks, resulted in the collapse of the AMOC (Alvarez-Solas et al., 2011, 2013; Marcott et al., 2013).

Compared to ice-core records in the GrIS, AIS ice-core records show a more gradual and symmetric sawtooth-like signal throughout the whole LGP. An increase in surface air temperature (SAT) is observed during Greenland stadials, most notably during Heinrich stadials, with cooling during interstadials. The amplitude of this signal can reach up to 2 K (Augustin et al., 2004; Petit et al., 1999; Ruth et al., 2007) and the peaks of the sawtooth signal are known as Antarctic isotope maxima (AIM). This bipolar seesaw behavior between Greenland and Antarctica is now well established (Blunier and Brook, 2001; EPICA Community Members, 2006). The paradigm to explain it is that intensifications of the AMOC translate into an increase in northward heat transport at the expense of the southernmost latitudes; conversely, a weakening of the AMOC reduces northward heat transport, thereby warming the south (Crowley, 1992; Stocker, 1998). The different timescale between northern and southern latitudes can be explained by the fact that the Southern Ocean (SO) acts as a heat reservoir that dampens and integrates in time the more rapid North Atlantic signal (Stocker and Johnsen, 2003). The occurrence of H events supports a high sensitivity of Northern Hemisphere (NH) ice sheets as well as their capability to react rapidly (Alvarez-Solas et al., 2013, 2017; Andrews and Voelker, 2018; Hemming, 2004). In the Southern Hemisphere (SH), data showing IRD deposition from the AIS are more scarce. There is evidence of ice discharges from the AIS (Kim et al., 2018; Weber et al., 2012, 2014), but neither a quantification of their contribution in terms of its sea-level equivalent (SLE) nor the identification of their triggering mechanism has yet been done, particularly for events during Marine Isotope Stage 3 (MIS-3). If a periodic deposition of IRD could be found in the SH analogous to the NH, it may hint at an Antarctic response to oceanic changes. This would consolidate the mechanism of the bipolar seesaw and the existence of the heat storage in the SO.

Finally, sea-level reconstructions show fast variations of more than 20 m at millennial timescales during MIS-3 (Frigola et al., 2012; Grant et al., 2012; Rohling et al., 2014) and rises of 4 m per century during meltwater pulse (MWP) 1A at ca. 14.5 ka (Liu et al., 2016) . However, the individual contribution of each paleo ice sheet remains unclear. Due to their location at lower latitudes compared to the AIS, NH ice sheets are more exposed to mass losing processes through atmospheric forcing (ablation). Therefore the majority of those rapid changes are thought to originate in the NH ice sheets (Arz et al., 2007; Ganopolski et al., 2010). However, during MIS-3 sea-level variations fluctuated on the Antarctic rhythm (Grant et al., 2012; Rohling et al., 2009; Siddall et al., 2008), suggesting that a considerable contribution from direct AIS waxing and waning cannot be excluded.

As far as we know, there have been no attempts to simulate Antarctic sea-level contributions at millennial timescales and their potential implications. The aim of this paper is thus to investigate the response of the AIS to millennial-scale variability during the LGP. In particular, we focus on the AIS advance and retreat and its potential sea-level contribution at these timescales. Some assumptions are made for the sake of simplicity, since our aim is to test if the AIS is likely to have responded at millennial timescales and to what extent. For this purpose we use a three-dimensional, thermomechanical, ice sheet–shelf model that is forced through a synthetic climatic forcing including both atmospheric and oceanic changes that evolve temporally through an index deduced from the Dome C deuterium ice-core record. To study the impact of ice–ocean interactions we use a basal melting parameterization that is a function of oceanic temperature anomalies.

The paper is structured as follows: first, the ice-sheet model, the forcing, and the experimental design are described (Sect. 2). Then the response of the AIS to the oceanic forcing is shown, focusing on the ice discharges and grounding line advances at millennial timescales (Sect. 3). Finally, the main results are discussed (Sect. 4) and conclusions summarized (Sect. 5).

2 Methods and experimental setup

2.1 Model

We use the three-dimensional, hybrid, thermomechanical ice-sheet model GRISLI-UCM based on the GRISLI model developed by Ritz et al. (2001) and further extended and tested at the Complutense University of Madrid (see Alvarez-Solas et al., 2017; Tabone et al., 2018). Important changes with respect to the original code include variations in boundary conditions (surface mass balance and basal melt), topography, and new auxiliary modules to calculate the basal drag. Simulations are run on a 40 km×40 km grid with 21 vertical layers corresponding to 157×147 grid points covering the whole Antarctic domain. Initial topographic conditions (ice thickness, surface, and bedrock elevation) are provided from the dataset RTopo-2 (Schaffer et al., 2016), which relies on Bedmap2 (Fretwell et al., 2013) with corrections for ice-shelf cavities. The grounded slow-moving ice, whose flow is dominated by shear processes, is computed by the non-sliding shallow ice approximation (SIA), whereas floating ice shelves, whose evolution is determined by stretching processes, are solved by the shallow shelf approximation (SSA) (Hutter, 1983; MacAyeal, 1989). Intermediate states, in which shearing and stretching regimes can appear simultaneously, are typical of fast-flowing ice streams and are evaluated by summing the velocities of the SIA and SSA. The SSA solution allows for basal sliding and thus includes basal drag depending on the topographic conditions. The model allows basal sliding when the ice base (land–ice interface) is at the melting point and the pressure of the basal water exceeds an imposed threshold.

The total mass balance is given by the difference between accumulation and ablation at the surface, melting at the base of the ice sheet, and ice discharge into the ocean via calving. The surface mass balance (SMB) is determined by atmospheric temperature and precipitation using the positive degree-day scheme (Reeh, 1989). The geothermal heat flux applied as a boundary condition to grounded ice is obtained from the field provided by Shapiro and Ritzwoller (2004). Submarine melt is determined through a linear equation, which transforms oceanic temperature anomalies into melting rates through a heat flux coefficient (details in Sect. 2.2). Calving occurs when the ice-shelf front grid point gets thin enough (200 m) and the incoming ice from upstream does not maintain the necessary ice thickness (Peyaud et al., 2007).

2.2 Forcing method and experimental design

GRISLI-UCM is forced through the same parameterization for atmospheric and oceanic forcing as in Banderas et al. (2018) and Tabone et al. (2018), who used it to specifically investigate the past evolution of the glacial NH and Greenland ice sheets, respectively, but here for the Antarctic domain. In the more general approach used in those studies, oceanic, atmospheric, and precipitation fields are scaled by two climatic indices, an orbital index α(t) (where α=0 represents the LGM state and α=1 the present day, PD), and a millennial index β(t) (β=0 at the LGM, β=1 at the AIM). Because our study focuses on millennial-scale variability, we fix α=0 to maintain constant glacial background conditions. The β index is extracted from the Dome C atmospheric temperature reconstruction (Jouzel and Masson-Delmotte, 2007) and is filtered between 1 and 19 ka to avoid both orbital and submillennial-scale variability. The time evolution of atmospheric temperature (T^atm(t)) and precipitation (P(t)) fields is given by the following equations:

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }{T}^{\mathrm{atm}}\left(t\right)=T_{\mathrm{LGM}}^{\mathrm{atm}}+\mathit{\beta }\left(t\right)\mathrm{\Delta }{T}_{\mathrm{mil}}^{\mathrm{atm}},\text{(2)}& {\displaystyle }& {\displaystyle }P\left(t\right)=P_{\mathrm{LGM}}\left[\left(\mathrm{1}-\mathit{\beta }\left(t\right)\right)\mathit{\delta }{P}_{\mathrm{orb}}+\mathit{\beta }\left(t\right)\mathit{\delta }{P}_{\mathrm{mil}}\right],\end{array}$$

where temperature and precipitation, $T_{\mathrm{LGM}}^{\mathrm{atm}}$ and P_LGM, respectively, are the LGM climatologies calculated from the ERA-Interim reanalysis (Dee et al., 2011) and corrected with orbital anomaly fields obtained from the climatic model of intermediate complexity CLIMBER 3-α (Montoya and Levermann, 2008). The millennial ($\mathrm{\Delta }{T}_{\mathrm{mil}}^{\mathrm{atm}}$, δP_mil) anomaly fields are obtained from the same climatic model.

The parameterization of the submarine melting rate under floating ice shelves follows a simple linear law based on Beckmann and Goosse (2003):

$$\begin{array}{}\text{(3)}& B=\mathit{\kappa }\left(T^{\mathrm{ocn}}-T_{\mathrm{f}}\right),\end{array}$$

where T^ocn is the oceanic temperature at the corresponding grid point, T_f the freezing point temperature at which the ice base is assumed to be, and κ the heat flux exchange coefficient between ocean and ice. Other possible choices are, for example, a quadratic approach (DeConto and Pollard, 2016; Pattyn, 2017; Pollard and DeConto, 2009). For the sake of simplicity, we assume a linear response between oceanic temperatures and melting rates, which was already tested previously (Alvarez-Solas et al., 2013; Golledge et al., 2015; Philippon et al., 2006; Tabone et al., 2018). The model distinguishes between basal melting at the grounding line (B_gl) and below the ice shelf (B_shlf).

$$\begin{array}{}\text{(4)}& B_{\mathrm{shlf}}=\mathit{\gamma }{B}_{\mathrm{gl}}\end{array}$$

Rignot and Jacobs (2002) have shown that melting rates at the ice shelves are about an order of magnitude lower than those close to the grounding line, and hence we set γ to 0.1. Following the same procedure as for the atmospheric forcing, the oceanic temperature can be rewritten as

$$\begin{array}{}\text{(5)}& B_{\mathrm{gl}}\left(t\right)=B_{\mathrm{LGM}}+\mathit{\kappa }\mathit{\beta }\left(t\right)\mathrm{\Delta }{T}_{\mathrm{mil}}^{\mathrm{ocn}},\end{array}$$

where B_LGM represents LGM melting rates and $\mathrm{\Delta }{T}_{\mathrm{mil}}^{\mathrm{ocn}}$ the millennial oceanic temperature anomaly. To avoid any accretion at the ice-shelf base, B_gl cannot become lower than 0 m a⁻¹.

To study the response of the AIS to millennial-scale variability alone, we spun up our model for 120 ka under fixed LGM conditions. Figure 1 illustrates the surface elevation and velocities after the spin-up procedure. We then impose the millennial-scale forcing. The oceanic temperature field and its resulting basal melt rates at the LGM, B_LGM, are complicated to obtain due to lack of proxy data. Moreover, B_LGM strongly determines the ice extent of the AIS during the LGM. Observations and reconstructions suggest that the ice sheet advanced to the continental shelf break at the LGM (Anderson et al., 2002; Bentley et al., 2014; Denton and Hughes, 2002; Hillenbrand et al., 2012; Kusahara et al., 2015; Whitehouse et al., 2012). Setting B_LGM= 0 m a⁻¹ (see Fig. 2a) allows for such an advance. In regions with ocean depths below 2000 m, an artificially large melting rate (50 m a⁻¹) is prescribed to avoid unrealistic ice-shelf growth beyond the continental slope, which would likely be subject to high melt rates in reality because of the intrusion of warm circumpolar deep waters into the ice-shelf cavities (Kusahara et al., 2015). The millennial-scale oceanic temperature anomaly is then obtained from the Dome C ice core (Jouzel and Masson-Delmotte, 2007): the LGM minus present atmospheric temperature at Dome C is estimated to be ca. −10 K and the maximum amplitude of AIM events ca. 2 K. Following Collins et al. (2013) and Golledge et al. (2015), the oceanic amplitude of temperature change is estimated to be up to one-fourth that of the air temperature change, and thus $\mathrm{\Delta }{T}_{\mathrm{orb}}^{\mathrm{ocn}}=-\mathrm{2.5}$ and $\mathrm{\Delta }{T}_{\mathrm{mil}}^{\mathrm{ocn}}=\mathrm{0.5}\mathrm{K}$. Oceanic temperature variations are applied uniformly in space. Figure 3a illustrates the index used for the perturbation. To assess the impact of the ice–ocean interaction we test different oceanic sensitivities. Thus, κ goes from no ice–ocean interaction (0 m a⁻¹ K¹) to a large sensitivity (15 m a⁻¹ K¹). All values of the tested parameters are provided in Table 1. Finally, sea-level variations are prescribed from Rohling et al. (2014).

Figure 1 Simulated ice-sheet (a) surface elevation (in km) and (b) ice velocities (in m a⁻¹) after the spin-up procedure. The thick black line indicates the simulated grounding line position. The thick grey line represents the continental shelf break (depth 2000 m).

Table 1 Summary of the studied parameter values used in each sensitivity test.

3 Results

In this section we present our main results focusing on the AIS response to oceanic changes (Fig. 3a) in terms of its SLE contributions (Fig. 3b) and grounding line migrations (Fig. 3c) at millennial timescales. When ignoring the interaction with the ocean (κ=0 m a⁻¹ K¹; dark blue curve), no SLE changes are observed, implying that the effect of the atmospheric forcing (temperature and precipitation variations) is negligible. When the oceanic forcing is considered, ice volume subsequently displays millennial-scale variations. The amplitude of these variations increases with increasing oceanic sensitivities (κ values). As long as the climatic index β stays positive, heat is transferred from the ocean to the AIS, ice is discharged from the ice sheet to the ocean, and the grounding line experiences migrations at millennial timescales. When the index becomes negative, the submarine melting is set to zero. In this way oceanic temperatures are assumed to remain close to the freezing point and no accretion is allowed; the ice-sheet volume grows through net accumulation and the ice sheet expands.

Figure 2 Mask used to evaluate grounding line migration. (a) Ice extent after glacial spin-up and (b) PD ice extension. Blue zones are model grid cells with grounded ice in marine zones. Grey zones are model grid cells without grounded ice in marine zones but the underlying bathymetry is shallow enough to potentially become grounded (i.e., marine zones with depths less than 2000 m). The thick black line indicates the grounding line position.

To quantify the grounding line migration we introduce a parameter called marine zone occupation (MZO), which is defined as

$$\begin{array}{}\text{(6)}& \mathrm{MZO}={\displaystyle \frac{N_{\mathrm{G}}}{N_{\mathrm{G}}+N_{\mathrm{P}}}},\end{array}$$

where N_G is the number of model grid cells with grounded ice in marine zones (i.e., zones in which the ice is grounded and its bedrock lies below sea level; see blue zones in Fig. 2a, b) and N_P is the number of grid cells of floating ice in marine zones that could potentially become grounded (i.e., zones in which the ice is not grounded but floating and where the underlying bathymetry is shallow enough to potentially become grounded; in practice, we identify these as marine zones with depths above −2000 m; see grey zones in Fig. 2a, b). Therefore, if MZO=1, the grounding line has advanced up to the continental shelf break, grounding all possible marine zones. If MZO is below 0.21, which corresponds to present-day (PD) values (Fig. 2b), the grounding line position has retreated beyond its PD limit. Finally, if MZO=0, the grounding line has entirely retreated up to the land with its bedrock fully above sea level (i.e., marine zones disappear). Figure 3c shows the evolution of the MZO for different oceanic sensitivities. After the spin-up, MZO=0.73 (Fig. 2a). The grounding line has thus advanced towards the continental shelf break but shelves like the Pine Island zone or George Land remain ungrounded (Fig. 1a). For κ=0 the position of the grounding line does not evolve away from the spin-up value. Only when the oceanic forcing is considered do grounding line migrations begin to be appreciable. When oceanic variability is considered, our modeled AIS reacts at millennial timescales.

Figure 3 (a) Millennial-scale forcing index (β). On the right-hand side the equivalent oceanic temperature anomaly is shown (in K). (b) Ice volume (in 10⁶ km³) and SLE contribution (in m). (c) MZO evolution for different oceanic sensitivities. Colors go from no ice–ocean interaction ($\mathit{\kappa }=\mathrm{0}\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$, dark blue) to large oceanic sensitivity ($\mathit{\kappa }=\mathrm{15}\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$, red). The solid grey line in (b) and (c) indicates the present-day value of the ice volume and MZO, respectively.

Figure 4 illustrates the surface elevation (a) and ice velocity (b) for three different oceanic sensitivities (κ=1, 5, and 10 $\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$) after a typical cold phase (at 61 ka). The configuration in the three cases is similar, with an advanced grounding line with grounded Ronne and Ross embayments. The grounding line retreat in the Ross shelf increases with increasing κ. Ice streams also penetrate further inland with increasing κ. Figure 5 illustrates the same fields after an AIM event (at 57 ka). While the lowest sensitivity case (κ=1 m a⁻¹ K¹) shows an extensive ice sheet close to the continental shelf break similar to the initial LGM state, as the sensitivity increases (κ=5 m a⁻¹ K¹) marine zones such as the Ronne ice shelf begin to retreat and velocities increase. For sufficiently high oceanic sensitivities (κ=10 m a⁻¹ K¹) the Ronne and Ross ice shelves experience a substantial retreat during AIM events. In addition, the ice velocity field shows ice streams penetrating further inland with increasing κ. The ice thickness difference between these two snapshots highlights the particular embayments for which the AIS is discharging for increasing ice–ocean sensitivities (see Fig. 6). The majority of the ice loss comes from the Ronne shelf as it is the most vulnerable zone to oceanic warming. The Ross shelf does not experience a substantial ice loss and grounding line retreat until κ>=10 m a⁻¹ K¹. The Pine Island zone responds in a similar manner to the oceanic warming but with less impact. Grounding line migrations and ice discharges are not restricted to the WAIS but also occur in the coastal zones of the EAIS, which goes all along the Amery shelf down to Wilkes Land.

Figure 4 Snapshots of the AIS simulations at a cold phase (61 ka) for three different oceanic sensitivities (κ=1, 5, and 10 $\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$). (a) Surface elevation (in km). The thick black line indicates the grounding line position, and the thick grey line is the continental shelf break. (b) Ice velocities (in m a⁻¹).

The longest ice regrowth periods, corresponding to cooling phases, happen between 70 and 60 ka and between 40 and 20 ka. During these periods, for medium to low sensitivities (up to κ=7 m a⁻¹ K¹), the grounding line position (as indicated by the MZO) advances close to its LGM value, whereas for high oceanic sensitivities the maximum MZO value reached decreases with increasing κ, indicating the irreversibility typical of hysteresis behavior (Fig. 3c). This suggests that the grounding line can readvance up to the continental shelf break if the oceanic forcing is suppressed long enough, which is not the case for large κ.

We further assess what determines the amplitude of ice discharges between 75 and 15 ka (Fig. 7a). During this time period we find six significant ice discharge events in response to enhanced submarine melting phases, marked with grey shading. Figure 7b shows the ice-volume loss and its corresponding sea-level contribution with respect to κ for every event. Again, for no ice–ocean interaction (κ=0 m a⁻¹ K¹) no ice discharges are found, implying that atmospheric millennial variability alone can not produce sea-level variations in the AIS. As the ice–ocean interaction increases with increasing κ, not only does the sea-level contribution of every event increase, but also a wider spread is found between the discharging events, meaning that the sea-level difference between the smallest and largest ice discharge increases. Finally, what determines the total amount of sea-level rise of an AIM event is the total heat exchange between ice and ocean (Fig. 8c). If the amplitude is large, generally major ice discharges will be likely, but if the time interval is too short, then this will not necessarily be true (Fig. 8a). The same is true for the AIM event duration: longer periods will have more potential time to discharge ice, but if the warming is smooth, less melting and ice retreat will happen (Fig. 8b).

Figure 5 Snapshots of the AIS simulations at the end of a warm phase (AIM) event (57 ka) for three different oceanic sensitivities (κ=1, 5, and 10 $\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$). (a) Surface elevation (in km). The thick black line indicates the grounding line position, and the thick grey line is the continental shelf break. (b) Ice velocities (in m a⁻¹).

Figure 6 Ice thickness difference between the AIM and the cold phase (AIM minus cold) for different values of oceanic sensitivity (κ=1, 5, and 10 $\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$). Zones with an intense red color illustrate a larger ice difference and hence a major ice loss. The thick blue line illustrates the grounding line position at the cold phase and the thick yellow line the grounding line position at the AIM phase. The thick grey line illustrates the position of the continental shelf break.

4 Discussion

Our experimental design follows the bipolar seesaw mechanism (Crowley, 1992; Stocker and Johnsen, 2003) according to which the SO acts as a heat reservoir during millennial-scale AMOC reorganizations. However, the extent to which the SO temperature increases during the slowdown of the AMOC is under debate. Pedro et al. (2018) have argued that the Antarctic Circumpolar Current (ACC) acts as a barrier for heat penetration into the SO and that the postulated heat reservoir is rather provided by the southern subtropical Atlantic and transferred to the AIS by the atmosphere; in addition, oceanic heat transport changes could be compensated for to a large extent by changes in heat transport by the atmosphere and the Pacific Ocean. Changes in SO overturning and/or convection can lead to much larger, albeit localized, warming (e.g., Martin et al., 2013, 2014). Positive feedbacks resulting from sea-ice and ice-shelf melting could further increase warming of the subsurface through enhanced stability of the water column (Weber et al., 2014).

For the sake of simplicity we also considered a spatially homogeneous oceanic warming in phase with the atmospheric temperature reconstruction of Dome C. We deduced the oceanic temperature anomaly from the atmospheric reconstruction of the Dome C ice core. This results in an oceanic temperature anomaly during AIM events of about 0.5 K. To our knowledge, there are no reconstructions available for the SO temperature of high enough temporal resolution. A lower amplitude for the oceanic temperature anomaly in our experimental setup would diminish the effect of the millennial-scale oceanic temperature variability. Nevertheless, our heat transfer coefficient κ can also be interpreted as a weighting parameter of the amount of heat transferred into the SO. However, Buizert et al. (2015) argue that the timing difference between the occurrence of DO events in Greenland ice cores and AIM events provides support for a slow (oceanic) versus a fast (atmospheric) propagation mechanism from north to south. Hence the main question of how much the SO warms up during AIM events is unclear and, again, requires oceanic temperature reconstructions that are yet not available.

Figure 7 (a) Simulated ice volume anomaly between 75 and 15 ka for different values of oceanic sensitivities. Anomalies are calculated relative to the state at 15 ka and detrended between 75 and 15 ka. Grey illustrates significant ice discharging events with increasingly darker grey colors for older events. (b) Scatterplot of the sea-level contribution of every discharging phase with respect to κ.

Figure 8 Ice-volume discharge and SLE contribution of every event against (a) the amplitude of the warming, (b) the duration of the warming phase, and (c) the integrated warming defined as the peak warming times the duration. Colors represent the different ice–ocean sensitivities.

We also found that if the heat flux transfer parameter between ice and ocean is larger than or equal to 10 m a⁻¹ K¹ then the ice sheet is not able to regrow to its initial state after spin-up, neither in volume nor in extent. This highlights the possibility that a heat flux parameter of 10 m a⁻¹ K¹ is maybe too large for our ice-sheet model as we know that during the LGM the AIS reached its maximum size from reconstructions.

Here we simulated the grounding line migration at millennial timescales for different oceanic sensitivities. We observed that at those relatively short timescales, the grounding line is capable of advancing to its initial state after retreating. Although here we mainly focus on ice-sheet dynamics, we think this variability could be relevant for brine rejection over the continental shelf as proposed by Paillard and Parrenin (2004). The underlying mechanism is that during grounding line advances, brine (salty water released during sea-ice formation) is pushed out of the continental shelf break. This salty water descends to the bottom of the ocean, having a strong impact on the carbon exchange. If a millennial oscillation of the grounding line took place, it could explain the rise of carbon into the atmosphere, which may be a potential explanation for DO events as well as glacial–interglacial shifts at orbital timescales.

Sea-level reconstructions during MIS-3 show millennial fluctuations that can reach more than 20 m SLE. These sea-level differences are generally attributed to paleo NH ice sheets (Arz et al., 2007). Our results highlight the possibility that a warming of the SO can have a strong impact on the AIS, producing substantial ice discharges. None of our results, including those with a high oceanic sensitivity, exceeded 20 m SLE. Low sensitivities (κ<5 m a⁻¹ K¹) do not produce discharging events of more than 5 m, which means that NH paleo ice sheets would still be the major contributors to millennial sea-level fluctuations. For κ>10 m a⁻¹ K¹, SLE contributions of more than 10 m occur, which would imply a significant Antarctic contribution as well. However, as discussed above, this contribution (and for larger oceanic sensitivities) seems unrealistic as our model does not support a regrowth of the AIS to the continental shelf break under LGM climate conditions. Intermediate values (κ=7 m a⁻¹ K¹) lead to discharges of around 6 m SLE. A non-negligible Antarctic contribution to sea-level changes at millennial timescales during the LGP will have an impact on reconstructing the size of other paleo ice sheets.

5 Conclusions

We have investigated the response of the AIS to millennial-scale climate variability and, in particular, its response to different oceanic sensitivities using a hybrid, three-dimensional, thermomechanical ice-sheet model. The model is forced using a method that has already been tested (Banderas et al., 2018) and is provided by an improved subglacial melting routine. Because SO temperature reconstructions are not available we assumed that oceanic temperatures covary with atmospheric temperature variations at millennial timescales based on Stocker and Johnsen (2003). Our simulations suggest that, contrary to the idea that the AIS is a slow reactive ice sheet, it could be more reactive to millennial-scale climate variabilities than previously thought. We found that whereas atmospheric millennial-scale variability had no appreciable impact on the AIS, SO warming could produce episodes of ice discharge, leading to substantial sea-level rise and grounding line migration. Although this timescale may seem short for such a large ice sheet, our simulations show, in the range of realistic values for oceanic sensitivities, that considerable grounding line retreat in the Ronne, Ross, and Wilkes Land embayment, as well as sea-level discharge of around 6 m SLE at millennial timescales, can occur. Our results highlight the possibility that, via the bipolar seesaw, a slowdown of the AMOC could have accumulated more heat in the Southern Ocean, resulting in significant sea-level rise produced by the AIS on millennial timescales.