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

There is growing evidence suggesting that ice-sheet models lack representation of important physical processes driving ice sheet retreat (for example, refs. ^{1,2,3,4,5,6,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⁷ 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¹⁶. During the Pliocene (∼3 million years before present), CO₂ levels were similar to present day¹⁷ and temperatures 2–3 °C above pre-industrial levels¹⁸, but the global mean sea level (GMSL) was 6–40 m above present-day levels⁴. Such sea-level rise is only possible with a sizeable contribution from the Antarctic ice sheet¹⁶. Whereas proxy records suggest that such Antarctic ice loss occurred during previous interglacials (for example, refs. ^19,20,21), ice-sheet models struggle to reproduce the corresponding ice-sheet retreat (for example, ref. ²²).

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²⁵ grounds. Ref. ²² 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.

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 salinity S_O. 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 along-grounding-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).

Full size image

Specifically, the model in ref. ²² 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. ²⁹), porous till layers (for example, ref. ³⁰) and channels on multiple lengthscales (for example, ref. ³¹) (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¹¹. 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³².

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. These models predict that warm water can intrude significant distances (up to kilometres) upstream of ice-sheet grounding lines⁵, 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⁵, 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⁵.

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.

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 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³³ (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; \({{\Delta}}T={T}_\mathrm{O}+{{\varGamma}}{{{\mathcal{S}}}}_\mathrm{O}-{T}_\mathrm{D}\) is the thermal forcing, with T_O the ocean temperature, \({{{{\mathcal{S}}}}}_\mathrm{O}\) the ocean salinity, Γ the freezing point slope with salinity and T_D the local freshwater freezing temperature; \({\mathcal{L}}\) is the latent heat of fusion of seawater; θ is the local grounding-zone slope, assumed constant (Supplementary Information); \({g}^{\,{\prime} }=g(\,{\rho }_{0}-{\rho }_\mathrm{w})/{\rho }_\mathrm{w}\) is the reduced gravity, with ρ₀ 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).

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). Snapshots are shown at times 1, 5, 10, 50 and 100 days after initialization with a flat channel (the same initial condition is used in both cases). y is the distance to the seabed, which forms the base of the channel. Filled points indicate the intrusion distance and correspond to points shown in a as indicated by Roman numerals. In b, the dashed black line indicates the steady state intrusion distance L. c,d,g,h, Profiles of thermal driving (c,g) and ice-adjacent flow velocity (d,h) corresponding to the snapshots shown in (b,f). Solutions here correspond to a flat bed (S = 0), a fairly inefficient drainage system (F = 0.25) and C = 0.1.

Full size image

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³⁴.

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.

Full size image

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³⁵, 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³⁶.

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. ⁵. 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).

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.

Full size image

Marine ice sheets grounded on retrograde bedslopes may be susceptible to the marine ice-sheet instability (MISI) (for example, refs. ^37,38,39). Retreat of the West Antarctic ice sheet, many areas of which have grounding zones on retrograde bed slopes⁴⁰, may be strongly controlled by MISI⁴¹. 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 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⁴³, basal slope⁴⁰, ocean thermal forcing⁴⁴ 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⁴⁵ and thinning⁴⁶ 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⁴⁰.

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 grounding-zone melting may represent another factor, in addition to a highly damaged ice shelf⁴⁷ 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¹⁴ 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⁵⁰.

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⁵. Tidal currents may also modulate near-grounding-zone ocean circulation², 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 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⁵⁴ 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.

Layered intrusion model

The layered intrusion model is identical to that of ref. ⁵ in the hard bed limit (γ → 0 in the nomenclature of ref. ⁵). This model builds upon that of ref. ¹¹ 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 \(\dot{m}\) is the melt rate, \({\mathcal{L}}={3.35\times {10}^{5}}\;{\mathrm{J}\;{\rm{Kg}}^{-1}}\) is the latent heat of fusion of seawater, c = 3.974 × 10³ J Kg⁻¹ °C⁻¹ 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⁻⁴; this value is standard in the literature and was obtained from a fit to data obtained from beneath the Ronne ice shelf⁵⁷.

Equation (3) results from the so-called ‘two-equation formulation’ of melting⁵⁵ 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⁵⁶) (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. ⁵⁷) and numerical (for example, ref. ⁵⁹) 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}_\mathrm{f}={T}_{\rm{ref}}+\lambda z-{{\varGamma}}{\mathcal{S}}\) the local freezing temperature, with T_ref = 8.32 × 10⁻² °C a reference temperature, λ = 7.61 × 10⁻⁴ °C m⁻¹ the liquidus slope with depth, Γ = 5.73 × 10⁻² the liquidus slope with salinity, \({\mathcal{S}}\) the local salinity and z the depth below sea level (more negative z corresponds to a greater depth)⁶⁰.

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 \({\left[(1+\epsilon )F\right]}^{2/3}\), 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⁵. In those results shown here, we use ϵ = 10⁻⁴ 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⁶².

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²/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. ⁵; 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⁵, 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⁴⁰ 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/ρ_ib, the ice thickness at which hydrostatic equilibrium is achieved, where ρ_w = 1,028.0 kg m⁻³ is the ocean density, ρ_i = 918.0 kg m⁻³ is the ice density and b is the bed elevation determined from Bedmachine V3 bed data⁴⁰. 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 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⁴³. 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₁, v₂) 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⁴⁰ onto the 500 m grid. We then compute gradients of basal elevation, ∇ b = (∂_xb, ∂_yb), 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 \theta\) 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. ⁴⁴. 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 \theta\) from observations, M and S are computed using standard values from the literature. We take C_d = 10⁻², which is consistent with a range of observations (for example, ref. ⁶³), modelling (for example, ref. ⁴⁸) and experiments (for example, ref. ⁶⁴) of ice–ocean interactions, alongside \({\mathcal{L}}={3.35\times {10}^{5}}\;{\mathrm{J}\;{\rm{Kg}}^{-1}}\), c = 3.974 × 10³ J Kg⁻¹ °C⁻¹ and St = 5.9 × 10⁻⁴ as discussed above. The mean upstream flow velocity U_∞ is taken to be 1 cm s⁻¹, which is consistent with observations^65,66 and modelling^67,68 of subglacial flow beneath ice sheets.