Sensitivity of the future evolution of the Wilkes Subglacial Basin ice sheet to grounding-line melt parameterizations

Abstract

Projections of Antarctic Ice Sheet mass loss and therefore global sea level rise are hugely uncertain, partly due to how mass loss of the ice sheet occurs at the grounding line. The Wilkes Subglacial Basin (WSB), a vast region of the East Antarctic Ice Sheet, is thought to be particularly vulnerable to deglaciation under future climate warming scenarios. However, future projections of ice loss, driven by grounding-line migration, are known to be sensitive to the parameterization of ocean-induced basal melt of the floating ice shelves and, specifically, to the adjacent grounding line – termed grounding-line melt parameterizations (GLMPs). This study investigates future ice sheet dynamics in the WSB with respect to four GLMPs under both the upper and lower bounds of climate warming scenarios from the present to 2500, with different model resolutions, ice shelf melt parameterizations (ISMPs) and choices of sliding relationships. The variation in these GLMPs determines the distribution and the amount of melt applied in the finite-element assembly procedure on partially grounded elements (i.e. elements containing the grounding line). Our findings indicate that the GLMPs significantly affect both the trigger timings of tipping points and the overall magnitude of ice mass loss. We conclude that applying full melting to the partially grounded elements, which causes melting on the grounded side of the grounding line, should be avoided under all circumstances due to its poor numerical convergence and substantial overestimation of ice mass loss. We recommend preferring options that depend on the specific model context, by either (1) not applying any melt immediately adjacent to the grounding line or (2) employing a sub-element parameterization.

1 Introduction

Melting beneath ice shelves and iceberg calving are recognized as equally important contributors to the current mass loss of the Antarctic Ice Sheet , accounting for a total contribution of approximately 5.2 $\mathrm{mm}$ to global sea level rise since 2003 . Basal melting plays a crucial role in the contemporary amplification of ice discharge in Antarctica . Variations in basal melt rates exert significant influence on ice shelf thickness, with thinning leading to a diminished ice shelf buttressing effect. The reduction in buttressing subsequently results in the acceleration of ice streams that supply the ice shelf. Such acceleration contributes to dynamic thinning of the ice upstream of the grounding line, inducing grounding-line retreat. The associated loss of basal resistance may, in turn, provoke a positive feedback if the subglacial topography deepens towards the interior of the continent. This unstable behaviour is known as the marine ice sheet instability (MISI) .

Figure 1

Bed topography of the Wilkes Subglacial Basin (WSB) and the designated catchment used as the model domain. The three primary outlet glaciers of the WSB (Cook, Ninnis and David glaciers) are marked. The orange contour delineates the model domain in this study.

[Figure omitted. See PDF]

The Wilkes Subglacial Basin (WSB; Fig. ), located west of the Transantarctic Mountains in East Antarctica, spans approximately 400 000 ${\mathrm{km}}^{\mathrm{2}}$, with depths extending as far as 2000 $\mathrm{m}$ below sea level in a deep marine-based setting. Ice flow predominantly occurs along two deep troughs extending subglacially towards the Cook and Ninnis ice shelves, which currently discharge 40.6 and 23.0 $\mathrm{Gt}{\mathrm{a}}^{-\mathrm{1}}$ of ice into the ocean, respectively Fig. . The WSB is notable for its extensive ice reserves and vulnerability to marine ice sheet instability . Tipping-point behaviour (onset and continuation of MISI) has been shown to occur in simulations , yet there is a paucity of observations and modelling efforts to inform this potentially unstable behaviour . As such, the WSB may be particularly sensitive to melting beneath the ice shelf and the grounding-line dynamics, thereby rendering the disparities among our sensitivity experiments more pronounced. These factors motivated us to select the WSB as the focus of our study.

Recent studies indicate that the migration of the grounding line is extremely sensitive to how basal melt occurs adjacent to the grounding line . However, due to constrained observations, our understanding of the actual melt rates at the grounding line and their underlying mechanisms remains in its infancy . Traditional plume and ocean models generally predict that the basal melt rates tend to approach zero near the grounding line e.g., with the peak melt occurring about 10 to 15 $\mathrm{km}$ away from it . In a detailed study, applied the ocean model NEMO to simulate Antarctic ice shelf melt rates, finding more than half of the ice shelves show melt rates approximating zero at the grounding line, with an average rate of 0.45 $\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}$ across all of them. Nevertheless, other studies challenge this traditional understanding represented by the plume model. discussed the possibility of high melting at and even glaciologically upstream of the grounding line caused by the intrusion of layered warm salty water. In their theoretical model experiments, seawater intrudes as far as several kilometres upstream of the grounding line, potentially doubling ice mass loss . validated the seawater intrusion theory by analysing satellite radar interferometry, revealing melt rates of up to 80 $\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}$ in the tidally influenced grounding zone of Petermann Glacier in Greenland. From another perspective, the Antarctic basal melt rates computed by , based on remote sensing observations and ice flux divergence calculation, do not show a pattern of melt rates approaching zero at the grounding line. In this study, both zero melt and high melt near the grounding line are examined through different ice shelf melt parameterizations (ISMPs).

Modelling studies suggest that ice sheet models may be more sensitive to melt rates near the grounding line than to cavity-integrated melt rates beneath ice shelves e.g.. As such, accurately simulating melt patterns, particularly near the grounding line, might be at least as important as simulating realistic integrated melt . Accurate representation of basal melt at the grounding line is crucial for ice flow models to reduce uncertainties in forecasting ice sheet dynamics and future mass loss . However, due to the discretization of the general fixed-grid ice sheet model, there inevitably exist grid cells or elements at the grounding line where ice is partially grounded and partially floating. How to represent basal melting within these cells remains a challenging and unresolved issue, which is further explored here.

In the past decade, various parameterization schemes for handling sub-grid scale features at the grounding line in basal friction and melt have been explored e.g.. The initial motivation to explore grounding-line parameterization was to optimize the treatment of basal friction at the grounding line, given its high impact on grounding-line dynamics . Sub-element parameterizations for the representation of basal friction generally over partially grounded elements provide improved convergence of model behaviour with a finer mesh resolution , and they are widely used in subsequent research on ice sheet modelling e.g.. pioneered a comprehensive study on representation of basal melt under partially floating cells, based on the MISMIP model configuration . They recommend for models to avoid the application of melt rates over entire partially floating cells, as this gives worse convergence with resolution and overestimates grounding-line retreat at typically used resolutions . Following this, a majority of subsequent ice sheet modelling efforts adopt melt parameterizations assuming zero melt at the grounding line . In ice sheet model intercomparisons, such as initMIP-Antarctica , it was found that marine ice sheet models using sub-element melt (SEM) parameterizations are consistently more sensitive to ocean forcing than those without melt applied to these elements (increasing the Antarctic contribution to sea level rise by 50 %–100 %; ). However, recent studies suggest that, in their finite-difference-based model experiments, models applying melt at the grounding line on the partially floating cells overall outperform those not applying melt in terms of convergence with resolution.

This study seeks to delve deeper into various parameterization solutions for basal melt at the grounding line applied to the domain of the Wilkes Subglacial Basin through a series of sensitivity experiments. We detail the methods in Sect. , including model configurations and inversions for ice viscosity and basal friction, as well as the experimental design of transient simulations. The results of a series of sensitivity experiments are presented in Sect. , with a subsequent discussion in Sect. . Conclusions are provided in Sect. .

Figure 2

Overview of the experimental workflow in this study. Marked in brackets is the resolution of the model grid. The results obtained from the inversion, including basal-friction parameter $\mathit{\beta }$ and viscosity enhancement factor $E_{\mathit{\eta }}$, are interpolated onto four grids to initialize the subsequent historical runs.

[Figure omitted. See PDF]

Table 1

Summary of the simulation naming convention. SSP: Shared Socioeconomic Pathway. WCS: water column scaling. NMP: no-melt parameterization. FMP: full-melt parameterization. SEM: sub-element melt.

Name part	Meaning	Possible values
FL	Basal-friction law	Weertman or Coulomb
SSP	Emission scenario of thermal forcing	SSP126 or SSP585
RES	Characteristic mesh resolution	250 $\mathrm{m}$, 500 $\mathrm{m}$, 1 $\mathrm{km}$, 2 $\mathrm{km}$
ISMP	Ice shelf melt parameterization	NoWCS, WCS75
GLMP	Grounding-line melt parameterization	NMP, FMP, SEM1, SEM3

2 Methods

We use Elmer/Ice to conduct a series of ice sheet simulations for the WSB. Elmer/Ice is an open-source finite-element ice sheet–ice shelf model capable of solving the full-Stokes equations but also allows for various simplifications, such as the shallow-shelf approximation we use here SSA;. We conduct a series of sensitivity experiments of the WSB with SSA to investigate the sensitivity of grounding-line movement and ice mass loss to different grounding-line melt parameterizations (GLMPs). The workflow is illustrated in Fig. . The sensitivity experiments encompass a range of model choices, including two basal-friction laws; two climate forcing scenarios; four characteristic mesh resolutions; two ice shelf melt parameterizations (ISMPs); and, as the focus of the study, four GLMPs for the partially floating elements. Each simulation is designated by the naming convention FL_SSP_RES_ISMP_GLMP, with the specific components detailed in Table . The model components involved will be introduced in detail in the following subsections, which include the model setup and inversions (Sect. ) and transient simulations (Sect. ), covering both historical and future runs sequentially.

Figure 3

Refined grid for the Wilkes Subglacial Basin with a 1 $\mathrm{km}$ characteristic mesh resolution (a); grid details at the Ninnis Ice Shelf, marked by a red box in (a), with 2 $\mathrm{km}$ (b), 1 $\mathrm{km}$ (c), 500 $\mathrm{m}$ (d) and 250 $\mathrm{m}$ (e) characteristic mesh resolutions.

[Figure omitted. See PDF]

2.1 Model setup and inversions

The two-dimensional (2-D) mesh used for the WSB domain is constructed using Gmsh . It features a quasi-uniform, unstructured triangular grid at a 1 $\mathrm{km}$ resolution. The inland domain boundary defining the glacier basin of the WSB model is sourced from MEaSUREs Antarctic Boundaries, Version 2 . The coastline boundary, initial ice geometry and bed topography are taken from MEaSUREs BedMachine Antarctica, Version 3 . The locations of the calving front and inland boundary are held fixed throughout the simulations. A minimum ice thickness of 15 $\mathrm{m}$ is maintained to preserve a thin ice shelf as it retreats. We then conduct mesh refinement using Mmg to optimize computational efficiency without compromising accuracy. We estimate the location of the grounding line in the year 2300, based on the projected grounding-line movement under the most severe ice loss scenario from the Antarctic model in the ISMIP6-2300 project . For the area downstream of this line, the grid is refined to characteristic resolutions of 250 $\mathrm{m}$, 500 $\mathrm{m}$, 1 $\mathrm{km}$ and 2 $\mathrm{km}$, respectively (Fig. ), in preparation for subsequent sensitivity experiments. Conversely, for its upstream inland region, the mesh resolution is progressively transitioned to coarser scales. The four grids maintain a very similar mesh resolution in the far inland area, characterized by elements with a horizontal extent of approximately 17 $\mathrm{km}$. This refinement strategy is designed to prevent the grounding line from retreating into areas of coarser resolution during centennial-scale transient runs. Besides, the local refinement metric draws upon both ice surface velocity observations and ice thickness , allocating slightly a finer resolution to regions with pronounced gradients in velocity and thickness. The statistics of the four grids are shown in Table 2.

Table 2

Summary of the four grids.

Mesh resolution	Nodes	Triangular elements
2 $\mathrm{km}$	54 771	94 894
1 $\mathrm{km}$	172 389	316 170
500 $\mathrm{m}$	612 204	1 142 726
250 $\mathrm{m}$	2 317 821	4 270 368

In this study, we solve the 2-D vertically integrated SSA equations for the stress balance. We consider two friction laws for the basal shear stress ${\mathit{\tau }}_{\mathrm{b}}$, the linear Weertman law and regularized Coulomb law : $$\begin{array}{}\text{1}& {\displaystyle }& {\displaystyle }{\mathit{\tau }}_{\mathrm{b}}=-C_{\mathrm{W}}{\mathit{u}}_{\mathrm{b}},\text{2}& {\displaystyle }& {\displaystyle }{\mathit{\tau }}_{\mathrm{b}}=-\mathit{\lambda }{C}_{\mathrm{C}}\cdot {\left({\displaystyle \frac{\Vert {\mathit{u}}_{\mathrm{b}}\Vert }{\Vert {\mathit{u}}_{\mathrm{b}}\Vert +u_{\mathrm{0}}}}\right)}^{\frac{\mathrm{1}}{m}}{\displaystyle \frac{{\mathit{u}}_{\mathrm{b}}}{\Vert {\mathit{u}}_{\mathrm{b}}\Vert }},\end{array}$$ where $C_{\mathrm{W}}$ and $C_{\mathrm{C}}$ are friction coefficients and ${\mathit{u}}_{\mathrm{b}}$ is the basal velocity field. This form of the regularized Coulomb law, Eq. (), follows , which subsumed the potentially non-linear dependence of effective pressure $N$ into the friction coefficient $C_{\mathrm{C}}$. $\mathit{\lambda }$ is used as a scaling factor: 3$$\mathit{\lambda }=\left\{\begin{array}{ll}\mathrm{1},& \text{for }{h}_{\text{af}}\ge h_{\mathrm{T}}\\ \frac{h_{\text{af}}}{h_{\mathrm{T}}},& \text{otherwise},\end{array}\right.$$ where $h_{\text{af}}$ is the height of ice above flotation and $h_{\mathrm{T}}$ is a threshold height. demonstrate that the Coulomb friction field has relatively low sensitivity to the choice of parameter $u_{\mathrm{0}}$ and suggest that their parameter setting can be transferred for use with general glaciers well. We set $u_{\mathrm{0}}$ $=$ 300 $\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}$ and $h_{\mathrm{T}}$ $=$ 75 $\mathrm{m}$ for all experiments that use the regularized Coulomb law, following the settings by . $m$ is a positive exponent corresponding to the creep exponent in Glen's law . Here, we use $m=\mathrm{3}$ following and . We assume a non-linear isotropic rheology following Glen's flow law . For the viscosity $\mathit{\eta }$, we use 4$$\mathit{\eta }=E_{\mathit{\eta }}^{\mathrm{2}}{\mathit{\eta }}_{\mathrm{0}},$$ where ${\mathit{\eta }}_{\mathrm{0}}$ represents the reference field for $\mathit{\eta }$. It is calculated from a 2-D temperature field, which is obtained by vertically averaging a three-dimensional (3-D) field. The 3-D field is derived from a multi-millennial spin-up of all of Antarctica, utilizing the ice sheet model SICOPOLIS . Furthermore, the values for activation energies and prefactors, essential for computing the temperature-dependent rate factor $A$ in accordance with Glen's flow law, are adopted from . The term $E_{\mathit{\eta }}$ in the equation stands for the viscosity enhancement factor, the determination of which will be achieved through inversion processes.

In this study, we invert the basal shear stress and ice viscosity using the refined 1 $\mathrm{km}$ resolution mesh (Fig. a), with ice velocity observations as the optimization target. We employ the linear Weertman law to compute the basal shear stress in the inversion process. More specifically, we utilize the adjoint inverse method with Tikhonov regularization, as described in and to invert the friction parameter $\mathit{\beta }$ and viscosity enhancement factor $E_{\mathit{\eta }}$ simultaneously. $\mathit{\beta }$ is given by $C_{\mathrm{W}}={\mathrm{10}}^{\mathit{\beta }}$. The inversion criterion is twofold: to minimize the velocity misfit and to avoid over-fitting of the inversion solution to non-physical noise in the velocity observation. We introduce three regularization terms in the total cost function: 5$$J_{\text{tot}}=J_{\mathrm{0}}+{\mathit{\lambda }}_{\mathit{\beta }}{J}_{\text{reg}\mathit{\beta }}+{\mathit{\lambda }}_{E_{\mathit{\eta }\mathrm{1}}}{J}_{\text{reg}{E}_{\mathit{\eta }\mathrm{1}}}+{\mathit{\lambda }}_{E_{\mathit{\eta }\mathrm{2}}}{J}_{\text{reg}{E}_{\mathit{\eta }\mathrm{2}}}.$$

The misfit between the magnitudes of simulated ($\mathit{u}$) and observed (${\mathit{u}}_{\text{obs}}$) surface velocity is encapsulated in the first cost term $J_{\mathrm{0}}$, which is a discrete sum evaluated directly at every grid node: 6$$J_{\mathrm{0}}=\sum _{\mathrm{1}}^{N_{\text{obs}}}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left|\right|\mathit{u}-{\mathit{u}}_{\text{obs}}|{|}^{\mathrm{2}},$$ where $N_{\text{obs}}$ is the total number of grid nodes with observations. The terms $J_{\text{reg}\mathit{\beta }}$ and $J_{{\text{reg}}_{E_{\mathit{\eta }\mathrm{1}}}}$ are implemented to penalize the first spatial derivatives of $\mathit{\beta }$ and $E_{\mathit{\eta }}$, respectively. Meanwhile, $J_{{\text{reg}}_{E_{\mathit{\eta }\mathrm{2}}}}$ penalizes the deviations from the prior (i.e. Glen's flow law; $E_{\mathit{\eta }}=\mathrm{1}$). The coefficients ${\mathit{\lambda }}_{\mathit{\beta }}$, ${\mathit{\lambda }}_{E_{\mathit{\eta }\mathrm{1}}}$ and ${\mathit{\lambda }}_{E_{\mathit{\eta }\mathrm{2}}}$ are positive regularization weighting parameters. We determine the optimal combination of these three parameters by conducting an “L-surface” analysis, resulting in ${\mathit{\lambda }}_{\mathit{\beta }}$ $=$ 20 000, ${\mathit{\lambda }}_{E_{\mathit{\eta }\mathrm{1}}}$ $=$ 10 000 and ${\mathit{\lambda }}_{E_{\mathit{\eta }\mathrm{2}}}$ $=$ 0.02. This L-surface analysis represents an innovative aspect of this study and is elaborated upon in Appendix .

Figure 4

The optimized basal resistance parameter $\mathit{\beta }$ (a), viscosity enhancement factor $E_{\mathit{\eta }}$ (b) and relative surface horizontal velocity discrepancy (c) for the WSB. The relative surface velocity discrepancy is the magnitude of the surface horizontal velocity difference between observations and simulations as a fraction of the observations. The three contours (yellow, orange and red) represent the observed surface velocities of 200, 700 and 1000 $\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}$, respectively. The white contour in (c) represents the observed surface velocity of 5 $\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}$. The black line represents the grounding line from BedMachine Antarctica V3 .

[Figure omitted. See PDF]

The spatial distributions of the two parameters are shown in Fig. a and b, respectively. As shown in Fig. c, the velocity difference between the inversion result and observations was assessed in terms of relative difference. The results indicated that the simulated velocities from the inversions were in good agreement with the observed velocities, especially in the fast-flow areas where velocities exceed 200 $\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}$ (Fig. c). In these fast-flow regions, relative differences are predominantly below 5 %. In Fig. c, the blue area indicates a high relative-velocity discrepancy and corresponds to regions with very slow flow (mostly below 5 $\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}$). Therefore, it does not present a concern. Such findings underscore that the inversion results can effectively serve as a reliable starting point for subsequent transient experiments. We interpolate the simulated basal-friction coefficient $\mathit{\beta }$ and viscosity enhancement factor $E_{\mathit{\eta }}$ from the 1 $\mathrm{km}$ resolution grid onto the 250 $\mathrm{m}$, 500 $\mathrm{m}$ and 2 $\mathrm{km}$ resolution grids. These interpolations serve as the starting points for the subsequent historical runs on the four distinct grids (Fig. ).

2.2 Transient simulations

We explore the sensitivity of ice dynamics to the four different GLMPs by conducting a series of transient simulations. After the inversions, we initiate historical runs to smoothly transition the model past an initial adjustment phase in the forward transient simulations (Fig. ). The historical runs span 20 years, from 1995 to 2015. Then we conduct future runs from 2015 to 2500 (Fig. ). Each future run is directly paired with a corresponding historical run, maintaining a consistent model configuration throughout.

Figure 5

Grounding-line discretization. Grounding line's exact location (a), no-melt parameterization (NMP; b), full-melt parameterization (FMP; c), sub-element melt parameterization 1 (SEM1; d) and sub-element melt parameterization 3 (SEM3; e). This figure is adapted from .

[Figure omitted. See PDF]

As the primary focus of this study, we test four GLMPs for partially floating elements, as shown in Fig. . We essentially adopt the parameterization schemes outlined by in an idealized domain. The “full-melt parameterization” (FMP) applies melt across all partially floating elements, irrespective of the grounding line's exact position. Conversely, for the “no-melt parameterization” (NMP), there is no melt applied to any part of these elements. The remaining two schemes employ sub-element parameterizations. In “sub-element melt 1” (SEM1), melt is applied to the entire area of partially floating elements, but its magnitude is reduced based on the fraction of the area of floating ice in the element. This ensures that the total melt over the element is proportionate to the floating ice area. In “sub-element melt 3” (SEM3), an increased number of 20 integration points are used during the finite-element assembly procedure within any partially floating element. We determine the float/ground status for each point and calculate the basal melt rate for the floating points based on its specific coordinates. It is named SEM3 to differentiate it from SEM2 in . In essence, our SEM3 aligns with the principles of sub-element parameterization 3 (SEP3) from , which indicate that with a sufficient number of integration points, the functionality of SEP3 closely mirrors that of sub-element parameterization 2 (SEP2). Thus, we anticipated that SEM3 in this study will perform similarly to SEM2 as described by . For basal friction on the partially floating elements, we consistently adopt SEP3 with 20 integration points for all transient experiments, following the methods discussed by .

We impose surface mass balance (SMB) and basal mass balance (BMB) data sourced from the ISMIP6-2300 project , based on CMIP6 (Coupled Model Intercomparison Project) climate model data, as the forcing. More specifically, the SMB consists of an average value for the reference period ${\mathrm{SMB}}_{\mathrm{ref}}$ and yearly SMB anomalies $\mathrm{aSMB}$: 7$$\mathrm{SMB}(x,y,t)={\mathrm{SMB}}_{\mathrm{ref}}(x,y)+\mathrm{aSMB}(x,y,t).$$

Figure 6

The basal mass balance distributions under the Cook Ice Shelf at a 1 $\mathrm{km}$ characteristic mesh resolution. Panel (a) shows the standard quadratic local parameterization from ISMIP6-2300. Panel (b) depicts the modified version using water column scaling with a threshold thickness of 75 $\mathrm{m}$.

[Figure omitted. See PDF]

In this equation, ${\mathrm{SMB}}_{\mathrm{ref}}$ represents the temporal average spanning 1995 to 2300 and is derived from MAR simulation products . $\mathrm{aSMB}$ is calculated based on thermal forcing from climate models, detailed below. Following the ISMIP6-2300 standard melting parameterization , the BMB is calculated using a quadratic function of thermal forcing as described by , complemented by a thermal forcing correction suggested by . Building upon this, we produce a revised version, whereby the basal melt rate smoothly transitions to zero as it approaches the grounding line: 8$$m_s(x,y)=m(x,y)\tanh \left({\displaystyle \frac{H_{\mathrm{c}}}{H_{\mathrm{c}\mathrm{0}}}}\right),$$ where $H_{\mathrm{c}}$ is the water column thickness beneath the ice shelf and $H_{\mathrm{c}\mathrm{0}}$ is a threshold thickness. An empirical value of $H_{\mathrm{c}\mathrm{0}}$ $=$ 75 $\mathrm{m}$ is adopted here, with the justification for this choice detailed in . This water-column-thickness-based scaling is inspired by prior research e.g. and serves as a comparison to Experiment 1 in . In the naming convention (Table ), this modified ISMP is designated WCS75 (water column scaling with a threshold thickness of 75 $\mathrm{m}$), while the original version is labelled NoWCS (no water column scaling). The comparison of the two ISMPs used is shown in Fig. . It is important to distinguish between the roles of ISMPs and GLMPs. ISMPs essentially represent two distinct physical assumptions regarding the melt rate around the grounding line, whereas GLMPs represent different parameterized implementations of the model.

We utilize the thermal forcing provided by the ISMIP6-2300 project to determine the BMB and $\mathrm{aSMB}$ applied during the transient simulations. Two emission scenarios are adopted in the two CMIP6 models for generating the thermal forcing: one sourced from the CESM2 climate model under SSP5-8.5 and the other from the UKESM1 model under SSP1-2.6. The original forcing data from the ISMIP6-2300 project span the period from 1995 to 2300. Beyond 2300, we extrapolate the forcing to the year 2500 by randomly sampling values from the 2280 to 2300.

Two basal sliding laws are employed in the sensitivity experiments, the linear Weertman law and the regularized Coulomb law, Eqs. () and (). The basal-friction parameter $C_{\mathrm{W}}$ for the linear Weertman law is derived directly from inversions. To derive the basal-friction parameter $C_{\mathrm{C}}$ for the regularized Coulomb law, we transform the inverted basal-friction parameter $\mathit{\beta }$ into $C_{\mathrm{C}}$ by substituting Eq. () into Eq. (): 9$$C_{\mathrm{C}}={\displaystyle \frac{{\mathrm{10}}^{\mathit{\beta }}}{\mathit{\lambda }}}{\left({\displaystyle \frac{\Vert {\mathit{u}}_{\mathrm{b}}\Vert +u_{\mathrm{0}}}{\Vert {\mathit{u}}_{\mathrm{b}}\Vert }}\right)}^{\frac{\mathrm{1}}{m}}\Vert {\mathit{u}}_{\mathrm{b}}\Vert .$$

This ensures that the basal shear stress remains consistent throughout the conversion process.

Figure 7

Total grounding-line flux simulated from 2015 to 2500 under the high-emission scenario (SSP5-8.5) using the linear Weertman sliding law. The figures represent NMP (a, e), FMP (b, f), SEM1 (c, g) and SEM3 (d, h), along with two ISMPs, NoWCS (a–d) and WCS75 (f–h). Each plot represents ice flux for the four mesh resolutions: 2 $\mathrm{km}$ (blue), 1 $\mathrm{km}$ (red), 500 $\mathrm{m}$ (yellow) and 250 $\mathrm{m}$ (purple).

[Figure omitted. See PDF]

Figure 8

Total ice mass simulated from 2015 to 2500 under the high-emission scenario (SSP5-8.5) using the linear Weertman sliding law. The figures are represent NMP (a, e), FMP (b, f), SEM1 (c, g) and SEM3 (d, h), along with two ISMPs, NoWCS (a–d) and WCS75 (f–h). Each plot represents the ice mass change for the four mesh resolutions: 2 $\mathrm{km}$ (blue), 1 $\mathrm{km}$ (red), 500 $\mathrm{m}$ (yellow) and 250 $\mathrm{m}$ (purple).

[Figure omitted. See PDF]

Figure 9

Total grounding-line flux simulated from 2015 to 2500 under the high-emission scenario (SSP5-8.5) using the regularized Coulomb sliding law. The figures represent NMP (a, e), FMP (b, f), SEM1 (c, g) and SEM3 (d, h), along with two ISMPs, NoWCS (a–d) and WCS75 (f–h). Each plot represents ice flux for the four mesh resolutions: 2 $\mathrm{km}$ (blue), 1 $\mathrm{km}$ (red), 500 $\mathrm{m}$ (yellow) and 250 $\mathrm{m}$ (purple).

[Figure omitted. See PDF]

Figure 10

Total ice mass simulated from 2015 to 2500 under the high-emission scenario (SSP5-8.5) using the regularized Coulomb sliding law. The figures represent NMP (a, e), FMP (b, f), SEM1 (c, g) and SEM3 (d, h), along with two ISMPs, NoWCS (a–d) and WCS75 (f–h). Each plot represents the ice mass change for the four mesh resolutions: 2 $\mathrm{km}$ (blue), 1 $\mathrm{km}$ (red), 500 $\mathrm{m}$ (yellow) and 250 $\mathrm{m}$ (purple).

[Figure omitted. See PDF]

Figure 11

The evolution of ice thickness along a characteristic flowline on Cook Glacier, as projected in the future run Coulomb_SSP585_500m_WCS75_SEM3 for illustration. The rainbow-coloured outlines represent the time series progression of ice thickness in the future run. The inset shows the location of the flowline in red. For a better visual presentation, ice at the front with a thickness of less than 20 $\mathrm{m}$ is not shown.

[Figure omitted. See PDF]

3 Results

This section presents the results of the future simulations from 2015 to 2500, featuring a comprehensive comparative analysis based upon the time series of two quantitative metrics: total ice mass and total grounding-line flux of the model. The analysis focuses on the high-emission scenario experiments because we can evaluate the effect of GLMPs best when the grounding line migrates. We also include results from simulations under a low-emission scenario in order for comparison. Figures  and  represent the evolution of total ice mass and the total grounding-line flux, respectively, under a high-emission scenario (SSP5-8.5) with the application of the linear Weertman sliding law. Figures  and  showcase these variables under the same emission scenario but using the regularized Coulomb sliding law. Figure  illustrates the evolution of ice thickness and grounding-line retreat in the future run. Although we have not demonstrated grounding-line hysteresis or irreversibility as discussed by , our projections of rapid grounding-line retreat across the retrograde section of the bedrock, compared to the retreat rates across the upsloping bed, strongly indicate that MISI can occur in this region.

In the linear Weertman experiments, a majority of the model configurations exhibit a relatively stable grounding-line flux over the initial 200-year span (Fig. ). During this period, the grounding line undergoes a retreat across the comparatively shallow and flat bed topography, as shown in Fig. , with persistent ice shelf thinning mainly caused by the basal melt. This phase is characterized by a stable total ice mass, as shown in Fig. . The onset of a surge in grounding-line flux (Fig. ) signals the tipping point of MISI, marked by an accelerated retreat of the grounding line into retrograde deep troughs (Fig. ; after the year 2200), subsequently manifesting itself as rapid ice mass loss in Fig. . The peak of grounding-line flux corresponds to a major rapid retreat of the grounding line within the troughs upstream of Cook Glacier (Fig. ). The tipping point of MISI, indicative of a critical transition in ice sheet dynamics, is generally attained around the year 2300 in experiments with water column scaling (Fig. e and f), while for the experiments without water column scaling (Fig. a–d), the timing of tipping point is significantly advanced. NoWCS_NMP reaches the tipping point around 2250 (Fig. a), NoWCS_FMP reaches it around 2150 (Fig. b), and both NoWCS_SEM1 and NoWCS_SEM3 attain it around 2200 (Fig. c and d), yielding very similar predictions. Notably, Weertman_SSP585_2km_NoWCS_FMP predicts the highest ice mass loss, at 1.04 $\times$ 10⁵ $\mathrm{Gt}$, doubling that of Weertman_SSP585_250m_NoWCS_FMP. This highlights the substantial dependency of the FMP scheme on grid resolution.

Figure 12

The simulated grounding lines in the year 2500 with NMP (red), FMP (purple), SEM1 (yellow) and SEM3 (blue) are presented in the bed topography map. They all are under the high-emission scenario (SSP5-8.5), without water column scaling, using the regularized Coulomb law at 1 $\mathrm{km}$ grid resolution. The grounding lines of SEM1 and SEM3 largely overlap. The grounding lines of all four GLMPs overlap around ice rises, covered by blue grounding lines. The position of the Cook and Ninnis glaciers is marked.

[Figure omitted. See PDF]

In the regularized Coulomb experiments, the system is relatively stable for the initial 100 years with NoWCS (Fig. a–d) and for around 150 years with WCS75 (Fig. e–h), after which MISI is triggered. A distinguishing feature of the Coulomb experiments is the earlier triggering of the tipping point, compared to the Weertman experiments, and the manifestation of two distinct peaks in grounding-line flux. The two peaks are dominated by the two major instances of rapid retreat of the grounding line in troughs upstream of the Cook (Fig. ) and Ninnis glaciers, respectively. The two peaks are experienced in all experiments without water column scaling (Fig. a–d), while the experiments with water column scaling only experienced the first peak in the last 100 years (Fig. e and f), due to its slower response. The overall magnitude of grounding-line flux and ice mass loss of regularized Coulomb experiments are 3 times greater than those of the linear Weertman experiments. Figure  provides a visual representation of the grounding-line position in the year 2500, comparing the four GLMPs under a specific model configuration, Coulomb_SSP585_1km_WCS75. The distance between the NMP and FMP grounding lines ranges from 20 to 70 $\mathrm{km}$, while the grounding-line locations are consistent between SEM1 and SEM3.

Tables B1 and B2 provide detailed data on the total ice mass change from 2015 to 2500 under the linear Weertman and regularized Coulomb laws, respectively. Among the four GLMPs, NMP consistently yields the lowest predictions of ice mass loss and FMP predicts the highest, with SEM1 and SEM3 being in between. Notably, the Weertman and Coulomb experiments reveal different yet internally consistent patterns of total grounding-line flux. The resolution dependence on the different parameterizations for partially grounded elements is comparable for both linear and the regularized Coulomb sliding laws, with the exception that a coarse resolution underestimates mass loss only in the case of WCS75 NMP Coulomb sliding. The choice of GLMPs exerts a significant impact on both the timing of the tipping point triggered and the cumulative magnitude of ice mass loss at a coarse resolution, while the incorporation of water column scaling can significantly diminish the discrepancies caused by different GLMPs and mesh resolutions.

Figure 13

The evolution of total ice mass (a, c) and total grounding-line flux (b, d) from 2015 to 2500 with four GLMPs. The results represent the experiments using the regularized Coulomb sliding law under the low-emission scenario (SSP1-2.6) at 1 $\mathrm{km}$ mesh resolution without (a, b) and with (c, d) water column scaling.

[Figure omitted. See PDF]

Regarding the low-emission experiments, we have opted to only present the results at 1 $\mathrm{km}$ grid resolution and using only the regularized Coulomb sliding law (Fig. ), as it did not exhibit notable distinctiveness as compared with the results of high-emission experiments. Without water column scaling, the system exhibits a continuous, albeit slight, loss of ice during the entire future simulation (Fig. a), and there is a substantial discrepancy in the total ice mass change (Fig. a) and total grounding-line flux (Fig. b) across different GLMPs. With water column scaling, the system experiences a discernible ice mass loss in the first 50 years; however, it subsequently stabilizes (Fig. c). The discrepancy is substantially reduced when water column scaling is applied (Fig. c and d), indicating a mitigation of the impact of melt scheme selections. In general, under a low-emission scenario, the predicted ice mass loss is less sensitive to the choice of GLMPs and mesh resolution in comparison to high-emission scenarios.

Figure 14

Convergence of total ice mass loss from 2015 to 2500 as a function of mesh resolution with the four GLMPs, NMP (blue), FMP (green), SEM1 (orange) and SEM2 (red). The results represent the experiments under the high-emission scenario (SSP5-8.5) for the Weertman (a, b) and Coulomb (c, d) sliding law with (b, d) and without (a, c) water column scaling. The coordinate axes are displayed on a dual logarithmic scale.

[Figure omitted. See PDF]

4 Discussion

In Fig. , we show the convergence of simulated ice mass loss with mesh resolution for different ISMPs and sliding laws. Specific data are presented in Tables  and  in Appendix . Our model, which simulates the real-world domain of the WSB, demonstrates a consistent convergence pattern with the idealized glacier model study by , showcasing a commendable level of agreement between the two ice sheet models, Elmer/Ice and ISSM (Ice-sheet and Sea-level System Model).

Among the four GLMPs, NMP tends to converge more rapidly with resolution in most cases, which is consistent with the findings of and . Our model results reveal a trend across all scenarios where ice mass loss diminishes as mesh resolution increases, except for the NMP scheme with the Coulomb law and water column scaling (Fig. d; Coulomb_SSP585_WCS75_NMP). In this scenario, the simulated ice mass loss actually increases with finer mesh resolutions. This result aligns with the simulation results from previous studies . A plausible explanation lies in the methodology of NMP, which, by definition, underestimates melt in partially grounded elements. As resolution becomes finer and elements become smaller, the area of no melting decreases, resulting in an increase in the area of melting close to the grounding line. However, this does not explain why NMP still overestimates mass loss in other cases, as resolution dependence exists not only due to the choice of GLMPs but also due to the sub-element parameterization of basal drag near the grounding line. The current study does not investigate impacts of basal drag on convergence with resolution, which has been more extensively studied, but the effects are present and not easily separated from the effects of melt parameterization. The cumulative impact of parameterizations on both basal drag and grounding-line melt is likely what determines convergence. Caution must be exercised regarding the potential for NMP to systematically under-represent melt at the grounding line and thus underestimate ice mass loss at coarse grid resolutions.

Conversely, FMP, by definition, overestimates melt in partially grounded elements, and our simulations using FMP always overestimate ice mass loss. In the experiments without water column scaling, the total ice mass loss simulated at a 2 $\mathrm{km}$ resolution is approximately double that simulated at a 250 $\mathrm{m}$ resolution (Fig. a and c). We notice that the ice sheet modelling community has largely moved away from the FMP scheme. We align with this perspective and concur with prior recommendations that the FMP scheme should be avoided under all circumstances.

Whilst FMP and NMP by definition always overestimate and underestimate melt in partially grounded elements, SEM1 and SEM3 are expected to fall in between and therefore give a more accurate estimation of melt in partially grounded elements. However, this does not translate into better convergence with resolution, with most simulations from both the current study and the work of showing significant overestimation of mass loss and grounding-line retreat when using SEM1 or SEM3 at coarse resolutions. This issue likely stems from fundamentally under-resolving the problem (i.e. the model's spatial resolution is insufficient to accurately capture and represent the dynamics at the grounding line). Although SEM1 and SEM3 provide a more viable average melt rate over partially grounded elements, the fact that they can cause thinning directly at “grounded” nodes (Fig. ) leads to ice detachment that would not occur with a fully resolved model (i.e. one with infinitely small elements, which is unachievable in practice). Consequently, this results in an overestimation of mass loss and grounding-line retreat. A more thorough handling of the partially grounded elements might be to implement runtime adaptivity with a specific focus on the grounding line itself, either by splitting partially grounded elements or by implementing a moving mesh that tracks grounding-line movement , but these approaches are beyond the scope of the current study.

The results of SEM1 and SEM3 consistently fall in between FMP and NMP results (Figs. –). The two sub-element GLMPs give almost identical results without water column scaling, which is similar to findings of the basal-friction parameterizations at the grounding line . Yet, with water column scaling, SEM1 and SEM3 diverge slightly, with SEM1 showing better convergence with resolution than SEM3. The SEM1 scheme shows the best convergence in the scenario with the Coulomb law and water column scaling. This appears contrary to the recommendation by against the use of SEM due to its overestimation of retreat of the grounding line. While NMP usually shows better convergence, SEM1 appears to outperform in specific scenarios, offering superior convergence.

In the vicinity of the grounding line, high melt rates essentially worsen the convergence with resolution and exacerbate the result discrepancies observed across all four GLMPs. This phenomenon is reflected in different aspects of the experimental results. Firstly, water column scaling significantly improves the convergence and reduces the disparities among the GLMPs (Figs.  and ). This is because when water column scaling is applied, the melt rates are significantly reduced near the grounding line, thereby minimizing the divergences represented by different GLMPs. Secondly, under a high-emission scenario, the predicted ice mass loss is more sensitive to the choice of GLMPs and mesh resolution in comparison to low-emission scenarios. In other words, the difference in simulated ice mass loss caused by the various GLMPs is significantly amplified under the high-emission scenario, as has been demonstrated by using a model of the Pine Island and Thwaites glaciers.

Numerical simulation methods and grid type significantly influence the performance of GLMPs. Consistent with previous model studies e.g., our research employs the finite-element method with an unstructured triangular grid and concludes that NMP always demonstrates superior convergence with resolution compared to FMP and usually compared to SEM. Notably, a model study employing a finite-volume method and a wavelet-based adaptive grid demonstrated the significant impact of sub-grid GLMPs on ice mass loss predictions, echoing our findings. However, respective studies by and , utilizing finite differences and a regular square grid, report contrary findings. Due to the distinct mechanism of the model implementation, the GLMPs they used differ from the four explored in our study. In addition to the NMP scheme (identical to ours), they incorporate a partial-melt parameterization (PMP; comparable to our SEM1) and a flotation condition melt parameterization (FCMP). Both FCMP and PMP outperform NMP in terms of convergence with resolution . Their agreement implies that for models using finite differences and regular grids, applying melt parameterizations to partially floating cells could be more advantageous. Thus, it is crucial to re-evaluate the performance of various GLMPs within specific model contexts.

Modelling studies emphasize the necessity of including significant melting processes within the grounding zone to replicate the observed retreat patterns . Further, satellite observations indicate pronounced melt rates at the grounding lines in both West Antarctica and Greenland . Drawing on the observations, recommend that ice sheet models adopt GLMPs that include melting at and upstream of the grounding line. We acknowledge the scientific rationale behind this suggestion; however, it may not directly translate to the parameterization strategies for the partially floating elements in ice sheet models. It is crucial to distinguish between the role of ISMPs and the specific function of GLMPs. We suggest that the ISMPs should reflect our current best understanding of ice–ocean interactions near the grounding line. Meanwhile, the design of GLMPs ought to prioritize model self-consistency and minimal resolution dependency.

The melting mechanism and the precise rates at the grounding line are still not well understood . Our NoWCS and WCS75 schemes encapsulate the divergent perspectives currently debated: one posits that the maximum melt rate occurs right at the grounding line, where the ice draft is deepest e.g., whereas ocean modelling studies suggest that the melt rate peaks downstream of the grounding line and diminish to zero towards the grounding line e.g.. A recent theoretical study suggests the possibility of melting extending upstream of the grounding line if warm salty seawater intrudes under the ice sheet e.g.. However, observations at Thwaites Glacier show melting at the grounding line there is suppressed by low current speeds and strong density stratification in the ice–ocean boundary layer . Moreover, melting processes and ice–ocean interactions at the grounding line likely differ among ice shelves due to variations in bathymetry, water mass characteristics and circulation beneath ice shelf cavities. To elucidate the melting mechanisms at play, there is a critical need for more extensive observations of melt rates in the vicinity of grounding lines.

The extensive exploration of model settings in this study underscores the significant uncertainties inherent in ice sheet modelling predictions. Utilizing the Coulomb sliding law, which is broadly considered superior, our analysis under the high-emission scenario of SSP5-8.5 suggests that the tipping point (onset of MISI; marked by a rapid increase in grounding-line flux, as shown in Fig. ) is anticipated to occur in the WSB between 2200 and 2300. After the tipping point, the grounding line retreats 110 $\mathrm{km}$ across the unstable retrograde bedrock within 100 years (as illustrated in Fig. ). The grounding-line flux consequently increases by a factor of 2.5, from 200 to 500 $\mathrm{Gt}{\mathrm{a}}^{-\mathrm{1}}$ (Fig. ). In this context, our simulations project an ice mass loss within the WSB in this scenario ranging from 0.26 to 0.42 $\times$ 10⁵ $\mathrm{Gt}$ by the year 2300. This corresponds to a mass above flotation of 0.21 to 0.33 $\times$ 10⁵ $\mathrm{Gt}$, equivalent to 0.06 to 0.09 $\mathrm{m}$ of global sea level rise. By 2500, the projected ice mass loss extends from 1.05 to 1.57 $\times$ 10⁵ $\mathrm{Gt}$, corresponding to a mass above flotation of 0.84 to 1.25 $\times$ 10⁵ $\mathrm{Gt}$, equivalent to a global sea level rise of 0.23 to 0.34 $\mathrm{m}$, assuming an extension in the final 2 decades of forcing from 2300. At a mesh resolution of 1 $\mathrm{km}$, which is commonly employed in ice sheet modelling, our model shows a change from NMP to SEM would induce a 15 % to 20 % increase in projected ice mass loss. Moreover, at a 1 $\mathrm{km}$ resolution, SEM could overestimate mass loss by up to 40 % compared to our finest mesh resolution of 250 m, whereas NMP might overestimate it by up to 25 % relative to the 250 $\mathrm{m}$ mesh, with specific overestimations dependent on the model configurations (Fig. ). These results provide a foundation for further detailed quantitative predictions and the examination of ice sheet dynamics in future stages of our ongoing research.

In our comparative analysis, both SEM and NMP schemes outperform FMP. As discussed earlier, SEM and NMP exhibit distinct advantages, each conducive to certain modelling contexts. The suitability of GLMPs is contingent upon the specific model and circumstances in question. The alignment between the results from idealized-model simulations and our comprehensive real-domain model experiments support the validity of a two-phased experiment process: one could firstly evaluate the performance of various GLMPs based on a cost-effective, idealized small ice flow model (e.g. MISMIP+; ) and then inform subsequent applications regarding more complex real-world domains. For the future explorations, mesh adaptation and re-segmentation at the sub-element scale during runtime would be a promising direction for more accurately representing basal friction and melting at the grounding line.

5 Conclusions

In this study, we explored the sensitivity of the future evolution of the Wilkes Subglacial Basin (WSB) ice sheet to grounding-line melt parameterizations (GLMPs) for the partially floating elements separating the grounded ice and the floating shelf. The study is conducted through a series of model simulations for the WSB spanning from 2015 to 2500. These simulations test the performance of four GLMPs under various model configurations, incorporating two basal-friction laws, two thermal forcing scenarios, four mesh resolutions and two ice shelf melt parameterizations (ISMPs). Drawing on our best model results, the tipping point, the onset of MISI, is projected to likely occur between 2200 and 2300 in the WSB under the high-emission scenario of SSP5-8.5, while the ice sheet system is expected to remain in a quasi-steady state under the low-emission scenario of SSP1-2.6. Under SSP5-8.5, our simulations suggest that the loss of ice from the WSB could contribute between 0.06 and 0.09 $\mathrm{m}$ to global sea level rise by 2300, while following the onset of MISI, this contribution could increase to between 0.23 and 0.34 $\mathrm{m}$ by 2500.

Our findings indicate that the GLMPs significantly affect both the timing of the tipping point triggered and the overall magnitude of ice mass loss. At a resolution considered high and commonly employed in ice sheet models (i.e. 1 $\mathrm{km}$), numerical errors due to inadequate convergence can lead to an overestimation of mass loss by up to 40 % when compared to our finest mesh resolution of 250 $\mathrm{m}$. This magnitude of overestimation is comparable to the impact of variations in basal-friction parameterizations at the grounding line. In the vicinity of the grounding line, high melt rates notably impair convergence with resolution and amplify the result discrepancies among the four GLMPs. This underscores the critical importance of not only knowing what melt rates are from an observational perspective but also choosing the appropriate melt parameterization in such scenarios.

Overall, both the SEM and NMP schemes outperform FMP in terms of mesh resolution convergence, with each exhibiting varying degrees of superiority over the other. The NMP scheme, in most scenarios, yields superior convergence of results but may systematically underestimate grounding-line retreat and ice mass loss. Conversely, the SEM scheme exhibited better convergence in the scenario with the Coulomb sliding law and water column scaling. The SEM scheme technically can more accurately represent the amount of melting in partially grounded elements and may capture some grounding-zone-like transitional behaviours, but it risks overestimating ice mass loss. As in prior studies , we advise against using the FMP scheme under all circumstances, due to its poor convergence and substantial overestimation of ice mass loss.

Our research suggests that there is currently no universally optimal melt scheme that suits all circumstances; the choice between NMP and SEM should be re-evaluated in their specific model contexts. Looking ahead, mesh adaptation and re-segmentation at the sub-element scale during runtime emerge as promising avenues for more accurately representing basal friction and melting at the grounding line. Idealized models, such as MISMIP+ , provide valuable insights for selecting GLMPs in more complex real-world domains. These improvements are critical to enhancing the accuracy of future predictions of ice mass loss and global sea level rise.

Appendix A L-surface analysis

In our cost function (Eq. ), we introduce three undetermined regularization parameters. Consequently, the conventional L-curve analysis is insufficient for our purposes, leading us to propose a more comprehensive L-surface analysis.

Throughout our analytical experiments, we adopt an empirical value of 0.02 for ${\mathit{\lambda }}_{E_{\mathit{\eta }}\mathrm{2}}$, as sensitivity experiments indicate that the inversion results are relatively insensitive to variations of ${\mathit{\lambda }}_{E_{\mathit{\eta }}\mathrm{2}}$ (as corroborated through personal communication with Fabien Gillet-Chaulet in 2020). To optimize the remaining regularization parameters ${\mathit{\lambda }}_{\mathit{\beta }}$ and ${\mathit{\lambda }}_{E_{\mathit{\eta }}\mathrm{1}}$, we undertake a systematic exploration of their feasible value combinations. As an initial step in our L-surface analysis, we conduct preliminary experiments to identify appropriate alternative values for these parameters. Specifically, we select 9 test values for ${\mathit{\lambda }}_{\mathit{\beta }}$ and 10 for ${\mathit{\lambda }}_{E_{\mathit{\eta }}\mathrm{1}}$. Pairwise combinations of these test values yield 90 distinct parameter sets for subsequent inversion experiments. The results of these experiments are presented in a 3-D visualization, as depicted in Fig. .

To identify the optimal combination of ${\mathit{\lambda }}_{\mathit{\beta }}$ and ${\mathit{\lambda }}_{E_{\mathit{\eta }}\mathrm{1}}$, we employ a metric defined as the relative distance $D_{\text{rel}}$, from each point to the origin in the 3-D coordinate system:

A1$$D_{\text{rel}}=\sqrt{{\left({\displaystyle \frac{J_{\mathrm{0}}}{max\left(J_{\mathrm{0}}\right)}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{J_{\text{reg}\mathit{\beta }}}{max\left(J_{\text{reg}\mathit{\beta }}\right)}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{J_{\text{reg}{E}_{\mathit{\eta }\mathrm{1}}}}{max\left(J_{\text{reg}{E}_{\mathit{\eta }\mathrm{1}}}\right)}}\right)}^{\mathrm{2}}}.$$

The point corresponding to the smallest $D_{\text{rel}}$ value is deemed to represent the most favourable combination of ${\mathit{\lambda }}_{\mathit{\beta }}$ and ${\mathit{\lambda }}_{E_{\mathit{\eta }}\mathrm{1}}$, marked as a red star in Fig. . Through the L-surface analysis, we determine the optimal parameter set to be ${\mathit{\lambda }}_{\mathit{\beta }}$ $=$ 20 000, ${\mathit{\lambda }}_{E_{\mathit{\eta }}\mathrm{1}}$ $=$ 10 000 and ${\mathit{\lambda }}_{E_{\mathit{\eta }}\mathrm{2}}$ $=$ 0.02.

Figure A1

L-surface. Black points represent the results from 90 parameter sets. Points connected by the red line correspond to the same ${\mathit{\lambda }}_{\mathit{\beta }}$, and points connected by the blue lines correspond the same ${\mathit{\lambda }}_{E_{\mathit{\eta }}\mathrm{1}}$. The 9 alternative values for ${\mathit{\lambda }}_{\mathit{\beta }}$ are 2000, 5000, 10 000, 20 000, 50 000, 100 000, 200 000, 500 000 and 1 000 000. The 10 alternative values for ${\mathit{\lambda }}_{E_{\mathit{\eta }}\mathrm{1}}$ are 10, 100, 1000, 2000, 5000, 10 000, 20 000, 50 000, 100 000 and 1 000 000.

[Figure omitted. See PDF]

Appendix B Simulated total ice mass loss

Table B1

Total ice mass loss ($\mathrm{Gt}$) under the high-emission scenario (SSP5-8.5) from 2015 to 2500 with the linear Weertman friction law for NoWCS (a) and WCS75 (b).

(a)	Melt parameterization (NoWCS)
Resolution	NMP	FMP	SEM1	SEM3
2 $\mathrm{km}$	57 799	103 715	73 319	72 807
1 $\mathrm{km}$	50 636	78 809	60 042	59 737
500 $\mathrm{m}$	45 245	62 448	51 152	50 857
250 $\mathrm{m}$	41 479	51 510	44 329	43 204
(b)	Melt parameterization (WCS75)
Resolution	NMP	FMP	SEM1	SEM3
2 $\mathrm{km}$	30 194	37 396	34 794	36 052
1 $\mathrm{km}$	29 637	33 530	32 338	32 858
500 $\mathrm{m}$	29 497	31 396	30 949	31 130
250 $\mathrm{m}$	29 179	30 055	29 851	29 946

Table B2

Total ice mass loss ($\mathrm{Gt}$) under the high-emission scenario (SSP5-8.5) from 2015 to 2500 with the regularized Coulomb friction law for NoWCS (a) and WCS75 (b).

(a)	Melt parameterization (NoWCS)
Resolution	NMP	FMP	SEM1	SEM3
2 $\mathrm{km}$	175 276	299 921	238 787	235 084
1 $\mathrm{km}$	171 745	234 377	206 355	204 773
500 $\mathrm{m}$	159 045	197 891	179 357	178 103
250 $\mathrm{m}$	145 699	168 594	156 673	156 217
(b)	Melt parameterization (WCS75)
Resolution	NMP	FMP	SEM1	SEM3
2 $\mathrm{km}$	90 755	126 957	115 825	$-$120 388
1 $\mathrm{km}$	98 431	114 097	110 440	$-$112 088
500 $\mathrm{m}$	103 380	110 595	109 121	$-$109 743
250 $\mathrm{m}$	105 871	109 585	108 194	$-$108 442

Code and data availability

All model simulations are implemented using Elmer/Ice 9.0 (rev bf10af7; 10.5281/zenodo.7892181, ) with the code available at https://github.com/ElmerCSC/elmerfem.git (last access: 7 November 2024) . Mesh and implementation scripts for the model are available at https://github.com/yuwang115/WSB-melt.git . Detailed output data for the model are available upon request to Yu Wang.

Author contributions

YW, CZ and RG designed the experiments together. YW, CZ, RG and TZ implemented the model simulations. YW processed, analysed and visualized the simulation results. YW drafted the paper. All authors contributed to the refinement of the experiments, the interpretation of the results and the final paper.

Competing interests

At least one of the (co-)authors is a member of the editorial board of The Cryosphere. The peer-review process was guided by an independent editor, and the authors also have no other competing interests to declare.

Disclaimer

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

Acknowledgements

Yu Wang, Chen Zhao and Benjamin K. Galton-Fenzi received grant funding from the Australian Government as part of the Antarctic Science Collaboration Initiative programme (ASCI000002). Chen Zhao is the recipient of an Australian Research Council Discovery Early Career Researcher Award (project no. DE240100267) funded by the Australian Government. This research project was undertaken with the assistance of resources and services from the National Computational Infrastructure (NCI), which is supported by the Australian Government. Rupert Gladstone and Thomas Zwinger were supported by Academy of Finland (grant nos. 322430 and 355572) and wish to acknowledge CSC – IT Center for Science, Finland, for computational resources. Rupert Gladstone was also supported by the Finnish Ministry of Education and Culture and CSC – IT Center for Science (decision diary no. OKM/10/524/2022).

Financial support

This research has been supported by the Australian Government as part of the Antarctic Science Collaboration Initiative programme (grant no. ASCI000002), the Academy of Finland (grant nos. 322430 and 355572), the Finnish Ministry of Education and Culture and CSC – IT Center for Science (decision diary no. OKM/10/524/2022), and the Australian Government (Australian Research Council Discovery Early Career Researcher Award (project no. DE240100267)).

Review statement

This paper was edited by Jan De Rydt and reviewed by Tijn Berends and one anonymous referee.