Ocean cavity regime shift reversed West Antarctic grounding line retreat in the late Holocene

Recent geologic and modeled evidence suggests that the grounding line of the Siple Coast of the West Antarctic Ice Sheet (WAIS) retreated hundreds of kilometers beyond its present position in the middle to late Holocene and readvanced within the past 1.7 ka. This grounding line reversal has been attributed to both changing rates of isostatic rebound and regional climate change. Here, we test these two hypotheses using a proxy-informed ensemble of ice sheet model simulations with varying ocean thermal forcing, global glacioisostatic adjustment (GIA) model simulations, and coupled ice sheet-GIA simulations that consider the interactions between these processes. Our results indicate that a warm to cold ocean cavity regime shift is the most likely cause of this grounding line reversal, but that GIA influences the rate of ice sheet response to oceanic changes. This implies that the grounding line here is sensitive to future changes in sub-ice shelf ocean circulation.


Ice sheet modelling with simple viscoelastic Earth deformation
The environmental forcings we use to force the model in the deglacial experiments are shown in Fig S1 (20 to 0 ka BP).The global mean sea level forcing is based on a compilation of global sea level proxy records from the last deglaciation (1)(2)(3)(4).Surface temperature (°C) and precipitation (%) anomalies are derived from the West Antarctic Ice Sheet (WAIS) Divide ice core reconstruction (5), and applied as anomalies to the modern climatology (6).The ocean temperature (°C) anomalies are derived from the TraCE-21ka deglacial climate model simulation (7,8), and are applied to a spatial field of ocean temperatures optimised for modern basal ice shelf melt rates (9).

Our model benchmarking focuses on the deglacial ice sheet evolution (Fig S2).
Parameter values for our reference simulations are shown in Table S1.Last Glacial Maximum ice volume of our reference simulations ranges from 10.9 to 11.5 m sea level equivalent (s.l.e.) above the present-day ice sheet (Fig S2a), within the range of estimates of various studies, i.e. 7.3 to 13.6 m s.l.e.(10)(11)(12).Deglacial ice thickness changes are within the ranges indicated by cosmogenic geochronology in the Transantarctic Mountains (13)(14)(15)(16), with rapid deglacial ice thinning occurring in the models from the early to middle Holocene (Fig S2b,c,d).
We also run a present-day reference simulation, in which we map ice temperature fields, internal velocities and bed conditions from our model spin-up to the present-day ice sheet configuration (17), and run forward in time using historical climate forcing from NorESM1-M (18).This produces reasonable fit to observations of grounding line position, ice thickness (Fig S3a,b) and ice surface velocity (Fig S3c,d).Bias with respect to ice surface velocity does occur at the Siple Coast of Antarctica, with lower-than-observed surface velocity of West Antarctic ice streams and regions of higher-than-observed ice shelf velocity.Modern basal melt rates are within the ranges estimated by Adusumilli et al. (2020) (19), with modelled melt rates as high as 18 m yr -1 in the Amundsen Sea sector ice shelves, but generally < 1 m yr -1 for the Ross Ice shelf, which overlies a cold ocean cavity (Fig S3e ).We extend this run for 500 years to demonstrate the ice sheet response to present-day climate over a longer duration.Sea level equivalent ice volume shows a small decrease (FigS3f), as expected given modern Antarctic mass loss, but grounding line position remains relatively unchanged at the Siple Coast (gold line in Fig S3b,d).Using the standard ocean forcing from TraCE-21ka (i.e.no mid-Holocene warming), we explore the parameter space of three solid earth parameters in the simple two-layer viscoelastic Earth deformation model from Bueler et al., (2007) (20): mantle viscosity, lithosphere flexural rigidity and mantle density.Experiments are run with 8 values of mantle viscosity: 1e19, 5e19, 1e20, 5e20, 7.5e20, 1e21, 2.5e21, and 5e21 Pa s; 5 values of lithosphere flexural rigidity: 1e23, 1e24, 5e24, 1e25, and 5e25 N m; and 2 values of mantle density: 3300 and 4500 kg m -3 .Mantle viscosity is an important control on the rate of isostatic rebound, with grounding-line retreat and readvance only occurring within the range of 5e20 to 1e21 Pa s (Fig S4).The relationship between grounding-line migration and lithosphere flexural rigidity is non-linear, with especially high and low values limiting grounding-line retreat.At intermediate values, in some cases (e.g.5e24 N m in Fig S4), pinning points on the outer continental shelf limit collapse of the outer ice shelf as well as limiting grounding line retreat and readvance.Of the two mantle density values tested, the higher value delays and limits grounding line retreat and limits readvance (Fig S5).Because large uncertainties exists with respect to the timing of WAIS retreat and advance of the Siple Coast grounding line, and because the radiocarbon input and decay modelling indicates differences in ice stream behaviour (21), we run additional ocean forcing experiments in which we vary the timing of onset and removal of anomalous ocean warming (Fig S7).The timing of retreat and advance of WIS and BIS corresponds to these forcings.Under anomalous ocean warming of +0.6°C from the standard TraCE-21ka forcing, the simulated ice shelf collapses between 6 and 4 ka BP in our simulations.We also explore the impact of model resolution by running an experiment at 10 km spatial resolution (as compared to our standard 20 km resolution experiments).Because increasing model resolution increases ice sheet and ice shelf velocity and produces thinner ice shelves, we ran this experiment using lower enhancement factors for the Shallow Ice and Shallow Shelf Approximations (i.e.1.5 and 0.4, respectively).The grounding line retreats and advances in response to the anomalous ocean forcing application and removal, demonstrating that the grounding-line reversibility in response to the change in ocean forcing is robust with different model resolutions.

Fig S1 .Surface pdd standard deviation value 2 K
Fig S1.Environmental forcings applied in the ice sheet model simulations.(a) Global mean sea level forcing (m), (b) surface temperature anomaly (°C), (c) snow accumulation (m yr -1 , but applied as a % anomaly in the ice sheet model), (d) ocean temperature anomaly in the Ross Sea region.The standard TraCE-21ka ocean forcing (black line) is the anomaly at ~500m depth relative to 0 ka, and examples of modified ocean forcings (+0.3, +0.4 and +0.5°C) are shown in color (purple, blue and orange, respectively).

Fig
Fig S3.Present-day ice sheet model reference simulation.This reference simulation was run for the period 1950 to 2450 using the model parameters listed in Table S1 and climate forcing derived from NorESM1-M.Following 2015, the climate forcing is held constant using a mean present-day climate.Shown are a) modelled ice thickness at 2015 and (b) the differences to Morlighem et al. (2020) (17); c) modelled ice surface velocity at 2015 and d) difference to Mouginot et al. (2017) (22); e) Modelled basal melt at 2015.The black lines indicate the modelled ice sheet grounding and calving line positions at 2015, the dark gold lines indicate the modelled grounding line position at 2450, and the purple lines in the difference plots indicate the observed grounding and calving line positions.f) Transient ice volume change (sea level equivalent in m) relative to 1950 of the reference simulation.

Fig S4 .
Fig S4.Predicted Holocene grounding line migration with varying lithosphere flexural rigidity and mantle viscosity.Bed topography (m above sea level) at 0 ka BP, and grounding line position every 1000 years from 9 ka BP (white) to 0ka BP (black), darkening 110 in scale, for ice sheet simulations with varying mantle viscosity (MV; units of Pa s) and lithosphere flexural rigidity (LFR; units of N m): (a-c) MV of 1e20 Pa s, (d-f) MV of 5e20 Pa s, (g-i) MV of 1e21 Pa s, and (j-l) MV of 2.5e21, with LFR of 1e24, 5e24, and 1e25, respectively.The orange line indicates the glacial grounding line position at 18 ka BP.Magenta circles indicate the West Antarctic subglacial sediment sites.Stars indicate that the 115 grounding retreats and readvances in the simulation.

Fig S5 .
Fig S5.Predicted Holocene grounding line migration with varying mantle density, lithosphere flexural rigidity and mantle viscosity.Bed topography (m above sea level) at 0 ka BP, and grounding line position every 1000 years from 9 ka BP (white) to 0 ka BP (black), darkening in scale, for ice sheet simulations with varying mantle viscosity (MV; units of Pa s), mantle density (MD; units of kg m -3 ) and lithosphere flexural rigidity (LFR; units of N m): (a) MV of 5e20 Pa s, MD of 3300, LFR of 1e24; (b) MV of 5e20 Pa s, MD of 4500, LFR of 1e24; (c) MV of 5e20 Pa s, MD of 4500, LFR of 1e25; (d) MV of 7.5e20 Pa s, MD of 3300, LFR of 1e24; (e) MV of 7.5e20 Pa s, MD of 4500, LFR of 1e24; (f) MV of 7.5e20 Pa s, MD of 4500, LFR of 1e25; (g) MV of 1e21 Pa s, MD of 3300, LFR of 1e24; (h) MV of 1e21 Pa s, MD of 4500, LFR of 1e24; (i) MV of 1e21 Pa s, MD of 4500, LFR of 1e25.The orange line indicates the glacial grounding line position at 18 ka BP.Magenta circles indicate the West Antarctic subglacial sediment sites.Stars indicate that the grounding retreats and readvances in the simulation.

Fig S6 .
Fig S6.Modelled grounding line position with alternative forcings and resolution.(a) Modified TraCE-21ka ocean forcings for the Ross Sea sector, with colors corresponding to the timing and magnitude of the anomalous warming.Distance of ice sheet grounding from the modern grounding line for the (b) Whillans and (c) Bindschadler ice streams for simulations with a mantle viscosity of 7.5e20 Pa s.The dotted orange line is for a simulation using higher spatial resolution and reduced SIA and SSA enhancement factors.Horizontal lines correspond to the proxy sites.The vertical blue line shows the radiocarbon age for retreat at Whillans Grounding Zone (WGZ) from Venturelli et al. (2020), with blue shading indicating age uncertainty.The vertical yellow and pink lines respectively indicate the modelled age for retreat and advance at the WIS and BIS sites from Neuhaus et al. (2021) with colored shading indicating age uncertainty.The asterisk for the +0.6°C simulation indicates that the ice shelf collapses in the mid-Holocene.

Fig S7 .
Fig S7.Late Holocene bed topography differences predicted with varying mantle viscosity.Left column (a, d, g): Bed topography (m above sea level) of the MV5e20 OT+0.5 simulation at 4, 2 and 0 ka BP, respectively.Middle column (b, e, h): Relative bed topography (m) of the MV7.5e20 simulation to the MV5e20 OT+0.5 simulation at 4, 2 and 0 150 ka BP, respectively.Right column (c, f, i): Relative bed topography (m) of the MV1e21 simulation to the MV5e20 OT+0.5 simulation at 4, 2 and 0 ka BP, respectively.The grounding line position is indicated by the black line.Magenta circles indicate the West Antarctic subglacial sediment sites.

Fig S9 .
Fig S9.Holocene bed topography changes predicted by the ice sheet and GIA model 195 with weak Earth structure.Left column: Bed topography changes (m) of the MV7.5e20OT+0.5 ice sheet model simulation (ISM) from (a) 8 to 6 ka BP, (d) 6 to 4 ka BP, (g) 4 to 2 ka BP, and (j) 2 to 0 ka BP.Middle column: Bed topography changes (m) of the weak earth GIA model (GIAWE) using ice thickness from the MV7.5e20OT+0.5 simulation from (b) 8 to 6 ka BP, (e) 6 to 4 ka BP, (h) 4 to 2 ka BP, and (k) 2 to 0 ka BP.Right column: Difference 200 between the bed topography changes of the ISM and the GIAWE from (c) 8 to 6 ka BP, (f) 6 to 4 ka BP, (i) 4 to 2 ka BP, and (l) 2 to 0 ka BP.The gray line indicates the grounding line position at the beginning of each time period, whereas the black line indicates the position at the end.

Fig S10 .
Fig S10.Holocene bed topography changes predicted by the ice sheet and GIA model with strong Earth structure.Left column: Bed topography changes (m) of the MV7.5e20OT+0.5 ice sheet model simulation (ISM) from (a) 8 to 6 ka BP, (d) 6 to 4 ka BP, (g) 4 to 2 ka BP, and (j) 2 to 0 ka BP.Middle column: Bed topography changes (m) of the strong earth GIA model (GIASE) using ice thickness from the MV7.5e20OT+0.5 simulation from (b) 8 to 6 ka BP, (e) 6 to 4 ka BP, (h) 4 to 2 ka BP, and (k) 2 to 0 ka BP.Right column: Difference between the bed topography changes of the ISM and the GIASE from (c) 8 to 6 ka BP, (f) 6 to 4 ka BP, (i) 4 to 2 ka BP, and (l) 2 to 0 ka BP.The gray line indicates the grounding line position at the beginning of each time period, whereas the black line indicates the position at the end.