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 marinebased 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. 5 Because their centre of action was in the North Atlantic, where their climate impacts were largest, modelling 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 to 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 10 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 almost 15 m of SLE involving substantial grounding line migrations at millennial timescales.

that the amount of floating ice is considerably smaller than in the WAIS, 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) 5 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, Weertman (1974); Schoof (2007)), 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 meters (Bamber et al., 2009;Feldmann and Levermann, 2015;10 Sutter et al., 2016), with major implications for coastal zones (Nicholls and Cazenave, 2010). From a modelling perspective, projections differ considerably in future sea-level contributions depending on the model used and process parameterisations 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, thousand years before present, ka BP). Thus, one way to gain insight in this respect is to assess its response 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 intriguing 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 which can last from several centuries to thousands of years (e.g. Dansgaard et al. (1993)). 20 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 reorganisations 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) 25 (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, 2013Marcott et al., 2013).
As 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 30 during Heinrich stadials, with cooling during interstadials. The amplitude of this signal can reach up to 2 K (Augustin et 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 behaviour 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 35 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., 2013Andrews and Voelker, 2018;Hemming, 2004). In the Southern Hemisphere (SH), data showing IRD deposition from 5 the AIS is more scarce. There is evidence of ice discharges from the AIS (Kim et al., 2018;Weber et al., 2012Weber et al., , 2014, but neither a quantification of their contribution in terms of its sea-level equivalent (SLE) nor the identification of their triggering mechanism have 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 to 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. 10 Finally, sea-level reconstructions show fast variations of more than 20 m at millennial time-scales 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 ca.
14.5 ka BP (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;15 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 time-scales and their potential implications. The aim of this paper is thus to investigate the response of the AIS to millennial-scale variability 20 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 time-scales 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 that is deduced from the Dome-C deuterium ice-core record. To study the impact of ice-ocean interactions we use a basal 25 melting parameterisation 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 (Section 2). Then the response of the AIS to the oceanic forcing is shown, focusing on the ice discharges and grounding line advances at millennial time-scales (Section 3). Finally the main results are discussed (Section 4) and conclusions summarized (Section 5).

Model
We use the three-dimensional, hybrid, thermo-mechanical, 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. (2018); Tabone et al. (2018)). Important changes with respect to the original code include variations in boundary conditions 5 (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 10 floating ice shelves, whose evolution is determined by stretching processes, are solved by the Shallow Shelf Approximation (SSA) (Hutter, 1983;MacAyeal, 1989). Intermediate states where shearing and stretching regimes can appear simultaneously are proper 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. 15 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 20 (details in Section 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).

Forcing method and experimental design
GRISLI-UCM is forced through the same parameterisation for atmospheric and oceanic forcing as in Banderas et al. (2018) and Tabone et al. (2018), who used it to investigate specifically the past evolution of the glacial NH and Greenland ice sheets, 25 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 alpha=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 ka and 30 19 ka to avoid both orbital and submillennial-scale variability. The time-evolution of atmospheric temperatures (T atm (t)) and precipitation (P (t)) fields are given by the following equations: LGM + β(t)∆T atm mil (1) where temperature and precipitation, T atm LGM 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 interme-5 diate complexity CLIMBER 3-α (Montoya and Levermann, 2008). The millennial (∆T atm mil , δ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): 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, 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 15 grounding line (B gl ) and below the ice shelf (B shlf ): 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, hence we set γ to 0.1. Following the same procedure as for the atmospheric forcing, the oceanic temperature can be rewritten as: where B LGM represents LGM melting rates and ∆T ocn mil the millennial oceanic temperature anomaly. To avoid any accretion at the ice-shelf base, B gl cannot become lower than 0 m a −1 .
To study the response of the AIS to millennial-scale variability alone, we spun up our model for 120 ka under fixed LGM conditions. Fig. 1 illustrates the surface elevation and velocities after the spin-up procedure. We then impose the millennial-25 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 −1 (see Fig. 2a) allows for such an advance. In regions with ocean depths below 2000m, an artificially large 30 melting rate (50 m a −1 ) 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 of 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, thus ∆T ocn orb =-2.5 K 5 and ∆T ocn mil =0.5 K. Oceanic temperature variations are applied uniformly in space. Fig. 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 −1 K −1 ), to a large sensitivity (15 m a −1 K −1 ). All values of the tested parameters are provided in Table 1. Finally, sea-level variations are prescribed from Rohling et al. (2014).

10
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 −1 K −1 ; 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). 15 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.
To quantify the grounding line migration we introduce a parameter called Marine Zone Occupation (MZO) which is defined 20 as: where N G is the number of model grid cells with grounded ice in marine zones (i.e., zones where 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 where the ice is not grounded but floating and where the underlying 25 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 30 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 Georges 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 modelled AIS reacts at millennial timescales. after a typical cold phase (at 61 ka BP). The configuration in the three cases is similar, with an advanced grounding line with 5 grounded Ronne and Ross embayments. The grounding-line retreat in the Ross shelf increases with increasing κ. Ice streams also penetrate further inland with increasing κ. Fig. 5 illustrates the same fields after an AIM event (at 57 ka BP). While the lowest sensitivity case (κ =1 m a −1 K −1 ) shows an extensive ice sheet close to the continental-shelf break similar to the initial LGM state, as the sensitivity increases (κ =5 m a −1 K −1 ), marine zones such as the Ronne ice shelf begin to retreat and velocities increase. For sufficiently high oceanic sensitivities (κ =10 m a −1 K −1 ) the Ronne and Ross ice shelves experience 10 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 where the AIS is discharging for increasing ice-ocean sensitivities (see Fig. 6). The majority of the ice loss comes from the Ronne shelf. 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 −1 K −1 . The Pine Island zone responds in a similar manner to the oceanic warming but with less impact.

15
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.
The longest ice regrowth periods, corresponding to cooling phases, happens between 70 and 60 ka BP and between 40 and 20 ka BP. During these periods, for medium-low sensitivities (up to κ =7 m a −1 K −1 ), the grounding-line position (as indicated by the MZO) advances close to its LGM value, whereas for high oceanic sensitivities the maximum MZO value 20 reached decreases with increasing κ, indicating irreversibility, typical of hysteresis behaviour (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 80-15 ka BP (Fig. 7a). During this time period we find six significant ice discharge events in response to enhanced submarine melting phases, marked with grey shaded colors.
25 Fig. 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 −1 K −1 ) 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 the sea-level contribution of every event increases but also a wider spread is found between the discharging events, meaning that the sealevel difference between the smallest and largest ice discharge increases. Finally, what determines the total amount of sea-level 30 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 to happen but if the time interval is too short, then this will not be necessarily 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).
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 reorganisations. However, the extent to which the SO temperature increases during the slowdown of the AMOC is under debate. dePedro 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 5 provided by the Southern subtropical Atlantic and transferred to the AIS by the atmosphere; in addition, oceanic heat transport changes could be compensated 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. (2013Martin et al. ( , 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). 10 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 temperatures 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 15 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 which are yet not available. 20 We also found that if the heat-flux transfer parameter between ice and ocean is strictly larger than 10 m a −1 K −1 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 −1 K −1 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 25 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 30 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 which can reach more than 20m of SLE. These sealevel 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 of SLE. Low sensitivities (κ < 5 m a −1 K −1 ) 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 −1 K −1 SLE contributions of more than 10 m occur, which would imply a significant Antarctic contribution as well. This could have an impact on the size of paleo ice sheets, as well as the origin of ice storage in sea-level 5 reconstructions.

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 through a forcing method, which has already been tested  and is provided by an improved subglacial melting 10 routine. Because SO temperature reconstructions are not available we assumed 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 discharges, leading to substantial sea-level rise and grounding-line migration. Although 15 this timescale may seem short for such a large ice sheet, some extreme cases showed considerable grounding line retreat in the Ronne and Ross embayments and sea-level discharge of more than 10 m at millennial timescales. 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.
Code and data availability 20 GRISLI-UCM code and the analysed data are available from the authors upon request.
Author contributions. JB carried out the simulations, analysed the results and wrote the paper. All other authors contributed to designing the simulations, analysing the results and writing the paper. Heat flux coefficient κ m a −1 K −1 0,1,3,5,7,10,15