Main

The West Antarctic ice sheet (WAIS) has been losing mass since the start of the satellite era1 and has contributed almost 90% of the overall Antarctic ice sheet (AIS) mass loss since 19922. In the Amundsen Sea Embayment (ASE), in particular, there has been widespread thinning3, accelerated ice flow1 and grounding-line retreat4, which has prompted questions about the future stability of the region5,6. Modelling studies have predicted further retreat under current and future climate conditions7,8,9,10,11 and there is a possibility of a complete collapse of the WAIS if local destabilization occurs12. Owing to the retrograde bed (sloping downwards in the inland direction) beneath its largest glaciers, the ASE is vulnerable to marine ice sheet instability13,14, where a perturbation in grounding-line position could result in irreversible mass loss and grounding-line retreat15,16. The floating extensions of glaciers, known as ice shelves, provide buttressing of upstream grounded ice and can be sufficient to restore stability to unstable grounding-line retreat17,18, particularly through the aid of pinning points such as ice rises and ice rumples19,20. However, ice shelves are susceptible to ocean-induced melting21,22, which can lead to thinning, weakened buttressing and accelerated ice flow23,24.

One of the largest ASE glaciers is Pine Island Glacier (PIG), which has contributed more to global mean sea-level rise in recent decades than any other glacier in Antarctica25. Thinning of the present-day ice shelf and grounding-line retreat can be traced back to the 1940s when an ocean cavity first started to form upstream of a subglacial ridge26. There was further grounding-line retreat and increased ice discharge in the 1970s with the ungrounding of an ice rumple over the highest part of the ridge1,27. These events in the 1940s and 1970s coincide with notable climate anomalies in the central tropical Pacific, which has been shown to have a teleconnection with the Amundsen Sea28. It is possible that tropically forced wind anomalies over the continental shelf break29 caused a shallowing of the thermocline, allowing more warm circumpolar deep water to access the cavity underneath the ice shelf, leading to higher melt and enhanced thinning30,31,32,33. Previous ice-flow modelling studies have shown that a shallower thermocline can cause irreversible retreat of an idealized representation of PIG and this happens when there is a sufficient gap between the subglacial ridge and ice shelf34,35.

Here, we investigate the retreat of PIG from the ridge and whether the marine ice sheet instability played a role in that retreat. To do this, we use the finite-element, vertically integrated ice-flow model Úa36 to solve the ice dynamics equations in the shallow ice-stream approximation. We first advance a present-day PIG configuration to a steady-state position on the subglacial ridge. This is then perturbed with control forcing that represents mean ocean conditions in the Amundsen Sea. Following this, a warm forced perturbation is applied, which has a shallower thermocline, to represent conditions during a warm period. We use a depth-dependent melt-rate parameterization with a piecewise-linear profile in both scenarios. The final experiment explores the stability regime of the glacier by incrementally changing the basal melting in retreat and advance steady-state simulations.

Pre-1940s Pine Island Glacier

The model starts from a present-day representation of PIG, with the grounding line of the main central trunk sitting on a 1,200 m deep section of bedrock, 47 km upstream of the subglacial ridge crest (Fig. 1). It is then run for 500 years, with no basal melting, allowing a new steady state to be reached. During this period, the ice stream thickens and advances forward, reduces in speed and fully grounds on the ridge. Steady state is reached within 150 years, with no further change in the central grounding-line position in the remaining 350 years of the run (Supplementary Fig. 4). The final ice flux, which is calculated along the present-day grounding-line position (dotted purple in Fig. 1), is 67 Gt yr−1, which is almost within the error range of the earliest observed ice flux, when PIG was still grounded on the ridge1. It is also similar to the overall surface mass balance of the PIG basin25, showing that the glacier is close to a balanced state.

Fig. 1: Pine Island Glacier subjected to different basal melt forcing.
figure 1

a,b, Bedrock elevation with overlain grounding lines (a) and flowline profiles (b) for the initial model setup, control and warm simulations. The flowline position is shown in dashed cyan in a. The location of PIG is shown in the inset map of Antarctica in a, which is displayed using the polar stereographic projection (xps, yps). In both panels, the present-day geometry is shown in dotted purple and the steady-state geometry after no basal melting for 500 years is shown in dash-dotted purple. The black solid line shows the geometry after 100 years of control forcing, with a 700 m thermocline depth and the red solid line shows the geometry following another 50 years of warm forcing, with a 600 m thermocline. The zero position along the flowline in b corresponds to the present-day grounding-line position.

Full size image

For the following experiments, we apply a simple depth-dependent melt-rate parameterization, similar to an approach in a previous PIG study8. The parameterization represents a two-layer ocean, typically used for conditions in the Amundsen Sea35,37, with zero melting at shallow depths and maximum melting in the deeper areas (Supplementary Fig. 5). Between the two layers is a linearly varying melt rate, which represents the ocean thermocline.

We first run a 100 year simulation with control forcing, which represents average conditions in the Amundsen Sea37. This has a maximum melt rate of 100 m yr−1 below a depth of 700 m (refs. 35,38,39), decreasing to zero melt at 300 m. The highest melt rate in this model run is at the depth of the ridge crest, hence, much of the ice shelf is initially exposed to high melting (Supplementary Fig. 6).

At the start of the run, the integrated melt rate across the ice shelf is 144 Gt yr−1, the mean melt rate is 60 m yr−1 and ice flux across the upstream gate (dotted purple in Fig. 1) is 67 Gt yr−1. After 100 years, the grounded ice has thinned by an average of 24 m, floating ice has thinned 200–300 m (Fig. 1) and the ice-stream central trunk has sped up by 20%. The ice shelf rapidly thins in response to the high melting, transforming the profile of the ice-shelf lower surface from convex to concave. The thinning causes grounding-line retreat across the ridge crest, with the slowest retreat occurring from the north end of the ridge, where the bedrock is shallow and wide and the fastest retreat occurs from the south, along the deep bedrock channel (Fig. 1). Upstream of the grounding line, the thinned grounded ice causes two isolated cavities to form. Despite the grounding line retreating between 5 and 20 km, it remains grounded along the ridge (Fig. 1). By the end of the simulation, the mean melt rate decreases to 18 m yr−1 and the integrated melt rate decreases to 48 Gt yr−1, which agrees with observations of average melt rates beneath PIG32,37. Owing to faster flowing ice, the ice flux increases to 79 Gt yr−1, which compares well with the earliest recorded observation1. The final configuration of the control case is an estimation of how PIG was situated before the 1940s.

Rapid retreat from a subglacial ridge

Following the control forcing experiment, we simulate the response of PIG to warmer ocean conditions by raising the melt-rate profile by 100 m so that the maximum melt rate is below a depth of 600 m and decreases to zero at 200 m. This is representative of the warmest temperature profiles that were observed in 200940 and this step change in forcing is a similar method to other studies8,35,41. In this experiment, the highest melt rate is above the depth of the ridge crest, which compared to the end of the control case results in more than a doubling of the starting mean melt rate (40 m yr−1) and integrated melt rate (120 Gt yr−1) across the ice shelf (Supplementary Fig. 6).

After 50 years of warm forcing (Fig. 2), there is an average of 25 m further thinning of grounded ice, 100–200 m thinning of floating ice and a speed up of almost 30% along the central trunk. During this warm simulation there is a further 10–20 km grounding-line retreat, which, in contrast to the control run, causes an ungrounding from the ridge crest and a new grounding-line position located at the next raised section of bedrock. By the end of the experiment, the mean melt rate decreases to 20 m yr−1 and the integrated melt rate decreases to 74 Gt yr−1. The ice flux at the end of the warm forcing simulation is 96 Gt yr−1, which is comparable to the 1996–2000 observations, when PIG was grounded in a similar position1.

Fig. 2: Warm forced retreat of Pine Island Glacier.
figure 2

a,b, Bedrock elevation with overlain grounding lines (a) and flowline profiles (b) during the warm forcing experiment. The zero position along the flowline in b corresponds to the present-day grounding line. c, Grounding-line position during the model simulation along the dashed cyan flowline that is shown in a. d,e, Total integrated melt rate over the entire ice shelf (d) and grounding-line flux and calving flux (e) during the experiment. The grounding-line flux is calculated along the present-day grounding-line position (dotted purple in Fig. 1) for all timesteps. The colour of grounding lines, profiles and plot markers in all panels show the model year during the experiment (increment of 2 yr). Shaded and unshaded regions in c–e indicate the different stages of retreat. Open markers in c–e show the steady-state grounding-line position, integrated melt rate and ice fluxes, respectively.

Full size image

The temporal changes in total melt and ice fluxes reveal different stages of the retreat during the warm experiment (Fig. 2), similar to a previous idealized study of PIG35. During the first stage, for approximately 8 years, there is a little thinning of floating and grounded ice as the ice shelf experiences higher melt rates. This causes a gradual increase in grounding-line ice flux, a small retreat across the ridge and a decrease in integrated melt rate as the ice shelf thins. During this period, the two isolated cavities start to enlarge and then merge with each other upstream of the ridge but they remain disconnected from the main outer ice-shelf cavity, so they do not experience any ocean-induced melting.

The next stage of retreat, between 8 and 17 years, is illustrated by rapid grounding-line retreat across several areas of retrograde bed (Fig. 2 and Supplementary Figs. 7 and 8) and the upstream cavities merge with the main outer cavity through the deep southern channel. This creates an ice rumple over the ridge in the North and then the ice shelf ungrounds completely around 17 yr. This stage of retreat is illustrated by a sharp increase in integrated melt rate as the grounding line enters a deeper section of the bedrock and experiences higher melting, causing a notable increase in grounding-line ice flux.

For the final retreat stage, from 18 years until the end of the simulation, there is gradual grounding-line retreat onto the next prominent section of bedrock and a slow decrease in integrated melt rate and ice flux as the ice shelf continues to thin (Fig. 2). The final grounding-line position, melt rate and ice fluxes all approach their steady-state values as the ice stream stabilizes in its new upstream position.

Hysteresis behaviour of Pine Island Glacier

To assess whether the warm forced retreat is reversible we perform a reversibility analysis, which consists of 38 separate steady-state simulations, with 19 comprising a retreat group and 19 a subsequent advance group. The retreat simulations all start from the no-melt steady-state solution at the ridge, approximately 47 km downstream of the present-day grounding line (Fig. 1). The advance simulations all start from the final steady-state solution of the last model run in the retreat group, which is approximately 11 km from the present-day grounding line. All model simulations in the two groups have a different thermocline depth in the melt forcing, which ranges from 1,300 to 400 m (Fig. 3) and each is run to a steady state, which indicates how far the grounding line can move under each forcing.

Fig. 3: Reversibility experiments.
figure 3

a–d, Final steady-state grounding lines (a,c) and flowline profiles (b,d) for model simulations with different thermocline depths. e, Flowline grounding-line position as a function of thermocline depth for each model simulation. Upward pointing triangles indicate the final grounding-line position for simulations that start at the subglacial ridge in the retreat group. Downward pointing triangles indicate the final grounding-line position for simulations that start at the upstream ice plain position in the advance group. The advance simulations all start from the final grounding-line position from the last retreat simulation (11 km along flowline). The open black and red triangles indicate the grounding-line position from the control (700 m) and warm (600 m) transient experiments, respectively, that were shown in Fig. 1. Note the solid black lines between markers in e are not results of model runs but are for visual purposes only. The grey shaded region in e corresponds to steep retrograde sections of bed, which are also indicated by black dash-dotted lines in all panels. The location of the flowline is shown as a dash-dotted line in a and c.

Full size image

The first six retreat simulations, with thermocline depths between 1,300 and 1,050 m, do not cause any thinning of the ice shelf because the highest melt rates are deeper than the ridge crest and lower ice surface. As the thermocline is raised above 1,050 m, the steady-state solutions show a gradual, continual thinning of the ice shelf and migration of the grounding line from the front of the ridge to the back (Fig. 3). Once the thermocline is raised above 700 m, the steady-state grounding-line retreats a further 20 km from the ridge crest to the next prominent high point in the bed. For thermocline depths above 650 m there is only a further 5 km of retreat, with the final steady-state grounding line stabilizing on the upstream ice plain.

The large migration in grounding-line position in response to the small change in thermocline depth above 700 m shows that the grounding line is highly sensitive to changes in the melt forcing but does not necessarily mean that a stability threshold has been crossed. Therefore, we reverse the forcing to explore the response of the grounding line. As the thermocline is lowered from a depth of 400 to 1,000 m, there is a gradual thickening of the ice shelf and 8 km grounding-line advance from the upstream bed rise. The thermocline must be lowered below 1,000 m for the melt rates to become small enough to allow for sufficient thickening of the ice shelf and regrounding on the ridge. There is no change in ice-shelf thickness or grounding line once the thermocline is lowered beneath a depth of 1,050 m and the steady-state position coincides with the original starting position, 47 km from the present-day grounding line.

It is evident from this experiment that a hysteretic behaviour exists when PIG is forced with a changing thermocline depth in the melt parameterization. There are non-unique steady states for the same forcing, whereby the final grounding-line position depends on the history of forcing applied, whether the glacier has been retreating or advancing. The stable steady-state positions are generally situated on the prograde slopes of the ridge and ice plain and unstable regions on the large retrograde sections (Fig. 3 and Supplementary Figs. 9 and 10). However, this does not hold everywhere as there are some local differences, which are possibly due to ice shelf buttressing36. Hence, we cannot make a general statement about bed slope and ice sheet stability as can be done for the one-dimensional example13,14. There are two threshold thermocline depths at 700 and 1,000 m, that when crossed, lead to irreversible grounding-line motion. These are irreversible transitions because the thermocline depth must be changed more than the reverse forcing to achieve the same grounding-line position. These results imply that PIG experienced a marine ice sheet instability retreat as it began to lose contact with the subglacial ridge after the 1940s climate anomaly.

Additional experiments were also carried out to test the dependency of our results on the selected model parameters and the choice of bedrock state. The first experiment used a smaller slipperiness coefficient in the Weertman sliding law (equation (7)) to account for a different sediment profile beneath the glacier42. The second experiment used a modified power law for the basal traction (equation (8)), which has been shown to affect grounding-line retreat and the rate of mass loss43,44,45. The third experiment was run on a lower bed to test the effect of solid-earth feedbacks. For this simulation the bed was lowered by 10 m at the start of the run, where we had assumed a high uplift rate of 20 cm yr−1 for our entire model period46. The final experiment used a different melt-rate parameterization, which has been used in a previous model intercomparison project for an idealized representation of the main trunk of PIG47. In all four experiments, a hysteresis was present in response to the changing thermocline depth in the melt forcing (Supplementary Figs. 11–14).

Discussion

Before the 1940s, it is likely that PIG had been grounded in a stable position on a subglacial ridge 47 km downstream from its present-day position26. Then, following notable climate anomalies and probably warmer basal conditions, in the 1940s and 1970s32,48, a pre-existing cavity beneath the ice shelf became connected with the open ocean and the glacier started to retreat from the ridge crest26,27. In the subsequent decades, PIG failed to recover its original position on the ridge, despite periods of cool ocean conditions that should have caused less melting and more thickening32. A decadal variability in local ocean conditions, largely driven by changes in the tropical Pacific Ocean32, is not reflected in the near monotonic increase in ice discharge that has been observed since the start of the satellite period in the 1970s1. By the early 1990s, the PIG grounding line had completely retreated off the ridge, across the retrograde bed, stabilizing at an ice plain 30 km upstream49 (Supplementary Fig. 15). This raises the question of whether its retreat from the ridge was an induced instability in response to the initial perturbation.

Using a vertically integrated ice-flow model and a depth-dependent melt-rate parameterization, we investigated the aspect of the retreat from the subglacial ridge that was due to internal dynamics of the system rather than changes in external forcing. The ocean forcing in this experiment is therefore simplified as we focus solely on whether the marine ice sheet instability played a role in the retreat of PIG from the ridge. Before the control simulation, the grounding line is in a stable position at the ridge crest. When basal melting is applied, to represent average ocean conditions in Amundsen Sea, the ice stream thins and the grounding line retreats but it remains grounded on the ridge. Therefore, before the 1940s, PIG probably experienced temporary periods of migration back and forth on the ridge in response to variable ocean conditions32.

When higher melt rates are applied for an extended period of time, to represent what may have happened during the 1940s El Niño event26,32, there is a rapid retreat down the retrograde slope facilitated by the merging of upstream cavities. Although we used a simple melt-rate parameterization, the initial behaviour of retreat, the speed at which it progresses and the final ungrounding of a pinning point above the ridge are all comparable with satellite observations and sediment records from the 1940s and 1970s26,27.

Our stability analysis suggests that by the early 1970s, when PIG had already started retreating from the ridge27, a threshold had been crossed, whereby its previous position could not be restored during subsequent cooler periods32. This irreversible phase came to a halt as the grounding line reached a new steady state on the next bed high point (Fig. 3). This location coincides with its early 1990s position, when PIG was grounded at an ice plain and had experienced a decrease in grounding-line ice flux1,49,50. During the suspected period of rapid retreat from the 1970s to the early 1990s, PIG was responsible for a third of the mass loss from West Antarctica and almost 13% of the overall AIS mass loss25. Despite its basin comprising of only 1.5% of the entire ice sheet area, PIG was the largest contributor to sea-level rise from the AIS during those years, adding 0.34 mm in total25.

Climate change is likely to cause further upstream migration of grounding lines of WAIS. In the Amundsen Sea, as local wind trends change in response to internal and external forcing29,51, this may deliver more warm water to the continental shelf30,31, leading to increased basal melt52 and ice-shelf thinning. Previous modelling studies of the behaviour of Amundsen Sea glaciers have suggested the existence of stability thresholds, which when crossed lead to irreversible mass loss at some point in the future9,16. This marine ice sheet instability is theoretically well understood13,14 and robustly replicated in numerical models8,9,16,53; however, the hypothesis has hitherto had little direct observational support. This is in part due to the long timescales involved and the sparsity of relevant observations.

Here, we have now shown that the recent observed grounding-line retreat of PIG, in the period from the 1940s to 1990s, was irreversible and thereby provided an observational validated example of the marine ice sheet instability. Our ice-flow model is based on the same physical assumptions used in previous future simulations9,16 and therefore this greatly strengthens our confidence in the capability of ice sheet models and their ability to simulate and predict highly nonlinear behaviours of large ice sheets. Furthermore, the results presented here are robust and insensitive to our choice of model parameters. These simulations suggest that the recent retreat phase of PIG may have been primarily internally driven, as opposed to external forced. While ocean-induced melt may have been the initial trigger, the retreat phase was driven by internal ice-dynamical processes leading to irrevocable loss of ice that could not be recovered by a reversal in external climatic condition. The implications for the future are clear: what has happened in the recent past, can happen again and, as predicted by ice-flow models, future ice loss from the WAIS may become self-sustaining, amplified and irreversible as the ice sheets enter unstable phases of retreat.

Methods

Ice-flow model

In this study we used the finite-element, vertically integrated ice-flow model Úa36,54 to solve the ice dynamics equations in the shallow ice-stream approximation (SSTREAM or SSA)55. The model has previously been used to study tipping points and drivers of retreat of PIG16,56, grounding-line stability and ice-shelf buttressing36,57,58 and in several intercomparison projects59,60,61.

The vertically integrated, or two horizontal dimension, momentum equations can be written in compact form as

$${\nabla }_{{xy}}\cdot \left({h\bf{R}}\right)-{{\bf{t}}}_{{\rm{bh}}}={{\rho}}_{{\rm{i}}}{gh}{\nabla }_{{xy}}{s}+\frac{1}{2}{g}{{h}}^{2}{\nabla }_{{xy}}{{\rho }}_{{\rm{i}}}$$
(1)

where h is the ice thickness, tbh is the horizontal component of the bed-tangential basal traction tb, \({{\rho }}_{{\rm{i}}}\) is the vertically averaged ice density, g is gravitational acceleration, s is the ice upper surface elevation and R is the resistive stress tensor defined as

$${\bf{R}}=\left(\begin{array}{cc}2{{\tau }}_{{xx}}+{\tau}_{{yy}} & {{\tau }}_{{xy}}\\ {{\tau }}_{{xy}} & 2{{\tau }}_{{yy}}+{{\tau }}_{xx}\end{array}\right)$$
(2)

and

$${\nabla }_{{xy}}={\left({\partial }_{{x}},{\partial }_{{y}}\right)}^{{T}}.$$
(3)

Here, \({{\tau }}_{{ij}}\) are the components of the deviatoric stress tensor. The relationship between deviatoric stresses \({\tau }_{ij}\) and strain rates \({{\epsilon }}_{{ij}}\) is given by Glen’s flow law

$${\dot{\epsilon}}_{ij}={A}{\tau}^{{n}-1}{\tau}_{ij},$$
(4)

where \({\rm{\tau }}\) is the second invariant of the deviatoric stress tensor

$${\tau }=\sqrt{{{\tau }}_{{ij}}{{\tau }}_{{ij}}/2},$$
(5)

A is a spatially varying ice rate factor determined using inverse methods and n = 3 is a creep exponent. In our main set of experiments, the basal traction is given by Weertman’s sliding law

$${{\bf{t}}}_{{\rm{b}}}={G}{{\beta }}^{2}{{\bf{v}}}_{{\rm{b}}},$$
(6)

where G is a floating mask, with G = 1 for grounded ice and G = 0 otherwise and vb is the horizontal component of the bed-tangential ice velocity. In equation (6), \({{\beta }}^{2}\) is given by

$${{\beta }}^{2}={{C}}^{-1/{m}}{\rm{|}}{{\bf{v}}}_{{\rm{b}}}{{\rm{|}}}^{1/{m}-1},$$
(7)

where C is a spatially varying slipperiness coefficient, determined using inverse methods, and m = 3, which gives a nonlinear viscous relationship. Downstream of the grounding line the slipperiness coefficient is set to a constant of C = 0.03 m yr−1 kPa−3, which allows the ice stream to advance forward. This constant is representative of upstream slipperiness values along the fast-flowing tributaries.

In two additional experiments, a different basal sliding setup was used. First, a downstream slipperiness coefficient of C = 0.01 m yr−1 kPa−3, representing a ‘stickier’ bed, was tested. Whilst in the second experiment, a modified power law was used for the basal traction47. This is given by

$${{\bf{t}}}_{{\rm{b}}}=\frac{{G}{{\beta }}^{2}{\rm{|}}{{\bf{v}}}_{{\rm{b}}}{\rm{|}}{{\mu }}_{{\rm{k}}}{N}}{{\left({\left({{\mu }}_{{\rm{k}}}{N}\right)}^{{m}}+{\left({G}{{\beta }}^{2}{\rm{|}}{{\bf{v}}}_{{\rm{b}}}{\rm{|}}\right)}^{{m}}\right)}^{1/{m}}}\frac{{{\bf{v}}}_{{\rm{b}}}}{{\rm{|}}{{\bf{v}}}_{{\rm{b}}}{\rm{|}}}$$
(8)

where μk is the coefficient of kinetic friction and is set to μk = 0.5.

Model domain and mesh

The model domain includes the grounded catchment of PIG (182,000 km2) and its floating ice shelf62 (Supplementary Fig. 1). The calving front is fixed throughout the study and corresponds approximately to the 2008/09 ice front, which is not far from its 1940s position63,64. For all experiments in this study, a Dirichlet boundary condition is imposed on the grounded portion of the boundary to set the velocity to zero along the ice divides and a Neumann boundary condition arising from ocean pressure is imposed along the ice front.

An irregular, triangular mesh was generated using MESH2D65 for the entire domain and consisted of 58,777 linear elements and 29,797 nodes. The mesh was refined for ice-shelf elements (1 km) and in areas of high strain rate and high strain rate gradients (0.7–1.5 km), whereas larger elements (10 km) were used for the slowest moving ice inland away from the main tributaries (Supplementary Fig. 2). This gave a mesh with minimum, median and maximum element sizes of 563, 1,311 and 11,330 m, respectively. For the control and warm experiments, a further grounding-line mesh adaption was applied to ensure that fine element sizes were used in a crucial transition area. Owing to computational and time limitations, no mesh adaption was used for the reversibility experiments.

Input data

This study aims to simulate the response of a 1940s PIG to a change in external forcing; however, with very little data available for that period we set up our model using present-day observations and then let the model evolve in time to get an approximate configuration for 1940. The bedrock topography, ice thickness, surface elevation and ice density were taken from BedMachine Antarctica v.2 (ref. 66). These datasets have a resolution of 500 m and nominal data of 2015. Some local adjustments were made to the ice-shelf thickness near the grounding line to ensure that the hydrostatic floating condition was met. As the BedMachine data represent a recent bed geometry, we also ran an additional experiment with a lower bed to test the effect of solid-earth feedbacks. The upper surface accumulation was given by the RACMO2.3p2 dataset67 and was averaged between 1979 and 2016.

Inversion

To initialize the model, we used present-day velocities from the MEaSUREs Annual Antarctic Ice Velocity Maps68,69 dataset to invert for the slipperiness parameter and the ice rate factor (Supplementary Fig. 3). For the inversion process, Úa minimizes a cost function containing a misfit and a regularization term, using the adjoint method and Tikhonov regularization, as has been done in previous studies70,71,72.

Melt-rate parameterization

The basal melt rate is given by a depth-dependent parameterization (Supplementary Figs. 5 and 6), similar to a previous study on PIG retreat8. Although this is a simple parameterization, it allows for conclusions to be made about the direct effect of basal melting. We also repeated our stability analysis using a different melt-rate parameterization that has been used in a previous model intercomparison project47. To ensure that the grounding-line retreat was not overestimated, we applied basal melting on mesh elements that are strictly downstream of the grounding line73. For the stability analysis, model simulations were run for hundreds of years until a steady state was reached. During these runs, to avoid unrealistic retreat along the southwest tributary, close to the model domain boundary, the basal melting was set to zero for elements in this region.