Abstract
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.
Similar content being viewed by others
Main
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 observations1,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 loss7 and struggle to reproduce observationally constrained sea-level highstands from previous interglacials3,4. Palaeoclimate ice-sheet reconstructions have largely been able to reproduce low-end estimates only when mechanisms to boost sensitivity to climatic forcing are invoked8,9,10. Recently, evidence from diverse sources has emerged suggesting that relatively warm ocean water can intrude long distances upstream of ice-shelf grounding lines5,6,11,12,13,14; such long intrusions have dramatic consequences for sea-level-rise contributions from ice sheets5,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 influence16. During the Pliocene (∼3 million years before present), CO2 levels were similar to present day17 and temperatures 2–3 °C above pre-industrial levels18, but the global mean sea level (GMSL) was 6–40 m above present-day levels4. Such sea-level rise is only possible with a sizeable contribution from the Antarctic ice sheet16. 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. 22).
Several palaeoclimate simulations8,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 physical23,24 and statistical25 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 melting15,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 line11. There is growing evidence of warm water intrusions from diverse sources including surface observations12,13, satellite data6,14 and ice-shelf basal features32.
Recent near-grounding-zone observations1,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 models5,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 lines5, delivering excess heat for sub-ice melting. Ice-sheet models including the intrusion mechanism show that ice loss is highly sensitive to this intrusion length5, 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 periods5.
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 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 interface33 (Supplementary Information); c is the specific heat capacity of water; cd 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 TO the ocean temperature, \({{{{\mathcal{S}}}}}_\mathrm{O}\) the ocean salinity, Γ the freezing point slope with salinity and TD 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 ρ0 a reference density and ρw the local water density; and ci 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 sheet5,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 Fc, 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 Fc 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 friction34.
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 indicators35, 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 projections36.
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, Sc, above which unbounded intrusion occurs, reduces with M (Fig. 4a). Over large regions of parameter space, Sc 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,38,39). Retreat of the West Antarctic ice sheet, many areas of which have grounding zones on retrograde bed slopes40, may be strongly controlled by MISI41. 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 stable38,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, Fc, which is shown as a function of M and S in Fig. 4c (as before, unbounded intrusion occurs for F < Fc). Assuming a uniform hydrological network efficiency for all ice shelves, Fc 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 velocity43, basal slope40, ocean thermal forcing44 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 accelerating45 and thinning46 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 retreat40.
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 shelf47 and predicted future increases in ice-shelf melting48,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 line14 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 Thwaites50.
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 cycles51,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 environment5. Tidal currents may also modulate near-grounding-zone ocean circulation2, 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 Greenland54 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.
Methods
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 \(\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 × 103 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 interface55,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 shelf57.
Equation (3) results from the so-called ‘two-equation formulation’ of melting55 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 another56) (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 − Tf, 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 Tref = 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, \({\mathcal{S}}\) 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), TD = 0 and SD = 0 are the temperature and salinity of the subglacial discharge layer, respectively, and TO and SO 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 convection52,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 channel5. 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 scheme62.
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 Lp = 1 − F2/4 − 3F2/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 = Lp/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 = F2/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 ℓ = 105, and results are insensitive to this value.
To determine the critical slope Sc and critical Froude number Fc 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 V340 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−3 is the ocean density, ρi = 918.0 kg m−3 is the ice density and b is the bed elevation determined from Bedmachine V3 bed data40. 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 velocities43. We first put the data onto the same 500 m resolution grid used to determine grounding-line locations and then extract velocity components v = (v1, v2) 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 V340 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. 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 \theta\) from observations, M and S are computed using standard values from the literature. We take Cd = 10−2, which is consistent with a range of observations (for example, ref. 63), modelling (for example, ref. 48) and experiments (for example, ref. 64) of ice–ocean interactions, alongside \({\mathcal{L}}={3.35\times {10}^{5}}\;{\mathrm{J}\;{\rm{Kg}}^{-1}}\), c = 3.974 × 103 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 observations65,66 and modelling67,68 of subglacial flow beneath ice sheets.
Data availability
Data used to create the figures contained in this paper are available via Zenodo at https://doi.org/10.5281/zenodo.10895498 (ref. 69).
Code availability
Code to perform simulations and produce the figures contained in this paper are available via Zenodo at https://doi.org/10.5281/zenodo.10895498 (ref. 69).
References
Davis, P. E. et al. Suppressed basal melting in the eastern Thwaites Glacier grounding zone. Nature 614, 479–485 (2023).
Schmidt, B. E. et al. Heterogeneous melting near the Thwaites Glacier grounding line. Nature 614, 471–478 (2023).
Dutton, A. & Lambeck, K. Ice volume and sea level during the last interglacial. Science 337, 216–219 (2012).
Dumitru, O. A. et al. Constraints on global mean sea level during Pliocene warmth. Nature 574, 233–236 (2019).
Robel, A. A., Wilson, E. & Seroussi, H. Layered seawater intrusion and melt under grounded ice. Cryosphere 16, 451–469 (2022).
Cirací, E. et al. Melt rates in the kilometer-size grounding zone of Petermann Glacier, Greenland, before and during a retreat. Proc. Natl Acad. Sci. USA 120, e2220924120 (2023).
Aschwanden, A., Bartholomaus, T. C., Brinkerhoff, D. J. & Truffer, M. Brief communication: a roadmap towards credible projections of ice sheet contribution to sea level. Cryosphere 15, 5705–5715 (2021).
DeConto, R. M. & Pollard, D. Contribution of Antarctica to past and future sea-level rise. Nature 531, 591–597 (2016).
Golledge, N. R. et al. Global environmental consequences of twenty-first-century ice-sheet melt. Nature 566, 65–72 (2019).
DeConto, R. M. et al. The Paris Climate Agreement and future sea-level rise from Antarctica. Nature 593, 83–89 (2021).
Wilson, E. A., Wells, A. J., Hewitt, I. J. & Cenedese, C. The dynamics of a subglacial salt wedge. J. Fluid Mech. 895, A20 (2020).
MacGregor, J. A., Anandakrishnan, S., Catania, G. A. & Winebrenner, D. P. The grounding zone of the Ross Ice Shelf, West Antarctica, from ice-penetrating radar. J. Glaciol. 57, 917–928 (2011).
Horgan, H. J. et al. Estuaries beneath ice sheets. Geology 41, 1159–1162 (2013).
Milillo, P. et al. Heterogeneous retreat and ice melt of Thwaites Glacier, West Antarctica. Sci. Adv. 5, eaau3433 (2019).
Seroussi, H. & Morlighem, M. Representation of basal melting at the grounding line in ice flow models. Cryosphere 12, 3085–3096 (2018).
Dutton, A. et al. Sea-level rise due to polar ice-sheet mass loss during past warm periods. Science 349, aaa4019 (2015).
Seki, O. et al. Alkenone and boron-based pliocene pCO2 records. Earth Planet. Sci. Lett. 292, 201–211 (2010).
Haywood, A. M. et al. The Pliocene Model Intercomparison Project Phase 2: large-scale climate features and climate sensitivity. Clim 16, 2095–2123 (2020).
Scherer, R. P. et al. Pleistocene collapse of the West Antarctic Ice Sheet. Science 281, 82–85 (1998).
Hillenbrand, C.-D. et al. West Antarctic Ice Sheet retreat driven by Holocene warm water incursions. Nature 547, 43–48 (2017).
Wilson, D. J. et al. Ice loss from the East Antarctic Ice Sheet during late Pleistocene interglacials. Nature 561, 383–386 (2018).
Golledge, N. R. et al. Antarctic climate and ice-sheet configuration during the early Pliocene interglacial at 4.23 Ma. Clim 13, 959–975 (2017).
Clerc, F., Minchew, B. M. & Behn, M. D. Marine ice cliff instability mitigated by slow removal of ice shelves. Geophys. Res. Lett. 46, 12108–12116 (2019).
Bassis, J., Berg, B., Crawford, A. & Benn, D. Transition to marine ice cliff instability controlled by ice thickness gradients and velocity. Science 372, 1342–1344 (2021).
Edwards, T. L. et al. Revisiting Antarctic ice loss due to marine ice-cliff instability. Nature 566, 58–64 (2019).
Parizek, B. et al. Dynamic (in) stability of Thwaites Glacier, West Antarctica. J. Geophys. Res. Earth Surf. 118, 638–655 (2013).
Arthern, R. J. & Williams, C. R. The sensitivity of West Antarctica to the submarine melting feedback. Geophys. Res. Lett. 44, 2352–2359 (2017).
Reese, R., Gudmundsson, G. H., Levermann, A. & Winkelmann, R. The far reach of ice-shelf thinning in Antarctica. Nat. Clim. Change 8, 53–57 (2018).
Graham, A. G. et al. Rapid retreat of Thwaites Glacier in the pre-satellite era. Nat. Geosci. 15, 706–713 (2022).
Anandakrishnan, S., Catania, G. A., Alley, R. B. & Horgan, H. J. Discovery of till deposition at the grounding line of Whillans Ice Stream. Science 315, 1835–1838 (2007).
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).
Whiteford, A., Horgan, H., Leong, W. & Forbes, M. Melting and refreezing in an ice shelf basal channel at the grounding line of the Kamb Ice Stream, West Antarctica. J. Geophys. Res. Earth Surf. 127, e2021JF006532 (2022).
Jenkins, A. Convection-driven melting near the grounding lines of ice shelves and tidewater glaciers. J. Phys. Oceanogr. 41, 2279–2294 (2011).
Hewitt, I. Seasonal changes in ice sheet motion due to melt water lubrication. Earth Planet. Sci. Lett. 371, 16–25 (2013).
Lenton, T. M. Early warning of climate tipping points. Nat. Clim. Change 1, 201–209 (2011).
Robel, A. A., Seroussi, H. & Roe, G. H. Marine ice sheet instability amplifies and skews uncertainty in projections of future sea-level rise. Proc. Natl Acad. Sci. USA 116, 14887–14892 (2019).
Weertman, J. Stability of the junction of an ice sheet and an ice shelf. J. Glaciol. 13, 3–11 (1974).
Schoof, C. Ice sheet grounding line dynamics: steady states, stability, and hysteresis. J. Geophys. Res. Earth Surf. 112, F03S28 (2007).
Katz, R. F. & Worster, M. G. Stability of ice-sheet grounding lines. Philos. Trans. R. Soc. A 466, 1597–1620 (2010).
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).
Favier, L. et al. Retreat of Pine Island Glacier controlled by marine ice-sheet instability. Nat. Clim. Change 4, 117–121 (2014).
Schoof, C. Marine ice sheet stability. J. Fluid Mech. 698, 62–72 (2012).
Gardner, A. S. et al. Increased West Antarctic and unchanged East Antarctic ice discharge over the last 7 years. Cryosphere 12, 521–547 (2018).
Adusumilli, S., Fricker, H. A., Medley, B., Padman, L. & Siegfried, M. R. Interannual variations in meltwater input to the Southern Ocean from Antarctic ice shelves. Nat. Geosci. 13, 616–620 (2020).
Mouginot, J., Rignot, E. & Scheuchl, B. Sustained increase in ice discharge from the Amundsen Sea Embayment, West Antarctica, from 1973 to 2013. Geophys. Res. Lett. 41, 1576–1584 (2014).
Pritchard, H. et al. Antarctic ice-sheet loss driven by basal melting of ice shelves. Nature 484, 502–505 (2012).
Lhermitte, S. et al. Damage accelerates ice shelf instability and mass loss in Amundsen Sea Embayment. Proc. Natl Acad. Sci. USA 117, 24735–24741 (2020).
Bradley, A. T., Bett, D., Dutrieux, P., De Rydt, J. & Holland, P. R. The influence of Pine Island Ice Shelf calving on basal melting. J. Geophys. Res. Oceans 127, e2022JC018621 (2022).
Bradley, A. T., De Rydt, J., Bett, D. T., Dutrieux, P. & Holland, P. R. The ice dynamic and melting response of Pine Island Ice Shelf to calving. Ann. Glaciol. 63, 111–115 (2022).
Schroeder, D. M., Blankenship, D. D. & Young, D. A. Evidence for a water system transition beneath Thwaites Glacier, West Antarctica. Proc. Natl Acad. Sci. USA 110, 12225–12228 (2013).
Brunt, K. M., Fricker, H. A. & Padman, L. Analysis of ice plains of the Filchner–Ronne Ice Shelf, Antarctica, using ICESat laser altimetry. J. Glaciol. 57, 965–975 (2011).
Begeman, C. B. et al. Ocean stratification and low melt rates at the ross ice shelf grounding zone. J. Geophys. Res. Oceans 123, 7438–7452 (2018).
Milillo, P. et al. On the short-term grounding zone dynamics of Pine Island Glacier, West Antarctica, observed with cosmo-skymed interferometric data. Geophys. Res. Lett. 44, 10–436 (2017).
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).
McPhee, M. G. Turbulent heat flux in the upper ocean under sea ice. J. Geophys. Res. Oceans 97, 5365–5379 (1992).
Bradley, A. T., Rosie Williams, C., Jenkins, A. & Arthern, R. Asymptotic analysis of subglacial plumes in stratified environments. Proc. R. Soc. London 478, 20210846 (2022).
Jenkins, A., Nicholls, K. W. & Corr, H. F. Observation and parameterization of ablation at the base of Ronne Ice Shelf, Antarctica. J. Phys. Oceanogr. 40, 2298–2312 (2010).
Holland, D. M. & Jenkins, A. Modeling thermodynamic ice–ocean interactions at the base of an ice shelf. J. Phys. Oceanogr. 29, 1787–1800 (1999).
Lazeroms, W. M., Jenkins, A., Rienstra, S. W. & Van De Wal, R. S. An analytical derivation of ice-shelf basal melt based on the dynamics of meltwater plumes. J. Phys. Oceanogr. 49, 917–939 (2019).
Hewitt, I. J. Subglacial plumes. Annu. Rev. Fluid Mech. 52, 145–169 (2020).
Kimura, S., Nicholls, K. W. & Venables, E. Estimation of ice shelf melt rate in the presence of a thermohaline staircase. J. Phys. Oceanogr. 45, 133–148 (2015).
Patankar, S. V. Numerical Heat Transfer and Fluid Flow (CRC Press, 1980).
Lu, P., Li, Z., Cheng, B. & Leppäranta, M. A parameterization of the ice–ocean drag coefficient. J. Geophys. Res. Oceans 116, C07019 (2011).
Ezhova, E., Cenedese, C. & Brandt, L. Dynamics of three-dimensional turbulent wall plumes and implications for estimates of submarine glacier melting. J. Phys. Oceanogr. 48, 1941–1950 (2018).
Joughin, I. et al. Basal conditions for Pine Island and Thwaites glaciers, West Antarctica, determined using satellite and airborne data. J. Glaciol. 55, 245–257 (2009).
Pattyn, F. Antarctic subglacial conditions inferred from a hybrid ice sheet/ice stream model. Earth Planet. Sci. Lett. 295, 451–461 (2010).
Hewitt, I. J. Modelling distributed and channelized subglacial drainage: the spacing of channels. J. Glaciol. 57, 302–314 (2011).
Carter, S. & Fricker, H. The supply of subglacial meltwater to the grounding line of the Siple Coast, West Antarctica. Ann. Glaciol. 53, 267–280 (2012).
Bradley, A.T. & Hewitt, I.J. ‘Tipping-point in ice-sheet grounding-zone melting due to ocean water intrusion’ by Bradley and Hewitt. Zenodo https://doi.org/10.5281/zenodo.10895498 (2024).
Acknowledgements
A.T.B. is supported by the Natural Environmental Research Programme (NERC) grant NE/S010475/1 and the PROTECT project, which received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement number 869304.
Author information
Authors and Affiliations
Contributions
A.T.B. conceived of the study and performed the analysis. A.T.B. and I.J.H. contributed to scientific discussion, interpretation of the results and wrote the manuscript.
Corresponding author
Ethics declarations
Competing interests
The authors declare no competing interests.
Peer review
Peer review information
Nature Geoscience thanks Andy Aschwanden, Katarzyna Warburton and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Primary Handling Editor: Tom Richardson, in collaboration with the Nature Geoscience team.
Additional information
Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary information
Supplementary Information
Supplementary Information and Figs. 1–10.
Rights and permissions
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/.
About this article
Cite this article
Bradley, A.T., Hewitt, I.J. Tipping point in ice-sheet grounding-zone melting due to ocean water intrusion. Nat. Geosci. (2024). https://doi.org/10.1038/s41561-024-01465-7
Received:
Accepted:
Published:
DOI: https://doi.org/10.1038/s41561-024-01465-7