Articles | Volume 14, issue 5

https://doi.org/10.5194/gmd-14-2545-2021

Articles | Volume 14, issue 5

Development and technical paper

07 May 2021

Development and technical paper |

| 07 May 2021

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

Thiago Dias dos Santos, Mathieu Morlighem, and Hélène Seroussi

Abstract

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.

Download & links

Article (PDF, 9141 KB)

Supplement (8766 KB)

Download & links

How to cite.

Received: 25 Nov 2020 – Discussion started: 17 Dec 2020 – Revised: 12 Mar 2021 – Accepted: 28 Mar 2021 – Published: 07 May 2021

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, 1984 b). 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., 2014 a; 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., 2014 a; 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{matrix} (1) & \frac{\partial H}{\partial t} + \nabla \cdot (\bar{v} H) = {\dot{m}}_{s} - {\dot{m}}_{b}, \end{matrix}

where $\bar{v} = {v_{x}, v_{y}}$ is the depth-averaged ice velocity in the horizontal plane, ${\dot{m}}_{s}$ the surface mass balance (positive for accumulation), and ${\dot{m}}_{b}$ the basal melt (positive for ablation). The velocity field $\bar{v}$ is a function of the ice geometry and therefore of the ice thickness H. Note that both surface mass balance ${\dot{m}}_{s}$ and basal melt ${\dot{m}}_{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}}_{s}$ as a constant in all transient simulations. The description of ${\dot{m}}_{b}$ is given in Sect. 4.1 (see Eq. 41).

The weak formulation of Eq. (1) is

\begin{array}{l} \int_{Ω} (\frac{\partial H}{\partial t} + \nabla \cdot (\bar{v} H)) ψ d Ω = \int_{Ω} ({\dot{m}}_{s} - {\dot{m}}_{b}) ψ d Ω, \\ (2) & \forall ψ \in 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{matrix} (3) & \frac{\partial H}{\partial t} + \nabla \cdot (\bar{v} H) - \nabla \cdot (D \nabla H) = {\dot{m}}_{s} - {\dot{m}}_{b}, \end{matrix}

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{matrix} (4) & D = \frac{h}{2} (\begin{array}{cc} |v_{x}| & 0 \\ 0 & |v_{y}| \end{array}), \end{matrix}

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{matrix} (5) & D = \frac{h}{2 ∥ \bar{v} ∥} \bar{v} \otimes \bar{v} = \frac{h}{2 ∥ \bar{v} ∥} (\begin{array}{cc} v_{x}^{2} & v_{x} v_{y} \\ v_{y} v_{x} & v_{y}^{2} \end{array}), \end{matrix}

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

\begin{array}{l} \int_{Ω} (\frac{\partial H}{\partial t} + \nabla \cdot (\bar{v} H)) ψ d Ω + \int_{Ω} \nabla ψ \cdot D \nabla H d Ω \\ = \int_{Ω} ({\dot{m}}_{s} - {\dot{m}}_{b}) ψ d Ω, \\ (6) & \forall ψ \in 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 $\hat{ψ}$ are generally defined as (Brooks and Hughes, 1982)

\begin{matrix} (7) & \hat{ψ} = ψ + τ \nabla \cdot (\bar{v} ψ), \end{matrix}

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

\begin{array}{l} \int_{Ω} (\frac{\partial H}{\partial t} + \nabla \cdot (\bar{v} H)) ψ d Ω \\ + \int_{Ω} (\frac{\partial H}{\partial t} + \nabla \cdot (\bar{v} H)) τ \nabla \cdot (\bar{v} ψ) d Ω \\ = \int_{Ω} ({\dot{m}}_{s} - {\dot{m}}_{b}) ψ d Ω + \int_{Ω} ({\dot{m}}_{s} - {\dot{m}}_{b}) τ \nabla \cdot (\bar{v} ψ) d Ω, \\ (8) & \forall ψ \in H . \end{array}

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

\begin{matrix} (9) & τ = \frac{h}{2 ∥ \bar{v} ∥} (\coth (P_{e}) - \frac{1}{P_{e}}), \end{matrix}

with

\begin{matrix} (10) & P_{e} = \frac{h ∥ \bar{v} ∥}{2 κ}, \end{matrix}

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{matrix} (11) & τ \sim τ_{1} = \frac{h}{2 ∥ \bar{v} ∥} . \end{matrix}

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{matrix} (12) & τ = \frac{h}{2 ∥ \bar{v} ∥} (\coth (C_{e}) - \frac{1}{C_{e}}), \end{matrix}

with

\begin{matrix} (13) & C_{e} = \frac{∥ \bar{v} ∥ Δ t}{h}, \end{matrix}

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, 1984 a; 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} / 3$ (Gudmundsson, 2020, p. 20). Therefore, for low CFL numbers, we can approximate Eq. (12) by

\begin{matrix} (14) & τ \sim τ_{2} = \frac{Δ t}{6} . \end{matrix}

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}{l} \int_{Ω^{e}} (\frac{\partial H}{\partial t} ψ - (\bar{v} H) \cdot \nabla ψ) d Ω^{e} + \int_{S^{e}} (\bar{v} H ψ) \cdot n d S^{e} \\ = \int_{Ω^{e}} ({\dot{m}}_{s} - {\dot{m}}_{b}) ψ d Ω^{e}, \\ (15) & \forall ψ \in 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 ( $\bar{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}{l} \int_{Ω^{e}} (\frac{\partial H}{\partial t} ψ - (\bar{v} H) \cdot \nabla ψ) d Ω^{e} + \int_{S^{e}} {(\bar{v} H)}_{u} \cdot (ψ n) d S^{e} \\ = \int_{Ω^{e}} ({\dot{m}}_{s} - {\dot{m}}_{b}) ψ d Ω^{e}, \\ (16) & \forall ψ \in H, \end{array}

where $(\bar{v} H)_{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 $CFL < 1 / 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{matrix} (17) & (M_{C} + θ Δ t K) H^{n + 1} = [M_{C} - (1 - θ) Δ t K] H^{n} + Δ t F, \end{matrix}

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 $θ = 1 / 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}{l} M_{C} = diag \{m_{i}\}, \\ (18) & m_{i} = \sum_{j = 1}^{N} m_{i j}, \end{array}

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

\begin{matrix} (19) & L = K + D, \end{matrix}

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

\begin{array}{l} d_{i j} = d_{j i} = - max \{0, k_{i j}, k_{j i}\}, \forall i \neq j, \\ (20) & d_{i i} = - \sum_{j = 1, j \neq i}^{N} d_{i j}, \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{matrix} (21) & (M_{L} + θ Δ t L) H^{n + 1} = [M_{L} - (1 - θ) Δ t L] H^{n} + Δ t F . \end{matrix}

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{matrix} (22) & (M_{L} + θ Δ t L) H^{n + 1} = [M_{L} - (1 - θ) L] H^{n} + Δ t F + F^{*}, \end{matrix}

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

\begin{matrix} (23) & f_{i}^{*} = \sum_{j = 1}^{N} α_{i j} r_{i j}, i = 1, \dots, N . \end{matrix}

In Eq. (23), $α_{i j} \in [0, 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}{l} R = & (M_{L} - M_{C}) (H^{n + 1} - H^{n}) \\ (24) & + Δ t D [θ H^{n + 1} + (1 - θ) H^{n}] . \end{array}

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

\begin{matrix} (25) & r_{i} = \sum_{j = 1, j \neq i}^{N} r_{i j}, r_{j i} = - r_{i j}, \end{matrix}

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}{l} r_{i j} = & m_{i j} (H_{i}^{n + 1} - H_{j}^{n + 1}) - m_{i j} (H_{i}^{n} - H_{j}^{n}) \\ - Δ t θ d_{i j} (H_{i}^{n + 1} - H_{j}^{n + 1}) \\ (26) & - Δ t (1 - θ) d_{i j} (H_{i}^{n} - H_{j}^{n}) . \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{matrix} (27) & (M_{L} + θ Δ t L) H_{L} = [M_{L} - (1 - θ) Δ t L] H^{n} + Δ t F . \end{matrix}

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

\begin{matrix} (28) & r_{i j} = m_{i j} ({\dot{H}}_{L, i} - {\dot{H}}_{L, j}) + d_{i j} (H_{L, i} - H_{L, j}), \end{matrix}

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

\begin{matrix} (29) & {\dot{H}}_{L}^{m + 1} = {\dot{H}}_{L}^{m} + M_{L}^{- 1} (K H_{L} - M_{C} {\dot{H}}_{L}^{m}), m = 0, 1, \dots, \end{matrix}

with ${\dot{H}}_{L}^{0} = 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{matrix} (30) & M_{L} H^{n + 1} = M_{L} H_{L} + Δ t F^{*}, \end{matrix}

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}{l} \int_{Ω} (\frac{H^{n + 1} - H^{n}}{Δ t} + \nabla \cdot ({\bar{v}}^{n} H^{n + 1})) ψ d Ω \\ + \int_{Ω} \nabla ψ \cdot D^{n} \nabla H^{n + 1} d Ω \\ = \int_{Ω} ({\dot{m}}_{s}^{n} - {\dot{m}}_{b}^{n}) ψ d Ω, \\ (31) & \forall ψ \in 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., 2014 a). This level set can be located anywhere in an element, so it does not necessarily coincide with elements' edges:

\begin{matrix} (32) & φ = H + \frac{ρ_{w}}{ρ} r, \end{matrix}

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{matrix} (33) & s (x, y) = \{\begin{cases} H (x, y) + r (x, y), & if φ (x, y) > 0, \\ H (x, y) (1 - ρ / ρ_{w}), & if φ (x, y) \leq 0 . \end{cases} \end{matrix}

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: $Ω_{G}^{e}$ and $Ω_{F}^{e}$ that are the grounded and floating parts of the element, respectively (i.e., $Ω^{e} = Ω_{G}^{e} \cup Ω_{F}^{e}$ and $Ω_{G}^{e} \cap Ω_{F}^{e} = \emptyset$ ). We then perform the numerical integration of the driving stress on these two subdomains, i.e.,

\begin{matrix} (34) & F^{e} = \int_{Ω^{e}} A^{e} \nabla s d Ω^{e} = \int_{Ω_{G}^{e}} A^{e} \nabla s d Ω_{G}^{e} + \int_{Ω_{F}^{e}} A^{e} \nabla s d Ω_{F}^{e}, \end{matrix}

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. (2014 a). 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 (1 - ρ / ρ_{w})$ 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., 2014 a)(Seroussi et al., 2014 a)

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

Download Print Version | Download XLSX

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{matrix} (35) & τ_{b} = C ∥ v_{b} ∥^{m - 1} v_{b}, \end{matrix}

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{matrix} (36) & r (x, y) = 100 - x, \end{matrix}

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}{l} H (x, y, t = 0) \\ (37) & = \{\begin{cases} {[H_{gl}^{m + 2} + \frac{m + 2}{m + 1} \frac{C {\dot{m}}_{s}^{m}}{ρ g} (x_{gl}^{m + 1} - x^{m + 1})]}^{\frac{1}{m + 2}}, \\ if x < x_{gl}, \\ {[\frac{A}{{\dot{m}}_{s}} - \frac{v_{gl}^{n + 1} (\frac{A}{{\dot{m}}_{s}} H_{gl}^{n + 1} - 1)}{{[{\dot{m}}_{s} (x - x_{gl}) + v_{gl} H_{gl}]}^{n + 1}}]}^{- \frac{1}{n + 1}}, \\ if x \geq x_{gl}, \end{cases} \end{array}

where H is the ice thickness, ${\dot{m}}_{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}{l} H_{gl} = \frac{ρ_{w}}{ρ} |r (x_{gl}, y)|, \\ (38) & v_{gl} = \frac{{\dot{m}}_{s} x_{gl}}{H_{gl}}, \end{array}

and 𝒜 is defined by

\begin{matrix} (39) & A = A {(\frac{ρ g (ρ_{w} - ρ)}{4 ρ_{w}})}^{n}, \end{matrix}

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, 2007 a). 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, 2007 a) and numerical convergence analyses employing SSA (Seroussi et al., 2014 a), 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{matrix} (40) & φ_{gr} (x, y) = x_{gl} - x, \end{matrix}

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}}_{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{matrix} (41) & {\dot{m}}_{b} (z_{b}) = \{\begin{cases} 0, & if z_{b} \geq z_{u}, \\ {\bar{m}}_{b} \frac{z_{b} - z_{u}}{z_{d} - z_{u}}, & if z_{d} < z_{b} < z_{u}, \\ {\bar{m}}_{b}, & if z_{b} \leq z_{d}, \end{cases} \end{matrix}

where ${\dot{m}}_{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), ${\bar{m}}_{b}$ is the maximum melt rate for vertical coordinates equal to or lower than z_d. We use ${\bar{m}}_{b} = 30$ m yr⁻¹, $z_{u} = - 50$ m, and $z_{d} = - 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} = - (ρ_{w} / ρ) r$ . The grounding line position is defined where H=H_f.

Table 2Constants and parameters used along the numerical experiments.

Download Print Version | Download XLSX

Table 3Mesh resolution and associated number of elements

Download Print Version | Download XLSX

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, 2007 b; 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., 2019 a) 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., 2014 b; 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_{d} = - 500$ m. We run two basal melt scenarios based on different values of ${\bar{m}}_{b}$ . For the first scenario, we set ${\bar{m}}_{b} = 50$ m yr⁻¹ and, for the second one, ${\bar{m}}_{b} = 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}{l} (42) & ∥ v_{i}^{h} - v_{i}^{r} ∥_{L_{2}}^{h} = \sqrt{\sum_{i = 1}^{N} (v_{i}^{h} - v_{i}^{r})^{2}}, \\ (43) & ∥ v_{i}^{h} - v_{i}^{r} ∥_{L_{\infty}}^{h} = max_{i = 1, N} |v_{i}^{h} - v_{i}^{r}|, \end{array}

where $v = ∥ \bar{v} ∥$ 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.

https://gmd.copernicus.org/articles/14/2545/2021/gmd-14-2545-2021-f01

Figure 1Ice 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).

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

2.1 Mass transport equation

2.2 Artificial diffusion and streamline upwinding

2.3 Streamline upwind Petrov–Galerkin

2.4 Discontinuous Galerkin

2.5 Flux-corrected transport

2.6 Time discretization of the mass transport equation

4.1 MISMIP3d – numerical setup

4.2 Amundsen Sea Embayment – numerical setup

5.1 MISMIP3d – diagnostic analysis

5.2 MISMIP3d – prognostic analysis

5.2.1 No external forcing experiment

5.2.2 Basal melt experiment

5.2.3 Friction perturbation experiment

5.3 Amundsen Sea Embayment – prognostic analysis