In this study we used the finite-element, vertically integrated ice-flow model Úa31 to solve the ice dynamics equations in the shallow ice-stream approximation (SSTREAM or SSA)35. The model has previously been used to study tipping points and drivers of retreat of Pine Island Glacier15,36, grounding-line stability and ice shelf buttressing31,37,38 and in a number of intercomparison projects39–41.

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

where *h* is the ice thickness, tbh is the horizontal component of the bed-tangential basal traction tb, \({{\rho }}_{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

$${R}\text{=}\left( \begin{array}{cc}\text{2}{\tau }_{xx}\text{+}{\tau }_{yy}& {\tau }_{xy}\\ {\tau }_{xy}& {\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 \(\dot{{\text{ϵ}}_{ij}}\) is given by Glen's flow law

$$\dot{{ϵ}_{ij}}=A{\tau }^{n-1}{\tau }_{ij}$$

4

where \(\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. For grounded ice, the basal traction is given by Weertman's sliding law

where **v**b is the horizontal component of the bed-tangential ice-velocity, *C* is a spatially varying slipperiness parameter, determined using inverse methods, and \(m=3\)which gives a non-linear viscous relationship. Downstream of the grounding line the slipperiness parameter is set to a constant of 0.03, which allows the ice stream to advance forward. This constant is representative of upstream slipperiness values along the fast-flowing tributaries.

**Model domain and mesh**

The model domain includes the grounded catchment of PIG (182,000 km2) and its floating ice shelf42 (Supplementary Fig. 1). The calving front position is fixed throughout the study and corresponds approximately to the 2008/09 ice front. 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 MESH2D43 for the entire domain, and consisted of 58777 linear elements and 29797 nodes. The mesh was refined for ice shelf elements (1 km) and in areas of high strain rate and high strain rate gradients (700–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 m, 1311 m and 11330 m respectively. For the control and warm experiments, a further grounding line mesh adaption was applied to ensure fine element sizes were used in a crucial transition area. Due 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, v244. 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 the hydrostatic floating condition was met. Upper surface accumulation was from RACMO2.3p245 and averaged between 1979 to 2016.

**Inversion**

To initialise the model, we used present day velocities from the MEaSUREs Annual Antarctic Ice Velocity Maps46,47 dataset to invert for the slipperiness parameter and the ice rate factor. For the inversion process, Úa minimises a cost function containing a misfit and a regularisation term, using the adjoint method and Tikhonov regularisation, as has been done in previous studies48–50. The inversion results are shown in Supplement Fig. 3.

**Melt rate parameterisation**

The basal melt rate is given by a depth dependent parameterisation (see Supplement Figs. 4 and 5), similar to previous study on Pine Island Glacier retreat7. Although this is a simple parameterisation, it allows for conclusions to be made about the direct effect of basal melting. To ensure the grounding line retreat was not overestimated, we applied basal melting on mesh elements that are strictly downstream of the grounding line51. For the stability analysis, model simulations were run for 100s 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.