Introduction
Several studies (e.g. Bassis and Walker, Reference Bassis and Walker2011; DeConto and Pollard, Reference DeConto and Pollard2016; Scambos and others, Reference Scambos2017) have previously suggested that the West Antarctic Ice Sheet (WAIS) might be vulnerable to a marine ice-cliff instability (MICI). The proposed instability mechanism involves the formation of high subaerial ice cliffs following a sudden loss of abutting ice shelves. Provided calving is a sufficiently strongly increasing function of cliff height, this might result in the calving rate exceeding ice flow velocities, causing further frontal retreat and exposing even higher ice cliffs. Currently, high-end estimates of sea-level projection in the IPCC Sixth Assessment Report (Meredith and others, Reference Meredith, Pörtner, Roberts, Masson-Delmotte, Zhai, Tignor, Poloczanska, Mintenbeck, Alegría, Nicolai, Okem, Petzold, Rama and Weyer2022) are based on model simulations with MICI scenarios where calving is prescribed as an increasing function of cliff height (DeConto and Pollard, Reference DeConto and Pollard2016). Those low-confidence, high-impact estimates predict a possible sea-level contribution in excess of one meter from the Antarctic Ice Sheet by the year 2100 (Fig. 5 of DeConto and Pollard (Reference DeConto and Pollard2016) and Table 4.3 of Meredith and others (Reference Meredith, Pörtner, Roberts, Masson-Delmotte, Zhai, Tignor, Poloczanska, Mintenbeck, Alegría, Nicolai, Okem, Petzold, Rama and Weyer2022), Figure 9.27 of Fox-Kemper and others (Reference Fox-Kemper, Masson-Delmotte, Zhai, Pirani, Connors, Péan, Berger, Caud, Chen, Goldfarb, Gomis, Huang, Leitzell, Lonnoy, Matthews, Maycock, Waterfield, Yelekçi, Yu and Zhou2023)). When included, this process could contribute to a sea-level rise that exceeds 1 m by 2100 relative to 1900 levels (Lee and others, Reference Lee2023). These studies, thus, imply that the MICI may result in ice loss considerably higher than the long-term potential of Marine Ice Sheet Instability (MISI).
The Antarctic Buttressing Model Intercomparison Project (ABUMIP) explored MISI impacts under sustained, complete loss of all Antarctic ice shelves an extreme and unrealistic scenario (Sun and others, Reference Sun2020). ABUMIP simulations suggest a maximum sea-level contribution of 
$1.05 \pm 0.4$ m from the Amundsen Sea sector over 500 years. These simulations also highlight the sensitivity of glacier dynamical responses to model numerics, particularly spatial resolution. Although the ice terminus of the Amundsen Sea glaciers is mainly characterised by floating ice shelves with ice cliff height less than 
$80\,\mathrm{m}$, there are nevertheless some frontal regions with ice cliffs above this threshold (Fig. 1), making those regions therefore potentially susceptible to MICI.
(a) Surface elevation distribution for Amundsen Sea glaciers, or cliff height assuming instantaneous downstream glacier collapse. The grey solid line denotes grounding lines, while the dark blue and red dotted lines indicate calving front locations: dark blue dots represent sites with current cliff heights below the 80 m threshold, and red dots indicate locations where cliff heights exceed 80 m. (b) Velocity normal to the ice front, calving rate and retreat rate at present-day calving fronts. The inset shows the cumulative distance between the calving front and its westernmost position. The straight dash lines indicate the profiles presented in Fig. 3.
In contrast to MISI, which is well understood theoretically (Schoof, Reference Schoof2007a, Reference Schoof2007b) and numerically verified in all modern ice-flow models, the concept of the MICI and its relevance for modern ice sheets is less well confirmed. The presence of high ice cliffs and the large magnitude of predicted ice loss due to MICI in the near future justifies further study of the implication of cliff-height-dependent calving rate parametrisations and the numerical robustness of those results.
For a given calving law, numerical models of ice flow can be used to simulate frontal retreat rates and the evolution of cliff height due to ice deformation. Ideally, such models should accurately capture the interplay between these competing processes, thereby producing numerically robust and physically realistic estimates of frontal retreat rates. However, modelling these processes presents significant numerical challenges, particularly in geographically complex regions such as the WAIS. Specifically, the dynamics of grounding lines are highly sensitive to mesh resolution, and parametrisations of calving rates based on cliff height demand that ice thickness evolution just upstream of calving fronts is resolved with sufficiently high spatial detail. The required horizontal numerical resolution can be determined through convergence experiments, but is typically expected to be on the order of a single ice thickness (approximately one kilometre) or less (Gladstone and others, Reference Gladstone, Lee, Vieli and Payne2010; Durand and others, Reference Durand, Gagliardini, Favier, Zwinger and le Meur2011; Houriez and others, Reference Houriez2025). Accurately representing changes in buttressing, as well as basal and lateral drag associated with frontal retreat, similarly requires high spatial resolution. Previous numerical studies employing cliff-height-dependent calving rate parametrisations have often utilised coarse grid resolutions of around 10 km (e.g. DeConto and Pollard, Reference DeConto and Pollard2016; Kopp and others, Reference Kopp2017; DeConto and others, Reference DeConto2021) which may be insufficient for resolving the dynamics of migrating grounding lines and the evolution of thickness profiles near actively calving fronts. Consequently, it cannot be ruled out that these studies may not have achieved the spatial resolution necessary to draw robust conclusions regarding frontal retreat rates. Furthermore, several of these studies have employed a parametrisation of grounding line flux intended to reduce the sensitivity of model results to resolution (Pollard and DeConto, Reference Pollard and DeConto2012; DeConto and others, Reference DeConto2021). However, applying the flux formula for the grounding line as a boundary condition in 2D models violates many of the assumptions under which the approximation was developed (Schoof, Reference Schoof2007a), and has been demonstrated to produce severely inaccurate results when applied to current grounding lines of the Antarctic Ice Sheet (Reese and others, Reference Reese, Winkelmann and Gudmundsson2018).
Recently, Morlighem and others (Reference Morlighem2024) studied the implications of a cliff-height-dependent calving rate parametrisation proposed by Crawford and others (Reference Crawford, Benn, Todd, Åström, Bassis and Zwinger2021). Morlighem and others (Reference Morlighem2024) concluded that the proposed calving rate parametrisation was unlikely to cause significant additional ice loss from the Thwaites glacier, or to give rise to an unstable frontal retreat. The calving rate parametrisation proposed by Crawford and others (Reference Crawford, Benn, Todd, Åström, Bassis and Zwinger2021) predicts that the calving rates are effectively zero for subaerial ice cliffs less than 136 m high, and then rapidly increase with subaerial cliff height. The calving rate parametrisation used in DeConto and Pollard (Reference DeConto and Pollard2016) assumes zero calving rate for ice cliff heights below 80 m, then a linear increase in calving rate with cliff height, reaching 
$\textstyle3\,\mathrm{km}\,\mathrm a^{-1}$ for 100 meter high cliffs. For cliffs above 100 m the calving rate is 
$3\,\mathrm{km}\,\mathrm{a}^{-1}$. The different cliff-height thresholds might explain why DeConto and Pollard (Reference DeConto and Pollard2016) observed a rapid frontal retreat in their numerical simulations. However, recently Rosier and others (Reference Rosier, Gudmundsson, Jenkins and Naughten2025) modelled the evolution of the WAIS using the calving rate parametrisations previously proposed by Pollard and others (Reference Pollard, DeConto and Alley2015), including the cliff-height-dependent formulation, and did not observe periods of accelerated frontal retreat under RCP8.5 or Paris 2C scenarios. However, due to the inclusion of several different calving rate parametrisations in Rosier and others (Reference Rosier, Gudmundsson, Jenkins and Naughten2025), the role of each individual calving rate parametrisation in controlling frontal retreat is difficult to assess. Furthermore, in Rosier and others (Reference Rosier, Gudmundsson, Jenkins and Naughten2025), the conditions required for the cliff-height-dependent calving law (e.g. DeConto and Pollard, Reference DeConto and Pollard2016) to be activated, i.e. the initial formation of sufficiently high cliffs, depended on various other competing factors. Furthermore, these previous studies (Morlighem and others, Reference Morlighem2024; Rosier and others, Reference Rosier, Gudmundsson, Jenkins and Naughten2025) focused on fast-flowing glaciers (Pine Island Glacier, Thwaites glacier), and did not include the slow-flowing ice front with high ice cliffs, e.g. Ellsworth Land (Fig. 1).
Here, we model the evolution of the WAIS, using the cliff-height-dependent calving rate parametrisation of DeConto and Pollard (Reference DeConto and Pollard2016), at various spatial resolutions. We calculate the grounding line flux by solving the momentum equations instead of using a flux parametrisation for the grounding line ice discharge (Schoof, Reference Schoof2007a, Reference Schoof2007b). We show that using resolution coarser than 2 km results in inaccurate numerical results. Application of the same calving law as used in DeConto and Pollard (Reference DeConto and Pollard2016) in simulations with finer spatial resolution does not lead to a self-enhancing frontal retreat of fast-flowing marine-terminating glaciers. Contrary to empirical observations, the direct application of the DeConto and Pollard (Reference DeConto and Pollard2016) calving rate parametrisation induces frontal retreat in several present-day calving fronts within slow-flowing regions of the WAIS, where no such retreat is currently observed.
Methodology
In order to study the implication of the cliff-height-dependent calving rate parametrisation that is used by (e.g. DeConto and Pollard, Reference DeConto and Pollard2016) for the Amundsen Sea glaciers, we apply Eq. (1), in isolation, and conduct several numerical experiments with all ice-shelves removed at the beginning of the experiments.
Numerical experiments
The high-end calving rate parametrisations investigated in this study include: 1) sustained collapse of ice shelves (ABUMIP) (Sun and others, Reference Sun2020), and 2) cliff-height-dependent calving parametrisation of DeConto and Pollard (Reference DeConto and Pollard2016). The ABUMIP calving scenario can be thought of as a ‘position-based’ calving rate parametrisation, where the position of the calving fronts coincides at all times with the grounding-line position, and all floating ice shelves are continuously removed. The cliff-height-dependent calving rate parametrisation, as proposed by Pollard and others (Reference Pollard, DeConto and Alley2015) and DeConto and Pollard (Reference DeConto and Pollard2016), provides the calving rate, 
$c$, as a function of ice-cliff height, 
$h_c$, as follows:
\begin{equation}
c=
\begin{cases}
0
& h_c \le 80 \\
\frac{3000}{20} ( h_c - 80)
& 80 \lt h_c \lt 100\\
3000
& h_c \ge 100
\end{cases}
\end{equation}where 
$c$ is the calving rate with unit 
$\mathrm{m}\,\mathrm{a^{-1}}$ and 
$h_c$ is the cliff height in meters.
Combining the use of the calving rate parametrisation Eq. (1) with collapse of ice shelves (ABUMIP), gives rise to several possible numerical experiments. We explore the following scenarios:
- Reference:
A standard transient run with fixed present-day calving fronts, and basal melting parametrisation from the ISMIP6 protocol (Jourdain and others, Reference Jourdain2020), is defined as the reference simulation. Here, we use the ‘PIGL’ approach, where basal melting of the ice shelf is a quadratic function of observed temperature, and the tuning parameters are calibrated by observed basal melt rates of the PIG ice shelf.
- DP:
A cliff-height-dependent calving rate parametrisation given by Eq. (1), with no basal melting and ice shelves are allowed to form.
- ABUMIP:
All ice shelves are removed at all times, and DP calving rate parametrisation Eq. (1) is not applied. This scenario from the ABUMIP protocol could represent an extreme case of hydrofracture, resulting in the collapse of any floating ice shelf.
- DP+ABUMIP:
The cliff-height-dependent calving rate parametrisation, Eq. (1), is used and, additionally, all ice shelves are removed at all times. This scenario explores the potential for cliff failure when high cliffs are exposed following ice shelf collapse due to processes like hydrofracture.
- DP+Reference:
The DP calving rate parametrisation given by Eq. (1) is used, with basal melting based on the ISMIP6 protocol, and ice shelves are allowed to form (This is an ISMIP6 experiment, with the DP calving rate parametrisation (Eq. 1)).
- DP+ABUMIP0:
All ice shelves are removed only at the beginning of the experiment, and the DP calving rate parametrisation, Eq. (1), applied to all calving fronts.
Ice-flow model
We used the version 2024a of the ice-flow model Úa to conduct this study (Gudmundsson, Reference Gudmundsson2024). Úa is a finite element model that solves the shallow ice stream equations (SSA or SSTREAM) on an irregular triangular mesh. Various criteria can be applied to automatically adapt the mesh for resolution requirements at defined locations, such as near the grounding line, calving front or locations experiencing high strain rates. The ice-flow model solves the momentum and mass conservation equations simultaneously using implicit time integration with respect to both ice velocities and ice thickness. In all numerical results presented, the computational domain was discretised using linear 3-node triangular elements.
Calving implementation
In the ice-flow model Úa, calving is implemented using an augmented level-set method based on a variational principle (Touré and Soulaïmani, Reference Touré and Soulaïmani2016; Gudmundsson, Reference Gudmundsson2024). The level-set method is frequently used in computational fluid dynamics to describe the evolution of internal interfaces and moving boundaries, and has previously been used in a number of glaciological studies (e.g. Bondzio and others, Reference Bondzio2016; Hossain and others, Reference Hossain, Pimentel and Stockie2020; Cheng and others, Reference Cheng, Morlighem and Gudmundsson2024). The level-set function, 
$\varphi : \mathbb{R}^2 \times \mathbb{R} \to \mathbb{R}$, is defined by
\begin{equation}
\varphi(\boldsymbol{x}_c(t),t)=0 \; ,
\end{equation}where 
$\boldsymbol{x}_c$ is a set of all points forming the calving fronts within the computational domain. The level-set function 
$\varphi$ is defined to be positive on one side of the calving front, and negative on the other side. Since, by definition, the value of the level-set function is always equal to zero at the calving front over time, the corresponding total time derivative is also equal to zero, that is
\begin{equation}
\left . \frac{\mathrm{d}}{\mathrm{d}t} \varphi(\boldsymbol{x}(t),t) \right |_{\boldsymbol{x}=\boldsymbol{x}_c} = \partial_t \varphi(\boldsymbol{x}_c,t) + \boldsymbol{u}_c \cdot \nabla \varphi =0
\end{equation}where
\begin{equation*}
\boldsymbol{u}_c = \partial_t \boldsymbol{x}_c \; .
\end{equation*} The velocity, 
$\boldsymbol{u}_c$, of the calving front is equal to the difference between the material velocity, 
$\boldsymbol{v}$, of ice at the calving front and the calving rate 
$\boldsymbol{c}$, that is
\begin{equation*}\boldsymbol{u}_c=\boldsymbol{v}-\boldsymbol{c} \; .\end{equation*} We therefore arrive at the following evolutionary equation for 
$\varphi$,
\begin{equation}
\partial_t \varphi + \boldsymbol{v} \cdot \nabla \varphi = \boldsymbol{c} \cdot \nabla \varphi \; .
\end{equation} Here, the ice velocity, 
$\boldsymbol{v}$, is calculated by an ice-flow model, and the calving velocity, 
$\mathbf{c}$, is provided by a calving rate parametrisation. The calving velocity vector, 
$\boldsymbol{c}$, is always normal to the calving front itself, and we can write
\begin{equation}
\partial_t \varphi + \boldsymbol{v} \cdot \nabla \varphi = -c \, \|\nabla \varphi \| \; ,
\end{equation}where the scalar 
$c$ is the calving rate, with 
$\boldsymbol{c} =\hat{\boldsymbol{n}} c$, where the normal, 
$\hat{\boldsymbol{n}}$ to the calving front is
\begin{equation}
\hat{\boldsymbol{n}}=-\frac{\nabla \varphi}{\| \nabla \varphi\|}\; .
\end{equation} The sign convention used in the definition of the normal in Eq. (6) is introduced in the anticipation that 
$\varphi$ will be defined as a decreasing function of distance as we travel across the calving front, from the ice-covered region to the ice-free region, with the normal 
$\hat{\boldsymbol{n}}$ pointing outwards.
Due to the dependency of 
$\boldsymbol{c}$ on 
$\nabla \phi$, Eq. (4) is a non-linear hyperbolic equation. The ice and calving velocities, i.e. 
$\boldsymbol{v}$ and 
$\boldsymbol{c}$, are functions of 
$x$ and 
$y$. It is well known that non-linear hyperbolic equations such as Eq. (4) can give rise to discontinuities over time, whenever the velocity, i.e. 
$\boldsymbol{u}_c=\boldsymbol{v}-\boldsymbol{c}$, is spatially variable. Numerically, this might result in overshoots and undershoots around the developing region of discontinuity. In extreme cases, this could cause the sign of the level-set function, 
$\varphi$, to change, resulting in a shift of the calving front position and potentially erroneously to the formation of new calving fronts. Several methods have been proposed in the literature to deal with this problem (e.g. Gomes and Faugeras, Reference Gomes and Faugeras2000). These include geometrical re-initialisation whereby 
$\varphi$ is periodically reset to be approximately a signed distance function, in which case 
$\|\nabla \varphi \| \approx 1$. However, geometrical re-initialisation can cause a shift in the position of the zero level, and results in a solution algorithm that is difficult to describe fully in terms of the partial differential equation being solved. Another common approach is to modify the level-set equation by adding a judiciously selected term to the equation (Touré and Soulaïmani, Reference Touré and Soulaïmani2016). In Úa, an augmented level-set equation with an additional non-linear forward-and-backward (FAB) diffusion term, derived from a variational principle, is used to ensure that the position of the zero level is accurately modelled, while at the same time the development of any shocks and discontinuities is suppressed. The augmented level-set equation takes the form
\begin{equation}
\partial_t \varphi + (\boldsymbol{v}-\boldsymbol{c}) \cdot \nabla \varphi - \nabla \cdot ( \kappa \nabla \varphi) = 0 \; .
\end{equation} The additional FAB diffusion term, i.e. 
$\nabla \cdot ( \kappa \nabla \varphi)$, results from minimising the functional, 
$\mathcal{P}$,
\begin{equation}
\mathcal{P} = \frac{1}{p q} \int_{\mathcal{A}} \left ( \|\nabla \varphi \|^q-1 \right )^p \, d\mathcal{A}
\end{equation}with respect to 
$\varphi$, where 
$p$ and 
$q$ are integers, and 
$\mathcal{A}$ is the computational domain. The minimum of 
$\mathcal{P}$ is obtained when 
$\varphi$ is a signed distance function, and the addition of the FAB term ensures that the solution to the augmented level-set equation, Eq. (7), is approximately a signed distance function. Further details are provided in the Úa Compendium—a technical manual for the Úa ice-flow model forming a part of the Úa distribution (Gudmundsson, Reference Gudmundsson2024)—where it is shown that the resulting FAB diffusion coefficient, which can be both positive and negative, has for the functional 
$\mathcal{P}$, the form
\begin{equation*}
\kappa=\mu \, \left (\|\nabla \varphi \|^q-1 \right )^{p-1} \|\nabla \varphi \|^{q-2} \; .
\end{equation*} The augmented level-set equation is solved using the Newton-Raphson method implicitly with respect to time with the 
$\theta$ time integration method combined with the consistent streamline-upwind Petrow–Galerkin finite-element formulation. One of several advantages of this approach is that no geometrical re-initialisation is required. Furthermore, the initial level-set function at the start of a simulation can be found by solving
\begin{equation}
\nabla \cdot ( \kappa \nabla \varphi) = 0 \; .
\end{equation} The solution of Eq. (9) will be an optimal approximation to a signed distance function for a given set of boundary conditions. In the numerical experiments shown here, we set 
$p=2$, 
$q=2$ and 
$\mu=0.2$. It can be shown that those 
$p$ and 
$q$ values result in a bounded Hessian as 
$\| \nabla \varphi \| \to 0$ (Gudmundsson, Reference Gudmundsson2024). It can furthermore be shown that 
$\mu$ can be expected to be on the order of unity, and numerical experiments where modelled calving rates were compared with analytical solutions suggested 
$\mu \approx 0.2$ (Gudmundsson, Reference Gudmundsson2024). Verification of this approach and the comparison with analytical test cases can be found in the Úa Compendium, which is distributed as part of the Úa source code and is freely available, and in Cheng and others (Reference Cheng, Morlighem and Gudmundsson2024).
Model configurations
The simulations of the Amundsen Sea glaciers are initialised with present-day geometry from Bedmachine Antarctica v2 (Morlighem and others, Reference Morlighem2020) and surface velocities (Rignot and others, Reference Rignot, Mouginot and Scheuchl2011). We set the Glen’s flow law exponent, 
$n$, and the Weertman sliding law exponent, 
$m$, to 
$n=3$ and 
$m=3$, respectively. Using a Bayesian inversion approach, we minimise the conditional probability of the system state for the given surface velocity field, where the system state are the spatial distributions of the logarithms of the rate factor, i.e. 
$\log_{10}(A)$, in Glen’s flow law, and the basal slipperiness coefficient, 
$\log_{10}(C)$, in Weertman’s sliding law (Fig. A1). The prior covariance functions for 
$\log_{10}(A)$ and 
$\log_{10}(C)$ are modelled as Matérn covariances, and the hyperparameters of the covariance function are determined through an 
$L$-curve analysis. This is a standard approach that has been followed in several previous publications (e.g Barnes and others, Reference Barnes, Dias dos Santos, Goldberg, Gudmundsson, Morlighem and De Rydt2021; De Rydt and others, Reference De Rydt, Reese, Paolo and Gudmundsson2021; Barnes and Gudmundsson, Reference Barnes and Gudmundsson2022; Gudmundsson and others, Reference Gudmundsson, Barnes, Goldberg and Morlighem2023; Reed and others, Reference Reed, Green, Jenkins and Gudmundsson2024).
We conducted a mesh sensitivity analysis to explore the sensitivity of the results to mesh resolution and to ensure that the model accurately captures the migration of the grounding line and the calving front. The sensitivity to mesh resolution was examined by repeating the ABUMIP experiment with eight different meshes that ranged from 32 km resolution to 250 m resolution. The coarser meshes had uniform element sizes of 32 km, 16 and 8 km. Finer meshes were then created using local mesh refinement of the 8 km global mesh near the grounding lines and the calving fronts. The criteria for mesh refinement were distance to the grounding lines, as in the previous study of Cornford and others (Reference Cornford, Martin, Lee, Payne and Ng2016), where the mesh size 
$\Delta x^l=2^{-l}\times 8~\mathrm{km}$, where 
$l\in N$ and 
$0 \lt l \lt L$, with the maximum value 
$0 \lt L\leq 5$. The fine resolution covers the distance of 
$8\Delta x^l$ from the grounding line from both the grounding side and the floating side. The mesh sizes are adapted with grounding line and calving front migration in transient runs. For each mesh, we base the initial geometry on the Bedmachine Antarctica v2 data set. Thus, as the mesh is refined, we are able to resolve increasingly finer details of the bed geometry.
Results
Resolution dependence
Our initial numerical experiments explore the sensitivity of modelled ice-dynamics to the spatial resolution of the computational domain. As described above, we used a wide range of meshes with varying spatial resolution ranging from 32 km down to about 250 m around all migrating grounding lines and calving fronts.
Modelled grounding-line migration and change in ice volume above flotation (sea-level equivalence) of the ABUMIP experiment for different mesh resolutions are presented in Figure 2. With mesh resolutions coarser than 2 km, the dynamic response of Amundsen Sea glaciers to the loss of ice shelves is limited and no area-wide collapse of the region results within the modelling period of 250 years. However, simulations with mesh resolutions finer than 
$2~\mathrm{km}$ all show a persistent frontal retreat, which, over time, leads to a collapse of the whole glaciated basin. We did not conduct reversibility experiments, as done in Reed and others (Reference Reed, Green, Jenkins and Gudmundsson2024) and Hill and others (Reference Hill, Gudmundsson and Chandler2024), to establish if the retreat is related to an instability, but we surmise that the retreat is related to an unstable grounding line retreat, i.e. the MISI. The calving front migration of the DP experiment for different mesh resolutions is presented in Figure A2. While calving rate depends solely on ice cliff height, solving ice dynamics with different mesh resolutions results in different extent and routes of calving front retreat, and simulations with mesh resolutions finer than 
$2~\mathrm{km}$ show a persistent frontal retreat. We, thus, suggest that a local mesh resolution of at least 
$1~\mathrm{km}$ is required around evolving grounding lines and calving fronts for the numerical model to correctly capture the dynamic evolution of the region. Further mesh refinement around the migrating grounding lines using 500 and 250 m resolution had modest additional impact on the numerical results (see Fig. 2). In conclusion, a mesh where the resolution is everywhere at least 8 km, and down to about 500 m around all evolving grounding lines and calving fronts, is sufficient to capture the transient evolution of the Amundsen sea embayment. The requirement for mesh resolution on the order of one ice thickness to faithfully simulate marine ice sheet dynamics is well established (e.g. Pattyn and others, Reference Pattyn2012, Reference Pattyn2013; Cornford and others, Reference Cornford, Martin, Lee, Payne and Ng2016). Nevertheless, our results highlight particularly pronounced differences in transient modelled evolution: no collapse is predicted for mesh sizes of 2 km or greater, while a basin-wide collapse emerges as soon as the resolution is refined below 2 km in the vicinity of all migrating grounding lines (see Fig. 2).
Grounding line migration and change of volume above flotation (sea-level equivalence) for the ice shelf collapsing experiment (ABUMIP) for different mesh refinement strategies. The legend on the right panel applies to both plots.
Calving rate parametrisations
After having established the required mesh resolution, we now explore the implications of the cliff-height-dependent calving rate parametrisation given by Eq. (1).
We conduct transient experiments described in the method section on adaptive meshes for 100 years. An average surface mass balance over the period 1995–2014 is calculated from RACMO2.3 (van Wessem and others, Reference van Wessem2018) and is constant in time in our simulations.
The DP calving rate parametrisation, as expressed by Eq. (1), stipulates that calving rates are zero for subaerial cliff heights of 80 metres or less, and increase to 
$3\,\mathrm{km}\,\mathrm{a}^{-1}$ for cliffs exceeding 100 metres. Presently, subaerial cliffs taller than 80 metres are relatively common across the Antarctic Ice Sheet, with several instances depicted in Figure 1 (ice front positions highlighted in red). Strict adherence to the DP calving rate parametrisation would therefore predict calving rates in these regions of several hundred metres to a few kilometres per year. However, in reality, there is no observational evidence of ongoing frontal retreat in those regions. For example, at the grounded ice front north of the Pine Island Ice Shelf, cliff heights exceed 80 metres, yet ice flow speeds are approximately one metre per year with no evidence of ongoing frontal retreat (see Fig. 1). One possible explanation is that the Bedmachine ice-thickness dataset may contain substantial errors in these areas, particularly if such errors affect the determination of bed elevation relative to sea level. Nevertheless, the direct application of the DP calving rate parametrisation, in conjunction with the Bedmachine dataset—currently the most comprehensive and accurate representation of Antarctic Ice Sheet geometry—results in predictions that are inconsistent with observed calving front behaviour. In contrast, fast-flowing ice shelves such as PIG have cliff heights below 80 meters at the calving front. Without other calving mechanisms such as hydrofracturing, the calving rate parametrisation yields zero calving rates, resulting in negative retreat rates or ice front advance (Fig. 1b) of these glaciers, which also conflicts with observations (Joughin and others, Reference Joughin, Shapero, Smith, Dutrieux and Barham2021; Greene and others, Reference Greene, Gardner, Schlegel and Fraser2022). Consequently, we are compelled to conclude that the DP calving rate parametrisation, as currently formulated in Eq. (1), is not valid for the Antarctic Ice Sheet in its current state and should not be applied without modification.
Modelled differences in calving front and grounding line positions generally grow over time as seen in Figs. 3, A3, A4 and A5. In the Reference run, the calving fronts are, by definition, not allowed to migrate, and a modest shift in grounding lines is observed (Figs. 3 and A5). In the DP experiment, the DP calving rate parametrisation, Eq. (1), is applied without any further newly formed ice shelves being removed. Here, calving fronts of all ice shelves advance, while at the same time calving fronts of several other slow-flowing areas retreat (Figs. 3 and A3). The frontal advance in the DP experiment persists (Fig. 3), while areas north of Pine Island ice shelf and nearby coastlines experience steady deglaciation due to the combination of slow flow and the presence of high cliffs (Figs. A3, A4 and A5). The DP+Reference experiment, where we apply ocean-induced basal melting and the DP calving rate parametrisation, Eq. (1), over floating ice shelves is observed the same locations of ice front retreat with the DP experiment. Employing the DP calving rate parametrisation, Eq. (1), and additionally systematically removing all ice shelves at the initiation of the simulation and as they form (i.e. experiment DP+ABUMIP), results in pronounced frontal retreat (Figs. 3, A3, A4 and A5). This frontal retreat is evident in several regions characterised by slower ice flow, such as the grounding ice front north of Pine Island ice shelf, north of Cosgrove ice shelf, the coastal area between Pine Island and Thwaites glaciers, and the region south of Thwaites ice shelf extending towards Crosson ice shelf—areas where no significant frontal retreat is currently observed. Combining the DP calving rate parametrisation, Eq. 1, with the instantaneous removal of all ice shelves at all times, i.e. Experiment DP+ABUMIP, leads to complete deglaciation of the entire area within the computational domain in about 100 years (Figs. 3 and A5).
Flow line profile of ice shelf evolution over 100 years. The lines are annotated as dashed lines in Figure 1 for the ice rise near Venable (left panel) and Pine Island Glacier (right panel).
Experiment DP+ABUMIP0 gives rise to complex patterns of migrating calving fronts and grounding lines. In most cases, the initial removal of all ice shelves results in the formation of 80 m or higher ice cliffs. In several areas, the calving rates, as calculated by Eq. (1), are then found to exceed the modelled forward velocity of the ice, resulting in a frontal retreat. However, for the fast-flowing Pine Island and Thwaites glaciers, ice deformation rapidly reduces cliff heights to below 80 m (within a time span of a few years), after which the calving fronts start to advance, while the grounding lines continue to retreat (Figs. 3, A4 and A5).
We note that the modelled mass loss (Fig. 4) is not a prediction of future mass loss, but an illustration of the impact of different processes. In all experiments, the Amundsen Sea glaciers lose ice and sea-level rises over time, but the ice loss varies by orders of magnitude between different experiments. In the reference experiment, Reference, using present-day ocean forcing induced basal melting and constant calving front location, modelled mass loss of the Amundsen Sea glaciers equates to approximately 
$30~\mathrm{mm}$ to sea-level rise over 100 years. Although no particular attempt was made to ensure that the reference run replicated current observations in any detail, the average contribution to sea level 
$0.3\,\mathrm{mm}\,\mathrm a^{-1}$, or similar to the current estimate of sea-level rise from West Antarctica of about 
$0.24\,\mathrm{mm}\,\mathrm a^{-1}$ (Otosaka and others, Reference Otosaka2023). This agreement between initial modelled and observed mass loss is presumably in parts due to our initialisation procedure, where we ensure that modelled ice velocities are in good agreement with observations, but also possibly a somewhat fortuitous consequence of the application of the ISMIP6 protocol.
Sea-level contribution for different combinations of ice shelf removal and use of a cliff-height-dependent DP calving rate parametrisation Eq. (1) suggested in DeConto and Pollard (Reference DeConto and Pollard2016). The experiments are defined in the main text. Note that this is not a prediction for future sea-level rise from glaciers of the Amundsen Sea region, but rather a sensitivity study aimed at exploring the consequences of applying different calving rate parametrisations in a numerical study.
In Experiment ABUMIP, simulating an immediate collapse of ice shelves leads to a much higher sea-level rise of about 80 mm over 100 years (red line in Fig. 4) compared to the reference run. Applying the DP calving rate parametrisation Eq. (1), both with and without ocean-induced melt (Experiments DP and DP+Reference), results in an even greater mass loss of about 400 mm over 100 years (yellow and blue lines in Fig. 4). However, this approach causes frontal retreat in slow-flowing areas where none is currently observed, and we therefore judge this result to contradict observational evidence.
Combining the use of the DP calving rate parametrisation Eq. (1) with initial instantaneous and complete collapse of ice shelves, Experiment DP+ABUMIP0, causes about 
$630\,\mathrm{mm}$ rise in sea level over 100 years (light-red curve in Fig. 4). Compared to the reference run, Experiment Reference, this corresponds to about 20-fold increase in the mean annual rate of sea-level rise. Interestingly, despite the fact that most calving fronts of fast-flowing areas advance as the ice-cliff heights are reduced below 80 m due to ice deformation, the grounding lines mostly continue to retreat (see Fig. A5), as do calving fronts across slow-flowing areas. The mass loss in the experiment is, thus, not primarily driven by calving front retreat of the fast-flowing Pine Island and Thwaites glaciers, but by deglaciation of Ellsworth Land (Fig. 1).
Finally, applying the DP calving rate parametrisation, Eq. (1), to all ice fronts, and removing all ice shelves at all times, Experiment DP+ABUMIP, causes a complete deglaciation of the whole computational domain (see Fig. A5), resulting in about 
$1.3\,\mathrm{m}$ sea-level rise from the Amundsen Sea region.
In summary, the application of the DP calving rate parametrisation given by Eq. (1) has decisive effect on the modelled ice loss. However, contrary to remote-sensing observations, frontal retreat is mostly found over slow-flowing areas where no retreat is currently observed. Initially, some parts of the calving fronts of the fast-flowing Thwaites and Pine Island glaciers retreat following the removal of the ice shelves. However, this initial retreat quickly reverses to persistent advance for the remainder of the simulation as the ice fronts of those areas are reduced in height due to ice deformation.
Discussion
DeConto and Pollard (Reference DeConto and Pollard2016) argued that sudden loss of protective ice shelves around Antarctica, caused, for example, by hydrofracturing, and the resulting formation of high subaerial cliff faces, may lead to a large-scale and self-perpetuating frontal retreat of fast-flowing marine-terminating glaciers (DeConto and Pollard, Reference DeConto and Pollard2016). As mentioned above, this process—the MICI—is based on the following key idea: Provided calving rates are increasing functions of subaerial cliff height, initial formation of sufficiently high subaerial cliff fronts may lead to calving rates exceeding the forward flow speeds of the ice. Consequently, the calving fronts will retreat, potentially leading to the formation of even higher cliffs and further increases in calving rates. The result is a self-sustained and accelerated frontal retreat.
This conceptual picture of the MICI assumes that the initial imbalance between calving and ice velocity, with calving rate exceeding the ice velocity, will persist over time. However, there are several processes that might affect this situation. Ice deformation can reduce cliff height and consequently decrease calving rate, assuming that calving rates are increasing functions of subaerial cliff height. The rate of reduction in cliff height due to ice deformation is expected to be a strong function of cliff height itself. The stress boundary condition at a calving front implies that the deviatoric horizontal stresses are proportional to subaerial cliff height, causing the ice to thin at a rate proportional to the 
$n$-th power of ice thickness, where 
$n$ is generally thought to be around 3 to 4 (Weertman, Reference Weertman1957). The resulting frontal retreat rate is, thus, the result of a balance between two different processes, both dependent on subaerial cliff height, but acting in opposite directions. Additionally, the ice loss through calving can cause reduction in ice-shelf buttressing, and with it a further increase in ice velocities (Joughin and others, Reference Joughin, Shapero, Smith, Dutrieux and Barham2021; Sun and Gudmundsson, Reference Sun and Gudmundsson2023). It is likely that the balance between these difference processes will be situation dependent.
Here, we present a numerical study of the evolution of the Amundsen Sea glaciers using a cliff-height-dependent DP calving rate parametrisation, Eq. (1). This calving rate parametrisation was previously used in ice-flow studies that resulted in rapid retreat of calving fronts (DeConto and Pollard, Reference DeConto and Pollard2016). Our results show how termini heights of fast-flowing glaciers, such as Pine Island and Thwaites glaciers, are quickly reduced in height following an initial ice shelf removal. This results in the calving fronts of these glaciers advancing rather than retreating, despite the initial cliff heights being above 100 m. Therefore, for those glaciers, we do not observe persistent frontal retreat. However, persistent frontal retreat is found when all ice shelves are continuously removed while applying the calving rate parametrisation.
We find that numerically robust estimates of the transient evolution of calving fronts and grounding lines require mesh resolution of less than 2 km around all migrating grounding lines and calving fronts. Some of the differences between our results and those of previous studies might be related to the lack of sufficient mesh resolution. The study of DeConto and Pollard (Reference DeConto and Pollard2016); DeConto and others (Reference DeConto2021), for example, appears to have used a uniform mesh resolution of 10 km. In our simulation, such a resolution does not simulate the ice-sheet collapse observed with a resolution of less than 2 km. This finding on mesh resolution disagrees with the mesh convergence study of DeConto and others (Reference DeConto2021), where mesh sizes of 1–10 km produced similar sea-level contributions from Pine Island and Thwaites glaciers over 100 years. A further source of difference between our results and those of (Pollard and DeConto, Reference Pollard and DeConto2012) could be their use of a flux formula (Schoof, Reference Schoof2007a, Reference Schoof2007b) for calculating ice flux across grounding lines. While the flux formula of Schoof (Reference Schoof2007a, Reference Schoof2007b) provides an accurate estimation of ice flux for unbuttressed grounding lines, for which the formula was derived, applications of this same formula for buttressed ice shelves of Antarctica have been shown to produce inaccurate and sometimes physically implausible results (Reese and others, Reference Reese, Winkelmann and Gudmundsson2018).
The yield strength of ice has been proposed to be between 0.1-1 mPa, and varies with individual glaciers (150 kPa in Bassis and Ultee (Reference Bassis and Ultee2019); 50-300 kPa in Ultee and Bassis (Reference Ultee and Bassis2016); 500 kPa-1 mPa in Crawford and others (Reference Crawford, Benn, Todd, Åström, Bassis and Zwinger2021)). The yield strength is then often transformed to a threshold cliff height of 80-540 m with assumptions of physical parameters in flow-line models, or by observed calving front migration in certain periods (DeConto and Pollard, Reference DeConto and Pollard2016; Clerc and others, Reference Clerc, Minchew and Behn2019; Parizek and others, Reference Parizek2019; Crawford and others, Reference Crawford, Benn, Todd, Åström, Bassis and Zwinger2021; DeConto and others, Reference DeConto2021; Needell and Holschuh, Reference Needell and Holschuh2023). While the derived cliff-height dependent calving rate parametrisations have been used to reproduce historical retreat of individual glaciers (e.g. Columbia glacier, Alaska (Ultee and Bassis, Reference Ultee and Bassis2016)), our results indicate that this may not be a suitable calving law for simulating Amundsen Sea glaciers.
Recently Rosier and others (Reference Rosier, Gudmundsson, Jenkins and Naughten2025) performed an ice-flow model simulation using a calving rate parametrisation based on cliff height but found no evidence of rapid retreat. While this agrees with our results on fast-flowing regions, we observed significant frontal retreat at slow-flowing glacier terminals. The apparent discrepancy with some of our results is most certainly due to different model setups, in particular to differences in model domains. Rosier and others (Reference Rosier, Gudmundsson, Jenkins and Naughten2025) focused on fast-flowing Amundsen Sea glaciers and excluded northern Peninsulas and ice rises near PIG and west of Dotson, the exact areas where we observe rapid retreat using the same parametrisation. Additionally, their three-year model relaxation followed by a three-year spin-up, while standard for minimising initial shocks, likely caused ice thickness drifts that could lower cliff heights in some locations.
Here, we have limited our analysis to one particular calving rate parametrisation formulation. Various other cliff-height-dependent parametrisations have been proposed, e.g. Mercenier and others (Reference Mercenier, Lüthi and Vieli2018); Schlemm and Levermann (Reference Schlemm and Levermann2019); Crawford and others (Reference Crawford, Benn, Todd, Åström, Bassis and Zwinger2021). The recent study by Morlighem and others (Reference Morlighem2024) used the cliff-height-dependent calving rate parametrisation proposed by Crawford and others (Reference Crawford, Benn, Todd, Åström, Bassis and Zwinger2021). This calving rate parametrisation formulation predicts that calving rates will be close to zero for cliffs of height less than 136 m. Currently, surface elevations around all grounding lines of the WAIS are less than 136 m, although such surface elevations are to be found in some locations only a few km upstream from current grounding lines. Similarly to what we find here in the DP+ABUMIP0 experiment using the DeConto and Pollard (Reference DeConto and Pollard2016) calving rate parametrisation, they find that irreversible calving front retreat is also unlikely for the Crawford and others (Reference Crawford, Benn, Todd, Åström, Bassis and Zwinger2021) calving rate parametrisation. However, as mentioned above, we do find persistent frontal retreat provided all ice shelves are continuously removed, i.e. in the DP+ABUMIP experiment.
Summary and outlook
Ice elevations across the grounding lines of the fast-flowing areas of all major glaciers flowing into the Amundsen Sea Embayment are generally above 80 m and in many cases above 100 m (Fig. 1). Based on the DP calving rate parametrisation, Eq. (1), removal of all ice shelves would therefore give rise to calving rates of up to 
$3\,\mathrm{km}\,\mathrm a^{-1}$, which is sufficient to cause frontal retreat based on current observed ice velocities. However, our numerical experiments where we initially remove all ice shelves and apply calving rate parametrisation Eq. (1), i.e. Experiment DP+ABUMIP0, do not result in persistent frontal retreats of those fast-flowing areas. Rather, we find that the combination of the increase in ice flow as the ice shelves are removed, and the reduction in the height of the ice cliffs due to ice deformation, results within a few years in a calving front advance. Importantly, we do not find that removing all ice shelves triggers an unstoppable retreat of calving fronts at Pine Island and Thwaites glaciers.
Continuously removing all ice shelves as they form in combination with the use of the DP calving rate parametrisation, Eq. (1), i.e. Experiment DP+ABUMIP, does cause persistent frontal retreat. When DP calving rate parametrisation is applied, the retreat initially occurs across slow-flowing areas where currently cliff heights are above 80 m. However, remote sensing observations show that the location of the calving fronts in those areas has not changed significantly in recent decades (MacGregor and others, Reference MacGregor, Catania, Markowski and Andrews2012; Greene and others, Reference Greene, Gardner, Schlegel and Fraser2022). Thus, these results lead us to conclude that the cliff-height-dependent calving rate parametrisation proposed in DeConto and Pollard (Reference DeConto and Pollard2016) results in unrealistic retreat rate for the current geometrical configuration of glaciers in the area. Direct application of DP calving rate parametrisation, Eq. (1), in combination with automatic removal of all ice shelves leads to a deglaciation of the whole area within just 100 years. In this regard, our results are similar to those of DeConto and Pollard (Reference DeConto and Pollard2016). However, this loss is due to modelled retreat in slow-flowing areas. For fast-flowing areas, we do not model sustained or accelerating frontal retreat. Any initial frontal retreat in those areas is, in all our simulations, quickly reversed to a frontal advance, as ice deformation causes reduction in subaerial cliff heights.
Calving has the potential to greatly increase mass loss from the Antarctic Ice Sheet, yet calving rate parametrisations continue to be poorly constrained and there is no general consensus on which, if any, calving rate parametrisations are most suitable for the Antarctic Ice Sheet. We argue that the most suitable approach to validate calving rate parametrisations is modelling recent frontal variations of glaciers, comparing those with observations, and discarding those calving rate parametrisations that do not provide a sufficiently good match (Amaral and others, Reference Amaral, Bartholomaus and Enderlin2020). This approach is becoming an increasingly viable option for future study, thanks to the availability of new observational datasets. For example, advances in remote sensing have already resulted in the production of highly accurate data sets on the heights of frontal cliffs in Antarctica (Howat and others, Reference Howat2022). Data sets with high spatial and temporal resolution on surface velocities and frontal variations have been available for some time (ENVEO and others, Reference ENVEO, Hetzenecker, Nagler and Scheiblauer2021; Greene and others, Reference Greene, Gardner, Schlegel and Fraser2024). From these data compilations, constraints on calving rate parametrisations can be extracted (Parsons and others, Reference Parsons, Sun and Gudmundsson2025) and modelled temporal variations in frontal positions validated against observations.
Data Availability Statement
The ice sheet model Úa and its manual are publicly available at: https://github.com/GHilmarG/UaSource. Configuration files for the Úa simulations are openly available at https://github.com/SainanSun0/ASE-MICI-revisit. All datasets used in this study are freely accessible through their original sources as cited in the References section.
Acknowledgements
S. Sun received support through the project Interacting ice Sheet and Ocean Tipping - Indicators, Processes, Impacts and Challenges (ISOTIPIC), funded by the Natural Environmental Research Council, Grant NE/Z503344/1. G. H. Gudmundsson was supported by the project PRECISE - PREdiction of Ice Sheets on Earth, provided by the Novo Nordisk Fonden, Grant Nr. NNF23OC0081251.
Appendix A
Basal sliding coefficient and rate factor estimated by the inverse method, where the red contour represents the grounding line and the blue contour represents the calving front at the initial state.
Modelled calving front locations for different mesh resolutions at 50 years after the start of the Experiment DP. Mesh resolutions range from 32 km down to 250 m. As in other numerical results presented, the computational domain was discretised using linear 3-node triangular elements. The size of the elements was defined as the leg of an isosceles right triangle with the same area as the triangular element. Alternatively, the element size can be defined as the length of the side of a perfect square with the same area as the triangular element, in which case the numbers provided need to be divided by 
$\sqrt{2}$.
Ice thickness changes after 10 years for all experiments. The solid magenta lines mark the locations of the grounding line, and the dashed magenta lines represent the locations of the calving fronts. The initial grounding line and calving front locations are shown as black solid and black dashed lines, respectively.
Changes of ice thickness in 30 years in all experiments. The solid lines mark the locations of the grounding line, and the dashed lines represent the locations of the calving front. The initial locations are in black, and the year 30 locations are in magenta.
Ice thickness changes over 100 years across all experiments. Solid lines indicate grounding line positions; dashed lines indicate calving front positions. Black lines show initial positions (year 0); magenta lines show final positions (year 100). White regions represent ice-free areas resulting from frontal retreat.
Open access
