Tipping point in ice-sheet grounding-zone melting due to ocean water intrusion

Marine ice sheets are highly sensitive to submarine melting in their grounding zones, where they transition between grounded and floating ice. Recently published studies of the complex hydrography of grounding zones suggest that warm ocean water can intrude large distances beneath the ice sheet, with dramatic consequences for ice dynamics. Here we develop a model to capture the feedback between intruded ocean water, the melting it induces and the resulting changes in ice geometry. We reveal a sensitive dependence of the grounding-zone dynamics on this feedback: as the grounding zone widens in response to melting, both temperature and flow velocity in the region increase, further enhancing melting. We find that increases in ocean temperature can lead to a tipping point being passed, beyond which ocean water intrudes in an unbounded manner beneath the ice sheet, via a process of runaway melting. Additionally, this tipping point may not be easily detected with early warning indicators. Although completely unbounded intrusions are not expected in practice, this suggests a mechanism for dramatic changes in grounding-zone behaviour, which are not currently included in ice-sheet models. We consider the susceptibility of present-day Antarctic grounding zones to this process, finding that both warm and cold water cavity ice shelves may be vulnerable. Our results point towards a stronger sensitivity of ice-sheet melting, and thus higher sea-level-rise contribution in a warming climate, than has been previously understood.

Marine ice sheets are highly sensitive to submarine melting in their grounding zones, where they transition between grounded and floating ice.Recently published studies of the complex hydrography of grounding zones suggest that warm ocean water can intrude large distances beneath the ice sheet, with dramatic consequences for ice dynamics.Here we develop a model to capture the feedback between intruded ocean water, the melting it induces and the resulting changes in ice geometry.We reveal a sensitive dependence of the grounding-zone dynamics on this feedback: as the grounding zone widens in response to melting, both temperature and flow velocity in the region increase, further enhancing melting.We find that increases in ocean temperature can lead to a tipping point being passed, beyond which ocean water intrudes in an unbounded manner beneath the ice sheet, via a process of runaway melting.Additionally, this tipping point may not be easily detected with early warning indicators.Although completely unbounded intrusions are not expected in practice, this suggests a mechanism for dramatic changes in grounding-zone behaviour, which are not currently included in ice-sheet models.We consider the susceptibility of present-day Antarctic grounding zones to this process, finding that both warm and cold water cavity ice shelves may be vulnerable.Our results point towards a stronger sensitivity of ice-sheet melting, and thus higher sea-level-rise contribution in a warming climate, than has been previously understood.
There is growing evidence suggesting that ice-sheet models lack representation of important physical processes driving ice sheet retreat (for example, refs.1-7), rendering their projections of sea-level rise less sensitive to climatic changes than they should be.Point-wise observations 1,2 suggest that ice-shelf basal melt rates are considerably smaller than those typically required by models to reproduce observed retreat rates.Moreover, ice-sheet models systematically underestimate recent ice loss 7 and struggle to reproduce observationally constrained sea-level highstands from previous interglacials 3,4 .Palaeoclimate ice-sheet reconstructions have largely been able to reproduce low-end estimates only when mechanisms to boost sensitivity to climatic forcing are invoked [8][9][10] .Recently, evidence from diverse sources has emerged suggesting that relatively warm ocean water can intrude long distances upstream of ice-shelf grounding lines 5,6,[11][12][13][14] ; such long intrusions have dramatic consequences for sea-level-rise contributions from ice sheets 5,15 , making them a candidate mechanism to reconcile modelled and observed sea-level rise.Here we investigate how previously ignored feedbacks between melting and the confining ice geometry make this intrusion mechanism even more powerful.
Sea levels during previous interglacials can be considered analogues for future sea-level rise under anthropogenic influence 16 .During the Pliocene (∼3 million years before present), CO 2 levels were similar to present day 17 and temperatures 2-3 °C above pre-industrial levels 18 , but the global mean sea level (GMSL) was 6-40 m above present-day levels 4 .Such sea-level rise is only possible with a sizeable contribution from the Antarctic ice sheet 16 .Whereas proxy records suggest that such https://doi.org/10.1038/s41561-024-01465-7 These models predict that warm water can intrude significant distances (up to kilometres) upstream of ice-sheet grounding lines 5 , delivering excess heat for sub-ice melting.Ice-sheet models including the intrusion mechanism show that ice loss is highly sensitive to this intrusion length 5 , with kilometre length intrusions potentially doubling rates of ice loss.Grounding-zone melting via warm water intrusion has, therefore, been suggested as a physically based mechanism for explaining ice-sheet retreat during past warm periods 5 .
Crucially, however, existing models lack feedbacks between melting and the confining ice geometry.As grounding zones widen in response to melting, both primary factors controlling melt rates-the amount of warm water entering the grounding zone and flow speeds within it-will be affected.In particular, this feedback may strongly influence the distance warm water is able to intrude beneath ice sheets, with significant implications for ice dynamics.References 5,11 demonstrated that with no feedback between melting and geometry, the distance warm water can intrude depends on the slope of the ice base and can become infinite if the slope is sufficiently steep.Here we show that with the melt-geometry feedback included, the system exhibits a tipping-point behaviour as the parameters controlling the melt rate are varied, causing unbounded intrusions to develop even on a flat or down-sloping bed; this behaviour can be triggered by changes in external forcing, such as ocean temperatures.

Melting causes enhanced warm water intrusion
To understand how the melt feedback affects grounding-zone behaviour, and, in particular, the distance warm water is able to intrude, we coupled the layered intrusion model of refs.5,11 with a common melt-rate model accounting for the dependence of melting on both temperature and flow velocity adjacent to the ice (Methods).
The modelled grounding-zone behaviour depends on four fundamental, dimensionless parameters (Supplementary Information): Here U ∞ and H ∞ are a characteristic flow velocity and thickness of the upstream hydrological network (Fig. 1); V is the grounding-zone ice Antarctic ice loss occurred during previous interglacials (for example, refs.19-21), ice-sheet models struggle to reproduce the corresponding ice-sheet retreat (for example, ref. 22).Several palaeoclimate simulations 8,10,22 that attain low-end Pliocene GMSL estimates have gained interest because of their pessimistic Antarctic ice loss projections.To reconcile observational constraints, these models incorporate processes that increase their sensitivity to past climatic change, naturally rendering them more sensitive to future anthropogenic warming.Refs.8,10 achieve low-end Pliocene GMSL by introducing a cliff-collapse mechanism, whereby sufficiently tall ice cliffs collapse, prompting rapid inland ice front retreat.However, these simulations have been questioned on both physical 23,24 and statistical 25 grounds.Ref. 22 obtained similar Pliocene GMSL using a model whose increased sensitivity results from a parameterization of basal melting in grounding zones, where grounded ice transitions into a floating ice shelf (Fig. 1).Flow of grounded ice is particularly sensitive to grounding-zone melting 15,[26][27][28] , because melt-induced thinning there both reduces basal drag and provides a thinning perturbation that propagates through the shelf, reducing buttressing.
Specifically, the model in ref. 22 interpolates melting across model grid cells either side of the traditional 'grounding line', where the ice is at the hydrostatic floatation thickness.Although this respects the fact that areas between grid cells around the grounding line may be exposed to warm ocean water, it has little physical basis beyond; in practice, grounding zones are highly complex regions, with myriad features including local topographic highs (for example, ref. 29), porous till layers (for example, ref. 30) and channels on multiple lengthscales (for example, ref. 31) (Fig. 1).Freshwater is delivered to the ocean from the upstream grounded ice through this region; where freshwater meets the relatively warm, salty ocean water, the lower density freshwater rises, permitting the warm water to intrude upstream of the grounding line 11 .There is growing evidence of warm water intrusions from diverse sources including surface observations 12,13 , satellite data 6,14 and ice-shelf basal features 32 .
Recent near-grounding-zone observations 1,2 confirm that warm ocean water can reach cavity extremities.However, observations are unable to probe regions far upstream of the grounding line.A lack of direct observations beneath grounded ice sheets, combined with their importance in large ice-sheet models, means that models of the grounding zone are essential.Recently published grounding-zone models 5,11 treat the region as a porous, two-layer system with cold, fresh subglacial discharge overlaying warm, saline ocean water.Here V is the velocity of ice above the channel, H ∞ is the characteristic vertical length scale of the upstream subglacial network and θ is the local angle of the seabed.The two-dimensional model should be taken to represent an alonggrounding-zone average of the complex three-dimensional drainage system in the centre panel; in particular, there will be areas where the ice remains in contact with the bed (that is, the model does not assume the ice is floating everywhere). https://doi.org/10.1038/s41561-024-01465-7 velocity; St is the Stanton number, the ratio between the thermal flux into the ice-ocean interface and the thermal capacity of the water, which effectively parametrizes exchange across a boundary layer at the ice-ocean interface 33 (Supplementary Information); c is the specific heat capacity of water; c d is a cross-sectional average drag coefficient between the water and the channel; c is the specific heat capacity of ocean water; ΔT = T O + Γ  O − T D is the thermal forcing, with T O the ocean temperature,  O the ocean salinity, Γ the freezing point slope with salinity and T D the local freshwater freezing temperature; ℒ is the latent heat of fusion of seawater; θ is the local grounding-zone slope, assumed constant (Supplementary Information); g ′ = g( ρ 0 − ρ w )/ρ w is the reduced gravity, with ρ 0 a reference density and ρ w the local water density; and c i is a cross-sectional average drag coefficient between the two layers.A full model description, including a discussion of underlying assumptions, can be found in the Supplementary Information.These parameters capture the complex ice-ocean interactions that occur in grounding zones: M (equation (1a)) is a dimensionless melt rate, describing the competing effects of increasing ocean temperature (increasing ΔT) in promoting enhanced melting, and increased ice advection (increasing V) in replacing this ice; S (equation (1b)) is a dimensionless bedslope, with positive (negative, respectively) S corresponding to retrograde (prograde) bedslopes-upward (downward) sloping in the direction of ice flow; F (equation (1c)) is the upstream Froude number, which describes the upstream hydrological network efficiency: efficient networks, with fast flow (high U ∞ ) through narrow confines (low H ∞ ) correspond to large F, whereas in inefficient networks (low F), meltwater is transported slowly through wide channels; C (equation (1d)) is a dimensionless interfacial drag coefficient, describing the relative importance of drag between the two water layers and between the water and solid (ice/bed) boundaries.
Figure 2 shows how the intrusion distance, grounding-zone geometry, thermal driving and flow velocity change as the geometry evolves from an initially flat ice base and no seabed slope in two cases with similar ocean conditions: one with ΔT = 2.3 °C and the other with ΔT = 2.5°C (we consider the case of a variable bedslope later).Melting is concentrated at the channel entrance, where the flow velocity and thermal driving are highest (Fig. 2c,d,g,h) and reduces with distance into the channel.As melting proceeds, the grounding-zone widens (Fig. 2b,f), permitting more warm water to enter, increasing the average temperature (Fig. 2c,g), and reductions in drag result in higher flow velocities (Fig. 2d,h); these work in tandem to promote enhanced melting.When ice is replaced by advection sufficiently quickly (left panels of Fig. 2), the system reaches a steady state with melting balancing ice advection and drag balancing the gravitational force that results from a titled warm-cold interface.This equilibrium is reached on a timescale of days (Fig. 2a).However, for marginally higher ocean temperatures, ice advection cannot balance melting and the grounding zone continually widens (Fig. 2e,f).The feedback between geometry, hydrology and melting results in runaway warm water intrusion.This system displays a tipping-point-like behaviour: a small change in the ocean temperature (and thus parameter M) results in a threshold being passed, across which a dramatic change in the modelled final intrusion length L occurs, from being bounded (Fig. 2a) to being unbounded (Fig. 2e).The timescale on which the grounding zone responds to melting is much shorter than that on which the grounded ice thickness responds to perturbations in melting; it is therefore reasonable to consider the late time behaviour, and we henceforth refer to the final intrusion length L as the intrusion length.The large increase in L, and thus melting beneath a large section of the grounded ice sheet, would have dramatic implications for the dynamics of a marine-terminating ice sheet 5,15 .Note, however, that we do not expect intrusions to penetrate indefinitely in practice, because processes not included in our model will play a role on long lengthscales and potentially stabilize the intrusion (see below).

Grounding-zone melt as a generic tipping point
This tipping point is generic: for any hydrological network efficiency F, the intrusion length increases with the melt parameter M (Fig. 3) and there is a critical M above which the intrusion becomes unbounded (solid line in Fig. 3).Equivalently, for any M, there is a critical F, named F c , below which the intrusion becomes unbounded.Although less efficient networks (lower F) correspond to lower upstream flow velocities (lower U ∞ , see equations (1a,b) and (1c)), this is outweighed by reduced drag between the layers, resulting in higher flow speeds adjacent to the ice and a thicker warm water layer, both of which promote increased melting.The critical hydrological network efficiency F c is increasing with M (Fig. 3): higher melting is required to cross the tipping point for more efficient subglacial networks.This suggests that increases in the flow of water beneath ice sheets may act as a stabilizing control on their dynamics via reduced grounding-zone melting, contrasting the common belief that subglacial flow predominantly enhances ice loss via reduced basal friction 34 .
The location of the transition between bounded and unbounded intrusions is relatively insensitive to the dimensionless drag coefficient C (Fig. 3).This provides support for our use of a two-dimensional model, because C encodes heterogeneity in the along-grounding-zone direction.
It is interesting to note that when the intrusion is bounded, the intrusion length L is fairly insensitive to the melt parameter M (Fig. 3).This suggests that early warning indicators 35 , which might indicate that a marine ice sheet is approaching such a grounding-zone melt tipping point as (say) the ocean temperature increases, may be hard to detect: an increase in the intrusion length would not appear as a detectable signal in ice dynamical observations until the tipping point is passed, propagating uncertainty into sea-level-rise projections 36 .

Widespread susceptibility to the tipping point in melting
Intrusion of dense seawater is ultimately gravity driven and therefore strongly modified by the grounding-zone bedslope.Retrograde bedslopes (S > 0) result in enhanced intrusion, whereas prograde bedslopes (S < 0) result in reduced intrusion, compared to a flat bed (Fig. 4a and Supplementary Fig. 8).Unbounded intrusion can occur with no melt feedback (equivalent to M → 0), provided that the slope is sufficiently retrograde, under conditions described by ref. 5.However, with melt feedbacks, unbounded intrusion can occur for any bedslope, including prograde (Fig. 4a): when melting is sufficiently strong, the widening of the warm layer accompanying channel opening creates a gravitational driving force that overcomes the retarding gravitational effect from the bedslope, which is otherwise unfavourable for intrusion.The critical bedslope, S c , above which unbounded intrusion occurs, reduces with M (Fig. 4a).Over large regions of parameter space, S c is negative (Supplementary Fig. 7): that is, our model suggests that unbounded intrusion may occur even on prograde bedslopes and particularly so for inefficient hydrological networks (low F; Supplementary Fig. 7).
Marine ice sheets grounded on retrograde bedslopes may be susceptible to the marine ice-sheet instability (MISI) (for example, refs.37-39).Retreat of the West Antarctic ice sheet, many areas of which have grounding zones on retrograde bed slopes 40 , may be strongly controlled by MISI 41 .Our modelling suggests warm water intrusion is most likely on retrograde bedslopes, potentially enhancing MISI.Conversely, it is commonly assumed that prograde grounding lines are stable 38,42 ; our results suggest prograde grounding zones can also host substantial intrusions and have the possibility for a switch in behaviour as ocean temperatures change, potentially questioning their assumed stability.
In practice, hydrological networks beneath ice sheets are poorly constrained, making it difficult to determine their efficiency and thus F. However, M and S can be determined from observations (below).Therefore, to place our results in a present-day context, we consider https://doi.org/10.1038/s41561-024-01465-7 the critical hydrological network efficiency, F c , which is shown as a function of M and S in Fig. 4c (as before, unbounded intrusion occurs for F < F c ).Assuming a uniform hydrological network efficiency for all ice shelves, F c is a proxy for susceptibility to passing the tipping point: regions of parameters space with darker (lighter, respectively) colours in Fig. 4c are more (less) susceptible to unbounded intrusions.Locations of key Antarctic ice shelves on this map are determined using the median of observations of grounding-line velocity 43 , basal slope 40 , ocean thermal forcing 44 and literature standard values for other parameters (Methods).
Despite a large spread in the data, particularly in the grounding-line slope, we find that on average, the rapidly accelerating 45 and thinning 46 Thwaites Glacier is the least susceptible of those ice shelves considered.This is perhaps surprising, given its high ocean forcing (Supplementary Fig. 5); however, its high grounding-line velocities (Supplementary Fig. 6) overcompensate for this melting potential.This highlights the stabilizing potential of high grounding-line velocities to warm water intrusion.However, Thwaites may be particularly sensitive to future changes: ice shelves corresponding to smaller M are more sensitive to changes in grounding-line slope (Fig. 4b), which Thwaites may be exposed to as it maintains its present retreat 40 .
The Filchner and Amery ice shelves also have low susceptibility (Fig. 4c).Compared to Thwaites, however, their low susceptibility results from low ocean forcing and strongly prograde bedslopes (Supplementary Figs.4-6), respectively, rather than high grounding-line velocities.Pine Island, on the other hand, has high susceptibility; like Thwaites, it has high grounding-line velocities, but its grounding-line slope is higher (more retrograde), on average, and thus more favourable for intrusion (Supplementary Fig. 4).Pine Island is currently Antarctica's largest contributor to sea-level rise; its high susceptibility to  https://doi.org/10.1038/s41561-024-01465-7 grounding-zone melting may represent another factor, in addition to a highly damaged ice shelf 47 and predicted future increases in ice-shelf melting 48,49 , which promotes its ongoing retreat.Both Getz and Larsen have similar susceptibility to Pine Island; although Getz has similar thermal forcing, its lower grounding-line velocities (higher M) are balanced by reduced grounding-line slopes (lower S).Larsen has high susceptibility owing to its low flow speeds (large M): ice is not advected quickly enough to replace that removed by melting.More generally, these examples highlight the complex interplay between ice velocity, ocean forcing and bedslope in controlling melt in grounding zones of marine ice sheets.Our results do not provide a prediction of which Antarctic ice shelves currently experience warm water intrusions but rather indicates their potential for such and their relative susceptibility.We speculate, however, that properties of subglacial hydrological networks may be inferred from our results; for example, recent observations of high melt upstream of the Thwaites grounding line 14 suggests this area may be in the unbounded intrusion regime.According to our modelling, such behaviour would require an inefficient subglacial network (Fig. 4c), which is consistent with observations of ponding in distributed canals beneath Thwaites 50 .
We stress that completely unbounded intrusions are not expected in practice but rather the possibility of large, rapid, increases in intrusion distance and thus melting; on longer lengthscales, processes such as bedslope variations and melt feedbacks on channel temperature may suppress intrusion (Supplementary Information), and along-grounding-zone variations, not included in our model, may play an important role.Our model does not provide a prediction of the ice-dynamic response to unbounded intrusion, which requires a coupled ice-hydrology model to determine; increases in melting may lead to ice acceleration (increasing V), potentially stabilizing the intrusion (effectively reducing M).Finally, our model does not include tides.Grounding-zone characteristics may vary substantially over tidal cycles (for example, refs.51,52), potentially affecting intrusion.As an ice shelf is raised in response to tides, water is forced into the grounding zone and evacuated as the ice shelf lowers.Given that grounding lines can migrate long (up to kilometre) distances over diurnal tidal cycles 51,53 , this has the possibility to create rapid flow in the grounding zone.In addition, the associated tidal flexure may lead to a modification of the grounding-zone geometry, potentially feeding back on intrusion, which is sensitive to the characteristic thickness of the subglacial environment 5 .Tidal currents may also modulate near-grounding-zone ocean circulation 2 , potentially altering flow boundary conditions on grounding zones and thus intrusion.
A complete treatment of grounding-zone flow on tidal timescales is beyond the scope of this work.However, when supplemented with a simple parameterization of tidal flow (Supplementary Information), our model still displays the tipping-point behaviour, with the location of the tipping point (in parameter space) modulated by the tidal amplitude (Supplementary Fig. 10).We find that tidal fluctuations can significantly enhance intrusion.This provides further motivation to better understand tidal influences on grounding zones, and, more generally, to develop high-resolution models of grounding zones and better constrain their characteristics via improved observations.
We have shown that feedbacks between subglacial water flow, melting and the confining ice geometry can result in increases in warm water intrusion into marine ice-sheet grounding zones, which would have implications for ice dynamics.In particular, we have identified a fundamental switch between bounded and unbounded warm water intrusions, occurring across a critical parameter threshold.The tipping  https://doi.org/10.1038/s41561-024-01465-7 point is generic: it exists for any marine-terminating ice sheet exposed to sufficiently warm ocean water, has sufficiently low grounding-line velocity or basal slopes or a sufficiently weak hydrological network.
We have shown that the intrusion mechanism is stronger than previously understood, lending further credence to the theory that it is a physically based 'sensitivity-boosting mechanism' to reconcile the gap between observed and modelled sea-level rise in previous warm periods and the basal melt rates required to reproduce observed retreat.Current sea-level-rise projections for Antarctica and Greenland 54 are based on simulations that lack grounding-zone melting via intrusion and may therefore represent underestimates.Although our model is a simplification of the myriad complex processes occurring in grounding zones, the possibility of tipping points in grounding-zone melt and the universality of susceptible shelves warrants a continued research effort to better constrain grounding-zone processes both from models and observations.https://doi.org/10.1038/s41561-024-01465-7

Layered intrusion model
The layered intrusion model is identical to that of ref. 5 in the hard bed limit (γ → 0 in the nomenclature of ref. 5).This model builds upon that of ref. 11 and was verified experimentally therein.It is described in full detail in the Supplementary Information.

Melt model
We couple the layered intrusion model to a simple model of melting, where ṁ is the melt rate, ℒ = 3.35 × 10 5 J Kg −1 is the latent heat of fusion of seawater, c = 3.974 × 10 3 J Kg −1 °C −1 is the specific heat capacity of water and St is a combined Stanton number, which parametrizes combined exchange of salt and heat across a thermal boundary layer that forms on the ocean side of the ice-ocean interface 55,56 (u* and τ are defined below).The Stanton number is fairly poorly constrained in general.In the results presented here, we take value St = 5.9 × 10 −4 ; this value is standard in the literature and was obtained from a fit to data obtained from beneath the Ronne ice shelf 57 .Equation ( 3) results from the so-called 'two-equation formulation' of melting 55 in the limit of low diffusive heat flux (this is reasonable as freezing and internal temperatures of the ice are typically within a few degrees of one another 56 ) (Supplementary Information).The two-equation formulation is a simplification of the more detailed 'three equation formulation' 33,58 in which salt and heat exchange across the boundary layer at the ice-shelf base are considered separately rather than together as in the two-equation formulation.However, the two formulations have been shown to work equally well in several observational (for example, ref. 57) and numerical (for example, ref. 59) studies.
In equation ( 3), u* and τ are the velocity of the water adjacent to the ice-ocean interface and the thermal driving, respectively (the latter should not be confused with the thermal forcing ΔT).We take u* to be the velocity of the fresh layer, which is the layer adjacent to the ice-ocean interface.The thermal driving is τ = T − T f , where T is the temperature adjacent to the ice-ocean interface (below) and T f = T ref + λz − Γ  the local freezing temperature, with T ref = 8.32 × 10 −2 °C a reference temperature, λ = 7.61 × 10 −4 °C m −1 the liquidus slope with depth, Γ = 5.73 × 10 −2 the liquidus slope with salinity,  the local salinity and z the depth below sea level (more negative z corresponds to a greater depth) 60 .
We take a simple model for the channel temperature and salinity, assuming that the relevant temperature and salinity that drive melting are the depth-weighted average of the layer temperatures: where ϕ is the column-wise proportion of the channel occupied by the freshwater layer (Supplementary Information), T D = 0 and S D = 0 are the temperature and salinity of the subglacial discharge layer, respectively, and T O and S O are the temperature and salinity of the warm ocean layer.The relations (4) and ( 5) capture the fact that the temperature and salnity in the channel increase with a greater proportion of warm, salty ocean water within it.Although entrainment between the two layers is not explicitly resolved, the column-wise averaging can be considered a simple proxy for mixing of the two layers.Observations also indicate that basal melting can occur where a cold fresh layer exists adjacent to an ice-ocean interface through double diffusive convection 52,61 .

Channel shape evolution
The dimensionless model equations (equations ( 19) and ( 20) in Supplementary Information) are solved numerically in MATLAB.For a given channel shape, the layered intrusion equations are solved using the ODE15S routine.The equations are solved backwards from the downstream end of the channel, where we apply a perturbed boundary condition, setting the dimensionless freshwater layer thickness equal to , where ϵ ≪ 1 (Supplementary Information).This perturbed boundary condition ensures that a singularity in the interfacial gradient is avoided at the downstream end of the channel 5 .In those results shown here, we use ϵ = 10 −4 but verified that results are insensitive to this value, provided that the ϵ ≪ 1 condition holds.The intrusion equations are integrated backwards until either the freshwater layer occupies the entirety of the channel or the end of the numerical grid is reached (we use a sufficiently large numerical grid to ensure that the latter is only realized in the case of unbounded intrusion).
Having determined the interfacial shape, and thus velocity in the fresh layer and channel temperature, the melt rate is determined from equation ( 3).This melt rate is interpolated onto a regular grid with spacing dz (below) and the channel thickness timestepped according to the kinematic condition (equation ( 14) in the Supplementary Information) using a first-order upwinding scheme 62 .
The numerical grid is made up of m blocks of n grid cells (giving a total number of grid cells of m × n); each block is of length L p = 1 − F 2 /4 − 3F 2/3 /4, which is the intrusion distance in the limit of no interfacial drag (C = 0), a flat bed (S = 0) and no melting (M = 0) as described by ref. 5; the grid size is then dz = L p /n.In the results shown herein, we use n = 100 and m = 20, with the latter value being sufficiently large that the intrusion only reaches the end of the channel in the case of an unbounded intrusion.

Steady intrusion length
To determine the steady intrusion length L for a given set of parameters (F, C, S, M), we integrate the steady form of the coupled layered intrusion-melt equations (equations ( 24) and ( 25) in the Supplementary Information) downstream from the nose of the wedge (where the freshwater layer occupies the width of the channel) using the ODE15S routine in MATLAB.At the nose, the problem is singular; to avoid this singularity, we linearize the problem about this point to determine the appropriate initial conditions (Supplementary Information).Solutions of this steady problem, which obtain the downstream boundary condition (dimensionless freshwater layer thickness = F 2/3 ) in a finite distance correspond to true steady states; otherwise, no steady solution exists for this particular set of parameters, and the intrusion will be unbounded (Supplementary Fig. 2).In practice, we specify a finite end of the domain ℓ ≫ 1, and, if the solution does not obtain the downstream boundary condition before this point, we assume the intrusion is unbounded; in the results shown here, we take ℓ = 10 5 , and results are insensitive to this value.
To determine the critical slope S c and critical Froude number F c shown in Fig. 4, we apply a bisection method.The algorithm to do so is described fully in the Supplementary Information.

Parameter estimation
Values of parameters M and S for Antarctic ice shelves shown in Fig. 4 were determined from observations of ice velocity (for grounding-line velocity V), thermal forcing (for ΔT) and bedslopes (for θ).
To determine grounding-line locations, we first obtain ice-shelf boundaries from Bedmachine V3 40 masks of ice-shelf location at 500 m resolution.Grid points within this mask corresponding to grounding-line and ice-shelf front positions were differentiated based on a floatation condition, relating to the floatation thickness ρ w /ρ i b, the ice thickness at which hydrostatic equilibrium is achieved, where ρ w = 1,028.0kg m −3 is the ocean density, ρ i = 918.0kg m −3 is the ice density and b is the bed elevation determined from Bedmachine V3 bed data 40 .Ice-shelf boundary points above 95% of floatation thickness were designated as grounding-line points, whereas the remaining points were designated as ice front points.Supplementary Fig. 9 shows https://doi.org/10.1038/s41561-024-01465-7 grounding-line and ice front points for each of the ice shelves shown in Fig. 4, indicating that this criterion does a good job at correctly identifying, and differentiating between, grounding-line and ice front points.
Ice velocities at grounding-line points were determined from NASA ITS_LIVE mosaics of 1985-2019 ice velocities 43 .We first put the data onto the same 500 m resolution grid used to determine grounding-line locations and then extract velocity components v = (v 1 , v 2 ) at these points (inset in Supplementary Fig. 9a).Grounding-line ice velocities used to compute M are taken as the median of all grounding-line velocities within the individual ice shelf (Supplementary Fig. 6).
Grounding-line slopes were determined by first extracting bed data from Bedmachine V3 40 onto the 500 m grid.We then compute gradients of basal elevation, ∇ b = (∂ x b, ∂ y b), where ∂ indicates a partial derivative, using second order finite differences.We then take the basal slope as the directional derivative in the direction of ice flow, that is The value of tan θ used for each is shelf is then determined as the median of all grounding-line basal slopes for that particular shelf (Supplementary Fig. 4).Thermal forcing is computed from maps of maximum thermal forcing between 200 m and 800 m in the water column from ref. 44.For each ice front grid point on the 500 m grid, the thermal forcing associated with that point is computed as the mean of the maximum thermal forcing within a 1.5 km × 1.5 km square centred around the grid point; the thermal forcing of each ice shelf is then determined as the median over all well-defined ice front points in the shelf (Supplementary Fig. 5).
Having determined V, ΔT and tan θ from observations, M and S are computed using standard values from the literature.We take C d = 10 −2 , which is consistent with a range of observations (for example, ref.  1 , c = 3.974 × 10 3 J Kg −1 °C −1 and St = 5.9 × 10 −4 as discussed above.The mean upstream flow velocity U ∞ is taken to be 1 cm s −1 , which is consistent with observations 65,66 and modelling 67,68 of subglacial flow beneath ice sheets.

Fig. 1 |
Fig. 1 | Complex grounding zones of ice sheets.The grounding zone (dashed circle in top left) of marine-terminating ice sheets features networks of tunnels, channels and porous sediments through which water moves (centre).Layered intrusion models (cross section at right) consider this region as a two-layer system: freshwater, at the local freezing temperature T f , enters the zone with average velocity U ∞ , where it meets warm ocean water of temperature T O and

Fig. 2 |
Fig. 2 | Tipping points in grounding-zone melting.a,e, Temporal evolution of the intrusion distance-the greatest upstream extent of the warm layer-for a thermal forcing of ΔT = 2.3 °C (left) and ΔT = 2.5 °C (right).For ΔT = 2.3 °C, the intrusion distance tends towards a bounded value (L ≈ 110 m), indicated by the black dashed line, whereas for ΔT = 2.5 °C, the intrusion becomes unbounded (L → ∞) (note the different ordinate scales between the left and right panels).Translucent points in e show the data in a. b,f, Evolution of grounding-zone channel surface (solid curves) and warm-cold interface (dashed curves).

Fig. 3 |
Fig. 3 | Regimes of warm water intrusion beneath marine ice sheets.Regime diagram of dimensionless intrusion distance L/(H ∞ /c d ) as a function of upstream Froude number F and dimensionless melt rate M for C = 0.1 and S = 0, corresponding to a flat base.Coloured lines indicate the bounded-unbounded intrusion length transition for different values of C, as indicated.The red line indicates schematically the transition between the two values of M used in Fig. 2.

Fig. 4 |
Fig. 4 | Conditions for unbounded intrusions.a, Dimensionless intrusion distance L/(H ∞ /c d ) as a function of dimensionless bedslope S for different M, as indicated by the colour bar.Filled circles indicate the critical slope S c , beyond which the intrusion becomes unbounded.The black dashed line indicates S = 0, which delineates pro-(S < 0) and retrograde (S > 0) bed slopes.Data shown here correspond to F = 0.5 and C = 0.1.b, Curves of critical upstream Froude number F c , below which unbounded intrusion occurs, as a function of dimensionless slope S for different M, as indicated by the colour bar in a. c, Map of critical upstream Froude upstream number F c .Points indicate ice shelves in Antarctica (PIG = Pine Island Glacier, PSK = Pope-Smith-Kohler), determined from observations (see main text).Grey contours correspond to F c = 0.3 and F c = 0.7.Ellipses indicate errors in distributions used to generate values of S and M. Points on the inset indicate locations of those ice shelves in the main panel with corresponding colours.