# Subglacial water amplifies Antarctic contributions to sea-level rise

Chen Zhao, Rupert Gladstone, Thomas Zwinger, Fabien Gillet-Chaulet, Yu Wang, Justine Caillet, Pierre Mathiot, Leopekka Saraste, Eliot Jager, Benjamin K Galton-Fenzi, et al.

► To cite this version:

Chen Zhao, Rupert Gladstone, Thomas Zwinger, Fabien Gillet-Chaulet, Yu Wang, et al. Subglacial water amplifies Antarctic contributions to sea-level rise. Nature Communications, 2025, 16 (1), pp.3187. 10.1038/s41467-025-58375-4. hal-05450904

HAL Id: hal-05450904

https://hal.inrae.fr/hal-05450904v1

Submitted on 9 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

# Subglacial water amplifies Antarctic contributions to sea-level rise

Received: 11 September 2024

Accepted: 12 March 2025

Published online: 07 April 2025

<figure><img src="image_3.png" /></figure>

Check for updates

Chen Zhao Ⓔ<fnref n="1" />,²,⁸✉, Rupert Gladstone Ⓔ³, Thomas Zwinger Ⓔ⁴, Fabien Gillet-Chaulet Ⓔ⁵, Yu Wang<fnref n="1" />, Justine Caillet Ⓔ⁵, Pierre Mathiot Ⓔ⁵, Leopekka Saraste Ⓔ⁴, Eliot Jager Ⓔ⁶, Benjamin K. Galton-Fenzi Ⓔ<fnref n="1" />,²,⁷, Poul Christoffersen Ⓔ<fnref n="1" />,²,⁸ & Matt A. King Ⓔ²,⁹

Antarctica's contribution to global sea-level rise is deeply uncertain, with subglacial water suspected to play a critical role, yet its impact remains unclear. We demonstrate that water at the base of ice sheets influences sliding behaviour and that its exclusion from models can underestimate sea-level rise projections and delay the predicted onset of tipping points. Here we use an Antarctic Ice Sheet model (Elmer/Ice) to explore how different assumptions about water pressure at the ice base affect sea-level rise projections from 2015 to 2300. Our results indicate that incorporating subglacial water can amplify ice discharge across the Antarctic Ice Sheet by up to threefold above the standard approach, potentially contributing an additional 2.2 metres to sea-level rise by 2300. Notably, a smoothly decreasing basal drag near the grounding line more than doubles grounding line flux by 2300 relative to scenarios where effective pressure is simplified into a spatially constant coefficient. Basin-specific responses vary significantly, with some scenarios advancing tipping points by up to 40 years. These findings underscore the critical need to integrate evolving subglacial hydrology into ice sheet models.

The Antarctic Ice Sheet provides the greatest source of uncertainty in current sea-level rise projections<fnref n="1" />,². This deep uncertainty, and the risk of potential tipping points that can cause rapid ice loss, is driven by the lack of understanding about fundamental processes influencing ice-sheet dynamics³. In the Ice Sheet Model Intercomparison for CMIP6 (ISMIP6), more than half of the total uncertainty is due to discrepancies between different ice-sheet simulations⁴,⁵. While it is challenging to model the complex interactions of Antarctic Ice Sheet with climate⁴,⁶, a key process that remains largely unconstrained by observations is basal sliding.

The widespread production of subglacial water beneath the Antarctic Ice Sheet, and the potential intrusion of warm seawater

near the grounding line, plays a crucial role in modulating basal sliding and enabling rapid ice flow⁷⁻¹¹. Subglacial water drainage beneath the ice sheet occurs through two primary pathways: distributed flow and channelized flow¹². The existence of channelized drainage beneath the Antarctic Ice Sheet has been inferred from the alignment of predicted subglacial meltwater flux across the grounding line with observed 'subglacially-sourced' sub-ice-shelf channels¹³⁻¹⁶. These drainage channels are fed by a distributed drainage system, which consists of inter-linked cavities and serves as an initially inefficient pathway for subglacial meltwater transport¹². The source of the distributed and channelised subglacial hydrologic systems is meltwater production

<footnote n="1">¹Australian Antarctic Program Partnership, Institute for Marine and Antarctic Studies, University of Tasmania, Hobart, Australia. ²Australian Centre for Excellence in Antarctic Science, University of Tasmania, Hobart, Australia. ³Arctic Centre, University of Lapland, Rovaniemi, Finland. ⁴CSC-IT Center for Science, Espoo, Finland. ⁵Univ. Grenoble Alpes, CNRS, INRAE, IRD, Grenoble INP, Institut des Géosciences de l'Environnement, Grenoble, France. ⁶Institute for Atmospheric and Earth System Research / Physics, Faculty of Sciences, University of Helsinki, Helsinki, Finland. ⁷Australian Antarctic Division, Kingston, Australia. ⁸Institute for Marine and Antarctic Studies, University of Tasmania, Hobart, Australia. ⁹School of Geography, Planning, and Spatial Sciences, University of Tasmania, Hobart, Australia. ✉e-mail: chen.zhao@utas.edu.au</footnote>
from both frictional heating, due to basal sliding and internal ice deformation, and geothermal heating from the Earth's interior⁷.

The Antarctic Ice Sheet flows through two primary mechanisms: internal ice deformation and sliding over the bed, with the latter dominating the dynamics of all fast flowing outlet glaciers. Accurately representing basal sliding is essential for projecting future ice sheet change. The processes controlling sliding broadly fall into two regimes: The 'Weertman' regime, represented in models by expressing basal resistance as a power law function of sliding speed¹⁷, and the 'Coulomb' regime, represented in models by expressing basal resistance as a function of the effective normal force at the ice-bedrock interface¹⁸,¹⁹. For ice sheets, this normal force is given by the ice overburden pressure minus the basal water pressure, henceforth referred to as 'effective pressure', N. The Weertman regime is expected to apply for the case of hard bedrock, low sliding speeds, and high effective pressure. The Coulomb regime is expected to apply for the case of a soft (sediment) bed, fast sliding speeds, and low effective pressure. Recent laboratory experiments²⁰ have yielded compelling evidence that a 'regularised Coulomb' (hereafter 'RC') sliding parameterisation, capturing both Weertman and Coulomb regimes, provides a unified framework for modelling sliding across diverse bed types¹⁸,²¹⁻²³. When applying RC sliding, the regime shift between Weertman and Coulomb type sliding is governed by sliding speed, effective pressure, and parameter choices.

Because the actual distribution of effective pressure under the Antarctic Ice Sheet is unknown, model-based sea-level rise projections typically either apply Weertman sliding everywhere, which neglects effective pressure entirely (e.g., refs. 24–26) or make simple assumptions about effective pressure when applying RC sliding (e.g., refs. 27–29). A summary of sliding laws used in three major Model Intercomparison Projects (Table S1) reveals that the majority of Antarctic Ice Sheet models (e.g., 10 of 16 in ISMIP6-2100², 10 of 19 in ISMIP6-2300³, and 7 of 15 in ABUMIP³⁰) employed the Weertman sliding relation.

For models employing RC sliding, assumptions must be made about effective pressure. One such assumption is perfect hydrological connectivity between the subglacial system and the ocean²⁷⁻²⁹,³¹. Another one is to prescribe effective pressure as proportional to ice overburden pressure³²,³³. While practical, these approaches lack robust observational support and may lead to systematic over- or underestimation of basal water pressures across large regions²⁷,³¹,³⁴⁻³⁶. An alternative approach with a stronger physical justification is to simulate effective pressure using more advanced hydrology models. Models such as GlaDS (Glacier Drainage System³⁷) and SHAKTI³⁸ offer distributed and channelized drainage systems that can capture spatial variability in effective pressure. While these models improve physical realism, they are computationally expensive and sensitive to uncertainties in basal boundary conditions (e.g., bed roughness, geothermal heat flux)³⁵,³⁹.

A simpler version of RC sliding has been proposed that removes explicit dependence on effective pressure²¹,⁴⁰. In this approach, effective pressure is instead implicitly embedded within spatially optimized

parameters of the sliding relation, which are calibrated to reproduce observed ice flow speeds²¹. While computationally efficient and practical for large-scale ice sheet simulations, this method may overlook important spatial and temporal variability in effective pressure.

A further assumption about the effective pressure and its impact on sliding near the grounding line can be imposed using a 'height-above floatation' (HAF) scaling²¹. HAF is defined as the height of the upper ice surface above sea level minus the height required for the ice overburden pressure to match the ocean water pressure at the bed. This assumption scales the basal resistance or the effective pressure from a prescribed HAF value down to zero at the grounding line, reflecting the influence of warm seawater intrusion beneath grounded ice³¹. However, the inland extent of this effect is poorly constrained due to limited observations, and this approach neglects significant spatial and temporal variability driven by evolving glacier geometry⁴²,⁴³. Both the HAF-scaling and perfect ocean connection assumption enforce the floating condition at the grounding line. However, HAF-scaling, applied only near the grounding line, adjusts basal resistance based on ice elevation relative to flotation, potentially under- or overestimating effective pressure. In contrast, the perfect ocean connection assumes uniform hydraulic connectivity across the domain, systematically overestimating effective pressure by neglecting local hydraulic potential gradients.

Building on these approaches, we conduct six experiments (Table 1) to assess the impact of these differing assumptions about spatial and temporal variations in effective pressure on projections of future Antarctic Ice Sheet behaviour. These experiments include approaches with both implicit and explicit representations of effective pressure, including the assumption of perfect ocean connection, hydrology model-derived estimates, and height-above-flotation (HAF) scaling. For comparison against the more physically justifiable RC sliding parameterisation, we also apply the linear Weertman sliding relation (hereafter 'LW'), which, despite its tendency to underestimate mass loss²¹, is still widely used for its simplicity in both numerical and computational aspects²⁵,²⁶,³⁰. By isolating these parameterisations, we evaluate their influence on basal sliding and the resulting ice dynamics, providing insights into the role of hydrology in ice sheet modelling.

## Results

### Subglacial effective pressure must be included in sliding relations

Our simulations show large differences in projected total Antarctic ice-mass loss between experiments. The largest differences are between LW and the RC relation with implicit effective pressure and HAF-scaling (RC\_iN\_HAF; Fig. 1). By 2300, the grounding line flux using RC\_iN\_HAF (8453 Gt yr⁻¹) is almost four times larger than using LW (2245 Gt yr⁻¹) under a high emission scenario. While LW suggests a negative sea-level contribution, the RC simulations all show a positive sea-level contribution by 2300, with RC\_iN\_HAF producing the largest rise of ~2 m sea-level equivalent (SLE; Fig. 1b).

The LW relation, which ignores the role of subglacial water pressure, underestimates sea-level rise contributions relative to the RC

<table><caption>**Table 1 | Summary of sliding parameterisations and effective pressure treatments**</caption><thead><tr><td>Experiment</td><td>Effective Pressure, N</td><td>Description</td><td>HAF Scaling<sup>a</sup></td></tr></thead><tbody><tr><td>LW</td><td>None</td><td>Sliding independent of effective pressure.</td><td>No</td></tr><tr><td>RC\_iN</td><td>Implicit</td><td>N subsumed into a constant sliding coefficient.</td><td>No</td></tr><tr><td>RC\_iN\_HAF</td><td>Implicit</td><td>N subsumed into a constant sliding coefficient.</td><td>Yes</td></tr><tr><td>RC\_eN\_POC</td><td>Explicit</td><td>Perfect ocean connection assumption for N.</td><td>No</td></tr><tr><td>RC\_eN\_GlaDS</td><td>Explicit</td><td>N simulated using a distributed hydrology model.</td><td>No</td></tr><tr><td>RC\_eN\_GlaDS\_HAF</td><td>Explicit</td><td>N simulated using a distributed hydrology model.</td><td>Yes</td></tr></tbody></table>

<footnote><sup>a</sup>HAF Scaling Height-above-flotation scaling.</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 1 | Projected Antarctic Ice Sheet mass loss and sea-level contributions.**</figcaption></figure>

Time series (2015–2300) of a grounding line flux and b sea-level contributions based on changes in volume above flotation from the Antarctic Ice Sheet across different experiments under the high emission scenario using UKESM model. ‘iN’ means implicit representation of effective pressure through a spatially tuned

<figure><img src="image_5.png" /></figure>

parameter while ‘eN’ means explicit representation of effective pressure. ‘HAF’ means height-above flotation scaling. GlaDS means simulated effective pressure from the GlaDS hydrology model, while POC means effective pressure is calculated based on perfect ocean connection assumption.

<figure><img src="image_6.png" /><figcaption>**Fig. 2 | Time series (2015–2300) of sea-level contributions in experiment using linear Weertman law (dashed lines) and regularised Coulomb law (solid lines) under various climate models.** Thick and thin lines represent high and low emission scenarios, respectively. The sea-level contributions are determined by the changes in volume above flotation relative to 2015.</figcaption></figure>

relation under both low and high emission scenarios (Fig. 2). Under high emission scenarios, all RC experiments produce a positive sea-level contribution, whereas this is the case for only two LW experiments. Under the two low emission scenarios, the RC relation suggests near equilibrium by 2300, which aligns more closely with observed mass loss in recent decades¹⁴, as opposed to the substantial mass gain (UKESM) or loss (NorESM1-M) with the LW relation (see further analysis in the Supplementary Material).

Substantial variations of sea-level rise contributions are produced in each Antarctic sector depending on the choice of sliding relation used (Fig. 3). The West Antarctic Ice Sheet exhibits the largest differences. For West Antarctica, under a high emission scenario and using the RC relation, we obtain a net mass loss of 6–24 mm SLE by 2100 for all climate models, except CCSM4, which predicts a mass gain. Conversely, with the LW relation, three of five climate models produce a mass gain. Beyond 2100, both sliding relations produce positive sea-level rise contributions from West Antarctica. However, using the RC relation results in contributions up to five times higher than using the LW relation by 2300. For the East Antarctic Ice Sheet, all climate

models predict net mass gain for this century, with LW predicting larger numbers than RC. However, East Antarctica starts to lose more ice after 2100 using the RC relation, with sea-level rise reaching 17–93 mm by 2200 and 23–248 mm by 2300. In contrast, LW produces overall mass gain, with three of five the climate models producing negative sea-level rise over the period 2200–2300. The Antarctic Peninsula shows less sensitivity to the choice of sliding relations, due to the steep valley glaciers and low water pressure in this region, where the high driving stress means the ice flow is relatively insensitive to the choice of sliding relation.

The RC simulations, incorporating either an explicit or implicit dependence on effective pressure, consistently predict higher mass loss than LW under strong warming scenarios, and exhibit a wide spread in sea-level projections due to differing assumptions about effective pressure distribution. These results highlight the importance of representing effective pressure in Antarctic Ice Sheet projections and carefully evaluating the assumptions about its spatial distribution.

## Subglacial effective pressure controls ice fluxes

Our simulations reveal the treatment of effective pressure significantly affects the spatial tuning of sliding coefficients, which, in turn, impacts the evolution of basal shear stress. For example, at the continental-scale, the implicit representation of effective pressure produces higher grounding line flux and mass loss than its explicit treatment (Fig. 1). A negative sea-level rise contribution before 2100 is simulated in most RC experiments (Fig. 1b), which is primarily due to high surface mass inputs from the climate models (see Supplementary Material), however, the ice dynamic loss characterised by grounding line flux is still positive in this century (–2000 Gt yr⁻¹).

The grounding line flux in the implicit run (RC_iN) increases after 2100, peaking at 4984 Gt yr⁻¹ around 2250. However, when HAF-scaling is introduced (RC_iN_HAF), the results diverge dramatically, showing exponential growth in grounding line flux and reaching over 8000 Gt yr⁻¹ – the highest ice discharge across the grounding line–displaying a significant disparity in magnitude. This suggests that HAF-scaling enhances basal sliding due to the progressive reduction of basal shear stress over time as the ice approaches flotation. The corresponding sea-level rise contribution reaches nearly 2 m by 2300, almost three times the increase relative to the non-HAF-scaling case (RC_iN). The divergence between the two curves beyond 2100 highlights how poorly constrained basal shear stress near the grounding line can significantly alter long-term sea-level rise projections.
<figure><img src="image_7.png" /><figcaption>**Fig. 3** | Sea-level contributions from various sectors of Antarctica across different periods. Sea-level contributions (sea level equivalent, SLE) from West Antarctic (WAIS), East Antarctic (EAIS), and Antarctic Peninsula (AP) from regularised Coulomb with implicit effective pressure and HAF-scaling (RC\_iN\_HAF; filled) and</figcaption></figure>

linear Weertman (hachured) sliding relations over **a** 2015–2100, **b** 2100–2200, **c** 2200–2300, and **d** 2015–2300. Different colours represent various climate models with high emission scenarios.

Under the explicit representation of effective pressure, the perfect ocean connection assumption has been commonly used in previous regional and continental-scale studies²⁷⁻²⁹, due to limited observations and the less advanced state of hydrology models. Under such assumption, RC\_eN\_POC updates effective pressure over time, which is more realistic near the grounding line where the ice approaches flotation, although this approach systematically underestimates basal water pressure which is not only influenced by the sea level but also the subglacial hydrological system (Fig. S1). While this parameterisation ensures the Coulomb regime dominates near the grounding line, the calculated effective pressure based on ice geometries is consistently higher as compared with simulations using the hydrology model (see Supplementary Material), resulting in the second lowest mass loss (540 mm SLE). The grounding line flux increases from 2185 Gt yr⁻¹ to 3698 Gt yr⁻¹, with a peak flux of ~4245 Gt yr⁻¹ at the year 2240.

While we adopted the simulated effective pressure from the GlaDS hydrology model, the sole experiment enforcing a constant effective pressure (RC\_eN\_GlaDS) produces the least multi-century ice-mass loss (183 mm SLE by 2300 in Fig. 1b), with the grounding line flux increasing from 2012 Gt yr⁻¹ to 3209 Gt yr⁻¹, with a later and lower peak at ~3578 Gt yr⁻¹ around the year 2270. Similar to the comparison between RC\_iN and RC\_iN\_HAF, with a HAF-scaling applied to the constant effective pressure from GlaDS (RC\_eN\_GlaDS\_HAF), the grounding line flux matches RC\_eN\_POC by the year 2140, and subsequently undergoes a rapid increase, peaking at 5152 Gt yr⁻¹ around the year 2200. The sea-level rise contribution reaches 682 mm SLE, nearly four times of the case with a constant effective pressure.

The assumption made about effective pressure substantially impacts the timing of tipping points (Fig. 1a), whereby peaks in grounding line flux correspond to rapid grounding line retreat along a retrograde slope. Both experiments with HAF-scaling reach a maximum grounding line flux around 2200, at least 40 years earlier than other experiments. The flux in the explicit HAF-scaled experiment with GlaDS gradually decreases after 2200, whereas the implicit treatment of effective pressure in RC\_iN\_HAF produces an additional maximum around 2250 which continues increasing by 2300. The modification of effective pressure through HAF-scaling substantially amplifies the grounding line flux. This underscores the significance of the magnitude and spatial distribution of effective pressure near the

grounding line, and subsequent impact on ice sheet discharge and sea-level rise.

Almost all basins exhibit the highest mass loss from RC\_iN\_HAF and lowest from LW (Fig. 4). Four basins with reverse bed slopes, that is G-H (feeding Pine Island and Thwaites Glaciers), E-EP (Ross Ice Shelf), D-Dp (Cook Glacier) and Cp-D (Totten Glacier), are chosen for a closer analysis because these are especially vulnerable to marine ice sheet instability (MISI)⁴⁵. In West Antarctica, Basin G-H, feeding Pine Island Glacier (PIG) and Thwaites Glacier, has experienced significant grounding line retreat and is particularly sensitive to the MISI²⁵,⁴⁶. Most experiments suggest a gradual increase in ice flux to 275–370 Gt yr⁻¹ by 2300. With implicit effective pressure and HAF-scaling, the projection is nearly 10 times larger than other parameterisations. In this region, the grounding line at 2300 retreats more than 300 km inland compared to its initial position in 2015 (Fig. 5c), while other experiments show a retreat of about 60 km. Notably, the grounding line retreats over 100 km in just 10 years around 2280 (Fig. 5c). The timing of reaching a tipping point varies, with RC\_eN\_POC reaching the tipping point around 2095, almost 100 years earlier than RC\_eN\_GlaDS\_HAF (2185) and LW (2215). In contrast, for Basin Ep-F, which feeds the Ross Ice Shelf, the total grounding line flux and movement is relatively invariant across all experiments.

In East Antarctica, the Wilkes Subglacial Basin (Basin D-Dp) projection with implicit effective pressure and HAF-scaling reaches its first tipping point with a peak grounding line flux of 503 Gt/yr by 2190, a factor of 1.6 larger and nearly a century earlier than the experiment with effective pressure calculated based on evolving ice geometry (RC\_eN\_POC). The Recovery basin (Basin Jpp-K) reaches its first tipping point around 2200 across all experiments, but another tipping point is reached around 2260 in RC\_iN\_HAF that is two times larger than the other RC experiments. In contrast, the Aurora Subglacial Basin (Basin Cp-D) is less sensitive, with peaks between 2170 and 2190 and a lower rate of mass loss until 2300. Separate research on this basin shows that ice-hydrology interaction increases sea-level contribution by 30% at 2100, compared to an ice-only model with constant effective pressure from a hydrology model⁹. This aligns with the comparison presented here between the two experiments adopting effective pressure from GlaDS, where HAF-scaling enhances the grounding line flux by approximately 20%.
<figure><img src="image_8.png" /><figcaption>**Fig. 4 | Time series (2015–2300) of grounding line flux (Gt yr⁻¹) across various basins under different experiments.** The background image is the simulated effective pressure and channel area (cross-sectional area of each channel)</figcaption></figure>

generated from GlaDS. Y axis is the grounding line flux (unit: Gt yr⁻¹) and X axis is the time (unit: year). The Antarctic drainage basins are from E. Rignot and J. Mouginot⁶⁵.

## Increased sensitivity to effective pressure near the grounding line

To better understand how assumptions about effective pressure can lead to such a wide spread of sea-level rise projections for RC sliding, we examine the physical processes near the grounding line. As subglacial water crosses the grounding line into the ocean, it reaches pressure balance with the ocean, driving the effective pressure at the bed near the grounding line to nearly zero. This results in a Coulomb regime with low basal resistance, with an upstream regime shift toward Weertman sliding. As the grounding line retreats, this low-drag region migrates upstream, allowing further speed up of ice flow, providing a positive feedback enhancing grounding line retreat and ice sheet mass loss. Grounding line retreat also exposes more of the ice shelf to ocean-induced melting, further accelerating the process. In regions with low gradients in bedrock elevation and ice thickness, such as the Siple Coast feeding the Ross Ice Shelf, this low-drag Coulomb regime can extend far from the grounding line (Fig. S12c), forming an extensive ice plain⁴⁷.

All RC implementations in this study (see Methods), except for the later stages of the RC_eN_GlaDS simulation (after significant grounding line retreat), impose a low-drag Coulomb region near the grounding line due to assumptions about the spatial distribution of effective pressure, either explicitly (e.g., the perfect ocean connection assumption) or

implicitly (e.g., through the HAF-scaling). These differing assumptions lead to variations in the location of the Weertman/Coulomb regime shift and the extent of the low-drag region. For example, the threshold $h_T$ in HAF-scaling determines the extent of basal drag reduction upstream of the grounding line as a function of ice geometry⁴⁸. Separate tests (see Supplementary Material) using a lower $h_T$ value (46 m) underscore the sensitivity of the low-drag region's extent to this threshold, highlighting its importance in sea-level rise projections and raising concerns about the broader validity of this approach.

Taking Thwaites Glacier as an example, the RC_iN_HAF simulation exhibits large-scale collapse during the 23rd century, whereas the other simulations show a more modest retreat. While changes in basal drag before 2015 are relatively small (Fig. S14d), the area of low basal drag in 2270 for RC_iN_HAF extends more than 50 km inland (Fig. S15d), much farther than in the other simulations (typically 10–15 km upstream of the grounding line). This enlarged low-drag region, which represents the formation of a near-flotation ice plain, is conditioned by the $h_T$ parameter as well as the bedrock geometry of the region into which the grounding line has retreated. The existence of such a zone of low basal drag explains the significantly high mass loss for Thwaites Glacier in this simulation. A similar low-drag ice plain exists at the ice streams feeding the Ross Ice Shelf⁴⁷, which is also a sensitive region in our model (Fig. S4).
<figure><img src="image_9.png" /><figcaption>**Fig. 5 | Time series (2015–2300) of ice-geometries along a Thwaites Glacier flowline.** Ice-geometries along the flowline (cyan colour in the inset) of Thwaites Glacier in Basin G-H, indicating the grounding line movement under different experiments: a linear Weertman (LW), b RC_eN_POC, c RC_iN_HAF, and d RC_eN_GlaDS_HAF. The coloured dots in the inset of (a) show the grounding line position in 2300 from LW (light grey), RC_iN_HAF (purple), RC_eN_POC (blue), and RC_eN_GlaDS_HAF (red).</figcaption></figure>

The primary difference between implicit and explicit representations of effective pressure lies in the location where the shift from the Weertman regime (slower inland flow) to the Coulomb regime (faster outlet flow) occurs, which determines the extent of the low-drag region or ice plain. The Coulomb limit defines the maximum basal shear stress under a Coulomb sliding regime, where sliding is primarily controlled by effective pressure and the sliding coefficient. In simulations using implicit representations of effective pressure (RC_iN and RC_iN_HAF), fast-flowing regions where basal sliding velocity exceeds $300 \text{ m yr}^{-1}$ exhibit Coulomb-plastic behaviour (see Methods for more details), effectively marking the Coulomb limit. Conversely, simulations with explicit effective pressure representation determine the Coulomb limit via the term $C''N''A_s$, which varies spatially and temporally with effective pressure. As a result, experiments with implicit effective pressure produce much larger Coulomb-dominated regions compared to those with explicit effective pressure, where the Coulomb regime is confined to less than 10 km upstream of the moving grounding line (Fig. S12). This distinction explains why RC_iN_HAF exhibits greater sensitivity to HAF-scaling compared to simulations that incorporate effective pressure from a hydrology model (RC_eN_GlaDS_HAF).

## Discussion

Our findings emphasize the critical role of subglacial water pressure and how a relatively small region close to the grounding line drives multi-century ice sheet mass loss. This aligns with the findings by ref. 11, which demonstrate that decreased effective pressure near the grounding line amplifies the ice sheet's sensitivity to climatic forcing, particularly for a given sliding relation power. Approaching the Coulomb limit and zero friction near the grounding line allows a more mobile grounding line with substantial mass loss potential. When this transition happens over a narrow area, model resolution becomes important, as ice-sheet models often struggle with small-scale transitions from high to zero drag near the grounding line³¹. We also find a strong sensitivity of ice sheet projections to assumptions about the spatial distribution of effective pressure, particularly to threshold parameters that govern the Weertman to Coulomb regime shift.

Previous studies²¹,⁴⁹ show that model projections based on LW cannot reproduce observations as well as RC that incorporate

sub-hydrology changes as part of the parameterisation. Among these RC experiments, we place greater confidence in the results from experiment RC_eN_GlaDS_HAF, as it uses a more physically realistic subglacial hydrology system where HAF-scaling adjusts effective pressure as the grounding line retreats. However, while the HAF-scaling accounts for some of these time-varying effects, the entire hydraulic regime may have shifted, particularly in regions experiencing significant changes in ice flow or grounding line dynamics¹²,³⁹. Moreover, the GlaDS model, originally designed to simulate subglacial hydrology on hard bedrock, may predict more efficient drainage systems than those occurring on a soft bed, potentially leading to lower modelled subglacial water pressures¹²,³⁷,³⁹. This highlights the importance of accurately predicting subglacial water pressure and developing a more comprehensive ice sheet-hydrology model that can account for the temporal coevolution of ice sheet dynamics and the subglacial hydrological system.

Our simulations explore the sensitivity to basal sliding parameterisations under various climate scenarios. However, limitations in our approach highlight areas for improvement. The use of SSA assumes vertical hydrostatic equilibrium and neglects vertical shear stresses, which restricts its ability to fully capture the ice dynamics at the grounding line. This is particularly relevant in regions with pronounced bedrock slopes, where non-hydrostatic effects and deviatoric (bridging) stresses-ignored by SSA become significant. While a fully resolved full-Stokes model, combined with a contact problem formulation for the lower ice surface, would be needed to consistently resolve these dynamics⁵⁰, we argue that the sliding relation itself is less sensitive to the choice of stress representation in the ice-flow model. For instance, a similar sensitivity to driving stress has been found between an SSA model and one incorporating simple vertical shear, particularly under a regularized Coulomb or Weertman sliding parameterisation²¹.

Nevertheless, unresolved processes may further alter projected sea-level contributions. Given the importance of the grounding line region to our results, we expect further changes when considering three-dimensional shear fields, hydrology-ocean-bedrock-sediment interactions, iceberg calving, and uncertainties in predicted future surface mass balance. While these processes are unlikely to negate the effects we identify-given that effective pressure is largely isolated from external forcings-they underscore the complexity of grounding line
dynamics and the need for urgent model development to advance the fidelity of sea-level projections.

This study underscores the necessity of incorporating evolving hydrology systems in ice sheet models, showing that assumptions made when attempting to integrate realistic basal processes strongly influence sea-level rise projections. The interaction between ice dynamics and the hydrological system, as parameterised with HAF-scaling, leads to a much more mobile grounding line compared to a non-evolving hydrological state, even if produced with an hydrology model. While flotation at the grounding line with HAF-scaling is physically justified, the threshold used to determine the scaling area is a crude approximation that overlooks the complex role of the hydrological system, which influences ice dynamics based on the water volume it carries and thus varies between outlet glaciers and ice shelves.

# Methods

## Model setup

Our simulations do not capture a fully coupled ice-hydrology system, as such a model is not yet available for continental-scale applications. Instead, our simulations aim to determine whether a coupling is indeed necessary, rather than predetermining what might be considered a valid future projection of ice-mass change. We report mass loss estimates from drainage basins in both East and West Antarctica, showing how the parameterisation of the effective pressure exerts control on tipping points characterized by sharp increases in grounding line flux and rapid grounding line retreat in deep basins with retrograde bed slopes. The sea-level contributions in this study are computed based on changes in total ice volume above flotation by assuming zero or very small changes in the bedrock elevation due to isostatic adjustment and a constant ocean area of typically $3.625 \times 10^{14} m^{251}$.

All experiments are conducted using the Shallow Shelf Approximation (SSA$^{52}$) implementation of Elmer/Ice$^{53}$, an open-source finite-element ice dynamic model, to solve the vertically integrated ice dynamics. Upper and lower ice surface elevation and bedrock elevation are from the MEaSURES BedMachine Antarctica, Version 3$^{54}$. An unstructured finite-element mesh is created initially using Gmsh$^{55}$ for the Antarctic domain following coastline positions provided by Bed-Machine V3, and refined using MMG$^{56}$ based on ice thickness gradients and the Hessian matrix of the observed velocity fields provided by MEaSURES Phase-Based Antarctica Ice Velocity map (Version 1)$^{57}$. The low-computation cost of SSA model allows a high-resolution mesh of 1 km in fast-flowing regions and near the grounding line, with a coarser resolution (25 km) further inland. This represents the highest resolution used in continental Antarctic ice sheet models from ISMIP6 studies$^{2,5}$. Sensitivity analysis using the LW sliding relation (Fig. S7) shows consistency between model outputs generated with 1 km mesh and 500 m mesh. Given the projected ice-mass loss using the Weertman sliding relation is more grid-size sensitive than the Coulomb relation used here$^{28,58}$, the mesh resolution artefact in this study is acceptable.

An inverse method is used to simulate both the basal drag coefficient $\beta$ and enhancement factor $\eta$ simultaneously through reducing the mismatch between the simulated and observed velocity. The inversion is performed with respect to a linear sliding coefficient:

$$ \tau_b = C_w \mathbf{u}_b = 10^\beta \mathbf{u}_b \quad (1) $$

where $u_b$ is the basal sliding velocity, $\beta$ is the sliding coefficient. The viscosity is defined by:

$$ \eta = E_\eta^2 \eta_0 \quad (2) $$

where $\eta_0$ is a reference viscosity field calculated based on a prescribed and vertically integrated 3-D volume temperature distribution$^{59}$. Three regularisation terms are introduced in the cost function to improve the conditioning of the inverse problem and balance between fitting with

observations and smoothness in $\beta$ and $\eta$. An 'L-surface' analysis$^{58}$ is conducted to determine the optimal regularisation parameters ($\lambda_\beta = 1e5$, $\lambda_{E_\eta} = 1e6$ and $\lambda_{E_\eta} = 0.02$). With the optimal parameters, we simulate the basal drag (Fig. S3a) and ice viscosity (Fig. S3b). The inverted basal drag is largely influenced by the bed geometry through affecting the ice flow resistance. The RMSEs between the simulated initial state in 1995 and observed ice thickness and ice surface velocity are shown in Fig. S8, which are among the low range of the RMSEs for the participating ice-sheet models in ISMIP6-2300 projections$^5$, except for the RMSE of velocity from RC\_iN and RC\_iN\_HAF.

Taking the initial geometry and inverted basal drag and viscosity as the starting point at 1995, we conduct a historical run from 1995 to 2014 to relax the model as an initial adjustment. A projection run from 2015 to 2300 is followed with climate forcing derived from selected climate models under the ISMIP6-2300 protocol$^{60}$. In the transient simulations, the surface mass balance (SMB) is given by the combination of a reference SMB and SMB anomalies, in which the reference SMB is the average mean SMB over 1995-2014 from the MAR product$^{61}$. The basal melt rate underneath the ice shelf is parameterised following the simplified ISMIP6 standard basal melt parameterisation based on a prescribed relation between ocean thermal forcing and ice shelf melting rates$^2$. Both SMB anomalies and thermal forcing are provided through the forcing provided by the ISMIP6 community$^5$. We adopt the forcing using NorESM1-M/RCP8.5 model for the historical run and UKESM/ssp5-85 model for the projection run. Forcing data from other four climate models under high emission scenario (CCSM4/RCP8.5, HadGEM2/RCP8.5, CESM2/ssp5-85, NorESM1-M/RCP8.5) and two model under low emission scenario (UKESM/ssp1-26, NorESM1-M/RCP2.6) are applied in the projection run to explore the sensitivity to various climate forcing. For partially floating elements, the representation of basal friction follows the sub-element parameterisation 3 (SEP3) with 20 integration points per element$^{60}$, while the basal melt only is applied on fully floating cells recommended by ref. 28,58. The ice front is fixed with a minimum ice thickness of 40 m applied, which may result in a vast and thin future Ice Shelf and could potentially lead to an overestimated buttressing effect on the grounded ice, though the effect might be small.

## Experiment design

Two different sliding relations, linear Weertman and RC, are chosen in the transient simulations. For the RC relation, the friction coefficients are converted based on the inverted basal shear stress, $\tau_b = 10^\beta \mathbf{u}_b$, to ensure that all experiments begin with the same initial basal shear stress (see details below). Importantly, the diverging transient behaviours arise from differences in how sliding evolves under the respective assumptions. To address uncertainties associated with the effective pressure, two formats of the RC sliding equations are considered (Table 1): one with explicit representation of effective pressure (hereafter 'RC\_eN') and another with implicit representation of effective pressure through a spatially tuned parameter (hereafter 'RC\_iN'). The transition between the Weertman and Coulomb-plastic behaviours differs in each format, as detailed below. Notably, in the Weertman regime, basal sliding is independent of effective pressure, making the transition critical for determining the extent of areas sensitive to spatial and temporal variations in effective pressure.

## Implicit representation of effective pressure

* Experiment RC\_iN and RC\_iN\_HAF

Here we follow ref. 21's version of RC sliding relation and assumed the roles of effective pressure, N, into the friction coefficient $C_1$:

$$ \tau_b = C_1 \left( \frac{\mathbf{u}_b}{\mathbf{u}_b + u_0} \right)^{\frac{1}{m}} \quad (3) $$
where m is a flow exponent (set to 3). C₁ is converted based on Eq (1) and Eq (3):

$$C_1 = 10^\beta \mathbf{u}_b \left( \frac{\mathbf{u}_b}{\mathbf{u}_b + u_0} \right)^{-\frac{1}{m}} \quad (4)$$

Here we use $u_0 = 300 \text{ m yr}^{-1}$ as the velocity threshold above which a Coulomb-plastic behaviour will occur. The model performance is proved to be relatively insensitive to $u_0$ values$^{21,49}$. The Coulomb limit occurs when:

$$\mathbf{u}_b \gg u_0 \quad (5)$$

Then Eq (3) reduces to:

$$\tau_b \approx C_1 \quad (6)$$

In RC\_iN, we assume N, and therefore C₁ obtained from inversion, do not change with time. However, this assumption may hold well for the interior of ice sheet but will break down in the fast flowing regions where the ice approaches flotation. To take this effect into account, a reduction on the basal shear stress $\tau_b$ is applied in experiment RC\_iN\_HAF so that $\tau_b$ can decrease smoothly towards zero at the grounding line, as shown below:

$$\tau_b = \lambda C_2 \left( \frac{\mathbf{u}_b}{\mathbf{u}_b + u_0} \right)^{\frac{1}{m}} \quad (7)$$

$$\lambda = \begin{cases} 1 & \text{if } h_{af} \ge h_T \\ \frac{h_{af}}{h_T} & \text{if } h_{af} < h_T \end{cases} \quad (8)$$

where $h_{af}$ is the height-above flotation, $\lambda$ is the HAF-scaling factor, $h_T$ is the threshold. In such case, the basal shear stress near the grounding line will evolve over time with ice geometry. To choose a reasonable value of $h_T$, we calculate $\lambda C_2$ for each element in our model:

$$\lambda C_2 = \tau_b \left( \frac{\mathbf{u}_b}{\mathbf{u}_b + u_0} \right)^{-\frac{1}{m}} \quad (9)$$

where $\tau_b$ and $u_b$ are simulated from inversion. The correlation between $\lambda C_2$ and $h_{af}$ (Fig. S9) indicates that $\lambda C_2$ reaches a relatively constant value of 0.21 Mpa where $h_{af} = 75 \text{ m}$. When $h_{af} < 75 \text{ m}$, a roughly linear trend is found, which indicates the linear decrease in basal drag when $h_{af}$ approaches zero. So we set $h_T$ to 75 m here and a sensitivity analysis on its values is discussed. In such case, C₂ is converted based on Eq (1) and Eq (8):

$$C_2 = 10^\beta \mathbf{u}_b \left( \frac{\mathbf{u}_b}{\mathbf{u}_b + u_0} \right)^{\frac{1}{m}} \lambda^{-1} \quad (10)$$

### Explicit representation of effective pressure

The implemented Coulomb sliding relation can be written as$^{19}$:

$$\tau_b = C_3 N \left( \frac{\chi \mathbf{u}_b^{-m}}{1 + a \chi^q} \right)^{\frac{1}{m}} \mathbf{u}_b \quad (11)$$

$$a = \frac{(q-1)^{q-1}}{q^q} \quad (12)$$

$$\chi = \frac{\mathbf{u}_b}{C_3 N^m A_s} \quad (13)$$

$q$ is a 'post-peak' exponent (set to 1), $C_3$ and $A_s$ are two different friction coefficients that affect the Coulomb and Weertman behaviour correspondingly. This equation can be rearranged to:

$$\tau_b = C_3 N \left( \frac{\mathbf{u}_b}{\mathbf{u}_b + C_3^m N^m A_s} \right)^{\frac{1}{m}} \quad (14)$$

The Coulomb limit occurs when the basal sliding resistance $\tau_b$ is primarily determined by $C_3N$, with minimal influence from the basal sliding velocity. This happens in regions where the sliding relatin transitions from Weertman regime to Coulomb regime when:

$$\mathbf{u}_b \gg C_3^m N^m A_s \quad (15)$$

Under this condition, Eq (14) reduces to:

$$\tau_b \approx C_3 N \quad (16)$$

This is classic Coulomb friction law, where basal shear stress $\tau_b$ is directly proportional to the effective pressure N and the Coulomb coefficient $C_3$. We now need to convert the linear Weertman coefficient $C_W = 10^\beta$ to the unknown Coulomb coefficients $C_3$ and $A_s$. The sliding parameter $A_s$ determines the ease of basal sliding, with higher values representing more resistance to sliding, which is influenced by both physical properties of the ice-bed interface and the effective pressure$^{18}$. Here we interpret the inverted linear friction coefficient $\beta$ everywhere in terms of a non-linear Weertman coefficient, using the following relation:

$$A_s = |\mathbf{u}_b|^{1-m} 10^{-m\beta} \tanh(2N) \quad (17)$$

Eq (14) is rearranged in terms of the Coulomb C parameter:

$$C_3 = |\mathbf{u}_b| 10^\beta N^{-1} (1 - 10^{m\beta} |\mathbf{u}_b|^{m-1} A_s)^{-1/m} \quad (18)$$

* Experiment RC\_eN\_POC

In this experiment, we assume a perfect hydrological connection between the subglacial drainage system and the ocean so that N can be calculated based on ice thickness and bedrock elevation:

$$N = \rho_i g H + \rho_o g z_b \quad (19)$$

where H is ice thickness, g is acceleration due to gravity, $\rho_i$ is ice density, $\rho_o$ is ocean water density, $z_b$ is bedrock elevation. It allows N to evolve with ice geometry under the assumption that subglacial water pressure equals the ocean pressure at the grounding line, which is more realistic near the grounding line when ice approaches flotation.

* Experiments RC\_eN\_GlaDS and RC\_eN\_GlaDS\_HAF

In experiments RC\_eN\_GlaDS and RC\_eN\_GlaDS\_HAF, the basal drag is decided following Eq. (14) and N is simulated with the Glacier Drainage System model (GlaDS$^{37}$), implemented in Elmer/Ice as a separate module$^{62}$. The GlaDS model was selected for this study as it combines a distributed water sheet and channelized drainage system, providing a good representation of subglacial hydrology that is particularly suited to investigating how bed topography influences effective pressure and basal drag$^{9,35,37,62}$. It is fully implemented within Elmer/Ice, allowing for seamless spin-up of the hydrology model with a restart from the ice sheet model. GlaDS is run into a steady state using the initial ice geometry to calculate the spatial distribution of effective pressure at the start of the simulation$^{63}$. In this setup, frictional heating (calculated using the inversion simulations described above) is assumed to be the sole source of subglacial meltwater production, as it
dominates in areas of significant basal sliding, where high ice velocities generate far more heat than the relatively small and more uniform geothermal flux⁶⁴. To ensure that the hydrology system reaches a steady state, the total subglacial water budget is held constant throughout the spin-up period (over -100 model years). This assumption introduces negligible variability in the final results compared to alternative assumptions on subglacial melt supply.

Bedrock geometry and ice thickness strongly influences subglacial water flow and channel distribution, as lower-elevation areas, such as valleys, accumulate water due to hydraulic potential gradients shaped by bedrock topography and ice thickness, resulting in higher water pressures and lower effective pressures. In RC eN GlaDS, we assume the initially computed effective pressure to remain constant with time. To account for a temporally evolving effective pressure near the grounding line, the previously described HAF-scaling is applied on effective pressure in RC eN GlaDS HAF when $h_{af} < h_T$, defining it as N.

$$ \tau_b = C_4 N^* \left( \frac{\mathbf{u}_b}{\mathbf{u}_b + C_3^m N^*} \right)^{\frac{1}{m}} \quad (20) $$

$$ N^* = \lambda N \quad (21) $$

$$ \lambda = \begin{cases} 1 & \text{if } h_{af} \ge h_T \\ \frac{h_{af}}{h_T} & \text{if } h_{af} < h_T \end{cases} \quad (22) $$

$$ A_s = |\mathbf{u}_b|^{1-m} 10^{-m\beta} \tanh(2 \cdot N^*) \quad (23) $$

$C_4$ is converted based on Eq (1):

$$ C_4 = |\mathbf{u}_b| 10^{\beta} N^{-1} (1 - 10^{m\beta} |\mathbf{u}_b|^{m-1} A_s)^{-1/m} \quad (24) $$

## Data availability

Original forcings datasets from various climate models are available on Ghub https://theghub.org/dataset-listingunder 'ISMIP6 Antarctica 2300'. The MEaSURES BedMachine Antarctica, Version 3 is available at https://nsidc.org/data/nsidc-0756/versions/3. The MEaSURES Phase-Based Antarctica Ice Velocity map (Version 1) is available at https://nsidc.org/data/nsidc-0754/versions/1. The model configurations and all data and scripts used to generate the figures are available at https://doi.org/10.5281/zenodo.14874036.

## Code availability

The most recent developments of the Elmer/Ice code are downloadable from https://github.com/ElmerCSC/elmerfem(last visited 2024-08-10). The most recent release of Elmer/Ice is available from https://doi.org/10.5281/zenodo.7892181.

## References

1. Edwards, T. L. et al. Projected land ice contributions to twenty-first-century sea level rise. Nature 593, 74–82 (2021).

2. Seroussi, H. et al. ISMIP6 Antarctica: a multi-model ensemble of the Antarctic ice sheet evolution over the 21st century. Cryosphere 14, 3033–3070 (2020).

3. Fox-Kemper, B. et al. Ocean, cryosphere and sea level change. In Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change (eds. Masson-Delmotte, V. et al.) Ch. 9, 1211–1361. https://www.ipcc.ch/report/ar6/wg1/downloads/report/IPCC_AR6_WGI_Chapter09.pdf. (Cambridge University Press, 2021)

4. Seroussi, H. et al. Insights into the vulnerability of Antarctic glaciers from the ISMIP6 ice sheet model ensemble and associated uncertainty. Cryosphere 17, 5197–5217 (2023).

5. Seroussi, H. et al. Evolution of the Antarctic ice sheet over the next three centuries from an ISMIP6 model ensemble. Earth's Future 12, e2024EF004561 (2024).

6. Park, J.-Y. et al. Future sea-level projections with a coupled atmosphere-ocean-ice-sheet model. Nat. Commun. 14, 636 (2023).

7. Dow, C. F. The role of subglacial hydrology in Antarctic ice sheet dynamics and stability: a modelling perspective. Ann. Glaciol. 63, 49–54 (2022).

8. Bradley, A. T. & Hewitt, I. J. Tipping point in ice-sheet grounding-zone melting due to ocean water intrusion. Nat. Geosci. 17, 631–637 (2024).

9. Pelle, T., Greenbaum, J. S., Ehrenfeucht, S., Dow, C. F. & McCormack, F. S. Subglacial discharge accelerates dynamic retreat of aurora subglacial basin outlet glaciers, East Antarctica, over the 21st century. J. Geophys. Res. Earth Surf. 129, e2023JF007513 (2024).

10. Goeller, S., Thoma, M., Grosfeld, K. & Miller, H. A balanced water layer concept for subglacial hydrology in large-scale ice sheet models. Cryosphere 7, 1095–1106 (2013).

11. Kazmierczak, E., Sun, S., Coulon, V. & Pattyn, F. Subglacial hydrology modulates basal sliding response of the Antarctic ice sheet to climate forcing. Cryosphere Discuss. 2022, 1–24 (2022).

12. Schoof, C. Ice-sheet acceleration driven by melt supply variability. Nature 468, 803–806 (2010).

13. Le Brocq, A. M. et al. Evidence from ice shelves for channelized meltwater flow beneath the Antarctic Ice Sheet. Nat. Geosci. 6, 945–948 (2013).

14. Hager, A. O., Hoffman, M. J., Price, S. F. & Schroeder, D. M. Persistent, extensive channelized drainage modeled beneath Thwaites Glacier, West Antarctica. Cryosphere 16, 3575–3599 (2022).

15. Alley, K. E., Scambos, T. A., Siegfried, M. R. & Fricker, H. A. Impacts of warm water on Antarctic ice shelf stability through basal channel formation. Nat. Geosci. 9, 290–293 (2016).

16. Drews, R. et al. Actively evolving subglacial conduits and eskers initiate ice shelf channels at an Antarctic grounding line. Nat. Commun. 8, 15228 (2017).

17. Weertman, J. On the sliding of glaciers. J. Glaciol. 3, 33–38 (1957).

18. Schoof, C. The effect of cavitation on glacier sliding. Proc. R. Soc. A Math. Phys. Eng. Sci. 461, 609–627 (2005).

19. Gagliardini, O., Cohen, D., Ráback, P. & Zwinger, T. Finite-element modeling of subglacial cavities and related friction law. J. Geophys. Res. Earth Surface 112. https://doi.org/10.1029/2006JF000576 (2007).

20. Zoet, L. K. & Iverson, N. R. A slip law for glaciers on deformable beds. Science 368, 76–78 (2020).

21. Joughin, I., Smith, B. E. & Schoof, C. G. Regularized Coulomb friction laws for ice sheet sliding: application to Pine Island Glacier, Antarctica. Geophys. Res. Lett. 46, 4764–4771 (2019).

22. Gillet-Chaulet, F. et al. Assimilation of surface velocities acquired between 1996 and 2010 to constrain the form of the basal friction law under Pine Island Glacier. Geophys. Res. Lett. 43, 310–311 (2016).

23. Minchew, B. & Joughin, I. Toward a universal glacier slip law. Science 368, 29–30 (2020).

24. DeConto, R. M. et al. The Paris Climate Agreement and future sea-level rise from Antarctica. Nature 593, 83–89 (2021).

25. Reed, B., Green, J. A. M., Jenkins, A. & Gudmundsson, G. H. Recent irreversible retreat phase of Pine Island Glacier. Nat. Clim. Change 14, 75–81 (2024).

26. Dawson, E. J., Schroeder, D. M., Chu, W., Mantelli, E. & Seroussi, H. Ice mass loss sensitivity to the Antarctic ice sheet basal thermal state. Nat. Commun. 13, 4957 (2022).
27. Brondex, J., Gillet-Chaulet, F. & Gagliardini, O. Sensitivity of centennial mass loss projections of the Amundsen basin to the friction law. *Cryosphere* **13**, 177–195 (2019).

28. Seroussi, H. & Morlighem, M. Representation of basal melting at the grounding line in ice flow models. *Cryosphere* **12**, 3085–3096 (2018).

29. Asay-Davis, X. S. et al. Experimental design for three interrelated marine ice sheet and ocean model intercomparison projects: MISMIP v. 3 (MISMIP +), ISOMIP v. 2 (ISOMIP +) and MISOMIP v. 1 (MISOMIP1). *Geosci. Model Dev.* **9**, 2471–2497 (2016).

30. Sun, S. et al. Antarctic ice sheet response to sudden and sustained ice-shelf collapse (ABUMIP). *J. Glaciol.* **66**, 891–904 (2020).

31. Gladstone, R. M. et al. Marine ice sheet model performance depends on basal sliding physics and sub-shelf melting. *Cryosphere* **11**, 319–329 (2017).

32. Ma, Y. et al. Sensitivity of Totten Glacier dynamics to sliding parameterizations and ice shelf basal melt rates. *EGUsphere* **2024**, 1–27 (2024).

33. Shreve, R. L. Movement of water in Glaciers. *J. Glaciol.* **11**, 205–214 (1972).

34. Fricker, H. A., Scambos, T., Bindschadler, R. & Padman, L. An active subglacial water system in West Antarctica mapped from space. *Science* **315**, 1544–1548 (2007).

35. Dow, C. F., Ross, N., Jeofry, H., Siu, K. & Siegert, M. J. Antarctic basal environment shaped by high-pressure flow through a subglacial river system. *Nat. Geosci.* **15**, 892–898 (2022).

36. Siegfried, M. R. et al. The life and death of a subglacial lake in West Antarctica. *Geology* **51**, 434–438 (2023).

37. Werder, M. A., Hewitt, I. J., Schoof, C. G. & Flowers, G. E. Modeling channelized and distributed subglacial drainage in two dimensions. *J. Geophys. Res.: Earth Surf.* **118**, 2140–2158 (2013).

38. Sommers, A., Rajaram, H. & Morlighem, M. SHAKTI: subglacial hydrology and kinetic, transient interactions v1.0. *Geosci. Model Dev.* **11**, 2955–2974 (2018).

39. Flowers, G. E. Modelling water flow under glaciers and ice sheets. *Proc. R. Soc. A Math. Phys. Eng. Sci.* **471**, 20140907 (2015).

40. Choi, Y., Seroussi, H., Gardner, A. & Schlegel, N.-J. Uncovering basal friction in Northwest Greenland using an ice flow model and observations of the past decade. *J. Geophys. Res. Earth Surf.* **127**, e2022JF006710 (2022).

41. Rignot, E. et al. Widespread seawater intrusions beneath the grounded ice of Thwaites Glacier, West Antarctica. *Proc. Natl. Acad. Sci. USA* **121**, e2404766121 (2024).

42. Damsgaard, A. et al. Ice flow dynamics forced by water pressure variations in subglacial granular beds. *Geophys. Res. Lett.* **43**, 112–165 (2016).

43. Rosier, S. H. R., Gudmundsson, G. H. & Green, J. A. M. Temporal variations in the flow of a large Antarctic ice stream controlled by tidally induced changes in the subglacial water system. *Cryosphere* **9**, 1649–1661 (2015).

44. Otosaka, I. N. et al. Mass balance of the Greenland and Antarctic ice sheets from 1992 to 2020. *Earth Syst. Sci. Data* **15**, 1597–1616 (2023).

45. Schoof, C. Ice sheet grounding line dynamics: steady states, stability, and hysteresis. *J. Geophys. Res. Earth Surf.* **112**. https://doi.org/10.1029/2006JF000664 (2007).

46. Naughten, K. A., Holland, P. R. & De Rydt, J. Unavoidable future increase in West Antarctic ice-shelf melting over the twenty-first century. *Nat. Clim. Change* **13**, 1222–1228 (2023).

47. Bindschadler, R. Siple Coast Project research of Crary Ice Rise and the mouths of Ice Streams B and C, West Antarctica: review and new perspectives. *J. Glaciol.* **39**, 538–552 (1993).

48. Joughin, I., Shapero, D. & Dutrieux, P. Responses of the Pine Island and Thwaites glaciers to melt and sliding parameterizations. *Cryosphere* **18**, 2583–2601 (2024).

49. Jager, E., Gillet-Chaulet, F., Mouginot, J. & Millan, R. Validating ensemble historical simulations of Upernavik Isstrøm (1985-2019)

using observations of surface velocity and elevation. *J. Glaciol.* **70**, 1–18. https://doi.org/10.1017/jog.2024.10 (2024).

50. Zhang, T. et al. A comparison of two Stokes ice sheet models applied to the Marine Ice Sheet Model Intercomparison Project for plan view models (MISMIP3d). *Cryosphere* **11**, 179–190 (2017).

51. Gregory, J. M. et al. Concepts and terminology for sea level: mean, variability and change, both local and global. *Surv. Geophys.* **40**, 1251–1289 (2019).

52. MacAyeal, D. R. Large-scale ice flow over a viscous basal sediment: theory and application to ice stream B, Antarctica. *J. Geophys. Res.: Solid Earth* **94**, 4071–4087 (1989).

53. Gagliardini, O. et al. Capabilities and performance of Elmer/Ice, a new-generation ice sheet model. *Geosci. Model Dev.* **6**, 1299–1318 (2013).

54. Morlighem, M. et al. Deep glacial troughs and stabilizing ridges unveiled beneath the margins of the Antarctic ice sheet. *Nat. Geosci.* **13**, 132–137 (2020).

55. Geuzaine, C. & Remacle, J.-F. Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities. *Int. J. Numer. Methods Eng.* **79**, 1309–1331 (2009).

56. Dobrzynski, C. & Frey, P. Anisotropic delaunay mesh adaptation for unsteady simulations. In *Proc. 17th International Meshing Roundtable* (ed. Garimella, R. V.) 177–194 (Springer Berlin Heidelberg, 2008).

57. Mouginot, J., Rignot, E. & Scheuchl, B. Continent-wide, interferometric SAR phase, mapping of Antarctic ice velocity. *Geophys. Res. Lett.* **46**, 9710–9718 (2019).

58. Wang, Y. et al. Sensitivity of the future evolution of the Wilkes Subglacial Basin ice sheet to grounding-line melt parameterizations. *Cryosphere* **18**, 5117–5137 (2024).

59. Greve, R. et al. Future projections for the Antarctic ice sheet until the year 2300 with a climate-index method. *J. Glaciol.* 1-11 https://www.cambridge.org/core/product/identifier/S0022143023000412/type/journal_article (2023).

60. Seroussi, H., Morlighem, M., Larour, E., Rignot, E. & Khazendar, A. Hydrostatic grounding line parameterization in ice sheet models. *Cryosphere* **8**, 2075–2087 (2014).

61. Agosta, C. et al. Estimation of the Antarctic surface mass balance using the regional climate model MAR (1979-2015) and identification of dominant processes. *Cryosphere* **13**, 281–296 (2019).

62. GAGLIARDINI, O. & WERDER, M. A. Influence of increasing surface melt over decadal timescales on land-terminating Greenland-type outlet glaciers. *J. Glaciol.* **64**, 700–710 (2018).

63. Zhang, Y. et al. The role of hydraulic conductivity in the Pine Island Glacier's subglacial water distribution. *Sci. Total Environ.* **927**, 172144 (2024).

64. Sergienko, O. V. & Hindmarsh, R. C. A. Regular patterns in frictional resistance of ice-stream beds seen by surface data inversion. *Science* **342**, 1086–1089 (2013).

65. Rignot, E. et al. Four decades of Antarctic Ice Sheet mass balance from 1979-2017. *Proc. Natl. Acad. Sci. USA* **116**, 1095–1103 (2019).

## Acknowledgements

C.Z. is the recipient of an Australian Research Council Discovery Early Career Researcher Award (project number DE240100267) funded by the Australian Government. C.Z., Y.W., P.C., and B.K.G. received grant funding from the Australian Government as part of the Antarctic Science Collaboration Initiative program (ASCI000002). R.G. was supported by the Research Council of Finland (grant numbers 322430 and 355572), and by the Finnish Ministry of Education and Culture and CSC-IT Center for Science (Decision diary number OKM/10/524/2022). T.Z. was supported by the Research Council of Finland (grant number 322978). E.J. was supported by the Research Council of Finland (grant number 355783). P.C., B.K.G. and M.A.K. are supported by the Australian Research Council Special Research Initiative, Australian Centre for
Excellence in Antarctic Science (project number SR200100008). J.C. and F.G.-C. received funding from the French National Research Agency (ANR) under grants ANR-19-CE01-0015 (EIS). P.M. received funding from Agence Nationale de la Recherche—France 2030 as part of the PEPR TRACCS programme under grant number ANR-22-EXTR-0010. All simulations were enabled by computational resources provided by CSC-IT Center for Science Ltd.

**Author contributions**

C.Z. led the research, conducted the simulations, analysed the results, led the writing. R.G. conducted the GlaDS simulation. R.G. and T.Z. contributed to the experiment design. R.G., T.Z., F.G.-C., Y.W., J.C., P.M., L.S. and E.J. contributed to the model setup and analysis. R.G., T.Z., F.G., Y.W., J.C., P.M., L.S., E.J., B.K.G., P.C. and M.A.K. contributed to the manuscript writing.

**Competing interests**

The authors declare no competing interests.

**Additional information**

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41467-025-58375-4.

**Correspondence** and requests for materials should be addressed to Chen Zhao.

**Peer review information** *Nature Communications* thanks Sebastian Beyer, Duncan Young who coiewed with Shuai Yan and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. A peer review file is available.

**Reprints and permissions information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© Crown 2025