Limited Impact of Thwaites Ice Shelf on Future Ice Loss From Antarctica

Thwaites Ice Shelf (TWIS), the floating extension of Thwaites Glacier, West Antarctica, is changing rapidly and may completely disintegrate in the near future. Any buttressing that the ice shelf provides to the upstream grounded Thwaites glacier will then be lost. Previously, it has been argued that this could lead to onset of dynamical instability and the rapid demise of the entire glacier. Here we provide the first systematic quantitative assessment of how strongly the upstream ice is buttressed by TWIS and how its collapse affects future projections. By modeling the stresses acting along the current grounding line, we show that they deviate insignificantly from the stresses after ice shelf collapse. Using three ice‐flow models, we furthermore model the transient evolution of Thwaites Glacier and find that a complete disintegration of the ice shelf will not substantially impact future mass loss over the next 50 years.

The transition of the TWGT from a largely intact ice tongue to a loosely connected assembly of ice blocks took place as recently as 2012 (Miles et al., 2020). The currently structurally more intact TEIS is pinned by a seafloor ridge upstream of its current terminus position. The areal extent of the pinning point has gradually become smaller over the last decade Benn et al., 2022;Smith et al., 2020;Wild et al., 2022) and if the thinning of TEIS continues, the possibility of a total loss of contact between TEIS and the seafloor ridge in the near future has been suggested (Wild et al., 2022). These recent changes, together with the prospect of further disintegration and loss of pinning points should these processes continue, raise the question of the role of TWIS in modulating ice flow upstream of the grounding line.
The goal of this study is to quantify the buttressing provided by the ice shelf, and investigate the impact on upstream flow. We use ice-flow models to calculate the stresses within the ice shelf, and provide an assessment of the stresses acting along the grounding line. Transient runs are conducted to quantify the effect of an ice-shelf disintegration event on ice loss from the area over the next 100 years.

Ice Flow Modeling
We use three large-scale ice-flow models in this study: Úa, the Ice-sheet and Sea-level System Model (ISSM) and STREAMICE. All models are initialized to ensure that they produce the observed surface velocity for the observed current geometry of the ice sheet, including the positions of all grounding lines. The results of these three models, when applied to the West Antarctic Ice Sheet, have previously been compared in detail by Barnes et al. (2021) and found to produce similar projections, despite differences in initialization procedures and discretization. For that reason, most of the runs shown were conducted by Úa, and a comparison was made between these models by running a subset of the numerical experiments by all three models. As shown in the Supporting Information S1, the three models, when compared, produced almost identical results. The models solve vertically integrated equations of ice flow. Descriptions of each model and a comparison of the inversion procedures can be found in Barnes et al. (2021).
The ice sheet geometry was based on BedMachine Antarctica V2 (Morlighem et al., 2020). Basal sliding was simulated using a sliding law that combines Weertman and Coulomb law behavior as where ‖ ‖ is the norm of the basal tangential traction, and and are the Weertman and Coulomb basal tractions, defined as respectively. Here  is the flotation mask, equal to 0 where the ice is afloat and 1 otherwise, v b is the basal velocity, m is the Weertman sliding-law stress exponent, N is the effective pressure, β 2 a spatially variable basal drag parameter, and μ k the coefficient of kinetic friction. The resulting sliding law, that is, has been used in numerous ice-flow simulations, and was one of the sliding laws used as a part of the MISMIP+ protocol (Equation 11 in Asay- Davis et al., 2016). We estimated the basal drag parameter β 2 using surface-to-bed inversion of measured surface velocities, and set μ k = 1/2, and made the assumption of perfect subglacial hydrological connection. The inversion procedure is based on the use of the momentum equations and does not assume steady-state conditions at the beginning of the runs.
As Equation 1 shows, the sliding law (Equation 4) is based on the idea of forming a reciprocal weighting of the Coulomb and Weertman tractions, with each term raised to the power m. We therefore abbreviate this sliding law as rCWm. To assess the sensitivity of our results to the type of sliding law, we also ran several additional experiments using a Weertman sliding law (W).

Quantification of Ice-Shelf Buttressing
As originally introduced by Schoof (2007a), ice-shelf buttressing can be quantified using the stress buttressing number: where R n , the normal component of the resistive stress vector (i.e., traction) at grounding line, is given by and where ̂ is a horizontal unit vector normal to the grounding line, and R o is the normal component of the resistive stress vector at the grounding line in the absence of an ice shelf, that is, where ρ and ρ o are the ice and ocean densities, respectively, h is the ice thickness, and g the gravitational acceleration. Here R is the resistive stress tensor where τ xx , τ yy , and τ xy are the (horizontal) components of the deviatoric stress tensor. R o is the normal resistive stress exerted on a calving front through the action of the vertically integrated ocean pressure.
For Θ n = 1, the normal component, R n , of the resistive stress vector (i.e., normal to the grounding line, or some other specified line of interest cutting across the ice shelf), equals the vertically integrated ocean stress, R o . This is the situation of no buttressing. For Θ n < 1, the ice shelf provides greater, and for Θ n > 1 less, backstress than the ocean alone. For Θ n > 0, the deviatoric stress regime normal to the grounding line is tensile, and compressive for Θ n < 0. Note that ice-shelf buttressing, being only a function of the stress regime, can be calculated diagnostically, that is, no transient simulation is required.

Buttressing Numbers
The buttressing ratio, Θ n , for the Amundsen Sea ice shelves is shown in Figure 1 along the current grounding line positions at equal spatial intervals of 2 km. Both Pine Island Ice Shelf (PIIS), as well as the combined Crosson and Dotson ice shelves, provide significant buttressing with Θ n values typically below 0.5. The buttressing provided by Crosson and Dotson ice shelves to upstream flow is particularly large with the buttressing ratio being negative in places along their grounding lines. (Negative Θ n values imply compressive deviatoric stress regime in the direction normal to the grounding line). This high degree of buttressing is presumably caused by Bear Peninsula and the laterally confined geometry of these ice shelves. The buttressing numbers for TWIS are, generally, close to unity. (See Figure 1, and also Figures A5, A6, and A7 in Supporting Information S1).
Histograms of buttressing values, Θ n , for Pine Island, Thwaites, and Pope, Smith and Kohler glaciers are shown Figure 2a. (The exact boundaries of these three grounding-line sectors are shown in Figure A8 in Supporting Information S1.) The histograms are based on Θ n values sampled along the respective grounding lines at equally spaced 500 m intervals, and the histogram has be normalized to ensure that the sum of bar heights is equal to one. Comparing the histograms of buttressing numbers in Figure 2a for PIIS, TWIS, and Crosson and Dotson ice shelves shows that Crosson, Dotson, and Pine Island ice shelves all significantly buttress upstream flow, with Θ n mostly smaller than 0.75. In comparison, upstream flow is less buttressed by TWIS. Although TWIS generally does not strongly buttress upstream flow, Θ n is less than 0.75 along some sectors of the Thwaites Glacier grounding line. There are broadly speaking three different mechanisms that contribute to buttressing: lateral confinement, internal pinning points, and hoop-stresses that arise through divergence of ice-shelf flow. TWIS is not laterally confined, and hoop-stresses are generally not able to provide significant buttressing (Wearing et al., 2020). The ice-shelf, however, has a prominent pinning point, and it therefore may appear likely that the pinning point is the primary source of buttressing. However, as inspection of Figure 1 reveals (see also Figure A5 in Supporting Information S1 for more spatial detail), areas of high buttressing appear to be linked to localized curves and bends in the position of the Thwaites' grounding line. To investigate further the factors controlling the TWIS buttressing regime, buttressing numbers were calculated along a smoothed line running about 2 km downstream of the current grounding line (red dashed line in Supporting Information S1, Figure A8). The resulting histogram of Θ n values along this smoothed curve is shown in Figure 2b (red), together with the distribution of those values along the grounding line (blue). As the figure shows, Θ n values along the smooth down-stream grounding line are shifted toward unity. Hence, a significant part of the buttressing along the current grounding line is due to those small scale bends and curves in the grounding absent from the smoothed curve, with the remainder primarily related to the existence of the pinning point within TEIS.

Sea Level Rise
To quantify the importance of TWIS to future sea level rise (SLR), we conducted transient simulations of a large sector of the West Antarctic Ice Sheet including PIIS, TWIS, Crosson and Dotson ice shelves, and the catchment area of Thwaites Glacier ( Figure A1 in Supporting Information S1). The surface mass balance was taken from RACMO2.3p2 (van Wessem et al., 2018), and ocean-induced melt prescribed using a depth-dependent melt-rate parameterization. To test the sensitivity of our results to the melt parameterization, we used two different parametrizations labeled M0 and M4 that assume a maximum melt rate of 100 and 50 m/yr, respectively ( Figure A2 in Supporting Information S1).
We first ran a reference run, starting with the current geometry of the ice sheet. We then ran further runs, using the same initialization procedure and external forcings, but now with either the PIIS or TWIS ice shelves removed at the beginning of the simulations, again for both the M0 and M4 melt parametrizations. In addition, we tested the sensitivity of our results to mesh resolution, and a selection of regularization parameters. Most of the simulations were run with the ice sheet model Úa, but for model comparison, we also ran a subset of the experiments using the ice-flow models ISSM and STREAMICE.  Figure A8 in Supporting Information S1).
Modeled mass losses over 100 years, relative to the respective reference runs, are shown in Figure 3. The removal of TWIS causes an increase in mass loss equivalent to about 2 mm over 100 years. This value is almost the same for both the M0 and the M4 melt parametrizations. This insensitivity to the choice of melt parameterization is expected, as factors common to two runs tend to (approximately) cancel out in model differences. The same TWIS removal experiment was also conducted using ISSM and STREAMICE and resulted in a similar mass loss of a few mm over 100 years, even though somewhat higher SLR perturbations (∼5 mm) were obtained with STREAMICE. Initializing the ice-flow model Úa using inversion estimates of the basal sliding parameter (β 2 ) obtained by STREAMICE, resulted in somewhat higher SLR estimates, suggesting that these small differences between model outputs might be related to differences in inversion methodology. Modeled SLR perturbations for the Weertman and the regularized Coulomb laws were almost identical for the first 50 years, or within 0.2 mm SLR from each other, but started to deviate for the remainder of the 100 years runs ( Figure A4 in Supporting Information S1). Similar insensitivity of modeled Thwaites Glacier behavior over the first few decades to basal sliding and ice shelf melt scenarios was previously noted by Barnes and Gudmundsson (2022) and Yu et al. (2018). In all cases considered, we find that removing TWIS causes a modification in modeled SLR of 1-2 mm SLR equivalent over the first 50 years. Compared to ongoing mass loss, and ongoing rate of SLR of a few mm per year, this additional sea level contribution due to the loss of TWIS appears negligibly small.
When conducting a similar experiment for Pine Island Glacier, we find that removing the PIIS causes about 20 mm additional increase in global sea level. Hence, in relative terms removing PIIS has about an order-of-magnitude larger SLR contribution, than the removal of TWIS. We expect our estimate of the impact of removing PIIS to be a lower bound, as our model domain does not include all of Pine Island Glacier catchment area, and that the actual impact might be larger. However, this will not affect our conclusion that the removal of PIIS has a much larger impact on ice loss, than the removal of TWIS.
To illustrate further the impact of TWIS removal on upstream flow, examples of modeled relative changes in height above flotation (HAF) and velocities are provided in Figure 4, 20 years after start of the TWIS removal run. As above, the relative changes shown are with respect to the reference run where TWIS was not removed. Panel   Figure A2 in Supporting Information S1) using the sliding law rCWm defined by Equation 4. Each curve represents the difference in ice volume loss with respect to a reference runs where the ice shelves were not removed.
as a consequence of the buttressing regime in this sector of the grounding line (See Figure 1 and, for more local detail; Figure A5 in Supporting Information S1). Along several grounding-line segments of this eastern sector the buttressing numbers are larger than 1. For those segments, the ice shelf therefore pulls rather than restrains ice flow. Thus, when the ice shelf is removed, ice velocities are reduced (i.e., less pull), resulting in thickening and concomitant increase in HAF. This relative reduction in ice speed can be seen in Figure 4b, showing the perturbed velocities Δv = v(TWIS removed) − v(TWIS not removed).
For the western sector, on the other hand, we find that removing TWIS causes surface lowering and an increase in flow speed toward the grounding line ( Figure 4). Here, HAF values are negative and the perturbed velocity vectors point toward the grounding line. These results can also readily be understood as a consequence of the changes in buttressing following the removal of the ice shelf. In this westerly region of Thwaites' grounding line, and especially directly upstream of a small local embayment of the grounding line in this area, buttressing numbers are considerably smaller than 1 and even negative in places. Thus, the ice shelf restrains upstream flow, even to the degree of causing the normal deviatoric stresses to become negative, that is, compressive. Hence, when the ice shelf is removed, upstream ice flow becomes less restricted and ice velocities increase.
Overall, the reduction in HAF outweighs any gains, and with time the ice volume above flotation (VAF) decreases. The result is a small additional contribution to sea level change of a few mm over the first 100 years due to the removal of TWIS.

Discussion
Because of the importance of Thwaites Glacier for future SLR, a large number of modeling studies have previously considered various aspects of its flow regime (Barnes & Gudmundsson, 2022;Benn et al., 2022;Docquier et al., 2014;Joughin et al., 2014;Parizek et al., 2013;Santos et al., 2021;Seroussi et al., 2017;Waibel et al., 2018;Yu et al., 2017Yu et al., , 2018. While none of these studies have provided direct quantification of ice-shelf buttressing as done here, several published studies have assessed the importance of the pinning point of TEIS. In the following we provide a short overview of related studies and put our results in context with previous work. Early modeling studies, for example, Docquier et al. (2014) and Parizek et al. (2013), assessed the influence of buttressing parameterization using flow-band models of Thwaites. Parizek et al. (2013) concluded that the Thwaites Ice Tongue provides limited "stability." Here, the term stability was presumably used to indicate buttressing strength. Docquier et al. (2014) estimated the impact of ice-shelf pinning points by adjusting the boundary condition applied at the calving front. They concluded that TWGT exerts limited buttressing on inland ice. While arguably ill-suited for quantification of buttressing-an effect inherently related to side drag and transverse variations in flow conditions -, these early studies already suggested that Thwaites Glacier is not strongly buttressed by its abutting ice shelves.  Reese et al. (2018) provided an assessment of ice-shelf buttressing of all grounding lines of the Antarctic Ice Sheet, including those of Thwaites and Pine Island Glacier. Here, we followed their methodology and similarly quantified buttressing through the buttressing number Θ n . The focus in Reese et al. (2018) was on the applicability of an explicit formula for grounding-line flux, and they did not provide specific information on the buttressing strength of TWIS. Reese et al. (2018) however provides an assessment of the buttressing flux response number, (θ B ), which describes the sensitivity of grounding-line flux to perturbation in ice-shelf thickness. While we cannot directly compare the buttressing (stress) number Θ n , calculated here, with θ B as calculated by Reese et al. (2018), Figure 1d in Reese et al. (2018) shows that flux response of TWIS is much smaller than that of PIIS and the Crosson-Dotson Ice Shelves, which is in good general agreement with our findings. Gudmundsson et al. (2019) estimated the instantaneous change in velocities and grounding line flux due to observed thinning of Antarctic Ice Shelves. This instantaneous response is caused by ice-shelf buttressing and is absent, that is, numerically equal to zero, for unbuttressed ice shelves. They did not provide a separate estimate for the flux perturbation of Thwaites Glacier and lumped the response together with that of the neighboring Crosson and Dotson ice shelves. However, their modeled change in speed shows the same spatial pattern with decrease in velocities upstream of TEIS and increase for TWGT.
In a modeling study using a similar type of Shallow-Shelf approximation vertically integrated ice-flow model with initialization similar to that done here, Joughin et al. (2014)  The conclusion that TWIS provides very limited buttressing to upstream flow has several important consequences. It suggests, as already shown by Joughin et al. (2014), that changes in ocean conditions and in related ocean-induced ice-shelf melt have limited potential to impact the dynamics of Thwaites Glacier in its current configuration. It should be stressed, however, that although currently Thwaites Glacier is mostly unbuttressed by its abutting ice shelves, this situation may potentially change in the future. Should the grounding line of Thwaites Glacier retreat significantly, and in the process form a new laterally confined ice shelf upstream of its current grounding line-as for example, occurs in our transient simulations-the resulting extended ice shelf is likely to provide a greater degree of buttressing. Hence, although changes in the thermal ocean conditions are currently unlikely to have any significant immediate impact of Thwaites Glacier, this situation could change in the future. It is also possible that, at some point in the past, TWIS may have been pinned more extensively than today and provided more buttressing stresses to upstream grounded ice.
This raises the question as to what is responsible for the current mass loss of Thwaites glacier. The insensitivity of the grounding line to changes in downstream conditions suggests that these changes are a transient response to previous, yet unidentified, variations in external conditions. Joughin et al. (2014) suggested similar reasons for current changes and stated that "the retreat that does occur is largely driven by the non-steady-state fixed velocity imposed at the start of the simulation." No modeling evidence published to date suggests that ongoing changes are related to a dynamical instability, and Benn et al. (2022) state specifically that: "there is currently no evidence that the imminent loss of TEIS will hasten marine ice sheet instability or the demise of Thwaites Glacier" (page 2550, Section 5.3 in Benn et al., 2022), although the same authors also write, in an apparent contradiction, that "the ongoing acceleration of parts of the ice shelf suggests the crossing of a threshold from stable to unstable" in (page 2546, Section 1 in Benn et al., 2022).
Our conclusions do not support a previous statement by Pettit et al. (2021) that removal of TEIS has the potential to increase the contribution of Thwaites Glacier to SLR by up to 25%. To the contrary, we conclude that the removal of TEIS will have almost no discernible effect on future mass loss of the glacier. We also do not find support for the notion put forward by Davis et al. (2023) that the rates of ocean-induced melt over TWIS can be expected to force notable change to grounded ice. We found no information in either Pettit et al. (2021) or Davis et al. (2023) to support their implicit assumption that the laterally unconfined TEIS provides significant buttressing to upstream flow.
In a recent study Urruty et al. (2022) investigated the stability regime of all grounding lines of the Antarctic Ice Sheet and concluded that they are currently dynamically stable. This finding does not imply that the grounding lines of Thwaites cannot become unstable at some future point in time (e.g., Reese et al., 2022). It also does not contradict the findings of Favier et al. (2014) or Joughin et al. (2014) that Thwaites and Pine Island Glaciers both show onset of unstable and irreversible grounding line retreat in simulations of their future behavior. Dynamically unstable grounding line retreat of marine type ice sheets, that is, the marine ice sheet instability, is commonly seen in numerical simulations. Rosier et al. (2021) recently showed, using the concept of critical slowing, that the neighboring Pine Island Glacier possesses at least three distinct 'tipping points' marking the onset of irreversible self-enhancing retreat. However, no numerical modeling work has shown that Thwaites Glacier is currently undergoing an irreversible retreat. The limited buttressing capacity of the TWIS makes it unlikely for it to significantly impact the stability regime of the grounding line. In conclusion, currently Thwaites glacier is stable and unsteady.

Conclusions
We have quantified the buttressing provided by TWIS to upstream flow in several different ways. First, we estimated the normal traction to the grounding line of Thwaites glacier for the current configuration of the ice shelf, and for the ice shelf removed. For Thwaites Glacier the ratio between these two quantities, that is, the buttressing number Θ n , deviates by less than 20% from unity, although in some other areas it is as low as 0.5. Significant buttressing-here defined as grounding-line sections where Θ n is either smaller than 0.75 or larger than 1.25 -, is mostly found to be related to localized bends and curves of the grounding line, rather than to the presence of the pinning point currently located within the ice shelf.
Second, we compared the buttressing numbers of Thwaites' grounding line with those of the neighboring Pine Island and Crosson and Dotson ice shelves. The buttressing numbers of Crosson and Dotson ice shelves range from about −0.5 to 0.5, and those of PIIS are mostly within 0-0.5. Hence, in contrast to TWIS, those neighboring ice shelves significantly buttress upstream flow.
Third, we conducted transient simulations of Thwaites Glacier to assess the impact of a possible disintegration of the TWIS on future mass loss. The difference between those two scenarios is small, corresponding to a few mm of global sea level change over 100 years. Thus, the presence or absence of the ice shelf has almost no impact on predicted mass loss.
We conclude that TWIS is of limited relevance to upstream flow dynamics in its current configuration, and there appears to be no reason to expect a possible disintegration of the ice shelf to meaningfully impact SLR projections over the next 50 years.

Data Availability Statement
The three ice-sheet models are open source. Úa can be downloaded from https://github.com/GHilmarG/UaSource (https://doi.org/10.5281/zenodo.3706623). ISSM is open source and can be downloaded and installed from https:// issm.jpl.nasa.gov/download/ as either binaries or from the source code. STREAMICE is part of the Massachusetts Institute of Technology general circulation model (MITgcm) and the source is freely available for download (https://github.com/MITgcm/MITgcm). No new observational data was generated as a part of the study and all data set used in this study are freely available. Velocity data products were downloaded from https://its-live.jpl. nasa.gov/ and the bedrock data from https://nsidc.org/data/nsidc-0756/versions/2.