the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

# Numerical stabilization methods for level-set-based ice front migration

### Mathieu Morlighem

### G. Hilmar Gudmundsson

Numerical modeling of ice sheet dynamics is a critical tool for projecting future sea level rise. Among all the processes responsible for the loss of mass of the ice sheets, enhanced ice discharge triggered by the retreat of marine-terminating glaciers is one of the key drivers. Numerical models of ice sheet flow are therefore required to include ice front migration in order to reproduce today's mass loss and to be able to predict their future. However, the discontinuous nature of calving poses a significant numerical challenge for accurately capturing the motion of the ice front. In this study, we explore different stabilization techniques combined with varying reinitialization strategies to enhance the numerical stability and accuracy of solving the level-set function, which tracks the position of the ice front. Through rigorous testing on an idealized domain with a semicircular and a straight-line ice front, including scenarios with diverse front velocities, we assess the performance of these techniques. The findings contribute to advancing our ability to model ice sheet dynamics, specifically calving processes, and provide valuable insights into the most effective strategies for simulating and tracking the motion of the ice front.

- Article
(6619 KB) - Full-text XML
- BibTeX
- EndNote

Ice sheet numerical modeling is the best tool to make future sea level rise projections (e.g., Seroussi et al., 2020; Goelzer et al., 2020; IPCC, 2021). One key process that significantly contributes to mass loss is the retreat of marine-terminating glaciers (Mouginot et al., 2019; Choi et al., 2021; Pattyn and Morlighem, 2020). For example, in Greenland, the increased ice discharge is mainly driven by the retreat of glacier fronts (King et al., 2020), which is a direct consequence of calving and undercutting at the ice front (Wood et al., 2021; Mouginot et al., 2019), possibly intensified by increased runoff and ocean temperatures (Black and Joughin, 2023). As Greenland has very few ice shelves, ice front retreat predominantly comprises small yet frequent calving events (Black and Joughin, 2023; Cheng et al., 2021). Future projections emphasize that ice front retreat will continue to be a primary driver of Greenland's mass loss by 2100 (Choi et al., 2021). Incorporating moving boundaries into numerical ice sheet models is a vital step in advancing our understanding of ice loss mechanisms and improving the accuracy of future sea level rise projections (Crawford et al., 2021; Bondzio et al., 2017; Cheng et al., 2022).

Ice sheets are commonly modeled as incompressible fluids governed by conservation laws (e.g., Greve and Blatter, 2009), with empirical calving laws to predict calving rates at the ice front (Pollard et al., 2015; Morlighem et al., 2016). These calving laws are parameterizations developed based on physical principles and observations, which offer computationally efficient and relatively straightforward expressions for calving rates (Benn and Astrom, 2018; Choi et al., 2018). In these parameterizations, the boundary of the model, which is generally the ice front, needs to be adjusted dynamically during the transient simulation. The way ice front migration is typically handled is through a level-set function, which is a signed distance function defined over the entire computational domain with the zero level-set contour representing the ice front position (e.g., Bondzio et al., 2016; Morlighem et al., 2016). The motion of the level-set function is determined by solving an advection equation, where the difference between the ice velocity and the calving (and melting) rate at the zero contour governs the evolution (Morlighem et al., 2016).

However, numerically solving the level-set function is challenging, especially when using the finite-element method (FEM), as it can lead to instabilities due to the unbounded gradient of the solution (Larson and Bengzon, 2013). To address this issue, stabilization techniques are employed to enforce the boundedness of the solution. Additionally, the transient solution of the level-set function may not always maintain its signed distance property due to inhomogeneities in the velocity field and the accumulation of numerical errors over time, particularly through the diffusion introduced by the stabilization method. Therefore, reinitialization is generally necessary during transient simulations to restore the signed distance function property. However, as highlighted in Henri et al. (2022), reinitialization may introduce an artificial subsequent displacement of the zero level-set contour. Therefore, the selection of the reinitialization interval is critical for obtaining an accurate solution of the signed distance function, which remains inherently dependent on the specific application.

In this paper, we aim to investigate and compare various stabilization techniques in combination with different choices of reinitialization intervals implemented in the Ice-sheet and Sea-level System Model v4.23 (ISSM; Larour et al., 2012; ISSM Team, 2023) and Úa 2019b (Gudmundsson et al., 2019; Gudmundsson, 2020). We present different stabilization and reinitialization procedures and apply them all in ISSM to solve the level-set equation on an idealized domain featuring a semicircular ice front shape (and a straight-line ice front shape case in Appendix A) representative of typical Greenland outlet glaciers. To evaluate the effectiveness of the stabilization techniques and reinitialization strategies, we perform several tests on three different spatially varying rates of ice front migration, encompassing both low- and high-speed scenarios. By exploring these approaches, we seek to investigate which combination leads to the best stability and accuracy of simulating the level-set function and effectively tracks the motion of the ice front in ice sheet models.

The level-set function *ϕ*(** x**,

*t*) is a scalar field defined on a two-dimensional domain Ω with zero contours implicitly representing the ice front position at every given time

*t*. Conventionally, the level-set function is set to be negative in the ice-covered region and positive in the ice-free region (Morlighem et al., 2016) in order for the gradient of the level set to be normal outward pointing to the ice front. The absolute value of the level set is the closest distance from

**to the ice front contour**

*x**ϕ*=0. Given an initial condition $\mathit{\varphi}(\mathit{x},\mathrm{0})={\mathit{\varphi}}_{\mathrm{0}}$, the evolution of the level-set function

*ϕ*(

**,**

*x**t*) is governed by the advection equation

where *v*_{f} is the front velocity of the level set, which is the difference between the ice velocity ** v** and the calving rate of

*c*, which is generally oriented perpendicular to the ice front:

where ** n** is the outward unit normal vector of the level set (Bondzio et al., 2016; Morlighem et al., 2016).

In order to solve Eq. (1) with the FEM, we introduce a Hilbert space ℋ^{1}(Ω) and define the variational form as follows. Find *ϕ*∈ℋ^{1}(Ω) such that for all the test functions *ψ*∈ℋ^{1}(Ω) the following equation is satisfied:

After replacing the space ℋ^{1}(Ω) with a continuous piecewise linear space Φ_{h}, the solution of Eq. (3) is then the numerical solution of Eq. (1). However, it is well known that Eq. (3) gives spurious oscillatory solutions without stabilization (Larson and Bengzon, 2013; dos Santos et al., 2021).

## 2.1 Stabilization

We consider four stabilization schemes in this paper. The first three methods are classical methods only to stabilize Eq. (3), namely the artificial diffusion (MacAyeal, 1989, AD) method, streamline upwind (Eriksson et al., 1996, SU) method, and streamline upwind Petrov–Galerkin (Brooks and Hughes, 1982, SUPG) method. The last one is a modification of the SUPG stabilization, where an additional forward-and-backward (Li et al., 2005, FAB) diffusion term is added to the SUPG scheme.

Among these methods, the simplest way to stabilize an advection equation is to add an additional diffusion term in the variational form Eq. (3) such that

where in two dimensions the coefficient of the artificial diffusion term is a scalar

where *h*_{x} and *h*_{y} are the characteristic mesh sizes in the *x* and *y* directions and *v*_{x} and *v*_{y} are the *x* and *y* components of the front velocity *v*_{f}.

The streamline upwind stabilization follows the same variational form as the artificial diffusion in Eq. (4) but with a modified coefficient derived from Eq. (5). Specifically, this modification ensures the addition of diffusion solely along the direction of the velocity vector *v*_{f} by using

where $h=\sqrt{{h}_{x}^{\mathrm{2}}+{h}_{y}^{\mathrm{2}}}$ and ⊗ is the Kronecker product. Due to the large dissipation introduced by these two stabilization methods, they are extremely stable but only have first-order accuracy (dos Santos et al., 2021).

A more accurate stabilization method is the streamline upwind Petrov–Galerkin (SUPG, Brooks and Hughes, 1982) method, which modifies the test function to be $\widehat{\mathit{\psi}}=\mathit{\psi}+\mathit{\mu}{\mathit{v}}_{\mathrm{f}}\cdot \mathrm{\nabla}\mathit{\psi}$ in the variational form in Eq. (3) such that

where $\mathit{\mu}=\frac{h}{\mathrm{2}\Vert {\mathit{v}}_{\mathrm{f}}\Vert}$ is a mesh-dependent coefficient (dos Santos et al., 2021).

The FAB diffusion was first introduced in Úa (Gudmundsson et al., 2019; Gudmundsson, 2020). We follow the same formulation and implement it in ISSM. The FAB term added to the variational form in Eq. (7) is derived from the potential

for which the directional derivative is

This results in the addition of a nonlinear diffusion term to the level-set equation, with a diffusion coefficient

which is bounded for $\Vert \mathrm{\nabla}\mathit{\varphi}\Vert \to \mathrm{0}$ provided *q*≥2. For even values of *p*, the diffusion term defined by Eq. (10) can be both negative and positive and is an example of a FAB diffusion. Note that the minimum of the potential 𝒫 in Eq. (8) is found for $\Vert \mathrm{\nabla}\mathit{\varphi}\Vert =\mathrm{1}$, i.e., when *ϕ* is a distance function. This approach therefore encourages the level set to remain a signed distance function and relaxes the need to reinitialize the level set.

## 2.2 Reinitialization

The formulation of the advection equation Eq. (1) describes the evolution of the level-set function; however, it does not guarantee that the level-set function is always a signed distance function due to the inhomogeneity of the front velocity. Indeed, as *v*_{f} is generally higher at the ice front than the far field, the gradient of the level-set function close to the zero contours tends to decrease during the transient simulation.

To maintain the gradient of the level-set function, a common practice is to reset the level set by calculating the signed distance every *n*_{R} time steps. This is often called “reinitialization” (Bondzio et al., 2016; Morlighem et al., 2016), and the reinitialization interval *n*_{R} is the number of time steps between two consecutive reinitializations. One method of reinitialization involves solving an eikonal equation (Sussman et al., 1994; Sethian, 1996):

generally expressed as a time-dependent problem for which we seek a steady-state solution:

However, this approach (Eq. 12) contains control parameters, and it is not clear what the optimal value of these parameters should be in practical application (Gross and Reusken, 2011). Moreover, the eikonal equation constitutes a nonlinear hyperbolic partial differential equation (PDE), posing challenges in achieving accurate discretization. An alternative is to use the fast marching method (Sethian, 1996; Toure and Soulaimani, 2016). This method offers a general framework capable of handling various scenarios.

Here, we use a straightforward geometric reinitialization algorithm, which is similar to the one described in Toure and Soulaimani (2016). At any point in time, the zero contour of the level set is represented by a set of segments if it is discretized using linear Lagrange elements. At the reinitialization step, we create a loop over all elements and generate this set of segments, with one segment per element containing a change in the sign of the level-set function. This set of segments is then shared across all model partitions through a message-passing interface in order to recompute a signed distance. Subsequently, at each vertex of the mesh, we compute the distances to these segments and keep the minimum distance as the new magnitude of the level set at that vertex, while preserving the original sign. When this reinitialization algorithm is applied, it is expected to yield exact results in terms of signed distance. Hence, we do not expect that the proposed FAB diffusion algorithms would outperform this method in terms of accuracy. However, as we show later in the numerical experiments, numerical errors highly depend on the reinitialization frequency. Here, we investigate different reinitialization intervals combined with the four stabilization methods described in Sect. 2.1.

## 2.3 Error quantification

In order to quantify the difference between two ice front positions represented by the level-set functions *ϕ*_{1} and *ϕ*_{2}, we introduce a misfit metric d(*ϕ*_{1},*ϕ*_{2}) such that

where

converts a level-set function to a sign function with −1 on the ice-covered side of the zero contour and 1 on the ice-free side of the contour. Therefore, if *ϕ*_{1} is ahead of *ϕ*_{2} in terms of the ice front positions (more advance), the misfit area in d(*ϕ*_{1},*ϕ*_{2}) will be negative.

We integrate the absolute misfit over the whole domain, Ω, and get the following metric:

which is actually the absolute misfit area between the two level-set functions.

We investigate the influence of the four stabilization methods described in Sect. 2.1 combined with different choices of reinitialization interval (Sect. 2.2). Here we consider a semicircle-shaped initial ice front, as shown in Fig. 1, where the ice-covered region is shown in light blue and the ice-free region is shown in light red. We apply spatially and temporally varying analytical velocity fields to mimic typical ice flow.

We run all the simulations on a two-dimensional square domain $\mathrm{\Omega}(x,y)=[\mathrm{0},L]\times [\mathrm{0},L]$, with *L*=20 km as the size of the domain. We create an unstructured triangular mesh on Ω with the element size of 100 m. In Fig. 1, the calving front is represented by a semicircle (red) centered at $({c}_{x},{c}_{y})=(\frac{\mathrm{5}L}{\mathrm{8}},\frac{L}{\mathrm{2}})$ with a radius of $r=\frac{L}{\mathrm{4}}$, and the side walls of the fjord are in blue and connect the semicircle to the right boundary of the domain. By construction, the width of the fjord is 10 km. The initial zero level set is the red ice front together with the blue side walls, which has a closed form as $\mathit{\left\{}\right(x,y\left)\right|(x-{c}_{x}{)}^{\mathrm{2}}+(y-{c}_{y}{)}^{\mathrm{2}}={r}^{\mathrm{2}},x\le {c}_{x}\mathit{\}}\bigcup \mathit{\{}(x,y)|x\in [{c}_{x},L],y={c}_{y}+r\mathit{\}}\bigcup \mathit{\left\{}\right(x,y\left)\right|x\in [{c}_{x},L],y={c}_{y}-r\mathit{\}}$.

We apply three distinct velocity fields to control the migration of the ice front. For simplicity, we assume that there is no ice flux across the side walls of the fjord so that the velocity field only contains a horizontal component as follows: ${\mathit{v}}_{\mathrm{f}}=({v}_{x},\mathrm{0}{)}^{T}$. The *x* component of the velocity fields are given in Table 1. They represent zeroth- (uniform), first- (triangle), and second-order (parabola) polynomial shapes of the velocities.

Temporal variations are introduced by flipping the sign of *v*(*t*) (as in Table 1) every half-year to mimic the typical annual cycle of the advance and retreat of an ice front such that

where *T*=1 year, $n=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{\dots},N$, and *v*_{0} is a velocity constant. We examine two scenarios with high- (*v*_{0}=5000 m a^{−1}) and low-velocity (*v*_{0}=1000 m a^{−1}) constants, respectively. All the simulations are run for *N*=50 periods (or years), with a constant time step at Δ*t*=0.005 years to satisfy the Courant–Friedrichs–Lewy (CFL) condition for both the high- and low-velocity scenarios. We reinitialize the zero level-set contour with the interval *n*_{R}=1, 10, 100, 200, which corresponds to a reinitialization every 2 d, two-thirds of a month, half a year, and year. We also set a control run with no reinitialization (*n*_{R}=∞) throughout the whole simulation period.

By applying the velocity for $\frac{T}{\mathrm{2}}$ in one direction and then flipping the sign of *v*_{x} for another $\frac{T}{\mathrm{2}}$, the ice front is expected to return to its initial position *ϕ*_{0} after every period *T*. Furthermore, the analytical solution at any given time *t*+*n**T* should be identical to the solution at time *t*. Therefore, we use the numerical solution at $t\in [\mathrm{0},T)$ as the exact solution and calculate the numerical error at *t*+*n**T* according to Eq. (15), with ${\mathit{\varphi}}_{\mathrm{1}}=\mathit{\varphi}(\mathit{x},t+nT)$ and ${\mathit{\varphi}}_{\mathrm{2}}=\mathit{\varphi}(\mathit{x},t)$.

The misfit between the numerical and the exact solution under a uniform velocity field at the low-velocity setting (*v*_{0}=1000 m a^{−1}) after 1.5, 2, and 50 periods (or years) is shown in Fig. 2 with *n*_{R}=1 and in Fig. 3 with *n*_{R}=100 for the four stabilization methods considered in this paper. The misfit at every time point is calculated according to Eq. (13), where the area with negative values (blue in the figures) indicates the ice front from the numerical solution is downstream (i.e., further advanced) of the exact solution. The errors of all the cases in Figs. 2 and 3 are almost evenly distributed along the ice front, and the total misfit grows as time increases. Indeed, all the errors are first-order in time, as we show the time series of the errors in Appendix B for different stabilizations, reinitializations, and velocity constants. Figures 2 and 3 also indicate that using *n*_{R}=100 gives more accurate results compared to reinitializing every time step (*n*_{R}=1).

To facilitate a better comparison of the different stabilization, reinitialization, and velocity constant choices, we show the total absolute misfit in Fig. 4, which is calculated according to Eq. (15) at the final time step *t*=50 years for the uniform velocity field. The numerical errors tend to decrease as the reinitialization interval *n*_{R} increases. Specifically, in Fig. 4a, the four largest errors occur when the level-set function is reinitialized at every time step (*n*_{R}=1), resulting in errors of 16.87 km^{2} in AD, 16.67 km^{2} in SU, 11.32 km^{2} in SUPG, and 10.28 km^{2} in SUPG+FAB. The spatial distributions of the errors are shown in Fig. 2c, f, i, and l. Given that the width of the fjord is 10 km, these errors correspond to an average offset of the ice front of 1 to 2 km along the flow direction. After *n*_{R}>10, the numerical errors remain almost constant, comparable to the ones of *n*_{R}=∞, for all the stabilization methods employed. We find a similar pattern in the high-velocity (*v*_{0}=5000 m a^{−1}) cases in Fig. 4b, where most of the numerical errors are approximately 5 times larger than those in the low-velocity (*v*_{0}=1000 m a^{−1}) cases in Fig. 4a. However, the high-velocity cases are less sensitive to *n*_{R} than the low-velocity cases. For instance, reinitializing every time step does not introduce exceptionally large errors as we found in the low-velocity cases. Indeed, the largest numerical error (88.62 km^{2}) among all the experiments is achieved by the AD stabilization without reinitialization.

Although all four stabilization methods tend to overestimate the advance of the ice front, the choice of stabilization method has a significant impact on the misfit area, and SUPG+FAB exhibits the lowest numerical errors. In the low-velocity scenario, e.g., Fig. 4a, with *n*_{R}=100, the final misfit for SUPG+FAB is 0.46 km^{2}, whereas the errors for AD, SU, and SUPG are 11.51, 4.92, and 0.44 km^{2}, respectively. The spatial distributions of these errors are shown in Fig. 3c, f, i, and l, where the misfit achieved by SUPG is equivalent to an offset of the ice front by approximately 46 m, which is even less than half of the mesh size. Similarly, in the high-velocity scenario, the errors are scaled by the front velocity in all the choices of stabilizations with *n*_{R}>1. For instance, in Fig. 4b, with *n*_{R}=100, the errors are 49.99 km^{2} in AD, 27.07 km^{2} in SU, 2.94 km^{2} in SUPG, and 2.92 km^{2} in SUPG+FAB.

We present the numerical errors at the final time step for the parabolic and triangular shape of velocity in Fig. 5 for both low- and high-velocity constants. Apparently, the shape of the velocity profile has a limited impact on the numerical errors. Nevertheless, the triangular velocity cases yield the smallest errors, while the parabolic velocity cases yield larger errors, but these errors are still smaller than the uniform velocity field scenario depicted in Fig. 4.

## 5.1 Reinitialization interval

From a finite-element method point of view, the reinitialization procedure is an ℒ^{2} projection of the zero level-set contour onto the mesh (Larson and Bengzon, 2013). It can be shown that the numerical errors of the projection are proportional to the mesh sizes (shown in Fig. C1), and they accumulate as the number of reinitializations increases (i.e., as *n*_{R} decreases). Furthermore, these errors are not only introduced during the projection but also transported and amplified by the governing equation Eq. (1) throughout the transient simulation. In the case of frequent reinitializations, such as *n*_{R}=1, the dominant source of numerical error is the ℒ^{2} projection, particularly evident when the front velocity is low (*v*_{0}=1000 m a^{−1}), as depicted in Figs. 4 and 5. However, in the high-velocity scenario, the projection error becomes less significant compared to the numerical errors resulting from discretization and stabilization techniques, which then become the primary sources of error.

As *n*_{R} increases, the numerical error decreases until no reinitialization is performed (*n*_{R}=∞). However, in the absence of reinitialization, additional errors emerge due to the distortion of the gradient of the level-set function. The worst-case scenario observed in this study is the high-uniform-velocity case with AD at *n*_{R}=∞ in Fig. 4b, where the zero contour of the final level-set solution is nearly halfway into the fjord, resulting in a total misfit of 84.02 km^{2}. This instance emphasizes the necessity to reinitialize the level-set when solving level-set functions in transient simulations. It is worth noting that the numerical errors are not significantly affected by the interval of reinitialization as long as *n*_{R} is sufficiently larger than 1. Consequently, for the remainder of this paper, the focus will be on discussing the cases with ${n}_{\mathrm{R}}=\mathrm{10},\mathrm{100}$, and 200, while disregarding those with *n*_{R}=1 and *n*_{R}=∞.

As discussed above, FAB penalizes deviations from the eikonal equation, ensuring $\Vert \mathrm{\nabla}\mathit{\varphi}\Vert =\mathrm{1}$ when solving the level set (Hartmann et al., 2010). The reinitialization interval is crucial in determining how often the level set needs to be reset using the geometric reinitialization algorithm described in Sect. 2.2. A naïve approach would be to reinitialize the level set after each step of solving the advection equation in order to maintain its signed distance property. However, in practice, frequent reinitialization introduces interface displacements due to numerical errors, resulting in artificial mass gain or loss, which may also alter the geometrical characteristics of the interface, with potential implications for topological changes (Hartmann et al., 2010; Gibou et al., 2018; Henri et al., 2022).

## 5.2 Stabilization method

The numerical errors in AD and SU are 5 to 20 times greater than those using SUPG and FAB as long as *n*_{R} exceeds 1. The main source of the numerical errors in AD and SU is the diffusion term ∇*ϕ*⋅** κ**∇

*ψ*added in the advection equation Eq. (4), which smears out the oscillations in the numerical solution and disperses the solution. The coefficient

**controls the magnitude and direction of the additional diffusion.**

*κ*In the AD case, the coefficient ** κ** is a scalar, which applies the diffusion to all directions with the same magnitude. In contrast,

**contains an outer product of the front velocity in SU, which only adds diffusion along the flow direction of**

*κ*

*v*_{f}. Therefore, the errors in SU are less dispersive than those in AD. Notably, the coefficients

**in AD and SU are also controlled by the mesh size, such that the additional diffusion term vanishes as the mesh size becomes zero. In numerical ice sheet modeling, the mesh size is generally limited by data accuracy and computational capacity. Therefore, the weak solution of the stabilized equation Eq. (4) does not necessarily satisfy the variational formulation Eq. (3), and the corresponding errors are proportional to the mesh size (Larson and Bengzon, 2013).**

*κ*On the other hand, the SUPG stabilizes the advection equation by adding an additional term in the test function as in Eq. (7), whose solution satisfies the weak form Eq. (3) almost everywhere, except for the position where the test functions are equal to 0. In this sense, the numerical error is expected to be much smaller than the other two stabilization methods. We therefore recommend using SUPG for the stabilization technique, together with a reinitialization interval greater than 10.

## 5.3 Front velocity

We anticipate the numerical errors to be scaled by the velocity magnitude when solving the advection equation using the finite-element method (Biswas et al., 1994) but not influenced by the shape of the calving front. As we construct the velocities in Table 1, for instance, with *v*_{0}=1000 m a^{−1}, the mean frontal velocity during the advance phase $t\in [nT,(n+\frac{\mathrm{1}}{\mathrm{2}}\left)T\right]$ is 1000 m a^{−1} for the uniform shape, 916.7 m a^{−1} for the parabola, and 750.0 m a^{−1} for the triangular shape. The corresponding numerical errors at *n*_{R}=100 with SUPG stabilization are 0.44, 0.36, and 0.29 km^{2}, respectively. Furthermore, as shown in Figs. 4 and 5, this relationship is found in almost all the reinitialization intervals *n*_{R}>1, all stabilization techniques, and both the low- and high-velocity scenarios considered in this study.

Note that while we do not model calving explicitly in this paper, the definition of the frontal velocity in Eq. (2) relies on *v*_{f}, which implicitly incorporates the effects of calving or calving rate. It is important to distinguish that the velocity of the front (*v*_{f}) is not the same as the calving velocity. The frontal velocity *v*_{f} is a sum of the ice speed (which is not necessarily normal to the ice front) and the calving rate *c*, which is generally defined along the normal ** n**. Therefore, the ice front velocity is not necessarily orthogonal to the front in practice.

This study primarily focused on comparing different stabilization and reinitialization strategies for solving the level-set equation, assuming that *v*_{f} is known (i.e., both ice velocity and calving rates are known). The main purpose of this study is to demonstrate that even with a simple prescribed frontal velocity, stabilization and reinitialization can have a significant impact depending on the choices made. Incorporating a realistic calving term may not necessarily provide additional insights into our study, as it is already accounted for through *v*_{f} in the level-set equation. Moreover, introducing a calving law would preclude the availability of analytical solutions, complicating the interpretability of our results. As a future continuation of this study, as part of the CalvingMIP project (https://github.com/JRowanJordan/CalvingMIP/wiki, last access: 20 August 2024) the ice sheet modeling community is testing more realistic calving velocities on more complex geometries, including constant and time-dependent calving rates.

## 5.4 Different front shapes

In Appendix A, we show the results of another shape of the ice front, which is a straight line with side walls orthogonal to the front. The final errors of the straight front cases with different stabilization methods, reinitialization intervals, and velocity shapes are more or less the same as those with the semicircle front. However, the spatial distribution of the numerical error differs significantly between the two shapes. To further investigate the source of the numerical errors, we show the animations of the evolution of misfits in the “Video supplement”. In the straight front cases, the misfit is initiated at the two corners, where the ice front meets the side wall of the fjord, and then propagates to the center. In contrast, the semicircle case generates numerical errors that do not initiate from single sources but grow along the entire ice front. The main reason for these differences is that the finite-element method approximates the level-set function by projecting it onto a piecewise linear functional space. As a result, the sharp corners and the curved level-set contours are the places where most of the numerical errors occur. On average, these approximation errors are proportional to the mesh size, whereas the shape of the ice front actually has a negligible influence on the numerical errors.

We studied multiple stabilization methods implemented in ISSM and Úa for solving a level-set equation on an idealized geometry with a reinitialization interval that varies from once every time step up to no reinitialization. We found that SUPG and SUPG+FAB are considerably more accurate than the other two methods (AD and SU) for all choices of reinitialization interval regardless of the front velocity and ice front shape. Using other stabilization methods results in more than 10 times larger errors in ice front positions. An optimal choice for the reinitialization interval is *n*_{R}>10, corresponding to a time period exceeding 2.5 weeks in our experiments. Excessively frequent reinitialization can introduce additional numerical errors surpassing those from other sources. By identifying the most effective stabilization techniques and reinitialization intervals, we can improve the reliability and robustness of simulations, enabling more accurate predictions of ice sheet behavior and its influence on future sea level rise.

We introduce an alternative ice front shape, represented as a straight line, as depicted in Fig. A1. Similar to Fig. 1, the ice-covered region is denoted in light blue, while the ice-free region is in light red. The red line signifies the ice front, and the blue lines represent the side walls of the fjord, with a width of 10 km and a length of 20 km. The same set of experiments outlined in Sect. 3 is conducted, and the total misfit at the final time step is presented in Fig. A2.

The numerical errors exhibit a linear scaling in time, as illustrated in Figs. B1 and B2 across nearly all cases. As expected, the slopes are dictated by the velocity *v*_{0}. Consequently, for the sake of simplicity in comparison, we exclusively consider the numerical errors at the final time step *T*=50 in the main text of this paper.

We also conducted this study using different mesh resolutions, namely 200 and 400 m, and the corresponding numerical errors are depicted in Fig. C1. To facilitate comparison, we scaled the *y* axis by a factor of 2 and 4 for the two mesh resolutions, respectively. As anticipated, the comparison with results in Figs. 4 and A2 reveals a linear scaling of numerical errors with the mesh size. Notably, in Fig. C1d, *n*_{R}=1 for all four stabilization methods reaches the maximum possible error, equivalent to the area of the fjord in the straight ice front case, i.e., 75 km^{2}.

We investigate the effects of structured meshes on the level-set solutions by conducting two additional sets of experiments with different mesh configurations. One experiment employs a diagonally aligned triangular mesh, while the other relies on a symmetric triangular mesh. The mesh illustrations are presented in Fig. D1, where Fig. D1a represents the unstructured mesh used in this study. Figure D1b depicts a diagonally aligned mesh extending from the top left to the bottom right, while Fig. D1c shows a symmetric triangular mesh.

Figure D2 shows the misfit between the numerical and exact solutions for the diagonally aligned triangular mesh with a mesh resolution of 100 m under a uniform velocity field with *v*_{0}=1000 m a^{−1} observed after 1.5, 2, and 50 periods. Conversely, Fig. D3 presents the results of the same experiment but on a symmetric triangular mesh. There is a clear asymmetry in the results shown in Fig. D2 when using the diagonally aligned triangular mesh but not for the unstructured mesh (e.g., Fig. 2) or the symmetric triangular mesh in Fig. D3.

To further examine the diagonally aligned triangular mesh, we refine the mesh resolution to 50 and 25 m. Figure D4 shows the misfit between the numerical and exact solutions observed after 1.5, 2, and 50 periods with AD stabilization at *n*_{R}=1 under a uniform velocity field with *v*_{0}=1000 m a^{−1}. Although numerical errors decrease with mesh refinement, the asymmetric error patterns persist even at a 25 m resolution, which is nearly the finest mesh resolution used in real-world applications. This experiment highlights the importance of the mesh structure, particularly when geometric reinitialization is performed, as it may significantly depend on the organization and orientation of the elements of the mesh.

We introduce an additional numerical experiment aimed at further improving the realism of the frontal velocity representation. In this experiment, we modify the frontal velocity in Eq. (16) to *v*(*t*)=*v*_{0}sin (2*π**t*), simulating seasonal variations rather than abrupt transitions. This adjustment seeks to emulate the dynamic movement of ice fronts influenced by seasonal changes.

We set *v*_{0}=1000 m a^{−1} with the uniform shape of front velocity and conduct simulations at a mesh resolution of 200 m over a 50-year period. Figure E1 illustrates the evolution of the total absolute misfit and their final values at *T*=50 years.

This experiment exhibits results consistent with other experiments in our study, wherein the SUPG and SUPG+FAB methods with *n*_{R}>10 have the smallest misfit areas among all other methods. In terms of magnitude, as discussed in Sect. 5.3, the errors are scaled by the mean velocity at the front, calculated as $\mathrm{2000}{\int}_{\mathrm{0}}^{\mathrm{1}/\mathrm{2}}\mathrm{sin}\left(\mathrm{2}\mathit{\pi}t\right)\phantom{\rule{0.125em}{0ex}}\text{d}t=\frac{\mathrm{2000}}{\mathit{\pi}}\approx \mathrm{636.56}$ m a^{−1} during each advance phase and −636.56 m a^{−1} for the retreat phase. Consequently, these misfits are approximately 0.64 times of those depicted in Fig. C1a with the same stabilization and reinitialization intervals.

ISSM version 4.23 is open source and available at https://doi.org/10.5281/zenodo.7850841 (ISSM Team, 2023). Úa (v2019b) is open source and available at https://doi.org/10.5281/zenodo.3706624 (Gudmundsson, 2020). The code and data analyses used in this paper are available at https://doi.org/10.5281/zenodo.10454657 (Gong, 2024).

The animations that show the evolution of the misfit are available at https://doi.org/10.5281/zenodo.10454554 (Cheng, 2024).

GC, MM, and GHG designed the study. GC did the numerical computations. GC wrote the manuscript with input from MM and GHG.

The contact author has declared that none of the authors has any competing interests.

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.

This work is from the PROPHET project, a component of the International Thwaites Glacier Collaboration (ITGC), with contribution no. ITGC-130.

This research has been supported by the Heising-Simons Foundation (grant nos. 2019-1161 and 2021-3059), the National Science Foundation (grant no. 1739031), and the Natural Environment Research Council (grant nos. NE/S006745/1, NE/S006796/1, and NE/T001607/1).

This paper was edited by Ludovic Räss and reviewed by Matt Trevers and one anonymous referee.

Benn, D. I. and Astrom, J. A.: Calving glaciers and ice shelves, Adv. Phys.-X, 3, 1513819, https://doi.org/10.1080/23746149.2018.1513819, 2018. a

Biswas, R., Devine, K. D., and Flahert, J. E.: Parallel, adaptive finite element methods for conservation laws, Appl. Numer. Math., 14, 255–283, https://doi.org/10.1016/0168-9274(94)90029-9, 1994. a

Black, T. E. and Joughin, I.: Weekly to monthly terminus variability of Greenland's marine-terminating outlet glaciers, The Cryosphere, 17, 1–13, https://doi.org/10.5194/tc-17-1-2023, 2023. a, b

Bondzio, J., Morlighem, M., Seroussi, H., Kleiner, T., Ruckamp, M., Mouginot, J., Moon, T., Larour, E., and Humbert, A.: The mechanisms behind Jakobshavn Isbræ's acceleration and mass loss: A 3-D thermomechanical model study, Geophys. Res. Lett., 44, 6252–6260, https://doi.org/10.1002/2017GL073309, 2017. a

Bondzio, J. H., Seroussi, H., Morlighem, M., Kleiner, T., Rückamp, M., Humbert, A., and Larour, E. Y.: Modelling calving front dynamics using a level-set method: application to Jakobshavn Isbræ, West Greenland, The Cryosphere, 10, 497–510, https://doi.org/10.5194/tc-10-497-2016, 2016. a, b, c

Brooks, A. N. and Hughes, T. J. R.: Streamline upwind Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations, Comput. Methods Appl. Mech. Engrg., 32, 199–259, 1982. a, b

Cheng, D., Hayes, W., Larour, E., Mohajerani, Y., Wood, M., Velicogna, I., and Rignot, E.: Calving Front Machine (CALFIN): glacial termini dataset and automated deep learning extraction method for Greenland, 1972–2019, The Cryosphere, 15, 1663–1675, https://doi.org/10.5194/tc-15-1663-2021, 2021. a

Cheng, G.: Animations of the evolution of misfits, Zenodo [video], https://doi.org/10.5281/zenodo.10454554, 2024. a

Cheng, G., Morlighem, M., Mouginot, J., and Cheng, D.: Helheim Glacier's Terminus Position Controls Its Seasonal and Inter-Annual Ice Flow Variability, Geophys. Res. Lett., 49, e2021GL097085, https://doi.org/10.1029/2021GL097085, 2022. a

Choi, Y., Morlighem, M., Wood, M., and Bondzio, J. H.: Comparison of four calving laws to model Greenland outlet glaciers, The Cryosphere, 12, 3735–3746, https://doi.org/10.5194/tc-12-3735-2018, 2018. a

Choi, Y., Morlighem, M., Rignot, E., and Wood, M.: Ice dynamics will remain a primary driver of Greenland ice sheet mass loss over the next century, Nature Commun. Earth Environ., 2, 26, https://doi.org/10.1038/s43247-021-00092-z, 2021. a, b

Crawford, A. J., Benn, D. I., Todd, J., Astrom, J. A., Bassis, J. N., and Zwinger, T.: Marine ice-cliff instability modeling shows mixed- mode ice-cliff failure and yields calving rate parameterization, Nat. Commun., 12, 2701, https://doi.org/10.1038/s41467-021-23070-7, 2021. a

dos Santos, T. D., Morlighem, M., and Seroussi, H.: Assessment of numerical schemes for transient, finite-element ice flow models using ISSM v4.18, Geosci. Model Dev., 14, 2545–2573, https://doi.org/10.5194/gmd-14-2545-2021, 2021. a, b, c

Eriksson, K., Estep, D., Hansbo, P., and Johnson, C.: Computational differential equations, Cambridge University Press, ISBN 9780521567381, 1996. a

Gibou, F., Fedkiw, R., and Osher, S.: A review of level-set methods and some recent applications, J. Comput. Phys., 353, 82–109, https://doi.org/10.1016/j.jcp.2017.10.006, 2018. a

Goelzer, H., Nowicki, S., Payne, A., Larour, E., Seroussi, H., Lipscomb, W. H., Gregory, J., Abe-Ouchi, A., Shepherd, A., Simon, E., Agosta, C., Alexander, P., Aschwanden, A., Barthel, A., Calov, R., Chambers, C., Choi, Y., Cuzzone, J., Dumas, C., Edwards, T., Felikson, D., Fettweis, X., Golledge, N. R., Greve, R., Humbert, A., Huybrechts, P., Le clec'h, S., Lee, V., Leguy, G., Little, C., Lowry, D. P., Morlighem, M., Nias, I., Quiquet, A., Rückamp, M., Schlegel, N.-J., Slater, D. A., Smith, R. S., Straneo, F., Tarasov, L., van de Wal, R., and van den Broeke, M.: The future sea-level contribution of the Greenland ice sheet: a multi-model ensemble study of ISMIP6, The Cryosphere, 14, 3071–3096, https://doi.org/10.5194/tc-14-3071-2020, 2020. a

Gong, C.: enigne/Levelset: Level-set Stabilization (v1.1), Zenodo [code and data set], https://doi.org/10.5281/zenodo.10454657, 2024. a

Greve, R. and Blatter, H.: Dynamics of Ice Sheets and Glaciers, Advances in Geophysical and Environmental Mechanics and Mathematics, Springer Science & Business Media, https://doi.org/10.1007/978-3-642-03415-2, 2009. a

Gross, S. and Reusken, A.: Numerical Methods for Two-phase Incompressible Flows, vol. 40 of Springer Series in Computational Mathematics, 1–472, https://doi.org/10.1007/978-3-642-19686-7, 2011. a

Gudmundsson, G. H., Paolo, F. S., Adusumilli, S., and Fricker, H. A.: Instantaneous Antarctic ice sheet mass loss driven by thinning ice shelves, Geophys. Res. Lett., 46, 13903–13909, https://doi.org/10.1029/2019GL085027, 2019. a, b

Gudmundsson, H.: GHilmarG/UaSource: Ua2019b, Zenodo [code], https://doi.org/10.5281/zenodo.3706624, 2020. a, b, c

Hartmann, D., Meinke, M., and Schroeder, W.: The constrained reinitialization equation for level set methods, J. Comput. Phys., 229, 1514–1535, https://doi.org/10.1016/j.jcp.2009.10.042, 2010. a, b

Henri, F., Coquerelle, M., and Lubin, P.: Geometrical level set reinitialization using closest point method and kink detection for thin filaments, topology changes and two-phase flows, J. Comput. Phys., 448, 110704, https://doi.org/10.1016/j.jcp.2021.110704, 2022. a, b

IPCC: Climate Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change, vol. In Press, Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, https://doi.org/10.1017/9781009157896, 2021. a

ISSM Team: Ice-sheet and Sea-level System Model source code, v4.23 r27696, Zenodo [code], https://doi.org/10.5281/zenodo.7850841, 2023. a, b

King, M. D., Howat, I. M., Candela, S. G., Noh, M. J., Jeong, S., Noel, B. P. Y., van den Broeke, M. R., Wouters, B., and Negrete, A.: Dynamic ice loss from the Greenland Ice Sheet driven by sustained glacier retreat, Commun. Earth Environ., 1, 1, https://doi.org/10.1038/s43247-020-0001-2, 2020. a

Larour, E., Seroussi, H., Morlighem, M., and Rignot, E.: Continental scale, high order, high spatial resolution, ice sheet modeling using the Ice Sheet System Model (ISSM), J. Geophys. Res., 117, 1–20, https://doi.org/10.1029/2011JF002140, 2012. a

Larson, M. G. and Bengzon, F.: The Finite Element Method: Theory, Implementation, and Applications, Springer Publishing Company, Incorporated, https://doi.org/10.1007/978-3-642-33287-6, 2013. a, b, c, d

Li, C., Xu, C., Gui, C., and Fox, M. D.: Level set evolution without re-initialization: A new variational formulation, in: 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition, Vol 1, Proceedings, edited by: Schmid, C., Soatto, S., and Tomasi, C., PROCEEDINGS – IEEE Computer Society Conference On Computer Vision And Pattern Recognition, 430–436, IEEE Comp Soc, conference on Computer Vision and Pattern Recognition, 20–25 June 2005, San Diego, CA, 2005. a

MacAyeal, D. R.: Large-scale ice flow over a viscous basal sediment: Theory and application to Ice Stream B, Antarctica, J. Geophys. Res., 94, 4071–4087, 1989. a

Morlighem, M., Bondzio, J., Seroussi, H., Rignot, E., Larour, E., Humbert, A., and Rebuffi, S.-A.: Modeling of Store Gletscher's calving dynamics, West Greenland, in response to ocean thermal forcing, Geophys. Res. Lett., 43, 2659–2666, https://doi.org/10.1002/2016GL067695, 2016. a, b, c, d, e, f

Mouginot, J., Rignot, E., Bjørk, A. A., van den Broeke, M., Millan, R., Morlighem, M., Noël, B., Scheuchl, B., and Wood, M.: Forty-six years of Greenland Ice Sheet mass balance from 1972 to 2018, P. Natl. Acad. Sci. USA., 116, 9239–9244, https://doi.org/10.1073/pnas.1904242116, 2019. a, b

Pattyn, F. and Morlighem, M.: The uncertain future of the Antarctic Ice Sheet, Science, 367, 1331–1335, https://doi.org/10.1126/science.aaz5487, 2020. a

Pollard, D., DeConto, R. M., and Alley, R. B.: Potential Antarctic Ice Sheet retreat driven by hydrofracturing and ice cliff failure, Earth Planet Sc. Lett., 412, 112–121, https://doi.org/10.1016/j.epsl.2014.12.035, 2015. a

Seroussi, H., Nowicki, S., Payne, A. J., Goelzer, H., Lipscomb, W. H., Abe-Ouchi, A., Agosta, C., Albrecht, T., Asay-Davis, X., Barthel, A., Calov, R., Cullather, R., Dumas, C., Galton-Fenzi, B. K., Gladstone, R., Golledge, N. R., Gregory, J. M., Greve, R., Hattermann, T., Hoffman, M. J., Humbert, A., Huybrechts, P., Jourdain, N. C., Kleiner, T., Larour, E., Leguy, G. R., Lowry, D. P., Little, C. M., Morlighem, M., Pattyn, F., Pelle, T., Price, S. F., Quiquet, A., Reese, R., Schlegel, N.-J., Shepherd, A., Simon, E., Smith, R. S., Straneo, F., Sun, S., Trusel, L. D., Van Breedam, J., van de Wal, R. S. W., Winkelmann, R., Zhao, C., Zhang, T., and Zwinger, T.: ISMIP6 Antarctica: a multi-model ensemble of the Antarctic ice sheet evolution over the 21st century, The Cryosphere, 14, 3033–3070, https://doi.org/10.5194/tc-14-3033-2020, 2020. a

Sethian, J. A.: Theory, algorithms, and applications of level set methods for propagating interfaces, Acta Numer., 5, 309–395, 1996. a, b

Sussman, M., Smereka, P., and Osher, S.: A Level set Aspproach for Computing Solutions to Incompressible Two-Phase Flow, J. Comput. Phys., 114, 146–159, https://doi.org/10.1006/jcph.1994.1155, 1994. a

Toure, M. K. and Soulaimani, A.: Stabilized finite element methods for solving the level set equation without reinitialization, Comput. Math. Appl., 71, 1602–1623, https://doi.org/10.1016/j.camwa.2016.02.028, 2016. a, b

Wood, M., Rignot, E., Fenty, I., An, L., Bjørk, A., van den Broeke, M., Cai, C., Kane, E., Menemenlis, D., Millan, R., Morlighem, M., Mouginot, J., Noël, B., Scheuchl, B., Velicogna, I., Willis, J. K., and Zhang, H.: Ocean forcing drives glacier retreat in Greenland, Sci. Adv., 7, 1315–1332, https://doi.org/10.1126/sciadv.aba7282, 2021. a

- Abstract
- Introduction
- Methods
- Numerical experiments
- Results
- Discussion
- Conclusions
- Appendix A: A straight ice front case
- Appendix B: Errors during the transient simulation
- Appendix C: Mesh resolution
- Appendix D: Numerical errors influenced by the mesh structure
- Appendix E: Additional experiment with a more realistic frontal velocity
- Code and data availability
- Video supplement
- Author contributions
- Competing interests
- Disclaimer
- Acknowledgements
- Financial support
- Review statement
- References

- Abstract
- Introduction
- Methods
- Numerical experiments
- Results
- Discussion
- Conclusions
- Appendix A: A straight ice front case
- Appendix B: Errors during the transient simulation
- Appendix C: Mesh resolution
- Appendix D: Numerical errors influenced by the mesh structure
- Appendix E: Additional experiment with a more realistic frontal velocity
- Code and data availability
- Video supplement
- Author contributions
- Competing interests
- Disclaimer
- Acknowledgements
- Financial support
- Review statement
- References