Assessment of numerical schemes for transient, finite-element ice flow models using ISSM v4.18

Time-dependent simulations of ice sheets require two equations to be solved: the mass transport equation, derived from the conservation of mass, and the stress balance equation, derived from the conservation of momentum. The mass transport equation controls the advection of ice from the interior of the ice sheet towards its periphery, thereby changing its geometry. Because it is based on an advection equation, a stabilization scheme needs to be employed when solved using the finite-element method. Several stabilization schemes exist in the finite-element method framework, but their respective accuracy and robustness have not yet been systematically assessed for glaciological applications. Here, we compare classical schemes used in the context of the finite-element method: (i) artificial diffusion, (ii) streamline upwinding, (iii) streamline upwind Petrov–Galerkin, (iv) discontinuous Galerkin, and (v) flux-corrected transport. We also look at the stress balance equation, which is responsible for computing the ice velocity that “advects” the ice downstream. To improve the velocity computation accuracy, the ice-sheet modeling community employs several sub-element parameterizations of physical processes at the grounding line, the point where the grounded ice starts to float onto the ocean. Here, we introduce a new sub-element parameterization for the driving stress, the force that drives the ice-sheet flow. We analyze the response of each stabilization scheme by running transient simulations forced by ice-shelf basal melt. The simulations are based on an idealized ice-sheet geometry for which there is no influence of bedrock topography. We also perform transient simulations of the Amundsen Sea Embayment, West Antarctica, where real bedrock and surface elevations are employed. In both idealized and real ice-sheet experiments, stabilization schemes based on artificial diffusion lead systematically to a bias towards more mass loss in comparison to the other schemes and therefore should be avoided or employed with a sufficiently high mesh resolution in the vicinity of the grounding line. We also run diagnostic simulations to assess the accuracy of the driving stress parameterization, which, in combination with an adequate parameterization for basal stress, provides improved numerical convergence in ice speed computations and more accurate results.

1 Introduction

Numerical modeling is routinely used to understand the past and future behavior of the ice sheets in response to the evolution of the climate (e.g., Ritz et al., 2015; DeConto and Pollard, 2016; Aschwanden et al., 2019; Goelzer et al., 2020; Seroussi et al., 2020). As is always the case with numerical models, one needs to minimize biases, numerical artifacts, or poor numerical convergence due to the choice of numerical scheme. It is therefore critical for the numerical solution to converge to the “true solution” regardless of the numerical scheme employed, and that the model is not overly sensitive to the mesh resolution (e.g., Szabó and Babuška, 1991, p. 4). One of the governing equations in ice-sheet numerical modeling is the mass transport equation, an advection equation derived from the conservation of mass that prescribes the evolution of the ice-sheet geometry (e.g., Cuffey and Paterson, 2010, p. 333). Another governing equation is the stress balance equation, a set of equations that describes the ice velocity over the entire ice sheet. This velocity field is used to “advect” the ice mass over time. These governing equations are often solved using numerical methods such as the finite-element method (FEM), widely employed in the ice-sheet modeling community (e.g., Larour et al., 2012; Gagliardini et al., 2013; Gudmundsson, 2020). As with any numerical method, several schemes exist to solve these equations within FEM in order to achieve stability, accuracy, and low computational cost, the desired properties of any numerical method (see, for example, Szabó and Babuška, 1991, Chapter 1).

It is well known that the discretization of advection-dominated equations by the traditional finite-element method leads to numerical instabilities and spurious oscillations (e.g., John et al., 2018). The mass transport equation therefore needs to be stabilized when solved using FEM. The choice of an adequate stabilization scheme is crucial to simulate the main characteristics of ice-sheet dynamics without introducing numerical artifacts in the solution, such as oscillations, non-physical diffusion, or poor numerical convergence (see, e.g., John and Schmeyer, 2008). While some stabilization techniques may have good properties for some specific applications, they may not be appropriate for others.

The finite-element method's literature presents a large number of stabilization schemes, with different levels of complexity and accuracy (Codina, 1998; Franca et al., 2006; John and Schmeyer, 2008). The simplest schemes are based on adding an artificial diffusive-type term, turning the advection equation into an advection–diffusion equation (LeVeque, 1992, p. 118). This approach is equivalent to the upwind differencing employed in the finite-difference method (Kelly et al., 1980; Selmin, 1993). While this method leads to a more stable formulation, the downside is the reduced accuracy especially in regions where the solution is not smooth or presents steep gradients (Brooks and Hughes, 1982; Donea, 1984b). Other schemes are based on counterbalancing the terms in the finite-element formulation, such that the numerical diffusion vanishes (Brooks and Hughes, 1982). A popular method is the streamline upwind Petrov–Galerkin (SUPG), but the accuracy and stability of this scheme rely on the definition of the stabilization parameter, which is problem-dependent (Almeida and Silva, 1997; Codina, 2000; Knopp et al., 2002; Bochev et al., 2004; John and Schmeyer, 2008; Burman, 2010). Alternatively, a discontinuous Galerkin (DG) formulation, a finite-volume-inspired scheme, produces accurate results in advection-dominated flows (Cockburn, 1998, 1999). However, its implementation requires specific data structures to handle (computational) nodes and (geometric) vertices of the mesh (e.g., Calle et al., 2005; Devloo et al., 2007). Also, the increased number of degrees of freedom introduced by this method may reduce the computational performance. Most recently, finite-element-oriented flux-correction schemes have been proposed as a promising alternative (Kuzmin and Turek, 2002; Kuzmin et al., 2003). These schemes manipulate the discretized algebraic system in order to add anti-diffusive terms without compromising numerical stability. All the schemes described above have been applied to a number of physical problems (e.g., Reed and Hill, 1973; Jameson, 1995; John and Schmeyer, 2008; Ngo et al., 2015; Watanabe and Kolditz, 2015; Diddens, 2017; Hansen et al., 2019). However, the performance of these schemes has not been evaluated in a systematic way for ice-sheet simulations.

The stress balance is another critical component of transient models. For simplified stress balance equations, such as the shallow shelf approximation (MacAyeal, 1989) or Blatter–Pattyn's higher-order models (Blatter, 1995), the right-hand side is a function of the ice surface gradient (Cuffey and Paterson, 2010, p. 295). Usually, in finite-element-based ice-sheet models, gradients are assumed to be continuous within each element. This is a reasonable assumption for most of the ice-sheet domain, except at the grounding line, the point where ice detaches from the underlying bedrock and starts to float over the ocean. From a numerical simulation point of view, the grounding line represents a discontinuity of several physical processes (e.g., basal friction, basal melt), and the accuracy of its dynamics requires a fine mesh resolution (Durand et al., 2009; Pattyn et al., 2012, 2013; Cornford et al., 2013). Sub-element parameterizations of such physical processes are commonly employed to improve numerical convergence of the ice velocity computation (Seroussi et al., 2014a; Feldmann et al., 2014; Cornford et al., 2016; Seroussi and Morlighem, 2018). In most basal friction parameterizations, for instance, the grounding line is free to evolve within the elements. Generally, in such a situation, the models assume the ice thickness to be continuous at the grounding line¹. This assumption implies that, for the grounded part of the element crossed by the grounding line, the ice surface is a function of both bedrock elevation and ice thickness, while, for the floating part, the ice surface is obtained by the hydrostatic floatation only, which only depends on the ice thickness. This makes the gradient of the ice surface and the resulting driving stress discontinuous within the elements containing the grounding line. While there exist comparison studies for basal friction and basal melt parameterizations (e.g., Seroussi et al., 2014a; Seroussi and Morlighem, 2018), little attention has been given to the sub-element parameterization of the driving stress in the context of the finite-element method. To the best of our knowledge, only studies based on finite-volume and finite-difference methods use driving stress parameterizations (Cornford et al., 2013; Feldmann et al., 2014).

In this context, the present paper aims to (i) assess the response of different stabilization schemes in transient simulations subject to ice-shelf basal melt and changes in basal friction, and (ii) develop and assess a sub-element parameterization for the driving stress. The numerical experiments are based on the Marine Ice Sheet Model Intercomparison Project for plan-view models (MISMIP3d) setup (Pattyn et al., 2013), a simple idealized ice-sheet geometry. Additional experiments are presented for the Amundsen Sea Embayment (ASE) of the West Antarctic Ice Sheet (WAIS), which includes the Pine Island and Thwaites glaciers. We use the Ice-sheet and Sea-level System Model (ISSM) to perform all the numerical experiments. In Sects. 2 and 3, we describe the technical details of the stabilization schemes and driving stress parameterization, respectively, and in Sect. 4 we present the numerical setup of experiments used to test them. The results are shown in Sect. 5, followed by discussions in Sect. 6 and final remarks in Sect. 7.

2 Mass transport equation and stabilization schemes

2.1 Mass transport equation

The evolution of the ice thickness is described by an advection equation with source terms on the right-hand side:

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle \frac{\partial H}{\partial t}}+\mathrm{\nabla }\cdot \left(\overline{\mathit{v}}H\right)={\dot{m}}_{\mathrm{s}}-{\dot{m}}_{\mathrm{b}},\end{array}$$

where $\overline{\mathit{v}}=\mathit{\{}{v}_x,v_y\mathit{\}}$ is the depth-averaged ice velocity in the horizontal plane, ${\dot{m}}_{\mathrm{s}}$ the surface mass balance (positive for accumulation), and ${\dot{m}}_{\mathrm{b}}$ the basal melt (positive for ablation). The velocity field $\overline{\mathit{v}}$ is a function of the ice geometry and therefore of the ice thickness H. Note that both surface mass balance ${\dot{m}}_{\mathrm{s}}$ and basal melt ${\dot{m}}_{\mathrm{b}}$ may depend on the surface elevation and ice-shelf depth, respectively. All these dependencies make Eq. (1) a non-linear advection equation. For the sake of simplicity, we keep ${\dot{m}}_{\mathrm{s}}$ as a constant in all transient simulations. The description of ${\dot{m}}_{\mathrm{b}}$ is given in Sect. 4.1 (see Eq. 41).

The weak formulation of Eq. (1) is

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\underset{\mathrm{\Omega }}{\int }\left({\displaystyle \frac{\partial H}{\partial t}}+\mathrm{\nabla }\cdot \left(\overline{\mathit{v}}H\right)\right)\mathit{\psi }\mathrm{d}\mathrm{\Omega }=\underset{\mathrm{\Omega }}{\int }\left({\dot{m}}_{\mathrm{s}}-{\dot{m}}_{\mathrm{b}}\right)\mathit{\psi }\mathrm{d}\mathrm{\Omega },\\ \text{(2)}& {\displaystyle }& {\displaystyle }\forall \mathit{\psi }\in \mathcal{H},\end{array}$$

where ℋ=ℋ(Ω) is a space of admissible functions for the model domain Ω, and ψ is called a test or weight function. We seek a solution H∈ℋ such that the weak form (Eq. 2) is satisfied. In the traditional finite-element method, both H and ψ belong to the same set of functions ℋ. It is known that in this approach Eq. (2) generates potentially large, spurious oscillations if not properly stabilized or if the mesh size is not excessively small (see, for example, Brooks and Hughes, 1982).

The weak form (Eq. 2) (and its alternative stabilized forms) requires approximating functions in ℋ with non-trivial first-order derivatives. In this sense, we employ subspace ℋ¹⊂ℋ whose functions (and their first derivatives) are square integrable. For discretization purposes, we use P1 Lagrange functions and Delaunay-based triangulation.

2.2 Artificial diffusion and streamline upwinding

In general, stabilization schemes may be seen as a consistent way of adding terms to Eq. (1) (or Eq. 2) in order to transform it into a more stable formulation. In the artificial diffusion (ArtDiff; MacAyeal, 1997, p. 172) and streamline upwinding (Streamline; Hughes and Brooks, 1979; Kelly et al., 1980) schemes, the resulting mass transport equation is

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle \frac{\partial H}{\partial t}}+\mathrm{\nabla }\cdot \left(\overline{\mathit{v}}H\right)-\mathrm{\nabla }\cdot \left(\mathcal{D}\mathrm{\nabla }H\right)={\dot{m}}_{\mathrm{s}}-{\dot{m}}_{\mathrm{b}},\end{array}$$

where 𝓓 is a second-order tensor, known as the diffusive tensor. In the artificial diffusion scheme, the tensor is defined as (MacAyeal, 1997, p. 172)

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }\mathcal{D}={\displaystyle \frac{h}{\mathrm{2}}}\left(\begin{array}{cc}\left|v_x\right|& \mathrm{0}\\ \mathrm{0}& \left|v_y\right|\end{array}\right),\end{array}$$

where h is the characteristic size of the element, and v_x and v_y are the horizontal components of the (depth-averaged) ice velocity. In the streamline upwinding method, the tensor is defined in such a way that the artificial diffusion is added only along streamlines and not in cross flows. In this sense, the tensor 𝓓 is defined as (Hughes and Brooks, 1979; Kelly et al., 1980; Brooks and Hughes, 1982)

$$\begin{array}{}\text{(5)}& {\displaystyle }{\displaystyle }\mathcal{D}={\displaystyle \frac{h}{\mathrm{2}\parallel \overline{\mathit{v}}\parallel }}\overline{\mathit{v}}\otimes \overline{\mathit{v}}={\displaystyle \frac{h}{\mathrm{2}\parallel \overline{\mathit{v}}\parallel }}\left(\begin{array}{cc}{v}_x^{\mathrm{2}}& v_x{v}_y\\ v_y{v}_x& v_y^{\mathrm{2}}\end{array}\right),\end{array}$$

where $\parallel \cdot \parallel$ is the Euclidean norm. The resulting weak formulation of Eq. (3), after integrating the diffusive term by parts, is

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\underset{\mathrm{\Omega }}{\int }\left({\displaystyle \frac{\partial H}{\partial t}}+\mathrm{\nabla }\cdot \left(\overline{\mathit{v}}H\right)\right)\mathit{\psi }\mathrm{d}\mathrm{\Omega }+\underset{\mathrm{\Omega }}{\int }\mathrm{\nabla }\mathit{\psi }\cdot \mathcal{D}\mathrm{\nabla }H\mathrm{d}\mathrm{\Omega }\\ {\displaystyle }& {\displaystyle }=\underset{\mathrm{\Omega }}{\int }\left({\dot{m}}_{\mathrm{s}}-{\dot{m}}_{\mathrm{b}}\right)\mathit{\psi }\mathrm{d}\mathrm{\Omega },\\ \text{(6)}& {\displaystyle }& {\displaystyle }\forall \mathit{\psi }\in \mathcal{H}.\end{array}$$

Both schemes are interpreted as an upwind-equivalent scheme employed in the finite-difference method (Kelly et al., 1980; Selmin, 1993), sharing therefore similar characteristics such as a first-order accuracy and large numerical dissipation (de Vahl Davis and Mallinson, 1976; Gresho and Lee, 1979). However, the resulting formulation is very stable, which made them popular in glaciology (e.g., MacAyeal, 1997, p. 172). Artificial diffusion was the default scheme employed in ISSM. The low accuracy of such schemes leads to the development of alternative methods with higher accuracy, such as the streamline upwind Petrov–Galerkin.

2.3 Streamline upwind Petrov–Galerkin

In the SUPG scheme (Brooks and Hughes, 1982), the diffusive or upwind term is not added directly into Eq. (1). Instead, the upwind effect is achieved by adding the upwinding term into the test function ψ (Christie et al., 1976; Heinrich et al., 1977). This procedure, where the solution and test functions belong to different spaces, is commonly called Petrov–Galerkin method (see, e.g., Griffiths and Lorenz, 1978; Brooks and Hughes, 1982). The modified test functions $\widehat{\mathit{\psi }}$ are generally defined as (Brooks and Hughes, 1982)

$$\begin{array}{}\text{(7)}& {\displaystyle }{\displaystyle }\widehat{\mathit{\psi }}=\mathit{\psi }+\mathit{\tau }\mathrm{\nabla }\cdot \left(\overline{\mathit{v}}\mathit{\psi }\right),\end{array}$$

where τ is a stabilization coefficient, defined below. Using the modified test functions $\widehat{\mathit{\psi }}$, the resulting weak formulation of Eq. (1) is written as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\underset{\mathrm{\Omega }}{\int }\left({\displaystyle \frac{\partial H}{\partial t}}+\mathrm{\nabla }\cdot \left(\overline{\mathit{v}}H\right)\right)\mathit{\psi }\mathrm{d}\mathrm{\Omega }\\ {\displaystyle }& {\displaystyle }+\underset{\mathrm{\Omega }}{\int }\left({\displaystyle \frac{\partial H}{\partial t}}+\mathrm{\nabla }\cdot \left(\overline{\mathit{v}}H\right)\right)\mathit{\tau }\mathrm{\nabla }\cdot \left(\overline{\mathit{v}}\mathit{\psi }\right)\mathrm{d}\mathrm{\Omega }\\ {\displaystyle }& {\displaystyle }=\underset{\mathrm{\Omega }}{\int }\left({\dot{m}}_{\mathrm{s}}-{\dot{m}}_{\mathrm{b}}\right)\mathit{\psi }\mathrm{d}\mathrm{\Omega }+\underset{\mathrm{\Omega }}{\int }\left({\dot{m}}_{\mathrm{s}}-{\dot{m}}_{\mathrm{b}}\right)\mathit{\tau }\mathrm{\nabla }\cdot \left(\overline{\mathit{v}}\mathit{\psi }\right)\mathrm{d}\mathrm{\Omega },\\ \text{(8)}& {\displaystyle }& {\displaystyle }\forall \mathit{\psi }\in \mathcal{H}.\end{array}$$

The most common definition of the stabilization term τ is given by (Franca et al., 2006)

$$\begin{array}{}\text{(9)}& {\displaystyle }{\displaystyle }\mathit{\tau }={\displaystyle \frac{h}{\mathrm{2}\parallel \overline{\mathit{v}}\parallel }}\left(\coth \left(P_e\right)-{\displaystyle \frac{\mathrm{1}}{P_e}}\right),\end{array}$$

with

$$\begin{array}{}\text{(10)}& {\displaystyle }{\displaystyle }{P}_e={\displaystyle \frac{h\parallel \overline{\mathit{v}}\parallel }{\mathrm{2}\mathit{\kappa }}},\end{array}$$

where P_e is the Péclet number of element e, and κ is known as the diffusion coefficient. The diffusion coefficient κ is not explicitly defined in Eq. (1). For fast ice streams where the flow is dominated by basal sliding rather than internal deformation², an alternative is to assume the asymptotic limit of the term within parentheses in Eq. (9) when κ and consequently P_e go to infinity. With this assumption, the stabilization coefficient is approximated by

$$\begin{array}{}\text{(11)}& {\displaystyle }{\displaystyle }\mathit{\tau }\sim {\mathit{\tau }}_{\mathrm{1}}={\displaystyle \frac{h}{\mathrm{2}\parallel \overline{\mathit{v}}\parallel }}.\end{array}$$

An alternative to defining τ is to replace the Péclet number by the Courant–Friedrichs–Lewy (CFL) number C_e in Eq. (9) (Gudmundsson, 2020, p. 19, 20), i.e.,

$$\begin{array}{}\text{(12)}& {\displaystyle }{\displaystyle }\mathit{\tau }={\displaystyle \frac{h}{\mathrm{2}\parallel \overline{\mathit{v}}\parallel }}\left(\coth \left(C_e\right)-{\displaystyle \frac{\mathrm{1}}{C_e}}\right),\end{array}$$

with

$$\begin{array}{}\text{(13)}& {\displaystyle }{\displaystyle }{C}_e={\displaystyle \frac{\parallel \overline{\mathit{v}}\parallel \mathrm{\Delta }t}{h}},\end{array}$$

where Δt is the simulation time step. The motivation behind this approach is the equivalence between the SUPG and Taylor–Galerkin methods in transient problems (see, e.g., Donea, 1984a; Donea et al., 1984; Codina, 1998; Blank et al., 1999; Akin and Tezduyar, 2004; Kuzmin, 2010, p. 73). For C_e≪1, which is basically true in ice-sheet simulations³, the term within parentheses in Eq. (12) is approximately $C_e/\mathrm{3}$ (Gudmundsson, 2020, p. 20). Therefore, for low CFL numbers, we can approximate Eq. (12) by

$$\begin{array}{}\text{(14)}& {\displaystyle }{\displaystyle }\mathit{\tau }\sim {\mathit{\tau }}_{\mathrm{2}}={\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{6}}}.\end{array}$$

In ISSM, we implement τ as given by Eq. (11), i.e., τ=τ₁. This definition is also employed in Elmer/Ice, a popular FEM-based ice-sheet model (Gagliardini et al., 2013). However, for some simulations performed here, we use the approximation given by Eq. (14) (τ=τ₂), as implemented in Úa (Gudmundsson, 2020), another popular ice-sheet model.

2.4 Discontinuous Galerkin

Strictly speaking, DG is not exactly a stabilization scheme in the sense of adding upwinding terms to Eq. (1) or to the space of test functions ψ. It is a variant of the traditional continuous Galerkin method, in which the functions approximating the solution are discontinuous across elements' edges. This adds the advantage of local conservation with typical FEM characteristics such as the ease of dealing with non-structured meshes, complex domains, and boundary conditions, and hp adaptivity (Cockburn, 1998; Kuzmin, 2010, p. 84). With appropriate definition of the numerical inter-element fluxes, the resulting formulation is known to be stable in advection-dominated problems (e.g., Brezzi et al., 2004). It was originally developed for the neutron transport simulation (Reed and Hill, 1973), and since then DG has been applied to other fields, including elliptic-type equations (e.g., Babuška et al., 1999; Brezzi et al., 2000; Arnold et al., 2002). A more complete history and an overview of the vast application fields of DG can be found in Cockburn (1998, 2003).

In DG, the weak formulation is written in an element-wise fashion, and the advection operator is integrated by parts such that Eq. (2) is rewritten as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\underset{{\mathrm{\Omega }}^e}{\int }\left({\displaystyle \frac{\partial H}{\partial t}}\mathit{\psi }-\left(\overline{\mathit{v}}H\right)\cdot \mathrm{\nabla }\mathit{\psi }\right)\mathrm{d}{\mathrm{\Omega }}^e+\underset{S^e}{\int }\left(\overline{\mathit{v}}H\mathit{\psi }\right)\cdot \mathit{n}\mathrm{d}{S}^e\\ {\displaystyle }& {\displaystyle }=\underset{{\mathrm{\Omega }}^e}{\int }\left({\dot{m}}_{\mathrm{s}}-{\dot{m}}_{\mathrm{b}}\right)\mathit{\psi }\mathrm{d}{\mathrm{\Omega }}^e,\\ \text{(15)}& {\displaystyle }& {\displaystyle }\forall \mathit{\psi }\in \mathcal{H},\end{array}$$

where Ω^e∈Ω is the domain of element e, S^e is the boundary of Ω^e, ℋ is built on the element's domain (i.e., ℋ=ℋ(Ω^e)), and n is a unit vector pointing outward along the boundary S^e. The stabilized version of Eq. (15) comes from the definition of the numerical flux ($\overline{\mathit{v}}H$) in the second integral on the left-hand side. In ISSM, an upwinding numerical flux is employed such that Eq. (15) is rewritten as (Brezzi et al., 2004)

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\underset{{\mathrm{\Omega }}^e}{\int }\left({\displaystyle \frac{\partial H}{\partial t}}\mathit{\psi }-\left(\overline{\mathit{v}}H\right)\cdot \mathrm{\nabla }\mathit{\psi }\right)\mathrm{d}{\mathrm{\Omega }}^e+\underset{S^e}{\int }{\left(\overline{\mathit{v}}H\right)}_{\mathrm{u}}\cdot \left(\mathit{\psi }\mathit{n}\right)\mathrm{d}{S}^e\\ {\displaystyle }& {\displaystyle }=\underset{{\mathrm{\Omega }}^e}{\int }\left({\dot{m}}_{\mathrm{s}}-{\dot{m}}_{\mathrm{b}}\right)\mathit{\psi }\mathrm{d}{\mathrm{\Omega }}^e,\\ \text{(16)}& {\displaystyle }& {\displaystyle }\forall \mathit{\psi }\in \mathcal{H},\end{array}$$

where $(\overline{\mathit{v}}H{)}_{\mathrm{u}}$ is the upwinding numerical flux.

As seen in Sect. 2.6, the time-derivative discretization of Eq. (16) in ISSM relies on a backward Euler scheme, which differs from the popular Runge–Kutta discontinuous Galerkin (RKDG; Cockburn and Shu, 1991; Cockburn, 1998, 2003), an explicit time projection scheme that is known for its stability (as long as $\text{CFL}<\mathrm{1}/\mathrm{3}$) and allows full parallelization since the resulting mass matrix is block diagonal. The DG implementation in ISSM was conceived to be an alternative to other stabilization schemes without enforcing large changes in an existing code. Therefore, we do not expect to achieve the same benefits as RKDG (and similar schemes) in our simulations. The Elmer/Ice model (Gagliardini et al., 2013) adopts a time discretization similar to the implementation of ISSM.

2.5 Flux-corrected transport

The flux-corrected transport (FCT) scheme operates in the resulting algebraic system of the traditional Galerkin discretization (i.e., Eq. 2) instead of modifying its weak form or the approximation/trial spaces. The scheme was developed to solve the continuity equation for compressible fluids in a finite-difference framework (Boris and Book, 1973), and it was extended to FEM by Löhner et al. (1987) and most recently by Kuzmin and Turek (2002) and Kuzmin et al. (2003). The latter is also named as FEM-FCT in the literature. The FEM-FCT scheme seems to be stable even in the presence of steep gradients and discontinuities (e.g., John and Schmeyer, 2008). For simplicity, we refer to the FEM-FCT scheme as FCT. This scheme can be described as a hybrid method combining a higher-order (but potentially oscillatory) scheme and a lower-order diffusive (not oscillatory) scheme. In regions where the solution is smooth, the higher-order scheme is applied, while the lower-order scheme is applied where necessary (e.g., steep gradient regions). However, this scheme adds just enough flux into the “low-order region” from the high-order (a.k.a. anti-diffusive) region to maintain its accuracy without inducing oscillations.

The scheme modifies the discrete form of Eq. (2) by employing a generic finite-difference scheme for the time derivative, i.e.,

$$\begin{array}{}\text{(17)}& {\displaystyle }{\displaystyle }\left({\mathbf{M}}_{\mathrm{C}}+\mathit{\theta }\mathrm{\Delta }t\mathbf{K}\right){\mathit{H}}^{n+\mathrm{1}}=\left[{\mathbf{M}}_{\mathrm{C}}-\left(\mathrm{1}-\mathit{\theta }\right)\mathrm{\Delta }t\mathbf{K}\right]{\mathit{H}}^n+\mathrm{\Delta }t\mathit{F},\end{array}$$

where Δt is the time step, M_C is the mass matrix⁴, K is the advection matrix, and F is the forcing vector. The superscripts n+1 and n indicate the next and current steps, and the fractional weight θ stands for the backward Euler scheme if θ=1 and for Crank–Nicolson if $\mathit{\theta }=\mathrm{1}/\mathrm{2}$. As explained in Sect. 2.6, the load vector F as well as the velocity field in K are defined in step n. To simplify the notation, we drop the superscripts of these variables in this section.

The first step consists of turning Eq. (17) into a stable, low-order algebraic system. This is achieved by replacing the consistent mass matrix M_C by a lumped mass matrix M_L and the advection matrix K by a matrix L:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathbf{M}}_{\mathrm{C}}=\text{diag}\left\{m_i\right\},\\ \text{(18)}& {\displaystyle }& {\displaystyle }{m}_i=\sum _{j=\mathrm{1}}^N{m}_{ij},\end{array}$$

where N is the total number of degrees of freedom. The matrix L is defined as

$$\begin{array}{}\text{(19)}& {\displaystyle }{\displaystyle }\mathbf{L}=\mathbf{K}+\mathbf{D},\end{array}$$

where D={d_ij} represents an artificial diffusion with elements defined as (Kuzmin and Turek, 2002)

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{d}_{ij}=d_{ji}=-\text{max}\left\{\mathrm{0},k_{ij},k_{ji}\right\},\forall i\ne j,\\ \text{(20)}& {\displaystyle }& {\displaystyle }{d}_{ii}=-\sum _{j=\mathrm{1},j\ne i}^N{d}_{ij},\end{array}$$

where k_ij are the elements of matrix K (={k_ij}). By construction, the matrix L does not contain any positive off-diagonal elements. The resulting stable, low-order system of equations is

$$\begin{array}{}\text{(21)}& {\displaystyle }{\displaystyle }\left({\mathbf{M}}_{\mathrm{L}}+\mathit{\theta }\mathrm{\Delta }t\mathbf{L}\right){\mathit{H}}^{n+\mathrm{1}}=\left[{\mathbf{M}}_{\mathrm{L}}-\left(\mathrm{1}-\mathit{\theta }\right)\mathrm{\Delta }t\mathbf{L}\right]{\mathit{H}}^n+\mathrm{\Delta }t\mathit{F}.\end{array}$$

Compared to the original system (Eq. 17), the modified system (Eq. 21) creates large numerical diffusion that prevents spurious oscillations. By doing so, however, it also reduces the accuracy of solution Hⁿ⁺¹.

In order to improve the accuracy of the solution while still preventing spurious oscillations, the second step of the scheme consists in adding an anti-diffusive term to the right-hand side of Eq. (21):

$$\begin{array}{}\text{(22)}& {\displaystyle }{\displaystyle }\left({\mathbf{M}}_{\mathrm{L}}+\mathit{\theta }\mathrm{\Delta }t\mathbf{L}\right){\mathit{H}}^{n+\mathrm{1}}=\left[{\mathbf{M}}_{\mathrm{L}}-\left(\mathrm{1}-\mathit{\theta }\right)\mathbf{L}\right]{\mathit{H}}^n+\mathrm{\Delta }t\mathit{F}+{\mathit{F}}^{*},\end{array}$$

where F^* is a vector whose elements $f_i^{*}$ are defined as

$$\begin{array}{}\text{(23)}& {\displaystyle }{\displaystyle }{f}_i^{*}=\sum _{j=\mathrm{1}}^N{\mathit{\alpha }}_{ij}{r}_{ij},i=\mathrm{1},\mathrm{\dots },N.\end{array}$$

In Eq. (23), ${\mathit{\alpha }}_{ij}\in [\mathrm{0},\mathrm{1}]$ are weights to be defined appropriately (see below and Appendix A), and R={r_ij} is based on the residual vector R between Eqs. (21) and (17):

$$\begin{array}{ll}{\displaystyle }\mathit{R}=& {\displaystyle }\left({\mathbf{M}}_{\mathrm{L}}-{\mathbf{M}}_{\mathrm{C}}\right)\left({\mathit{H}}^{n+\mathrm{1}}-{\mathit{H}}^n\right)\\ \text{(24)}& {\displaystyle }& {\displaystyle }+\mathrm{\Delta }t\mathbf{D}\left[\mathit{\theta }{\mathit{H}}^{n+\mathrm{1}}+\left(\mathrm{1}-\mathit{\theta }\right){\mathit{H}}^n\right].\end{array}$$

The residual vector R can be decomposed as (see, e.g., Kuzmin, 2009)

$$\begin{array}{}\text{(25)}& {\displaystyle }{\displaystyle }{r}_i=\sum _{j=\mathrm{1},j\ne i}^N{r}_{ij},r_{ji}=-r_{ij},\end{array}$$

where r_ij represents the raw anti-diffusive flux from node j into node i. Using Eqs. (18) and (20), r_ij can be written as (John and Schmeyer, 2008)

$$\begin{array}{ll}{\displaystyle }{r}_{ij}=& {\displaystyle }{m}_{ij}\left(H_i^{n+\mathrm{1}}-H_j^{n+\mathrm{1}}\right)-m_{ij}\left(H_i^n-H_j^n\right)\\ {\displaystyle }& {\displaystyle }-\mathrm{\Delta }t\mathit{\theta }{d}_{ij}\left(H_i^{n+\mathrm{1}}-H_j^{n+\mathrm{1}}\right)\\ \text{(26)}& {\displaystyle }& {\displaystyle }-\mathrm{\Delta }t\left(\mathrm{1}-\mathit{\theta }\right)d_{ij}\left(H_i^n-H_j^n\right).\end{array}$$

Note that r_ij depends on the solution Hⁿ⁺¹. Roughly speaking, there are two approaches to proceed with Eq. (26): (i) a non-linear algorithm and (ii) a linear algorithm. Here, we describe the latter (Kuzmin, 2009), which is currently implemented in ISSM. Further details of both approaches can be found in John and Schmeyer (2008) or in Kuzmin (2009).

In the linear FCT algorithm, the solution Hⁿ⁺¹ in Eq. (26) is replaced by the solution H_L obtained in the low-order system (Eq. 21), i.e.,

$$\begin{array}{}\text{(27)}& {\displaystyle }{\displaystyle }\left({\mathbf{M}}_{\mathrm{L}}+\mathit{\theta }\mathrm{\Delta }t\mathbf{L}\right){\mathit{H}}_{\mathrm{L}}=\left[{\mathbf{M}}_{\mathrm{L}}-\left(\mathrm{1}-\mathit{\theta }\right)\mathrm{\Delta }t\mathbf{L}\right]{\mathit{H}}^n+\mathrm{\Delta }t\mathit{F}.\end{array}$$

In ISSM, the Crank–Nicolson scheme (i.e., $\mathit{\theta }=\mathrm{1}/\mathrm{2}$) is used in Eq. (27), and the raw anti-diffusive flux (Eq. 26) is replaced by an alternative form (Kuzmin, 2009):

$$\begin{array}{}\text{(28)}& {\displaystyle }{\displaystyle }{r}_{ij}=m_{ij}\left({\dot{H}}_{\mathrm{L},i}-{\dot{H}}_{\mathrm{L},j}\right)+d_{ij}\left(H_{\mathrm{L},i}-H_{\mathrm{L},j}\right),\end{array}$$

where ${\dot{\mathit{H}}}_{\mathrm{L}}$ is an approximation of the time derivative $\partial \mathit{H}/\partial t$. This approximation is computed using Richardson's iteration:

$$\begin{array}{}\text{(29)}& {\displaystyle }{\displaystyle }{\dot{\mathit{H}}}_{\mathrm{L}}^{m+\mathrm{1}}={\dot{\mathit{H}}}_{\mathrm{L}}^m+{\mathbf{M}}_{\mathrm{L}}^{-\mathrm{1}}\left(\mathbf{K}{\mathit{H}}_{\mathrm{L}}-{\mathbf{M}}_{\mathrm{C}}{\dot{\mathit{H}}}_{\mathrm{L}}^m\right),m=\mathrm{0},\mathrm{1},\mathrm{\dots },\end{array}$$

with ${\dot{\mathit{H}}}_{\mathrm{L}}^{\mathrm{0}}=\mathbf{0}$. The convergence of Eq. (29) usually takes one to five iterations (Kuzmin, 2009).

Once the anti-diffusive flux is obtained by Eq. (28) using Eqs. (27) and (29), the last step is to compute the solution Hⁿ⁺¹. In this linearized FCT version, the solution Hⁿ⁺¹ is explicitly obtained by solving (Kuzmin, 2009)

$$\begin{array}{}\text{(30)}& {\displaystyle }{\displaystyle }{\mathbf{M}}_{\mathrm{L}}{\mathit{H}}^{n+\mathrm{1}}={\mathbf{M}}_{\mathrm{L}}{\mathit{H}}_{\mathrm{L}}+\mathrm{\Delta }t{\mathit{F}}^{*},\end{array}$$

where F^* is obtained by Eq. (23) with r_ij computed by Eq. (28). The weights α_ij are obtained using the so-called Zalesak algorithm (Zalesak, 1979). Appendix A presents the algorithm as implemented in ISSM.

2.6 Time discretization of the mass transport equation

We employ a semi-implicit finite-difference time-stepping scheme to solve the temporal evolution of the ice thickness. This scheme involves a backward Euler method⁵ for the time derivative in Eq. (1), but the other variables (velocity, surface mass balance, and basal melting) are based on the previous time step. To illustrate this scheme, we apply this time-derivative discretization in Eq. (6), as follows:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\underset{\mathrm{\Omega }}{\int }\left({\displaystyle \frac{H^{n+\mathrm{1}}-H^n}{\mathrm{\Delta }t}}+\mathrm{\nabla }\cdot \left({\overline{\mathit{v}}}^n{H}^{n+\mathrm{1}}\right)\right)\mathit{\psi }\mathrm{d}\mathrm{\Omega }\\ {\displaystyle }& {\displaystyle }+\underset{\mathrm{\Omega }}{\int }\mathrm{\nabla }\mathit{\psi }\cdot {\mathcal{D}}^n\mathrm{\nabla }{H}^{n+\mathrm{1}}\mathrm{d}\mathrm{\Omega }\\ {\displaystyle }& {\displaystyle }=\underset{\mathrm{\Omega }}{\int }\left({\dot{m}}_{\mathrm{s}}^n-{\dot{m}}_{\mathrm{b}}^n\right)\mathit{\psi }\mathrm{d}\mathrm{\Omega },\\ \text{(31)}& {\displaystyle }& {\displaystyle }\forall \mathit{\psi }\in \mathcal{H},\end{array}$$

where Δt is the time step, and superscripts n and n+1 indicate the current and next simulation steps, respectively. Algorithm 1 presents the solution sequence employed for the mass transport computation.

3 Sub-element parameterization of driving stress

The position of the grounding line in non-full Stokes models is generally tracked with a level set condition based on a floatation criterion (Seroussi et al., 2014a). This level set can be located anywhere in an element, so it does not necessarily coincide with elements' edges:

$$\begin{array}{}\text{(32)}& {\displaystyle }{\displaystyle }\mathit{\phi }=H+{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{w}}}{\mathit{\rho }}}r,\end{array}$$

where H is the ice thickness, ρ and ρ_w are ice and ocean densities, respectively, and r is the bedrock elevation (negative if below sea level). The ice is grounded if φ>0; otherwise, it is floating. The grounding line is implicitly defined where φ=0. In the element e containing the grounding line, the ice surface, s, is recovered as follows:

$$\begin{array}{}\text{(33)}& {\displaystyle }{\displaystyle }s\left(x,y\right)=\left\{\begin{array}{ll}H\left(x,y\right)+r\left(x,y\right),& \text{if}\mathit{\phi }\left(x,y\right)>\mathrm{0},\\ H\left(x,y\right)\left(\mathrm{1}-\mathit{\rho }/{\mathit{\rho }}_{\mathrm{w}}\right),& \text{if}\mathit{\phi }\left(x,y\right)\le \mathrm{0}.\end{array}\right.\end{array}$$

The second condition in Eq. (33) guarantees the continuity of the ice surface at the grounding line. However, its gradient is discontinuous within the element: in the grounded part (i.e., φ>0), the surface gradient is a function of both thickness and bedrock elevation, whereas in the floating part (φ<0), it is proportional to the thickness gradient only. The driving stress is therefore also discontinuous in elements partially floating and partially grounded, and we propose to use a sub-element driving stress parameterization to account for this discontinuity. A similar approach to subgrid cells was proposed (Cornford et al., 2013; Feldmann et al., 2014) based on finite-volume or finite-difference methods by applying one-sided differences to compute surface gradients on each side of the grounding line.

The driving stress parameterization is based on recovering the ice surface, and consequently its gradient, on the element e containing the grounding line. We divide the element domain Ω^e in two subdomains: ${\mathrm{\Omega }}_{\mathrm{G}}^e$ and ${\mathrm{\Omega }}_{\mathrm{F}}^e$ that are the grounded and floating parts of the element, respectively (i.e., ${\mathrm{\Omega }}^e={\mathrm{\Omega }}_{\mathrm{G}}^e\cup {\mathrm{\Omega }}_{\mathrm{F}}^e$ and ${\mathrm{\Omega }}_{\mathrm{G}}^e\cap {\mathrm{\Omega }}_{\mathrm{F}}^e=\mathit{\varnothing }$). We then perform the numerical integration of the driving stress on these two subdomains, i.e.,

$$\begin{array}{}\text{(34)}& {\displaystyle }{\displaystyle }{\mathit{F}}^e=\underset{{\mathrm{\Omega }}^e}{\int }{\mathbf{A}}^e\mathrm{\nabla }s\mathrm{d}{\mathrm{\Omega }}^e=\underset{{\mathrm{\Omega }}_{\mathrm{G}}^e}{\int }{\mathbf{A}}^e\mathrm{\nabla }s\mathrm{d}{\mathrm{\Omega }}_{\mathrm{G}}^e+\underset{{\mathrm{\Omega }}_{\mathrm{F}}^e}{\int }{\mathbf{A}}^e\mathrm{\nabla }s\mathrm{d}{\mathrm{\Omega }}_{\mathrm{F}}^e,\end{array}$$

where F^e is the element load vector, A^e is a matrix representing the rest of components of the driving stress and element basis functions⁶, and ∇s is evaluated according to the recovered ice surface (Eq. 33). For comparison purposes, we named the proposed driving stress parameterization SED2 (sub-element driving stress 2), since the approach here is similar to the basal friction parameterization SEP2 developed by Seroussi et al. (2014a). The non-parameterized case, i.e., when the ice surface is evaluated on the elements' vertices and then linearly interpolated elsewhere within the elements, regardless of the grounding line position, is referred to as no driving stress parameterization (NSED). Mathematically, the ice surface in the NSED scheme is defined as $s=H\left(\mathrm{1}-\mathit{\rho }/{\mathit{\rho }}_{\mathrm{w}}\right)$ over Ω^e, and the resulting driving stress is proportional to ∇H.

4 Numerical experiments

4.1 MISMIP3d – numerical setup

In this section, we describe the idealized geometry experiments used to evaluate the stabilization schemes and the proposed driving stress parameterization. For the latter, we employ different parameterization schemes for basal friction. The list of all the schemes tested is summarized in Table 1.

(Seroussi et al., 2014a)(Seroussi et al., 2014a)

Table 1 List of the numerical schemes analyzed in this work.

The numerical experiments are based on the MISMIP3d setup (Pattyn et al., 2013). The ice-sheet flows along the x axis in an 800×50 km² rectangular domain entirely filled with ice. For the ice divide, x=0, we set v_x=0. A free-slip condition is applied to the lateral boundaries of the domain, i.e., v_y=0 for y=0 and y=50 km. The calving front is fixed in time and located at x=800 km, where we apply a Neumann boundary condition based on ocean water pressure. We employ a Weertman-type friction law given by

$$\begin{array}{}\text{(35)}& {\displaystyle }{\displaystyle }{\mathit{\tau }}_{\mathrm{b}}=C\parallel {\mathit{v}}_{\mathrm{b}}{\parallel }^{m-\mathrm{1}}{\mathit{v}}_{\mathrm{b}},\end{array}$$

where τ_b is the basal friction, C is the friction coefficient, v_b is the basal velocity, and m is the sliding law exponent. The bed elevation is defined as

$$\begin{array}{}\text{(36)}& {\displaystyle }{\displaystyle }r(x,y)=\mathrm{100}-x,\end{array}$$

where r(x,y) is the bedrock elevation (in meters; negative if below sea level), and x∈[0 800] is the x coordinate in kilometers. All experiments start from the same initial geometry, defined by the following ice thickness profile:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }H(x,y,t=\mathrm{0})\\ \text{(37)}& {\displaystyle }& {\displaystyle }=\left\{\begin{array}{l}{\left[H_{\text{gl}}^{m+\mathrm{2}}+\frac{m+\mathrm{2}}{m+\mathrm{1}}\frac{C{\dot{m}}_{\mathrm{s}}^m}{\mathit{\rho }g}\left(x_{\text{gl}}^{m+\mathrm{1}}-x^{m+\mathrm{1}}\right)\right]}^{\frac{\mathrm{1}}{m+\mathrm{2}}},\\ \text{if }x<x_{\text{gl}},\\ \\ {\left[\frac{\mathcal{A}}{{\dot{m}}_{\mathrm{s}}}-\frac{v_{\text{gl}}^{n+\mathrm{1}}\left(\frac{\mathcal{A}}{{\dot{m}}_{\mathrm{s}}}{H}_{\text{gl}}^{n+\mathrm{1}}-\mathrm{1}\right)}{{\left[{\dot{m}}_{\mathrm{s}}\left(x-x_{\text{gl}}\right)+v_{\text{gl}}{H}_{\text{gl}}\right]}^{n+\mathrm{1}}}\right]}^{-\frac{\mathrm{1}}{n+\mathrm{1}}},\\ \text{if }x\ge x_{\text{gl}},\end{array}\right.\end{array}$$

where H is the ice thickness, ${\dot{m}}_{\mathrm{s}}$ is the accumulation rate, x_gl is the grounding line position (x axis), n is the Glen law exponent, H_gl and v_gl are the ice thickness and the ice velocity (x direction) at the grounding line, respectively, defined as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{H}_{\text{gl}}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{w}}}{\mathit{\rho }}}\left|r(x_{\text{gl}},y)\right|,\\ \text{(38)}& {\displaystyle }& {\displaystyle }{v}_{\text{gl}}={\displaystyle \frac{{\dot{m}}_{\mathrm{s}}{x}_{\text{gl}}}{H_{\text{gl}}}},\end{array}$$

and 𝒜 is defined by

$$\begin{array}{}\text{(39)}& {\displaystyle }{\displaystyle }\mathcal{A}=A{\left({\displaystyle \frac{\mathit{\rho }g\left({\mathit{\rho }}_{\mathrm{w}}-\mathit{\rho }\right)}{\mathrm{4}{\mathit{\rho }}_{\mathrm{w}}}}\right)}^n,\end{array}$$

where A is the Glen law rate factor. The parameters used in all experiments are summarized in Table 2. In Eq. (37), the thickness expression for x<x_gl is the steady-state profile of a 1-D ice sheet considering uniform accumulation rate and negligible longitudinal stresses (Schoof, 2007a). The x≥x_gl case is the steady-state profile of an unconfined 1-D ice shelf under a uniform accumulation rate. Note that the initial thickness profile defined by Eq. (37) is a function of the grounding line position, x_gl.

We define the initial grounding line position to be close to its steady-state position. According to boundary layer (Schoof, 2007a) and numerical convergence analyses employing SSA (Seroussi et al., 2014a), the grounding line should be located at x≃600 km. We therefore set x_gl=600 km. Note that this steady-state grounding line position only applies to SSA models; other stress balance models (full Stokes, Blatter–Pattyn, L1L2, etc.) produce steady-state positions upstream of 600 km (Pattyn et al., 2013). The ice thickness as defined by Eq. (37) is not an exact steady-state profile but represents a slightly perturbed profile used to initialize all the experiments. The numerical models rely on unstructured triangular meshes (see some examples in Appendix B), and the mesh resolutions (and respective number of elements) chosen for the numerical convergence experiments are shown in Table 3. The analytical thickness (Eq. 37) is interpolated onto each vertex of the mesh, and the floatation criterion is applied to generate both surface and base profiles of the ice sheet. To enforce the same initial position of the grounding line (x_gl) in the models, we employ a friction level set defined according to a distance-based function:

$$\begin{array}{}\text{(40)}& {\displaystyle }{\displaystyle }{\mathit{\phi }}_{\text{gr}}(x,y)=x_{\text{gl}}-x,\end{array}$$

where φ_gr is the initial friction level set (positive if ice is grounded, and therefore friction is applied). The level set is evaluated at all elements' vertices.

The numerical experiments are divided in two sets of analyses: (i) diagnostic analysis and (ii) prognostic analysis. The diagnostic analysis consists of solving the stress balance equations under different sets of sub-element parameterization schemes (driving stress and basal friction; Table 1), with different mesh resolutions (Table 3). Here, we compare the ice speeds calculated by each set of sub-element parameterization schemes and mesh resolutions. We employ the two-dimensional SSA (MacAyeal, 1989) to compute the velocity field. The SSA equations are solved using Picard iterations and each linear system is solved using an iterative linear solver (conjugate gradient). The ice-sheet geometry is given by Eq. (37), with grounding line position defined by Eq. (40). The aim of the prognostic analysis is to solve the mass transport equation and compare the transient response using different stabilization schemes (Table 1). The transient simulations start from the same initial condition (Eq. 37) and grounding line position (Eq. 40) and run forward in time for 100 years under the same accumulation rate (${\dot{m}}_{\mathrm{s}}$; Table 2) under three different scenarios of external forcings: (i) no external forcing, (ii) basal melt under the ice shelf, and (iii) changes in basal friction. The first experiment (no external forcing) aims to analyze the models' adjustments under no external perturbation. In the second experiment, basal melt is applied at the base of the ice shelf in order to thin it, generating a large change in ice thickness close to the grounding line. Here, the basal melt is applied on (a) only fully floating elements and (b) fully and partially floating elements. It is important to note that applying melt on partly floating elements is not suitable for realistic marine ice-sheet simulations (see discussion in Seroussi and Morlighem, 2018), and we only assess the implications of using different stabilization schemes. In (a), no grounding line retreat is expected, since no buttressing effect is present in the initial condition (unconfined ice shelf). In (b), the grounding line is expected to retreat because part of the melt rate is applied on the first grounded vertices (Seroussi and Morlighem, 2018). The basal melt is defined below (see Eq. 41). The last experiment (friction perturbation) is based on the MISMIP3d perturbation phase, as a non-symmetric change in the basal friction coefficient is introduced. The goal of this experiment is not to assess grounding line reversibility, as proposed in the original MISMIP3d benchmark, but instead to assess the general migration of the grounding line (and resulting mass loss) for each stabilization scheme and different mesh resolution.

The basal melt applied in the second set of experiments is defined as follows (Seroussi and Morlighem, 2018):

$$\begin{array}{}\text{(41)}& {\displaystyle }{\displaystyle }{\dot{m}}_{\mathrm{b}}\left(z_{\mathrm{b}}\right)=\left\{\begin{array}{ll}\mathrm{0},& \text{if }{z}_{\mathrm{b}}\ge z_{\mathrm{u}},\\ {\overline{m}}_{\mathrm{b}}{\displaystyle \frac{z_{\mathrm{b}}-z_{\mathrm{u}}}{z_{\mathrm{d}}-z_{\mathrm{u}}}},& \text{if }{z}_{\mathrm{d}}<z_{\mathrm{b}}<z_{\mathrm{u}},\\ {\overline{m}}_{\mathrm{b}},& \text{if }{z}_{\mathrm{b}}\le z_{\mathrm{d}},\end{array}\right.\end{array}$$

where ${\dot{m}}_{\mathrm{b}}$ is the depth-dependent basal melting in m yr⁻¹ (positive if melting), z_b, in meters, is the vertical coordinate of the ice-shelf base (negative if below sea level), ${\overline{m}}_{\mathrm{b}}$ is the maximum melt rate for vertical coordinates equal to or lower than z_d. We use ${\overline{m}}_{\mathrm{b}}=\mathrm{30}$ m yr⁻¹, $z_{\mathrm{u}}=-\mathrm{50}$ m, and $z_{\mathrm{d}}=-\mathrm{200}$ m, unless otherwise specified.

In all prognostic experiments, the grounding line is free to migrate, and its position over the simulation time is updated according to a hydrostatic floatation criterion, following Eq. (32): ice is floating (grounded) if the thickness, H, is lower (higher) than the floatation height, $H_f=-({\mathit{\rho }}_{\mathrm{w}}/\mathit{\rho })r$. The grounding line position is defined where H=H_f.

Table 2 Constants and parameters used along the numerical experiments.

Table 3 Mesh resolution and associated number of elements

4.2 Amundsen Sea Embayment – numerical setup

In order to quantify the performance of the stabilization schemes with real ice-sheet geometries and numerical setups, we run transient simulations (prognostic analyses) of the ASE, which includes the fastest glaciers of WAIS. The glaciers in the ASE are subject to high ocean-induced melt rates and are prone to the marine ice-sheet instability (MISI), a positive feedback of grounding line retreat and increased ice discharge sustained by a retrograde bedrock slope (Weertman, 1974; Schoof, 2007b; Gudmundsson et al., 2012). Ice-shelf buttressing is present for most glaciers in the ASE, which helps stabilize the grounding line on retrograde slopes, reducing the possibility of MISI (Dupont and Alley, 2005; Goldberg et al., 2009; Docquier et al., 2014; Gudmundsson et al., 2019; Martin et al., 2019). Therefore, ice-shelf thinning plays an important role in the dynamics of this sector, and it is important to assess the impact of different stabilization schemes on the response of ASE due to ocean-induced melt and resulting sea-level-rise contribution.

Our ASE domain includes Pine Island, Thwaites, and neighboring glaciers (Haynes, Pope, Smith, and Kohler glaciers). We use BedMachine Antarctica v1 (Morlighem et al., 2020) to build the digital elevation model and interferometric synthetic aperture radar (InSAR)-derived surface velocities (Mouginot et al., 2019a) to infer basal friction coefficient and ice rheology using control algorithms (Morlighem et al., 2010, 2013). We generate the mesh based on an interpolation error estimate of the observed ice velocity and on the distance to the grounding line: a mesh resolution equal to 1 km is employed in the vicinity of the grounding lines, while coarser resolution (up to 16 km) is employed for the rest of the grounded ice. The mesh contains 261 375 elements and 131 087 vertices. Note that this is a typical numerical setup employed in many of ice-sheet studies (e.g., Favier et al., 2014; Joughin et al., 2014; Seroussi et al., 2014b; Cornford et al., 2015). Details of the model setup and initialization are described in Barnes et al. (2020).

For this setup, we perform only transient simulations with different stabilization schemes. All simulations start from the same initial condition and are forced by a constant surface mass balance obtained from the regional climate model (RACMO v2.3, Van Wessem et al., 2014) and the parameterized basal melt defined by Eq. (41). Here, the parameters in Eq. (41) are based on ice–ocean simulations (Seroussi et al., 2017; Nakayama et al., 2019): z_u=0 and $z_{\mathrm{d}}=-\mathrm{500}$ m. We run two basal melt scenarios based on different values of ${\overline{m}}_{\mathrm{b}}$. For the first scenario, we set ${\overline{m}}_{\mathrm{b}}=\mathrm{50}$ m yr⁻¹ and, for the second one, ${\overline{m}}_{\mathrm{b}}=\mathrm{200}$ m yr⁻¹. Melt is applied only on fully floating elements. We chose the combination SEP1 + NSED as a set of parameterization schemes for basal friction and driving stress, respectively. The densities of ice and ocean are set to 917 and 1027 kg m⁻³, respectively. We run forward in time for 50 years with a fixed time step equal to one-eighth of a year. Like the other experiments performed here, the ice flow is computed by the SSA equations (MacAyeal, 1989).

5 Results

5.1 MISMIP3d – diagnostic analysis

To compare the ice speed from different sets of sub-element parameterizations, we compute the speed from a reference model. The reference model is based on a triangular structured conforming mesh with resolution of 50 m. This structured conforming mesh is constructed in such a way that the elements' edges (around the grounding line) match perfectly the grounding line position (x_gl=600 km; see Sect. 4.1; see also examples of structured meshes in Appendix B), and therefore no error related to driving stress and basal friction modeling is introduced in the stress balance computation in this reference model. We compare the speeds using two norms, the L₂ norm and L_∞ norm, defined as follows:

$$\begin{array}{}\text{(42)}& {\displaystyle }& {\displaystyle }\parallel v_i^h-v_i^r{\parallel }_{L_{\mathrm{2}}}^h=\sqrt{\sum _{i=\mathrm{1}}^N(v_i^h-v_i^r{)}^{\mathrm{2}}},\text{(43)}& {\displaystyle }& {\displaystyle }\parallel v_i^h-v_i^r{\parallel }_{L_{\mathrm{\infty }}}^h=\underset{i=\mathrm{1},N}{\text{max}}\left|v_i^h-v_i^r\right|,\end{array}$$

where $v=\parallel \overline{\mathit{v}}\parallel$ is the absolute value of the ice speed, h and r refer to the mesh resolution and reference model, respectively. To evaluate the norms, the results from all models are interpolated onto a finer regular grid (25 m resolution), and the speed differences ($v_i^h-v_i^r$) are calculated on each vertex i of that grid. For the error convergence comparison, we also run models with structured meshes (conforming to x_gl=600 km) considering the mesh resolutions shown in Table 3 (see examples in Appendix B). We refer to the speeds obtained by these conforming-mesh-based models as the reference speeds.

Figure 1 Ice speeds obtained by the diagnostic analysis along a flow line (y=25 km) and around the grounding line (x=600 km, dashed line) for different mesh resolutions: 5000, 2000, 1000, and 500 m. All sets of sub-element parameterizations are shown: SEP1 + NSED, SEP2 + NSED, SEP1 + SED2, and SEP2 + SED2. The speed from the 50 m resolution structured conforming mesh (reference model) is also shown (reference, dotted line).

Figure 2 Convergence of the ice speed at the grounding line (x=600 km, a) and at the ice front (x=800 km, b) obtained in the diagnostic analysis. All sets of sub-element parameterizations are shown: SEP1 + NSED, SEP2 + NSED, SEP1 + SED2, and SEP2 + SED2. Reference speeds from structured conforming-mesh-based models are also displayed (reference, dotted line).

Figure 3 Error convergence of the ice speed for the diagnostic analysis in L₂ norm (a) and L_∞ norm (b). All sets of sub-element parameterizations are shown: SEP1 + NSED, SEP2 + NSED, SEP1 + SED2, and SEP2 + SED2. The error convergence from structured conforming-mesh-based models is also shown (reference, dotted line).

Upstream of the grounding line, all sets of parameterizations “approach” the reference speed ($v_i^r$) for mesh resolutions finer than or equal to 1 km (Fig. 1). In the downstream part (floating ice), a 500 m mesh resolution or finer is required for all models to reach speeds at the grounding line and ice front similar to the reference model (Figs. 1 and 2). Overall, the errors in speed from models employing SEP1 + NSED and SEP2 + SED2 are closer to the ones obtained from structured conforming meshes (dotted line in Fig. 3). In other words, the convergence of sets SEP1 + NSED and SEP2 + SED2 is similar to the convergence of conforming-mesh-based models. The other combinations, SEP1 + SED2 and SEP2 + NSED, present relatively higher error levels in both norms (Fig. 3) for mesh resolutions coarser than 500 m.

5.2 MISMIP3d – prognostic analysis

We compare the transient results using the volume above floatation changes (ΔVAF) generated by each model and mesh resolution. The changes in VAF over time, t, are calculated as follows:

$$\begin{array}{}\text{(44)}& {\displaystyle }{\displaystyle }\mathrm{\Delta }{\text{VAF}}_t^h={\text{VAF}}_t^h-{\text{VAF}}_{t_{\mathrm{0}}}^h,\end{array}$$

where h refers to mesh resolution and t₀ is the initial time of the transient simulations. For some simulations, we also compare the grounding line positions at the end of the experiments. Based on the diagnostic analysis (Sect. 5.1), we employ the following sets of parameterizations in all transient simulations: SEP1 + NSED and SEP2 + SED2. We also compute the relative error in ΔVAF^h obtained at the end of the experiments for a mesh resolution h, defined as

$$\begin{array}{}\text{(45)}& {\displaystyle }{\displaystyle }{\mathit{\epsilon }}_{\mathrm{\Delta }\text{VAF}}^h=\left|{\displaystyle \frac{\mathrm{\Delta }{\text{VAF}}^{\mathrm{2}h}-\mathrm{\Delta }{\text{VAF}}^h}{\mathrm{\Delta }{\text{VAF}}^{\mathrm{2}h}}}\right|,\end{array}$$

where 2h indicates the next coarser mesh resolution. We use similar quantities to compute the relative error in grounding line position at the end of the experiments, replacing ΔVAF by the grounding line position, GL_pos, in Eq. (45).

5.2.1 No external forcing experiment

Since the initial ice-sheet profile (Eq. 37) is not exactly in steady state, some changes in VAF are expected to occur along the transient simulation due to grounding line adjustments (Fig. 4). In this control experiment, for mesh resolutions equal to or finer than 2 km, all stabilization schemes produce similar evolution of ΔVAF for both sets of parameterizations (Fig. 5). At the end of the experiment (t=100 years), all models produce a VAF loss equal to 116±4 Gt (Fig. 6). As shown in Fig. 6, the stabilization schemes generate different convergence curves of relative errors (except for streamline upwinding and artificial diffusion), but all show a decrease in error as the resolution increases, as expected. For a mesh resolution equal to 250 m, the relative errors for all stabilization schemes and sub-element parameterizations are below 5 % (Fig. 6).

Figure 4 No external forcing experiment: ice surface and ice base at the end of the experiment (t=100 years) for different stabilization schemes (see legend). Two sets of parameterizations are employed: SEP1 + NSED (a) and SEP2 + SED2 (b). The dotted black line is the initial ice-sheet geometry, defined by Eq. (37). Here, the mesh resolution is equal to 500 m.

Figure 5 No external forcing experiment: evolution of volume above floatation change (ΔVAF) for different mesh resolutions and stabilization schemes (see legend). Two sets of parameterizations are employed: SEP1 + NSED and SEP2 + SED2.

Figure 6 No external forcing experiment: convergence of ΔVAF at the end of the experiment (t=100 years) for different stabilization schemes (see legend). Two sets of parameterizations are employed: SEP1 + NSED (left panels) and SEP2 + SED2 (right panels). The relative errors are computed using Eq. (45).

5.2.2 Basal melt experiment

In the setup where basal melt is applied only to fully floating elements (i.e., no melt on partly floating elements), models using artificial diffusion and streamline upwinding schemes produce almost 4 times the VAF losses observed in the control experiment (Fig. 7). At the end of the experiment, t=100 years, and for a mesh resolution equal to 500 m, the expected ΔVAF is $-\mathrm{116}\pm \mathrm{4}$ Gt (see Sect. 5.2.1). A comparable amount of mass loss is obtained with models employing SUPG and discontinuous Galerkin. However, using artificial diffusion and streamline upwinding, the resulting ΔVAF is $\sim -\mathrm{360}$ Gt, while with FCT, the ΔVAF reaches $\sim -\mathrm{200}$ Gt (Fig. 8). The grounding line positions at the end of this basal melt experiment are expected to be close to the ones obtained with the control experiment (no external forcing). This is virtually achieved by models running with SUPG and discontinuous Galerkin, as illustrated on the left panel of Fig. 9. The grounding lines obtained from models employing artificial diffusion and streamline upwinding and even FCT have retreated further inland, resulting in an overestimated mass loss in comparison to SUPG and discontinuous Galerkin (Fig. 9). Both sets of sub-element parameterizations (SEP1 + NSED and SEP2 + SED2) lead to similar ΔVAF “convergence” for all stabilization schemes, although the convergence errors differ among them (Fig. 8). Only SUPG shows a decrease in relative error with mesh resolution, reaching error levels smaller than 10 % for a 250 m mesh resolution. Note that, in Fig. 8, artificial diffusion and streamline upwinding present relatively smaller errors in comparison to the others schemes (mainly DG and FCT), but the errors produced by these schemes do not decrease with mesh resolution and seem to be far from convergence, even for a mesh resolution of 250 m.

Figure 7 Basal melt experiment (no melt on partly floating elements): evolution of ΔVAF for different mesh resolutions and stabilization schemes (see legend). Two sets of sub-element parameterizations are employed: SEP1 + NSED and SEP2 + SED2.

Figure 8 Basal melt experiment (no melt on partly floating elements): convergence of ΔVAF at the end of the experiment (t=100 years) for different stabilization schemes (see legend). Two sets of sub-element parameterizations are employed: SEP1 + NSED (left panels) and SEP2 + SED2 (right panels). The relative errors are computed using Eq. (45).

Figure 9 Basal melt experiment: ice surface and ice base at the end of the experiment (t=100 years) for different stabilization schemes (see legend). (a) No melt on partly floating elements. (b) Melt on fully and partly floating elements. The dotted black line is the profile obtained with no basal melt applied (no external forcing experiment; Sect. 5.2.1 and Fig. 4). Here, the set of parameterizations is SEP1 + NSED, and the mesh resolution is equal to 500 m.

Figure 10 Basal melt experiment (melt on partly floating elements): evolution of ΔVAF for different mesh resolutions and stabilization schemes (see legend). Two sets of sub-element parameterizations are employed: SEP1 + NSED and SEP2 + SED2.

Figure 11 Basal melt experiment (melt on partly floating elements): convergence of ΔVAF at the end of the experiment (t=100 years) for different stabilization schemes (see legend). Two sets of sub-element parameterizations are employed: SEP1 + NSED (left panels) and SEP2 + SED2 (right panels). The relative errors are computed using Eq. (45).

When some basal melt is also applied to partly floating elements, all models generate VAF losses higher than those generated with the previous basal melt setup (Fig. 10), as expected (Seroussi and Morlighem, 2018). Models employing artificial diffusion and streamline upwinding schemes produce almost twice the change in VAF compared to models using SUPG and discontinuous Galerkin. In this setup, FCT tends to perform as well as SUPG and discontinuous Galerkin as the mesh resolution becomes finer (Figs. 10 and 11). The grounding lines obtained with artificial diffusion and streamline upwinding schemes evolve upstream of the grounding lines computed with the SUPG, discontinuous Galerkin, and FCT schemes (Fig. 9b). The sensitivity to mesh resolution is higher for models using artificial diffusion and streamline upwinding: VAF losses vary from 400 to 1300 Gt for the range of mesh resolutions used in this study, 250 to 5000 m (Fig. 11). For models employing SUPG, discontinuous Galerkin, and FCT, VAF losses vary from about 200 to 800 Gt for the same range of mesh resolutions. In this melt setup, a mesh resolution equal to 250 m is not enough to achieve relatively small errors in VAF changes: the relative errors for this mesh resolution vary from 20 % to 30 %, depending on the stabilization scheme used, and they do not seem to have entered the asymptotic region (Fig. 11). Both sets of sub-element parameterizations, SEP1 + NSED and SEP2 + SED2, generate virtually the same ΔVAF convergence for all stabilization schemes (see Fig. 11).

5.2.3 Friction perturbation experiment

Virtually all stabilization schemes produce the same ΔVAF evolution for both sets of sub-element parameterizations (Fig. 12). The amount of mass loss at the end of the experiment varies with mesh resolution: from 1300 to 1600 Gt for the range of mesh resolutions considered here, 5000 to 250 m. Using the SUPG, discontinuous Galerkin, and FCT schemes, the models produce slightly less mass loss than artificial diffusion and streamline upwinding, ∼40 Gt (Fig. 13). Both sets of parameterizations generate similar ΔVAF convergence for all stabilization schemes (Fig. 13). For mesh resolutions equal to or finer than 1 km, the relative errors decrease to values below 5 %, reaching ∼0.5 % for a resolution of 250 m. The positions of the grounding line along the lateral boundaries (i.e., y=0 and y=50 km) reach 618 and 600.3 km, respectively, for mesh resolutions finer than 1 km (Fig. 14). For these mesh resolutions, the relative errors are smaller than 0.1 % (Fig. 15).

Figure 12 Basal friction perturbation experiment: evolution of ΔVAF for different mesh resolutions and stabilization schemes (see legend). Two sets of sub-element parameterizations are employed: SEP1 + NSED and SEP2 + SED2.

Figure 13 Basal friction perturbation experiment: convergence of ΔVAF at the end of the experiment (t=100 years) for different stabilization schemes (see legend). Two sets of sub-element parameterizations are employed: SEP1 + NSED (left panels) and SEP2 + SED2 (right panels). The relative errors are computed using Eq. (45).

Figure 14 Basal friction perturbation experiment: convergence of grounding line positions at the end of the experiment (t=100 years) for y=0 (a) and y=50 km (b). Different stabilization schemes and two sets of sub-element parameterizations are employed (see legend).

Figure 15 Basal friction perturbation experiment: convergence of grounding line positions (relative errors) at the end of the experiment (t=100 years) for y=0 (top panels) and y=50 km (bottom panels) and for different stabilization schemes (see legend). Two sets of sub-element parameterizations are employed: SEP1 + NSED (left panels) and SEP2 + SED2 (right panels). The relative errors are computed using Eq. (45) and replacing ΔVAF by the grounding line position.

5.3 Amundsen Sea Embayment – prognostic analysis

To evaluate the performance of the stabilization schemes in real ice-sheet simulations (i.e., the ASE setup), we compare the VAF changes obtained with transient simulations employing the five schemes considered in this work. For the SUPG scheme, we chose the stability parameter as defined by Eq. (14). Models running with the other definition (i.e., Eq. 11) present spurious oscillations in this experiment. The two major glaciers in ASE, i.e., the Thwaites and Pine Island (PIG) glaciers, may respond differently to ocean-induced melt: Pine Island presents a more confined ice shelf compared to Thwaites. Therefore, we also compute the changes in VAF for these two glaciers.

In the experiment forced by the first basal melt scenario (i.e., ${\overline{m}}_{\mathrm{b}}=\mathrm{50}$ m yr⁻¹), the VAF loss after 50 years varies from 3200 to 2800 Gt for Thwaites and from 1900 to 1700 for Pine Island (Fig. 16, upper panels). In both glaciers, the artificial diffusion scheme overestimates the amount of VAF loss up to 10 % in comparison to the SUPG scheme. Discontinuous Galerkin and SUPG produce a smaller change in VAF, while streamline upwinding and FCT lies in between.

Under a higher basal melt scenario (${\overline{m}}_{\mathrm{b}}=\mathrm{200}$ m yr⁻¹), the VAF losses vary from 4400 to 3900 Gt and from 5100 to 3500 Gt for Thwaites and Pine Island, respectively (Fig. 16, lower panels). For Thwaites, the streamline upwinding scheme produces less mass loss during the entire simulation time, while artificial diffusion and FCT generate the highest amount of mass losses. The effect of the artificial damping is more pronounced in Pine Island: at the end of the experiment, the artificial diffusion scheme leads to ∼50 % more VAF loss compared to SUPG, the scheme that produces the lowest change in VAF for PIG. Discontinuous Galerkin generates spurious oscillations in ice thickness for this experimental setup, and therefore it is not shown in Fig. 16.

Considering the entire ASE domain, in simulations forced by a low melt rate (${\overline{m}}_{\mathrm{b}}=\mathrm{50}$ m yr⁻¹), the model running with artificial diffusion overestimates by 10 % the VAF loss in comparison to the one employing SUPG (Fig. 17). Under a higher rate of basal melt (${\overline{m}}_{\mathrm{b}}=\mathrm{200}$ m yr⁻¹), the VAF loss of ASE is overestimated by about 20 % in the same artificial diffusion–SUPG rate comparison (see Fig. 17). Streamline upwinding and FCT present similar responses for both melt scenarios: these schemes generate less VAF losses in comparison to artificial diffusion. This difference is more pronounced in the experiment forced by 200 m yr⁻¹ of melt rate.

Figure 16 ΔVAF along transient simulations for Thwaites (left panels) and Pine Island (PIG, right panels) glaciers. The transient simulations are forced by two different basal melt rate scenarios: 50 m yr⁻¹ (upper panels) and 200 m yr⁻¹ (lower panels). The basal melt is applied only on fully floating elements, and the parameterization schemes for basal friction and driving stress are SEP1 and NSED, respectively.

Figure 17 ΔVAF along transient simulations for the Amundsen Sea Embayment (ASE) domain. The transient simulations are forced by two different basal melt rate scenarios: 50 m yr⁻¹ (a) and 200 m yr⁻¹ (b). The basal melt is applied only on fully floating elements, and the parameterization schemes for basal friction and driving stress are SEP1 and NSED, respectively.

6 Discussion

The diagnostic analysis using the analytical ice-sheet profile (Sect. 5.1) shows that the convergence of ice speeds depends on the set of parameterizations chosen. The convergence curve associated with the new driving stress parameterization (SED2) is similar to the ones of structured conforming meshes when combined with the parameterization SEP2 for basal friction. Specifically, considering an a priori error estimate of kh^β, where k and β are constants, k is independent of h, β depends on the polynomial order of elements and smoothness of the solution (Szabó and Babuška, 1991, p. 61, 62), employing SEP2 + SED2 reduces the value of k to a value closer to the ones produced by structured conforming meshes (see SEP2 + SED2 and SEP2 + NSED in Fig. 3). The same is observed for SEP1 + NSED in comparison to SEP1 + SED2. Overall, models employing SEP1 + NSED and SEP2 + SED2 achieve relatively low levels of errors with mesh resolutions at least 4 times coarser in comparison to the other schemes, SEP1 + SED2 and SEP2 + NSED (see Fig. 3). These results suggest that improved convergence (i.e., smaller k) is achieved by discretizing consistently the friction coefficient and driving stress, i.e., by employing SEP1 + NSED or SEP2 + SED2. This is consistent with previous studies based on finite-volume and finite-difference methods (Cornford et al., 2013; Feldmann et al., 2014), where one-sided differences were employed on each side of the grounding line to compute surface gradients and basal friction. The driving stress and basal friction should be equally balanced (i.e., discretized) at the grounding line to improve the accuracy of the ice velocity as well as the dynamics of the grounding line (Cornford et al., 2013; Feldmann et al., 2014). The error norms for fine mesh resolutions (500 m or finer) shown in Fig. 3 are probably impacted by our iterative solver. However, note that the errors for these mesh resolutions are below 1 m yr⁻¹ in the L_∞ norm, or $\mathrm{3.17}\times {\mathrm{10}}^{-\mathrm{8}}$ m s⁻¹, which is the same order of magnitude as our solver tolerance and machine precision.

Table 4 Steady-state GL positions for the MISMIP3d setup using different sub-element parameterization schemes. GL positions for SEP1 + NSED and SEP2 + NSED are extracted from Seroussi et al. (2014a). We employ the same numerical setup as described in Seroussi et al. (2014a). The estimated GL position from BLT is also included (Schoof, 2007b).

Figure 18 Convergence of steady-state grounding line (GL) positions (y=0 km) for the MISMIP3d setup. Different sets of sub-element parameterizations are tested: SEP1 + NSED, SEP2 + NSED, SEP1 + SED2, and SEP2 + SED2. The dotted line is the grounding position from the boundary layer theory (BLT) (Schoof, 2007b).

Employing an analytical expression of ice geometry based on a predefined grounding line position allows the setup of reference models (i.e., models whose mesh captures the exact position of the grounding line), in which no errors due to parameterization schemes are introduced during the stress balance solution (diagnostic analysis). Therefore, using the reference setup improves the confidence of this analysis. Comparing grounding line positions at steady state is another approach (Table 4 and Fig. 18) where the boundary layer theory (BLT) provides an estimated position of the steady grounding line. The steady-state comparison shows that 1 km mesh resolution is enough for models using SEP1 + NSED and SEP2 + SED2 to achieve the grounding line position predicted by the BLT, within a tolerance of 0.5 %, while models employing other schemes (SEP1 + SED2 and SEP2 + NSED) need finer mesh resolution (at least 16 times more elements to generate a 250 m resolution mesh) to achieve the BLT prediction with the same tolerance. This corroborates the conclusions obtained with the diagnostic analysis in Sect. 5.1.

The prognostic analysis performed with the MISMIP3d-type geometry shows that the numerical damping produced by the artificial diffusion and streamline upwinding schemes impacts the accuracy of grounding line dynamics mainly in simulations when large ΔH appear at its vicinity, such as the ice-shelf melt experiments. These two schemes generate the same stabilization term in the x-direction flow performed in this numerical setup (see Sect. 2.2). The numerical damping induces a positive feedback of mass loss: the grounded ice upstream of the grounding line thins due to the (artificial) damping and starts to unground; once it is floating, it is subject to basal melt, which induces further thinning inland. For melt-induced experiments, where some melt is also applied to partly floating elements, models running with artificial diffusion and streamline upwinding overestimate mass loss in comparison to models running with other schemes (SUPG, DG, and FCT). Note, however, that applying melt on partly floating elements leads to high levels of errors for all the stabilization schemes and mesh resolutions used in this study and therefore is not a suitable approach. The basal friction perturbation experiment (Sect. 5.2.3) shows that, for a given mesh resolution, the relative errors in ΔVAF are at least 1 order of magnitude higher than the errors in grounding line position (see Figs. 13 and 15). For instance, the grounding line position has an error of 1 % for a 2 km mesh, while the corresponding error in ΔVAF is higher than 10 %. Small changes in grounding line position can therefore lead to large differences in ΔVAF. This suggests that ΔVAF should be used as a metric of numerical convergence in ice-sheet model intercomparison projects. The MISMIP3d setup used in this analysis is suitable because (1) there is no buttressing effect involved in these basal melt experiments, and therefore the grounding line dynamics is expected to be independent of ice-shelf basal melt; (2) the bedrock is the same for all mesh resolutions, which eliminates the source of errors related to bedrock resolution (Durand et al., 2011); (3) it allows the definition of an analytical ice-sheet profile and guarantees the same initial condition for all models independently of the mesh resolution or stabilization schemes. These numerical characteristics therefore eliminate the influence of several sources of errors, allowing the analysis to focus only on the response to the stabilization schemes.

For the prognostic analysis performed with real glaciers in West Antarctica (Sect. 5.3), streamline upwinding performs as well as the FCT scheme, which may be explained by the “anisotropic balancing dissipation” of the stabilization term (Kelly et al., 1980; Brooks and Hughes, 1982) (Sect. 2.2) that has prevented numerical damping over transverse flows. Interestingly, streamline upwinding generates less mass loss in comparison to SUPG for the Thwaites Glacier in the high-melt-rate scenario (Fig. 16). We attribute this to the low performance of the stabilization parameter (τ) in this simulation. In this same melt scenario, the mass loss of PIG is clearly overestimated using artificial diffusion, which is likely associated with the positive feedback explained above: the grounding line of PIG retreats several kilometers more using ArtDiff compared to the other schemes. For the entire ASE domain considered here, SUPG, FCT, and streamline upwinding yield similar VAF evolution within a difference of ∼5 % in the high-melt experiment. Discontinuous Galerkin has similar performance to SUPG, at least for the lower-basal-melt scenario. However, producing spurious oscillations in ice thickness in the second melt experiment highlights that the DG scheme as implemented in ISSM may be not completely robust.

The choice of the stabilization scheme relies on the balance between stability, accuracy, computational cost, and implementation effort. Yet, the “best” choice could be never reached (e.g., John et al., 2018). Artificial diffusion provides stability to numerical solutions, even in the presence of strong discontinuities and shocks (which is not the case in ice-sheet dynamics). But at the same time, it excessively smooths the solution for large ΔH, which impacts the results. For simulations with large buttressing effect, as is the case for the Pine Island and Thwaites glaciers, numerical damping could artificially “enhance” the marine ice-sheet instability feedback existing for retrograde bedrock slopes and, consequently, overestimate mass loss in ocean-induced melt simulations. Note that numerical damping does not always lead to grounding line retreat, but it can also prevent its advance. High mesh resolution could be employed to decrease the diffusion effects, but it comes at a higher computational cost. Adaptive mesh refinement could be an alternative in this case, although it is not available in all ice-sheet models. Discontinuous Galerkin and SUPG may generate spurious oscillations in idealized experiments where discontinuities are present (e.g., John and Schmeyer, 2008), and in real ice-sheet simulations, as observed here for ASE experiments. The performance and stability of SUPG clearly rely on the definition of the stabilization parameter τ, and the numerical issues and results observed here for the ASE simulations indicate that the definitions of τ as given by Eqs. (11) and (14) are not totally robust and optimized for real ice-sheet simulations, at least for the SUPG as currently implemented in ISSM. FCT presents better results on idealized cases (John and Schmeyer, 2008), but some excessive VAF loss is observed in some experiments performed here (e.g., Fig. 7). Besides that, this scheme performs as well as the streamline upwinding in the two basal melt scenarios using the ASE setup.

Apart from differences observed in terms of accuracy (i.e., VAF change), the remaining differences between the stabilization schemes used here are their numerical implementations and computational costs. The implementation of the artificial diffusion, streamline upwinding, and SUPG is straightforward in most of ice-sheet FEM-based models. However, the definition of the stability coefficient for SUPG (Eq. 7) is problem dependent, and possibly, its optimum value may remain unclear in many real ice-sheet simulations. Discontinuous Galerkin requires specific coding of data structures, at the minimum requiring information on elements' neighbors, and significant implementation effort to compute the integrals along elements' edges. Also, the number of degrees of freedom in comparison to other schemes is considerably increased (up to a factor of 6 for triangular P1 Lagrange elements), which impacts the computational cost when the backward Euler approach is used, as is the case in ISSM. An alternative would be an explicit approach (Runge–Kutta discontinuous Galerkin): in this case, the solution using the DG scheme would be completely parallel (e.g., Cockburn, 2003). While this would improve the computational cost, there are more stringent restrictions on the CFL time step. The FCT scheme requires operations on global matrices and vectors. While this is straightforward in codes relying on shared memory (e.g., one process, multiple threads), these operations require additional CPU communications in codes based on distributed memory (e.g., Message Passing Interface), potentially translating into a larger computation time. Finally, we note that the current ISSM implementations of the stabilization schemes presented here are based on the classical literature of FEM, where the numerical analyses of such schemes are carried out in idealized problems. There is still room for development of stabilization schemes and improved numerical accuracy, stability, and computational performance in the specific field of ice-sheet modeling; most stabilization schemes were designed to handle the presence of shocks and strong discontinuities, which are mostly absent from this field.

7 Final remarks

The convergence error of ice speed depends on the combination of parameterizations chosen for basal friction and driving stress. Given that the a priori error estimate is kh^β (h is the mesh resolution), the sub-element parameterization for driving stress proposed here (SED2) presents smaller values of k when combined with a similar approach for basal friction (SEP2). In models employing the SEP1 basal friction parameterization, a smaller k is achieved with no driving stress parameterization (NSED). As already suggested by previous simulations based on finite-volume or finite-difference methods, in order to achieve improved numerical convergence (i.e., smaller values of k) for a given computational cost, the discretization of the basal friction should match the discretization of the driving stress; the following combinations of sub-element parameterizations therefore provide the best results: SEP1 + NSED and SEP2 + SED2. In all transient simulations performed here for both idealized and real ice-sheet configurations, models relying on artificial diffusion generate the largest amount of mass loss in comparison to models running with the other schemes. SUPG and discontinuous Galerkin produce the expected grounding line dynamics for the idealized case (e.g., MISMIP3d-type setup). For the ASE experiments, these two schemes produce less mass loss than the others ones, although some spurious oscillations are observed, which compromised the results. By design, the streamline upwinding has the same behavior as artificial diffusion in the idealized case (x-direction flow). However, in the real ice-sheet simulations (ASE setup), streamline upwinding performs as well as the flux-corrected transport scheme. Based on the numerical tests performed here and the ease of implementation, SUPG seems a preferred scheme, although careful attention shall be given to the definition of the stabilization parameter, which may be problem dependent. A second choice would be the streamline upwinding scheme, as long as a high-enough mesh resolution is employed around the discontinuities (e.g., grounding line). The development of new stabilization schemes and/or improvements of existing ones in FEM remains an active field of research. Nevertheless, since most theoretical studies and convergence analyses involve, in general, smooth data and regular boundaries, the conclusions drawn from these studies may not necessarily hold for real cases, such as the Amundsen Sea Embayment simulation performed here. This highlights the importance of testing future stabilization schemes with real geometries and external forcing. In general, most stabilization schemes were developed in the context of compressible flow, where shocks and strong discontinuities appear, which is not the case of ice-sheet modeling, opening new development opportunities for this specific field.

Appendix A: Zalesak's algorithm

We compute the weights α_ij in Eq. (23) using the anti-diffusive flux r_ij (Eq. 28) and Zalesak's algorithm (Zalesak, 1979) (see also Möller et al., 2004; John and Schmeyer, 2008, and Kuzmin, 2009):

⁷

Appendix B: Example of meshes

Six examples of structured and unstructured meshes are shown in Fig. B1. We use the package Triangle (Shewchuk, 1996) to generate the unstructured meshes.

Figure B1 Examples of meshes employed in this work. Left panels are structured conforming meshes, and right panels the unstructured meshes. Three mesh resolutions are shown: 2, 1, and 0.5 km. The color maps are the ice speeds obtained in the diagnostic analysis (Sect. 5.1) considering the grounding line at x=600 km. Note that the structured meshes are conforming to the grounding line; i.e., they are generated such that the elements' edges match the grounding line position (in this case, x_gl=600 km). Sub-element parameterizations SEP1 + NSED are employed in the unstructured meshes to computed the ice speeds presented here.