Strong Ocean Melting Feedback During the Recent Retreat of Thwaites Glacier

Accelerating ice loss from Thwaites Glacier is contributing approximately 5% of global sea‐level rise, and could add tens of centimeters to sea level over the coming centuries. We use an ocean model to calculate sub‐ice melting for a succession of Digital Elevation Models of the main trunk of Thwaites Glacier from 2011 to 2022. The ice evolution during this period induces a strong geometrical feedback onto melting. Ice thinning and retreat provides a larger melting area, thicker and better‐connected sub‐ice water column, and steeper ice base. This leads to stronger sub‐ice ocean currents, increasing melting by over 30% without any change in forcing from wider ocean conditions. This geometrical feedback over just 12 years is comparable to melting changes arising from plausible century‐scale changes in ocean conditions and subglacial meltwater inflow. These findings imply that ocean‐driven ice loss from Thwaites Glacier may only be weakly influenced by anthropogenic emissions mitigation.


Introduction
The West Antarctic Ice Sheet (WAIS) is losing ice at an ever-increasing rate Shepherd et al., 2019) and forms a major uncertainty in projections of global sea-level rise (Fox-Kemper et al., 2021). The most rapid ice loss is occurring in the Amundsen Sea sector, where thinning and retreat of the floating ice shelves is causing acceleration of their tributary glaciers (De Rydt et al., 2021;dos Santos et al., 2021). Thwaites Glacier is one of the largest contributors to sea-level rise Shepherd et al., 2019), and offers the largest, and least-certain, potential future contribution (Alevropoulos-Borrill et al., 2020;Arthern & Williams, 2017;Joughin et al., 2014;Seroussi et al., 2020;Yu et al., 2018). Thwaites Glacier has accelerated , thinned (Konrad et al., 2017), and experienced grounding-line retreat (Milillo et al., 2019;Rignot et al., 2014) throughout the satellite era. The Thwaites sector lost 76 ± 6 Gt/y during 2012-2016, contributing ∼5% of global-mean sea-level rise (Fox-Kemper et al., 2021).
The floating section of Thwaites Glacier can be divided into Thwaites East Ice Shelf and Thwaites West Ice Tongue (TEIS and TWIT; Figure 1). TEIS is grounded on a prominent seabed ridge, buttressing ice flow, but TWIT has only been grounded on a tiny pinning point on this ridge during the satellite era ( Figure 1a; Hogan et al., 2020;Jordan et al., 2020;Rignot et al., 2014;Tinto & Bell, 2011). The ongoing ice loss is focused on the fast-flowing trunk of Thwaites Glacier, which flows into TWIT (Figure 1b; Konrad et al., 2017;Mouginot et al., 2014). TWIT was a coherent ice tongue that calved giant icebergs, but has disaggregated completely in recent decades (Figure 1d), with only a ∼10 km contiguous ice shelf remaining in a small embayment that we refer to as Thwaites Inlet (TI; Figure 1a; Ferrigno et al., 1993;MacGregor et al., 2012;Miles et al., 2020).
Abstract Accelerating ice loss from Thwaites Glacier is contributing approximately 5% of global sea-level rise, and could add tens of centimeters to sea level over the coming centuries. We use an ocean model to calculate sub-ice melting for a succession of Digital Elevation Models of the main trunk of Thwaites Glacier from 2011 to 2022. The ice evolution during this period induces a strong geometrical feedback onto melting. Ice thinning and retreat provides a larger melting area, thicker and better-connected sub-ice water column, and steeper ice base. This leads to stronger sub-ice ocean currents, increasing melting by over 30% without any change in forcing from wider ocean conditions. This geometrical feedback over just 12 years is comparable to melting changes arising from plausible century-scale changes in ocean conditions and subglacial meltwater inflow. These findings imply that ocean-driven ice loss from Thwaites Glacier may only be weakly influenced by anthropogenic emissions mitigation.

Plain Language Summary
The West Antarctic Ice Sheet is losing ice, making a substantial contribution to global sea-level rise. This ice loss is known to be triggered by changes in ocean melting of the floating parts of the ice sheet. Computer predictions show that this ice loss could make a large contribution to global sea-level over the coming centuries, but the future trajectory is very uncertain. In this study we simulated the ocean melting of Thwaites Glacier during 2011-2022, a period when the glacier rapidly thinned and retreated. We show that the geometrical evolution of the glacier during this period led to a substantial increase in ocean melting, caused by the exposure of more ice base to warm ocean waters, and changing ocean currents beneath the ice. This change in melting is similar to what might be expected from 100 years of ocean warming under anthropogenic climate change. These results imply that the future melting of such glaciers is strongly controlled by the geometrical evolution of the ice through internal ice and ocean feedbacks, and will therefore only weakly be influenced by reductions in the emissions of greenhouse gases.
These changes must be caused by the rapid oceanic ice-shelf melting in this region (Shepherd et al., 2004), but it is not clear how this melting has changed, and whether anthropogenic forcing played a role. The ice loss may be driven by climatic trends (either anthropogenic or natural) that led to an increase in melting during the 20th century (Holland et al., 2019Naughten et al., 2022). It is also possible that the changes were triggered historically by climate variability, and have since been perpetuated by feedbacks in the ice/ocean system (Holland  Figure S3 in Supporting Information S1. Axis tick marks are every 10 km.  Figure S2 in Supporting Information S1. et al., J. A. Smith, Andersen, et al., 2017;Steig et al., 2012). Understanding the relative importance of forcing and feedbacks is needed to quantify the future sea-level rise from Thwaites Glacier and its sensitivity to greenhouse-gas emissions.
In this study we investigate the ocean melting feedback during the recent retreat of Thwaites Glacier. We use ocean simulations to model the ice-shelf melting for Digital Elevation Models (DEMs) of Thwaites Glacier geometry each year from 2011 to 2022 (Bevan et al., 2021), holding other forcings constant. This set of simulations allows us to quantify the geometrical feedback onto melting during the retreat. Next, we determine the relative importance of this feedback by changing ocean forcings and subglacial meltwater in the model. This provides insight into how far ocean melting of Thwaites Glacier is controlled by internal feedbacks rather than external climatic forcing.

Methods
A general lack of observations beneath TWIT leads to many sources of uncertainty in our modeling, including seabed and ice geometries, ocean conditions, the influx of subglacial water from beneath the ice sheet, and in model parameterizations of turbulent mixing and ice-shelf melting. As a result, our model experiments should be viewed as hypothesis tests, not realistic model hindcasts.

Model Setup
We use a hydrostatic implementation of the MITgcm ocean model in a small 120 km × 120 km domain focused on the ocean cavity beneath Thwaites Glacier ( Figure 1a). Wider ocean forcing is applied through restoring boundary conditions for temperature and salinity on the open-ocean boundaries. No ocean surface fluxes or sea-ice model are applied. The simulations rapidly attain a steady state in which heat and salt sources from the boundaries are balanced by the heat sink and freshwater source from ice-shelf melting. Simulations are run for 6 months and results averaged over the final month. The equations are solved on a 200 m (horizontal) by 10 m (vertical) grid, with a timestep of 30 s. Ice-shelf melting and subglacial inflow are implemented as virtual heat and salt fluxes to avoid a rising sea surface. A standard three-equation ice-shelf melting parameterization is used, with parameter choices following Jenkins et al. (2010) (see Text S1 in Supporting Information S1). Constant viscosities (1 m 2 /s horizontal, 5 × 10 −4 m 2 /s vertical) and diffusivities (0.1 m 2 /s horizontal, 5 × 10 −5 m 2 /s vertical) are used.
Text S2 in Supporting Information S1 describes in detail how the model seabed and ice draft fields are generated, which is summarized here. A simulation is performed for each year between 2011 and 2022, and the different simulations are compared to reveal the effect of the changing ice geometry. All simulations use the seabed geometry from BedMachine version 3 (Morlighem et al., 2020; Figure 1a), but each simulation has a different steady ice-shelf draft, constructed as follows. First, ice-shelf surface elevation is taken from the Reference Elevation Model of Antarctica (REMA) (Howat et al., 2018). Next, a series of TanDEM-X DEMs (Bevan et al., 2021) are used to replace REMA over TWIT. We select one DEM from each year during 2011-2022, yielding 12 yearly surface elevation datasets that differ over TWIT only (e.g., Figures 1c and 1d). These surface elevations are converted to ice draft assuming floatation, with a correction near grounded ice. Each ice draft field defines a grounding line wherever the floatation ice base is below the seabed. The resulting grounding-line evolution agrees closely with satellite radar interferometry (Figure 1a; Milillo et al., 2019;Rignot et al., 2014).

Model Experiments and Forcing
We use the term "experiment" to refer to a set of 12 simulations, one for each year from 2011 to 2022, under a particular choice of forcings. By comparing the 12 simulations within an experiment we derive the influence of changing ice geometry. By comparing different experiments to each other we derive the influence of changing forcings. In geophysical experiments, we test the influence of warm and cool ocean forcing and a subglacial inflow. In diagnostic experiments, we test the influence of changes in the grounding line and water-column thickness. The different experiments are summarized in Table S1 in Supporting Information S1 and described fully below.
To quantify the influence of ocean forcings, we follow De  in conducting two illustrative experiments using extreme warm and cool ocean boundary conditions guided by observations taken on the Amundsen Sea shelf. The thickness of the warm Circumpolar Deep Water (CDW) layer in this region varies substantially in response to decadal climatic variability Jenkins et al., 2018). To the best of our knowledge, 2009 was the warmest year on record, and 2012 the coldest , despite the availability of more recent observations Dotto et al., 2022;Wåhlin et al., 2021). Ocean forcing profiles in WARM and COOL model experiments are shown in Figure S1 in Supporting Information S1. In the WARM experiment, representative of 2009, we prescribe a warm (+1.2°C) and saline (34.7 PSU) deep CDW layer below 600 m depth, with a colder (−1°C) and fresher (34 PSU) Winter Water layer above 200 m depth, and thermocline in between . This WARM experiment is used as a baseline against which other experiments are compared. In the COOL experiment, representative of 2012, the thermocline is shifted downwards by 200 m, to provide a thinner CDW layer and cooler water column .
To quantify the influence of subglacial meltwater inflow, we introduce a 150 m 3 /s flux into TI (Hager et al., 2022) with a salinity of 0 PSU and temperature of −0.5°C (approximately the freshwater freezing point at depth). This influx is uniformly distributed (vertically and horizontally) over all cells in a selected volume that is within 2 km of the grounding line and is ungrounded in all of the 12 yearly geometries (Figure 1a). This steady inflow represents the mean meltwater production at the base of the ice sheet, and does not account for subglacial lake filling or drainage. The subglacial inflow is enabled in all experiments except the NOSG experiment, which has WARM ocean forcing and no subglacial inflow.
We explain the geometrical feedback using experiments in which the 12 yearly simulations have parts of their geometrical change disabled. In the FIXGL experiment we modify the ice draft as before, but set the grounding line in all 12 simulations equal to the maximum grounded extent in any of the 12 yearly geometries. This is practically identical to fixing the grounding line at its 2011 position. In the FIXWC experiment we fix the grounding line in the same way, but also preserve the water-column thickness field at its 2011 values in all 12 simulations, by artificially changing the seabed to follow the changes in ice-shelf draft. Both FIXGL and FIXWC experiments use WARM ocean forcing. Figure 1e shows the ice-shelf melt rate in the 2022-geometry simulation from the WARM experiment, which totals 71.1 Gt/y. The melting is strongly focused upon the thickest, steepest ice near the western grounding line and in TI (Nakayama et al., 2019). There are no ocean observations in this area and estimates of melting are highly uncertain. Satellite-derived estimates of overall Thwaites melting are similar to the model, 60-100 Gt/y (Adusumilli et al., 2020;Depoorter et al., 2013;Liu et al., 2015;Rignot et al., 2013), but these assume that the ice is a continuum, which is not the case (Figure 1d). The model agrees with observations of weak melting beneath TEIS Dotto et al., 2022), but this provides little constraint over the stronger melting to the west. These higher melt rates are plausible: ∼100 m/y melting is consistent with ice thickness changing by ∼300 m over ∼10 km in ∼3 km/y ice flow (Figure 1). Therefore we proceed with the standard melting parameter choices of Jenkins et al. (2010), rather than tuning the model further (Text S1 in Supporting Information S1). The remainder of this paper focusses on the rapidly melting region marked on Figure 1e. Figure 2 shows the ice cavity geometry and modeled ocean behavior in the 2011 and 2022 simulations from the WARM experiment. The first two columns show the dramatic geometry change that occurred during this period (Bevan et al., 2021;Milillo et al., 2019). The ice thinned in TI and along the grounding line to the east, increasing the water-column thickness (Figure 2j), causing grounding-line retreat, and ungrounding an ice rumple in TI (Figure 2i). The ice shelf thickened offshore of these regions, presumably in response to the changed ocean melting and ice flow patterns.

Response of Melting to Geometrical Changes
Modeled melting in the region increases from 34.0 Gt/y to 45.1 Gt/y in response to this geometry change, with the greatest increases in the newly ungrounded areas and in TI (Figure 2k). The changes in melting are driven by the exposure of newly ungrounded areas to warm ocean waters, and changes in ocean currents immediately beneath the ice (Figure 2l). Both of these factors control melting through their effect on turbulent ocean heat fluxes (see Text S1 in Supporting Information S1). The overall increase in melting arises through a complex pattern of positive and negative melting changes, as spatial shifts in the currents move the rapidly-melting regions (Figure 2l).
Strong sub-ice currents flow westward along the steep ice slopes near the grounding line (Figure 2), driven by buoyant melt water under the influence of Coriolis force (Holland & Feltham, 2006). In 2011 this current seems to disappear within TI, leading to low melting (Figure 2d). A vertical section through this region ( Figure S2 in Supporting Information S1) shows that the ice slopes upward toward the grounding line in 2011, trapping a layer of cold, fresh water. This displaces the westward current downwards, away from the ice, leading to weak melting (Figure 2c). This situation is analogous to recent observations beneath TEIS Schmidt et al., 2023). By 2022 the ice slope is reversed ( Figure S2 in Supporting Information S1) and the sub-ice current remains attached to the ice base throughout the sector, dominating melting (Figure 2h).  Figure 2. The fastest increase in melting occurs in response to a geometrical change during the Austral summer 2012/13. Closer inspection of the geometries shows that this is when the ice rumple in TI ungrounded and the ice basal slope changed. It is not clear that this ungrounding had any particular driver, since it preceded subglacial lake drainage events that occurred later in 2013 (Malczyk et al., 2020;B. E. Smith, Gourmelen, et al., 2017); it may have simply been caused by the general ice thinning in this region.

Factors Contributing to the Geometrical Response
Figure 3a also shows the FIXGL and FIXWC experiments, which have a much smaller geometrical melting response. By comparing these to the WARM experiment, we may gain insight into the underlying causes of the geometrical melting response. Text S3 in Supporting Information S1 presents additional investigation into these causes, which we summarize here.
Overall melting increases by 11.1 Gt/y between 2011 and 2022 geometries in the WARM experiment, but in FIXGL this melting increases by only 2.9 Gt/y (Figure 3a). This implies that 8.2 Gt/y (∼80%) of the overall melting increase is caused by grounding-line retreat (which is present in WARM but absent in FIXGL). This is easily explained because grounding-line retreat exposes new ice-base area to ocean waters. These new ice-shelf areas typically host high melt rates, as they feature steeply sloping ice bathed in deep warm water.
In FIXWC the melting actually decreases by 1.2 Gt/y between 2011 and 2022 (Figure 3a). This means that all of the overall melting increase in WARM is explained by the combination of grounding-line retreat and water-column thickness changes (which are present in WARM but absent in FIXWC). The thicker water column allows ocean currents to flow unimpeded into TI.
The small melting decrease in FIXWC must be caused by changes in ice draft, despite the absence of grounding-line retreat or water-column thickness changes. This must be caused by changes in ice-base slope. While these changes have a remarkably weak effect overall, they dominate the spatial pattern of the melting response ( Figure S3 in Supporting Information S1).

of 10
In summary, the FIXGL and FIXWC experiments demonstrate that the overall change in melting is dominated by grounding-line retreat and water-column thickening, while the regional pattern of change is dominated by changes in ice-shelf slope. All of these changes may be important to the evolution of Thwaites Glacier .

Comparison to Changes in Forcing
It is unsurprising that a change in geometry induces a change in melting, so the pertinent question is whether this change is large enough to matter. In Figure 3b we compare overall melting in the WARM experiment to the COOL and NOSG experiments. The geometric response is captured in the evolution of the curves over time, while the response to forcings is captured in the offset between curves. It is immediately clear that the geometric response (11.1 Gt/y increase between WARM 2011 and WARM 2022 simulations) is substantial relative to the ocean temperature response (16.3 Gt/y increase between COOL 2022 and WARM 2022 simulations) and the subglacial inflow response (6.1 Gt/y increase between NOSG 2022 and WARM 2022 simulations). Figure 4 shows the spatial patterns of melting changes caused by these different factors. The ocean temperature response (Figure 4g) is substantially larger than the geometric response (Figure 4f) overall, but is more evenly distributed with much weaker local changes in melting, particularly near the grounding line. This is because the geometric melting response is strongly influenced by changes in currents, while the ocean temperature response reflects the more spatially uniform warming. The subglacial inflow response is focused on its inflow region (Nakayama et al., 2021) but also induces a modest increase in melting to the east, accelerating currents throughout the region (Figure 4h). These results show that the geometric melting response is of leading importance to Thwaites Glacier melting in both a local and area-averaged sense.

Discussion
Our results show a remarkably strong oceanic melting response to the recent geometrical changes in Thwaites Glacier. These changes occurred over only 12 years, but their effect on melting is comparable to extreme changes in other forcings (the coldest-to-warmest oceanic conditions on record, and a complete removal of subglacial inflow). These strong melting changes are expected to further influence ice geometry, implying a geometric melting feedback.
This geometric feedback is particularly strong when considered in the context of centennial changes. The COOL and WARM experiments have a temperature difference peaking at 1.2°C in the thermocline. This can be compared to modeled Amundsen Sea warming of just ∼0.35°C over the 20th century , or ∼1.5°C over the 21st century under extreme anthropogenic forcing (Jourdain et al., 2022). Thus, the temperature increase between COOL and WARM is a reasonable analog for the warming expected from 100 years of strong anthropogenic forcing. Based on Figure 4, we conclude that the geometric feedback over the 12-year period studied here has a comparable influence on melting to a whole century of anthropogenic ocean warming.
The production of subglacial meltwater is controlled by frictional heating at the bed of Thwaites Glacier (Hager et al., 2022;Joughin et al., 2009). The glacier flow speed approximately doubled since the 1970s . If the subglacial inflow had also doubled, that suggests a melting increase of ∼3 Gt/y over 50 years (half of the 6.1 Gt/y total influence of subglacial inflow). This weak feedback is dwarfed by the 11.1 Gt/y geometrical feedback over 12 years. However, intermittent subglacial lake drainage events may strongly influence melting on shorter timescales (Malczyk et al., 2020;B. E. Smith, Gourmelen, et al., 2017), and the resulting changes in ice geometry may further trigger geometric melting feedbacks.
The role of ice and ocean feedbacks has been investigated in models (Arthern & Williams, 2017;De Rydt & Gudmundsson, 2016;Donat-Magnin et al., 2017;Seroussi et al., 2017), but we believe this is the first time a strong geometric ocean melting feedback has been demonstrated for the real ice retreat. It is unclear how prevalent these processes are in different regions and time periods, and whether the feedback might operate in reverse during ice advance. It also remains to be determined how such fine-scale processes may be incorporated into ice-sheet projections. This geometric ocean melting feedback is expected to operate in addition to other feedbacks, such as ice acceleration related to bed geometry (Favier et al., 2014;Joughin et al., 2014) and ice damage (Lhermitte et al., 2020;Surawy-Stepney et al., 2023).
Our findings have grave implications for the influence of greenhouse-gas emissions policy on sea-level rise from Thwaites Glacier. It might be hoped that lower emissions would lead to lower future ocean warming and lower rates of ice loss. However, our results show that ocean temperature changes exert a limited influence on the melting of Thwaites Glacier because it is also subjected to a strong geometric feedback. The feedback could continue to increase melting even in the absence of further ocean warming. This suggests that greenhouse-gas emissions mitigation will not prevent Thwaites Glacier from making a substantial sea-level contribution in the coming centuries.

Conclusions
This study examines ocean melting beneath the floating section of the main trunk of Thwaites Glacier, which connects the current grounding line to the deepest grounded ice inland, with the largest potential sea-level contribution (Joughin et al., 2014;Yu et al., 2018). The results show that a strong geometric feedback enhanced ocean melting during the recent ice retreat. This feedback is caused by changes in ocean currents and temperatures, driven by changes in grounding-line position, water-column thickness, and ice base slopes. The geometric feedback over just 12 years has a comparable effect on melting to the ocean warming expected from a century of anthropogenic forcing. This has important implications for the future ice loss from Thwaites Glacier, suggesting that it may only be weakly influenced by greenhouse-gas emissions mitigation.