Sensitivity of the Antarctic ice sheets to the warming of marine isotope substage 11c

Studying the response of the Antarctic ice sheets during periods when climate conditions were similar to the present can provide important insights into current observed changes and help identify natural drivers of ice sheet retreat. In this context, the marine isotope substage 11c (MIS11c) interglacial offers a suitable scenario, given that during its later portion orbital parameters were close to our current interglacial. Ice core data indicate that warmer-than-present temperatures lasted for longer than during other interglacials. However, the response of the Antarctic ice sheets and their contribution to sea level rise remain unclear. We explore the dynamics of the Antarctic ice sheets during this period using a numerical ice sheet model forced by MIS11c climate conditions derived from climate model outputs scaled by three glaciological and one sedimentary proxy records of ice volume. Our results indicate that the East and West Antarctic ice sheets contributed 4.0–8.2 m to the MIS11c sea level rise. In the case of a West Antarctic Ice Sheet collapse, which is the most probable scenario according to far-field sea level reconstructions, the range is reduced to 6.7–8.2 m independently of the choices of external sea level forcing and millennialscale climate variability. Within this latter range, the main source of uncertainty arises from the sensitivity of the East Antarctic Ice Sheet to a choice of initial ice sheet configuration. We found that the warmer regional climate signal captured by Antarctic ice cores during peak MIS11c is crucial to reproduce the contribution expected from Antarctica during the recorded global sea level highstand. This climate signal translates to a modest threshold of 0.4 C oceanic warming at intermediate depths, which leads to a collapse of the West Antarctic Ice Sheet if sustained for at least 4000 years.

As described in Sect. 2.1 (main text), melting rates at the base of ice shelves are computed through a parameterisation based on Beckmann and Goosse (2003), following the form presented in Martin et al. (2011), but with a quadratic dependence on thermal forcing (Holland 5 et al., 2008): where p 0 is the specific heat capacity of the ocean mixed layer (3974 J kg −1 K −1 ), the thermal change velocity T is 10 −4 m s −1 (as in Beckmann and Goosse, 2003;Martin et al., 2011), i is the density of ice (910 kg m −3 ), and i is the latent heat of fusion (3.35×10 5 J kg −1 ) of ice. Here, forc is given by where zb is the ocean temperature (obtained from the forcing data set) at the depth of the ice shelf base and melt is calculated from Eq. 2 10 of Beckmann and Goosse (2003): where b and sl are the modeled depths of the ice shelf base and the sea level respectively, and o is the average sea water salinity, set as a constant value of 35 since initial sensitivity tests using spatially variable salinity showed a negligibly small effect on the parameterised basal melting rates.
We optimize the choice of the melting coefficient melt by first iteratively adjusting the spatial distribution of basal melting rates, bm , 15 over all Antarctic ice shelves (i.e., at every ice shelf grid point) such that it keeps the thickness of each grid point as close as possible to the Bedmap2 data (cf. Bernales et al., 2017). We then derive spatially varying values for melt through Eq. 4: On the one hand, and because model simulations over long time scales need to consider the evolution of ice shelf distribution (including changes in the position of the grounding line), a spatially heterogeneous distribution of melt is impractical for the experiments presented in this study. On the other hand, the use of a single value for melt is not able to reproduce the stark difference in magnitude between the 20 melting rates at grounding lines and ice-shelf interiors (cf. Rignot and Jacobs, 2002). We thus redefine the value of melt as where SH is a value chosen for ice-shelf interiors (away from the grounding line) and melt GL is the value at grounding lines. These 1 Supplementary Material S1 Calibration of the ice shelf basal melting parameterisation two values are weighed by fg , which varies according to the flux of ice across the grounding line, given by where Φ ref = 2 × 10 5 m 2 a −1 and Φ gl is the flux across the grounding line also in m 2 a −1 . The value of Φ ref is simply hand-tuned so that the method identifies the majority of modern flux gates (according to Rignot and Jacobs, 2002;Rignot et al., 2013). 25 We then find minimum, average, and maximum fg and melt GL  these pairs.
All combinations tested yielded 40 steady-state simulations of 10 kyr under CCSM3 PI forcing. We chose the best pair of melt SH and melt GL values based on the run that revealed the least changes compared to the original PI configuration. We find the combination of values that met this criteria to be 0.5 K −1 for melt GL and 8.5 × 10 −3 K −1 for melt SH . It is important to stress that the former is very high as it is a hypothetical limit which is never attained due to the weighting applied by fg during the simulations (see Eq. 6). However, we 35 do identify three regions where the model was not able to sustain a grounding line that matches Present Day (PD) and PI configurations: the southeast part of the Filchner-Ronne Ice Shelf (related to Support Force, Foundation, and Academy ice streams), the Scott Glacier discharge area on the Ross Ice Shelf, and in the strait between George VI Island and the Antarctic Peninsula. The difficulty in accurately representing these sites is not exclusive to our model and has been reported before (e.g., Golledge et al., 2012). Part of this could also be due to inaccuracies in the used topography dataset (Pritchard, 2014;Kingslake et al., 2018). Higher-resolution runs and sub-grid 40 interpolation schemes at the grounding lines could potentially help to solve this problem, but would increase the computational costs of such long runs. These caveats are taken into account when interpreting the results. While this simplification is less realistic than having the melting rate distribution under the entire shelf, it allows to incorporate the calibration into a freely evolving run without having to adjust the basal melting parameters point by point, especially since ice shelf geometries change over a transient simulation. Figure S1a-c shows the resulting PI Surface Mass Balance, the distribution of basal melt rates post-optimization, and the PI topography from which we 45 do the model spinup, respectively.
The CCSM3 model results used seem to represent surface air temperature satisfactorily given its spatial resolution (T31 spectral resolution), capturing the pattern of colder surfaces in the highest regions in the ice sheet interior, and milder temperatures at its margins and over the Antarctic Peninsula (Figs. S2a,b). Regarding the ocean temperatures, since all datasets for PI or PD of continent-wide sub-shelf ocean 50 temperatures that could be considered as our reference are also highly subject to statistical and model biases, we chose to provide an assessment of the modelled ice-shelf basal melting against observation-inferred values instead (Fig. S2d). As reference we take the basal melt rates of Rignot et al. (2013). By comparing the sub ice-shelf basal melt rates (Figs. S2c,d), we can see that the parameterization S2 Brief assessment of CCSM3 forcings Figure S1: Ice sheet calibration using CCSM3 forcings. (a) Imposed Surface Mass Balance (m a −1 ) from CCSM3 over the ice sheet during the calibration period. (b) Optimized basal melting rates ( in m a −1 ) for PI, and (c) resulting PI topography (in meters) and geometry under a steady state from our calibration.
captures successfully the higher melting rates at the grounding lines of ice streams, but smooths the uneven pattern of melting under more interior regions of the ice shelves, and has little basal melting at the calving front (Figs. S2c,d). One limitation of this method is its inability 55 to represent sub-shelf refreezing. A different parameterization that accounts for this effect has already been developed (Lazeroms et al., 2018), but its implementation for use at paleo timescales has not been done yet. Nevertheless, the absence of this process does not seem to impair our results. As for precipitation (Figs. S2e,f), the model manages to capture the overall pattern, i.e., increased precipitation at the fringes of the ice sheet and less precipitation in the interior. It does not reproduce, however, higher precipitation induced by orography (e.g., west of the Antarctic Peninsula and north of the Transantarctic Mountains), mainly due to the model resolution. Finally, despite 60 the precipitation being significantly lower over the ice sheet interior compared to its margins as expected, the values are significantly overestimated over the interior and underestimated over the high precipitation areas, similar to what has been observed for more modern CMIP5 models (Palerme et al., 2017).
Initial simulations were carried using both CCSM3 and CESM1.2 as climate forcings. CCSM3 was chosen rather than CESM1.2 for two 65 important reasons: 1. With CESM1.2 forcing (in-house simulations, see Bakker et al., 2020), the ice-sheet model failed to match the geological constraints for MIS11c by Raymo and Mitrovica (2012), i.e., not showing a volume loss that would cause the expected contribution from the AIS for this period (Fig. S3, "CESM"). In order to understand what exactly caused this difference in performance between the two climate models, we ran a series of sensitivity experiments where we replaced one forcing field (air temperature, ocean temperature, or precipitation; 70 Fig. S3) or two forcing fields (e.g. air temperature + ocean temperatures; not shown) from CCSM3 by their CESM1.2 equivalent. We found that the evolution of the ice sheet throughout MIS11c was most sensitive to the differences in precipitation. The differences in precipitation fields show that CESM1.2 precipitation rates are, in general, lower than those of CCSM3, especially in key areas of the WAIS, such as East of Siple Coast. As a consequence, the calibration of the ice-sheet model to the CESM1.2 forcing fields resulted in the need for a higher basal friction (relative to CCSM3) to compensate for the reduced precipitation in order to match the modern reference 75 observational data sets. In turn, the combination of this higher basal friction and the higher precipitation rates during MIS11c results in 3 S3 Comparison of CCSM3 with other datasets for MIS11c a much reduced sensitivity of the ice sheet to MIS11c warming, thus not capturing the AIS contribution to sea level rise shown by the geological record (Raymo and Mitrovica, 2012).

Several studies have shown that CCSM3 satisfactorily simulates Southern Ocean conditions during glacials and interglacials, which
is important for the simulation of the AIS. It has been shown that CCSM3 correctly simulates characteristics of water masses produced 80 in the Southern Ocean (AABW, AAIW) for the LGM (Otto-Bliesner and Brady, 2008;Ronge et al., 2015) and the transition into the Holocene (Ronge et al., 2020). Moreover, Marzocchi and Jansen (2017) have shown that CCSM3 has a better skill in simulating glacial Antarctic sea ice and deep-ocean circulation than CCSM4. Figure S4 shows how CCSM3 ΔT (LGM minus PI) compare to those derived from ice cores (Werner et al., 2018). The CCSM3 results suggest a stronger cooling than the proxy data, which is likely related to a too thick Antarctic ice sheet prescribed in the LGM (PMIP2) 85 simulations, and hence does not substantially affect our ice-sheet forcing due to lapse-rate correction.
To assess how our forcing's biases impact our results, we performed an EDC-forced simulation using RACMO2 (Van Wessem et al., 2014) as our modern reference, with CCSM3 as the anomalies and our ocean (Fig. S5). We also include a similar experiment using RACMO2 modern fields and CESM1.2 anomalies and ocean forcing. In Fig. S5 we show how they differ from the simulations fully forced Figure S3: (a) Grounding lines (dashed for an easier comparison) at 405 ka (i.e., the MIS11c sea level highstand), and (b) ice volume [10 6 km 3 ] throughout MIS11 for the CESM1.2 and CCSM3 runs. 'airt','ocnt' and 'prec' denote atmospheric surface temperature, ocean temperatures, and precipitation rates respectively, and these refer to the runs where CCSM3 forcing variables were replaced by the above referred variables from CESM1.2 (e.g., 'CCSM3 + CESM airt' means that all variables were from CCSM3, except for atmospheric surface temperatures, which were from CESM1.2). with CCSM3 (which are presented in the manuscript's main text) and CESM1.2 (which we presented in Fig. S3). Differences do exist 90 between RACMO2 and CCSM3, but are relatively small and most evident during the time of sea level highstand at 405 ka. This difference in ice volume at 405 ka means that the RACMO2 runs contribute 1.1 m less to global mean sea level. The position of the grounding line shows that this difference seems to relate mainly to the fringes of the EAIS (particularly in Dronning Maud Land) and to the WAIS sector just south of the Peninsula. The runs in which CESM1.2 anomalies were applied yield a rather insensitive ice sheet. As discussed above, and considering that basal sliding conditions in these simulations were calibrated to RACMO2 climate, the 'RACMO+CESM' experiment 95 highlights the fact that its precipitation anomalies are what makes it not capture the sea level contribution constraints for MIS11c. Figure S5: (a) grounding lines (dashed for an easier comparison) at 405 ka (i.e., the MIS11c sea level highstand), and (b) Ice volume [10 6 km 3 ] throughout MIS11c for the bias assessment runs for CCSM3 and CESM1.2, where we test for the impact of using RACMO2 atmospheric fields as the PI reference, while keeping the anomalies from CCSM3 and CESM1.2.
We performed an ensemble forced by the EDC GI with three simulations: Two where we add a ΔT of -0.5 and -1.0, • C to the ocean, and one where we reduce the ocean GI amplitude by 50%. We include the original ensemble member for comparison. Figures S6 and S7 show respectively the resulting ice volume and grounding lines at times of interest. In short, the tipping point at ca. 412 ka still persists despite 100 the fact that the ocean is substantially colder, but a total collapse is not achieved. Considering that (i) there is no connection between the Ross and Weddell seas at 405 ka for the ΔT ocn = −0.5 run (Fig. S7d), and (ii) the passage is open between the Ross and Amundsen seas, a comparison of their ocean temperatures at intermediate depths (Fig. S8) further strengthens the fact that a 0.4 • C warming at intermediate depths (which we define in the main text to be between 400 and 1000 m) is the threshold for which WAIS collapse is possible. Figure S9 shows the applied thermal forcing under the ice shelves for the "EDC control" run shown in Fig. S7.

105
As mentioned in Sect. 2 of the main text, ice-shelf calving in our model is done by a simple thickness threshold of 50 m, since there are no observed ice shelves thinner than this threshold. This particular value allows the growth of new ice shelves, but avoiding growth beyond reasonable extents for present-day simulations. As a consequence, whenever the ocean shows a strong warming and ice shelves lose area, calving fluxes are also reduced. Similarly, fluxes increase again as the ice-shelf areas grow (Fig. S10). Finally, we have also tested for different lags applied to the ocean response to the atmospheric forcing (Fig. S11), which did not change our results.
110 Figure S6: Sensitivity of the AIS response expressed in total ice volume [10 6 km 3 ] to three simulations where we test for the ocean sensitivity for a collapse. In two runs we apply a ΔT of -0.5 and -1.0 • C, and in a third the ocean GI amplitude was reduced by 50%.
8 S4 Sensitivity of the WAIS collapse to ocean forcing   Figure S9: Spatial distribution of Thermal forcing (i.e., the difference between the ocean temperature and the ice shelf base temperature) at time steps of interest for each of the CFEN members. Figure S10: Calving fluxes (in Gt/a) integrated along the Filchner-Ronne and Ross ice shelves calving fronts. Figure S11: Sensitivity of the AIS response expressed in total ice volume [10 6 km 3 ] to different lags in the GI applied to the ocean between 420 and 394 ka. Figure S12: Relaxation period between 425 and 420 ka for all four CFEN members (cf. Sect. 2.2 in the main text). (a) shows the grounding line at 425 ka (solid line) and at 420 ka for each member (dashed lines); (b) shows the evolution of total ice volume [10 6 km 3 ] during this 5 kyr period for each member. Figure S13: (a-f) Forcing fields used to construct the climate forcings used in this study. Left fields show PI mean states, while right fields show applied anomalies (i.e., LGM − PI for temperatures, and LGM ÷ PI for precipitation). Anomalies are multiplied by the GI and then combined with the left-side fields as in Eqs. (4) and (5).

Supplementary figures
Figure S14: Comparison between surface atmospheric temperature anomalies (ΔT, • C) obtained by the ice cores GI scaling and those inferred from their respective D values. Figure S15: Comparison between EDC GI-forced runs at different horizontal resolutions: 20 km (orange, presented in the main text), 16 km (blue), 15 km (green), 12 km (red), and 10 km (black). (a) grounding lines at 405 ka; (b) total grounded ice area (in 10 6 km 2 ), and (c) sea level contribution (in m s.l.e.) for the entire AIS.