Increasing numerical stability of mountain valley glacier simulations: implementation and testing of free-surface stabilization in Elmer/Ice

This paper concerns a numerical stabilization method for free-surface ice flow called the free-surface stabilization algorithm (FSSA). In the current study, the FSSA is implemented into the numerical ice-flow software Elmer/Ice and tested on synthetic two-dimensional (2D) glaciers, as well as on the real-world glacier of Midtre Lovénbreen, Svalbard. For the synthetic 2D cases it is found that the FSSA method increases the largest stable time-step size at least by a factor of 5 for the case of a gently sloping ice surface (∼ 3°) and by at least a factor of 2 for cases of moderately to steeply inclined surfaces (∼ 6° to 12°) on a fine mesh. Compared with other means of stabilization, the FSSA is the only one in this study that increases largest stable time-step sizes when used alone. Furthermore, the FSSA method increases the overall accuracy for all surface slopes. The largest stable time-step size is found to be smallest for the case of a low sloping surface, despite having overall smaller velocities. For an Arctic-type glacier, Midtre Lovénbreen, the FSSA method doubles the largest stable time-step size; however, the accuracy is in this case slightly lowered in the deeper parts of the glacier, while it increases near edges. The implication is that the non-FSSA method might be more accurate at predicting glacier thinning, while the FSSA method is more suitable for predicting future glacier extent. A possible application of the larger time-step sizes allowed for by the FSSA is for spin-up simulations, where relatively fast-changing climate data can be incorporated on short timescales, while the slow-changing velocity field is updated over larger timescales.

1 Introduction

Ice sheets and glaciers are important constituents of the global climate system, and the mass loss from these is expected to be a main contributor to future sea-level rise (DeConto and Pollard, 2016; Hock et al., 2019; Meredith et al., 2019; Fox-Kemper et al., 2021). In order to reliably estimate future sea-level rise the accurate representation of ice-sheet and glacier dynamics is crucial, and higher-order physics models have proven to be instrumental in increasing the confidence in predictions (Hanna et al., 2013; Shepherd and Nowicki, 2017; Pattyn, 2018).

The most accurate description for the flow of ice, in the sense that all stress components are present in the Cauchy stress tensor, is the Stokes equations (Greve and Blatter, 2009). Approximations of the model are made by neglecting various components, with some of the most notable examples being the shallow-ice approximation (SIA) (Hutter, 1983; Morland, 1984) and the first-order Stokes approximation (FOS) (Blatter, 1995; Pattyn, 2003). Owing to its simplicity and computational efficiency, the SIA method has a long history of use (see, e.g., Blatter et al., 2010, for a historical overview); however, the SIA method has been found to be insufficient at reproducing the flow at regimes with a steep sloping bedrock (Meur et al., 2004; Dukowicz et al., 2011; Leng et al., 2012), as well as for smaller glaciers with a complex bedrock topography (Zwinger et al., 2007).

It has been demonstrated for lower-order physics models of shear-dominated flow (such as the SIA) that when coupled to the free-surface equation governing the evolution of ice sheets and glaciers, they are subject to a parabolic-type time-step size constraint that is highly dependent on the ice-domain thickness (Bueler et al., 2005; Gong et al., 2017; Bueler, 2022; Robinson et al., 2022). However, for the Stokes equations the same type of time-step size restriction does not necessarily hold true – even for setups where the SIA and the Stokes equations give qualitatively similar solutions (Löfgren et al., 2022). Still, for ice-sheet simulations using the Stokes equations, the restriction on the time-step size to have numerical stability, herein broadly defined as the unbounded growth of numerical errors, is typically found to be on the order of 0.1 to 10 years (Gong et al., 2017; Löfgren et al., 2022). This is considerably smaller than typical timescales at which ice sheets evolve (Hindmarsh and Payne, 1996), which can be as large as 10 000 years (Greve and Blatter, 2009). Computation times can thus be cut if time-step sizes can be increased beyond the largest stable time-step size (LST) without compromising the desired accuracy of the solution.

One way of stabilizing the problem is to use a fully implicit time-stepping scheme, which has been demonstrated by Bueler (2016) for the SIA in a frozen-bed setting. However, the Stokes equations are considerably more expensive to solve than the simpler SIA equations and, since the nonlinear Stokes equations would have to be solved multiple times in each time step, make such a scheme computationally infeasible for long-term simulations. Instead Kaus et al. (2010) propose the free-surface stabilization algorithm (FSSA), which modifies the weak formulation of the Stokes equations in order for the free-surface coupled system to mimic an implicit time-stepping scheme. The method was originally developed for mantle-convection simulations where a similar viscous-flow problem is solved, and multiple studies have indeed demonstrated that the method lengthens the LST substantially (Kaus et al., 2010; Duretz et al., 2011; Kramer et al., 2012; Andrés-Martínez et al., 2015; Rose et al., 2017).

From a glaciological perspective, a limitation of the original FSSA method is that linear rheologies are used on domains that are geometrically isotropic, meaning that they span equally in horizontal and vertical directions; i.e., domain aspect ratios are 1:1. A notable exception is Glerum et al. (2020), a study which considers both a similar shear-thinning nonlinear rheology and domains with aspect ratios on the order of 1:10. Still, the values of the physical parameters describing ice flow are different from those in mantle convection, and aspect ratios of ice sheets can be as small as 1:1000.

These issues were addressed by Löfgren et al. (2022), where the FSSA method was adapted to ice-flow modeling. It was concluded that the method works well in an ice-dynamical setting, and for the problems presented it showed the potential to increase the LST by an order of magnitude. Nevertheless, one of the shortcomings in this case is that the method was only applied to simple ice-sheet benchmark problems, and more complex glacier simulations, e.g., using variable bedrock topography and sliding conditions, were only touched on briefly in that paper's “Supplementary material” and not studied thoroughly.

This work focuses on addressing these issues and applying the FSSA method to the regime of glacier modeling, considering both slip conditions and steep bedrock and surface inclinations. The method is assessed with regards to stability and accuracy for a synthetic two-dimensional (2D) case with a randomly generated bedrock topography, using a novel method based on so-called Perlin noise (Perlin, 1985) and a real-world application to the glacier of Midtre Lovénbreen, Svalbard. The experiments are carried out using the ice-sheet solver Elmer/Ice (Gagliardini et al., 2013), in which the FSSA method has been implemented.

The rest of the paper is structured as follows: in Sect. 2 the equations governing the flow of ice are presented; Sect. 3 introduces the numerical methods, including a presentation of the FSSA method for ice-sheet and glacier simulations; in Sect. 4, the experiments are presented along with their results; and finally the paper is concluded in Sect. 5 with a discussion of the results and the general outlook of the FSSA method from the perspective of glacier modeling.

2 Governing equations

2.1 The Stokes equations

The dynamics of ice flow can be described as a very slow-moving gravity-driven highly viscous fluid and is as such governed by the Stokes equations (see, e.g., Greve and Blatter, 2009):

$$\begin{array}{}\text{(1)}& {\displaystyle }\mathrm{\nabla }\cdot \left(\mathrm{2}\mathit{\eta }\left(\mathit{u}\right)\dot{\mathit{\epsilon }}\left(\mathit{u}\right)\right)-\mathrm{\nabla }p& {\displaystyle }=\mathit{\rho }g\widehat{\mathit{z}},\mathit{x}\in \mathrm{\Omega },\text{(2)}& {\displaystyle }\mathrm{\nabla }\cdot \mathit{u}& {\displaystyle }=\mathrm{0},\mathit{x}\in \mathrm{\Omega },\end{array}$$

where Eq. (1) follows from conservation of momentum and Eq. (2) follows from the conservation of mass. Furthermore, $\dot{\mathit{\epsilon }}=\frac{\mathrm{1}}{\mathrm{2}}\left(\mathrm{\nabla }\mathit{u}+\mathrm{\nabla }{\mathit{u}}^T\right)$ is the strain-rate tensor, and u and p are the ice velocity and pressure, respectively, at spatial coordinate x in the domain Ω⊂ℝ^d (see Fig. 1), where $d\in \mathit{\{}\mathrm{2},\mathrm{3}\mathit{\}}$ is the geometrical dimension. Furthermore, ρ=910 kg m⁻³ is the ice density and g=9.8 m s⁻² is the acceleration due to gravity. Lastly, η is the effective viscosity, which for ice depends on the velocity and temperature through Glen's flow law (Glen, 1955; Nye, 1957):

$$\begin{array}{}\text{(3)}& \mathit{\eta }(\mathit{u},T^{\prime })=A(T^{\prime }{)}^{-\frac{\mathrm{1}}{n}}{\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}tr\left({\dot{\mathit{\epsilon }}}^{\mathrm{2}}\right)+{\dot{\mathit{\epsilon }}}_{\mathrm{0}}^{\mathrm{2}}\right)}^{\frac{\mathrm{1}-n}{\mathrm{2}n}},\end{array}$$

where n=3 (Cuffey and Paterson, 2010) is the Glen or power-law exponent and ${\dot{\mathit{\epsilon }}}_{\mathrm{0}}^{\mathrm{2}}={\mathrm{10}}^{-\mathrm{10}}$ yr⁻² is a small regularization term added in order to avoid an infinite viscosity at strain rates of 0. The rate factor A(T^′) depends on the ice temperature relative to the pressure melting point T^′ through the Arrhenius equation (Glen, 1955):

$$\begin{array}{}\text{(4)}& A\left(T^{\prime }\right)=A_{\mathrm{0}}\exp \left(-{\displaystyle \frac{Q}{R{T}^{\prime }}}\right),\end{array}$$

where $A_{\mathrm{0}}=\mathrm{2.89165}\times {\mathrm{10}}^{-\mathrm{13}}$ s⁻¹ Pa⁻³ is a pre-exponential factor, Q=60 kJ mol⁻¹ is the activation energy, and R=8.314462 J K⁻¹ mol⁻¹ is the ideal gas constant. The values stated here are the ones recommended by Paterson (1994) for ice at a temperature $T^{\prime }\ge -\mathrm{10}$ °C when n=3.

Figure 1 Cross section of a generic glacier domain Ω with its bedrock Γ_b marked in black and its surface Γ_s marked in gray. The west Γ_W and east Γ_E sides are marked in orange. The dark-red lines in the interior are the flow lines of the velocity field u, and a_s is the accumulation/ablation rate.

2.2 Boundary conditions

In order to specify the appropriate boundary condition (BC), the glacier boundary ∂Ω is divided into non-overlapping boundary parts Γ_s, Γ_W, Γ_b, and Γ_E (see Fig. 1). The ice surface Γ_s is the only non-stationary part of the domain, meaning the future of the glacier is determined purely by the evolution of the surface. For the different parts of the boundary, the following BCs are considered:

$$\begin{array}{}\text{(5)}& {\displaystyle }\mathit{\sigma }\widehat{\mathit{n}}& {\displaystyle }=\mathbf{0},\mathit{x}\in {\mathrm{\Gamma }}_{\mathrm{s}},\text{(6)}& {\displaystyle }\mathit{u}\cdot \widehat{\mathit{n}}& {\displaystyle }=\mathrm{0},\mathit{x}\in \partial \mathrm{\Omega }/{\mathrm{\Gamma }}_{\mathrm{s}},\text{(7)}& {\displaystyle }{\widehat{\mathit{t}}}_i\cdot \mathit{\sigma }\widehat{\mathit{n}}& {\displaystyle }=-\mathit{\beta }|\mathit{u}{|}^{m-\mathrm{1}}\mathit{u}\cdot {\widehat{\mathit{t}}}_i,\mathit{x}\in {\mathrm{\Gamma }}_{\mathrm{b}}^{\mathrm{s}},\text{(8)}& {\displaystyle }\mathit{u}& {\displaystyle }=\mathbf{0},\mathit{x}\in {\mathrm{\Gamma }}_{\mathrm{b}}^{\mathrm{f}},\end{array}$$

where $\mathit{\sigma }=\mathrm{2}\mathit{\eta }\dot{\mathit{\epsilon }}-p\mathbf{I}$ (I is the identity matrix) is the Cauchy stress tensor, $\widehat{\mathit{n}}$ is the unit normally outward pointing to the boundary, $\mathit{\{}{\widehat{\mathit{t}}}_i{\mathit{\}}}_{i=\mathrm{0}}^{d-\mathrm{1}}$ is tangent vectors spanning the plane defined by $\widehat{\mathit{n}}$, β≥0 is the drag coefficient, and m≥1 is an exponent.

The explanation of each BC is as follows: Eq. (5) is a stress-free condition on the glacier surface, following from the assumption that the stresses asserted on the surface due to, for instance, wind or the atmospheric pressure are negligible compared to the internal stresses (Greve and Blatter, 2009). The second BC, Eq. (6), is an impenetrability condition under which ice cannot flow into the bedrock, meaning its velocity in the direction normal to the bedrock must necessarily be 0. The third BC, Eq. (7), is a Weertman-type sliding law (Weertman, 1957), stating that the ice may slip along the bedrock, following a power-law relation between the slip velocity and the shear stress. This study focuses only on the case for which m=1 such that the relation is linear. Lastly, the fourth BC, Eq. (8), is a no-slip BC representing conditions where the ice is frozen to bedrock. The bedrock thus consists of the following parts: ${\mathrm{\Gamma }}_{\mathrm{b}}^{\mathrm{s}}$, where slip is present, and ${\mathrm{\Gamma }}_{\mathrm{b}}^{\mathrm{f}}$, where no slip occurs.

2.3 The free-surface equation

The time evolution of a glacier (or an ice sheet) is determined by its surface position $z_{\mathrm{s}}=z_{\mathrm{s}}(x,y,t)$ and is governed by a separate equation called the free-surface equation (Greve and Blatter, 2009):

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{\partial z_{\mathrm{s}}}{\partial t}}+u_x^{\mathrm{s}}{\displaystyle \frac{\partial z_{\mathrm{s}}}{\partial x}}+u_y^{\mathrm{s}}{\displaystyle \frac{\partial z_{\mathrm{s}}}{\partial y}}=u_z^{\mathrm{s}}+a_{\mathrm{s}},\end{array}$$

where a_s is the vertical rate of mass accumulation (or ablation) and ${\mathit{u}}^{\mathrm{s}}=(u_x^{\mathrm{s}},u_y^{\mathrm{s}},u_z^{\mathrm{s}})$ is the velocity field from the Stokes equations, Eqs. (1)–(2), evaluated on the surface boundary Γ_s (see Fig. 1). Furthermore, the bedrock z_b(x,y) is assumed to be impenetrable and rigid such that the following constraint is fulfilled at all times t:

$$\begin{array}{}\text{(10)}& z_{\mathrm{s}}(x,y,t)\ge z_{\mathrm{b}}(x,y)+H_{min},\end{array}$$

where H_min is the minimum ice thickness. Equation (10) together with the weak formulation of Eq. (9) forms a variational inequality, which is solved using a method of imposed Dirichlet conditions, as described in Gagliardini et al. (2013).

3 Computational aspects

3.1 Solution procedure

A first-order time-stepping approach for solving the Stokes equations coupled to the free-surface equation is shown in Fig. 2a and consists of first solving the Stokes equations, Eqs. (1)–(2), for the velocity field evaluated on the surface u^s, which then enters as coefficients into the free-surface equation, Eq. (9). The free-surface equation is then solved for a new height function $z_{\mathrm{s}}(x,y,t+\mathrm{\Delta }t)$ which in turn determines the new domain Ω(t+Δt). The mesh is updated based on an extruded-mesh principle, wherein nodes are vertically aligned in columns such that the mesh can be updated by simply displacing nodes vertically according to the new height function $z_{\mathrm{s}}(x,y,t+\mathrm{\Delta }t$); see, e.g., Löfgren et al. (2022) for implementation details. This process is repeated until the final simulation time is reached. This is the standard approach in ice-sheet modeling (used in, e.g., Elmer/Ice; Gagliardini et al., 2013) and is in this study referred to as an explicit time-stepping scheme in terms of velocity.

This explicit time-stepping scheme can be contrasted with a Picard linearized implicit time-stepping scheme, with respect to the coupled system, updating both velocities and geometry simultaneously (Bueler, 2016, 2022). An example of a first-order implicit scheme, available in Elmer/Ice, is shown in Fig. 2b, where an extra loop is needed in order to the resolve the velocity field u(t+Δt) over the next time step. This has the disadvantage that the computationally expensive nonlinear Stokes equations need to be solved repeatedly in each time step, by iterating back and forth between the domain at the old time step Ω(t) and the domain at the new time step Ω(t+Δt). The advantage is that it is numerically stable, allowing for large time-step sizes. The goal of this paper is to evaluate an approach, the FSSA, that finds a solution which is close to the solution yielded by the implicit time-stepping scheme without adding the extra computationally costly iteration. It is thus an approach that uses the explicit time-stepping scheme in a way that is stable and without substantial loss of accuracy.

Figure 2 Examples of first-order (a) explicit (b) and implicit (Picard linearization) time-stepping schemes. For each time step the explicit scheme solves the Stokes equations, Eqs. (1)–(2), for the velocity field u(t) on the domain of the current time step Ω(t); the solution u(t) then enters as a coefficient into the free-surface equation, Eq. (9), from which the domain at the next time step Ω(t+Δt) is obtained directly and the time is updated, i.e., $t=t+\mathrm{\Delta }t$. The implicit scheme on the other hand solves for the velocity at the next time step u(t+Δt). The procedure for obtaining u(t+Δt) is as follows: in each time step an initial guess u₀(t) is inserted into the free-surface equation to obtain an estimate of the new domain $\stackrel{\mathrm{̃}}{\mathrm{\Omega }}(t+\mathrm{\Delta }t)$. The Stokes equations are then solved on $\stackrel{\mathrm{̃}}{\mathrm{\Omega }}(t+\mathrm{\Delta }t)$ to obtain an estimate $\stackrel{\mathrm{̃}}{\mathit{u}}(t+\mathrm{\Delta }t)$. This estimate is then checked if it is close to the initial guess u₀(t). If the estimate is not close, the process is repeated, using $\stackrel{\mathrm{̃}}{\mathit{u}}(t+\mathrm{\Delta }t)$ as the new initial guess, until convergence is obtained. Finally, after convergence, the time is updated. For both cases the algorithm terminates when the final simulation time T is reached.

3.2 The Stokes weak formulation – the basis for the stabilization method

The Stokes equations, Eqs. (1)–(2), are discretized and solved numerically using the finite element method (FEM), which first requires recasting the problem in its weak form: find $(\mathit{u},p)\in \mathcal{X}\times \mathcal{Q}$ such that

$$\begin{array}{}\text{(11)}& \begin{array}{rl}\left(\dot{\mathit{\epsilon }}\right(\mathit{v})& :\mathrm{2}\mathit{\eta }\dot{\mathit{\epsilon }}\left(\mathit{u}\right){)}_{\mathrm{\Omega }}-\left(\right(\mathrm{\nabla }\cdot \mathit{v}),p{)}_{\mathrm{\Omega }}-(q,\mathrm{\nabla }\cdot \mathit{u}{)}_{\mathrm{\Omega }}\\ & +(\mathit{v},\mathit{\beta }\mathit{u}{)}_{{\mathrm{\Gamma }}_{\mathrm{b}}}=-(\mathit{v},\mathit{\rho }g\widehat{\mathit{z}}{)}_{\mathrm{\Omega }},\end{array}\end{array}$$

for all $(\mathit{v},q)\in \mathcal{X}\times \mathcal{Q}$. Here the colon operator A:B between matrices A and B denotes their Frobenius inner product. Furthermore, 𝒳 and 𝒬 are appropriate function spaces satisfying the so-called inf-sup stability condition (Ladyzhenskaya, 1969; Babuška, 1971; Brezzi, 1974). The fact that the forcing term is constant inside the integrals of the inner products opens up for the construction of the FSSA method, as will be described in the next section.

The nonlinear nature of Eq. (11) requires linearizing the viscous term with Picard or Newton iterations. Convergence issues sometimes prohibit using Newton solvers for glaciological problems. To overcome these issues relaxation methods can be employed.

3.3 Free-surface time discretization and stabilization

The weak formulation of the Stokes-coupled free-surface equations, Eqs. (9)–(11), is discretized in time by first evaluating all integrals in Eq. (11) at time $t=t^{k+\mathit{\theta }}$ so that the weak formulation reads as follows: find $({\mathit{u}}^{k+\mathit{\theta }},p^{k+\mathit{\theta }})\in \mathcal{X}\times \mathcal{Q}$ such that

$$\begin{array}{}\text{(12)}& \begin{array}{rl}\left(\dot{\mathit{\epsilon }}\right(\mathit{v})& :\mathrm{2}\mathit{\eta }\dot{\mathit{\epsilon }}\left({\mathit{u}}^{k+\mathit{\theta }}\right){)}_{{\mathrm{\Omega }}^{k+\mathit{\theta }}}-(\mathrm{\nabla }\cdot \mathit{v},p^{k+\mathit{\theta }}{)}_{{\mathrm{\Omega }}^{k+\mathit{\theta }}}\\ & -(q,\mathrm{\nabla }\cdot {\mathit{u}}^{k+\mathit{\theta }}{)}_{{\mathrm{\Omega }}^{k+\mathit{\theta }}}+(\mathit{v},\mathit{\beta }{\mathit{u}}^{k+\mathit{\theta }}{)}_{{\mathrm{\Gamma }}_{\mathrm{b}}^{k+\mathit{\theta }}}\\ & =-(\mathit{v},\mathit{\rho }g\widehat{\mathit{z}}{)}_{{\mathrm{\Omega }}^{k+\mathit{\theta }}},\end{array}\end{array}$$

for all $(\mathit{v},q)\in \mathcal{X}\times \mathcal{Q}$. Here k denotes the time step and θ∈ℝ is an implicitness parameter for which θ=0 yields an explicit solver and θ=1 yields an implicit solver. Secondly, a semi-implicit Euler discretization is employed for the free-surface equation, Eq. (9), such that

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{z}_{\mathrm{s}}^{k+\mathrm{1}}& +\mathrm{\Delta }t\left((u_x^{\mathrm{s}}{)}^{k+\mathit{\theta }}{\displaystyle \frac{\partial z_{\mathrm{s}}^{k+\mathrm{1}}}{\partial x}}+(u_y^{\mathrm{s}}{)}^{k+\mathit{\theta }}{\displaystyle \frac{\partial z_{\mathrm{s}}^{k+\mathrm{1}}}{\partial y}}\right)\\ & =z_{\mathrm{s}}^k+\mathrm{\Delta }t\left((u_z^{\mathrm{s}}{)}^{k+\mathit{\theta }}+a_{\mathrm{s}}^{k+\mathrm{1}}\right),\end{array}\end{array}$$

where ${\mathit{u}}_s^{k+\mathit{\theta }}=\left(\right(u_x^{\mathrm{s}}{)}^{k+\mathit{\theta }},(u_y^{\mathrm{s}}{)}^{k+\mathit{\theta }},(u_z^{\mathrm{s}}{)}^{k+\mathit{\theta }})$ is the surface velocity obtained from Eq. (12), Δt is the time-step size, and $z_{\mathrm{s}}^{k+\mathrm{1}}$ is the unknown surface at time step k+1 to be solved for. The scheme is called semi-implicit since it is implicit in terms of the surface z_s and explicit in terms of velocity u when θ=0.

The explicit nature, when θ=0, of the time-stepping scheme in Eq. (13) makes it prone to numerical instabilities. A fully implicit scheme (θ=1) would, on the other hand, involve a highly expensive computation in order to evaluate u^k+1 (see Fig. 2b). To circumvent this problem, the FSSA was introduced by Kaus et al. (2010) for mantle-convection problems and adapted to glaciological problems in Löfgren et al. (2022). The FSSA mimics a fully implicit scheme by approximating all integrals on the left-hand side in Eq. (12) by ${\int }_{{\mathrm{\Omega }}^{k+\mathit{\theta }}}\cdot \mathrm{d}\mathrm{\Omega }\approx {\int }_{{\mathrm{\Omega }}^k}\cdot \mathrm{d}\mathrm{\Omega }$ and estimating the force term over the new domain by applying the Reynolds transport theorem (see, e.g., Greve and Blatter, 2009):

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\underset{\mathrm{\Omega }\left(t\right)}{\int }f\mathrm{d}\mathrm{\Omega }=\underset{\mathrm{\Omega }\left(t\right)}{\int }{\displaystyle \frac{\partial f}{\partial t}}\mathrm{d}\mathrm{\Omega }+\underset{\partial \mathrm{\Omega }\left(t\right)}{\int }({\mathit{u}}_{\mathrm{b}}\cdot \widehat{\mathit{n}})f\mathrm{d}\mathrm{\Gamma }.\end{array}$$

In ice-sheet modeling ${\mathit{u}}_{\mathrm{b}}=\mathit{u}+a_{\mathrm{s}}\widehat{\mathit{z}}$ is the velocity of the moving boundary and $f=-\mathit{\rho }g\widehat{\mathit{z}}\cdot \mathit{v}$. Since f is constant in time the first term on the right-hand side of Eq. (14) becomes 0. Taking this into account gives, using a forward Euler scheme, the estimate of the forcing on the time step k+θ as

$$\begin{array}{}\text{(15)}& \begin{array}{rl}{\left(\mathit{v},\mathit{\rho }g\widehat{\mathit{z}}\right)}_{{\mathrm{\Omega }}^{k+\mathit{\theta }}}& \approx {\left(\mathit{v},\mathit{\rho }g\widehat{\mathit{z}}\right)}_{{\mathrm{\Omega }}^k}\\ & +{\left(\mathit{\theta }\mathit{\rho }g\mathrm{\Delta }t(\mathit{v}\cdot \widehat{\mathit{z}}),(\mathit{u}+a_{\mathrm{s}}\widehat{\mathit{z}})\cdot \widehat{\mathit{n}}\right)}_{{\mathrm{\Gamma }}_{\mathrm{s}}^k}.\end{array}\end{array}$$

Note how the integrals of the weak form and the fact that ρg is constant are features that open up the use of the Reynolds transport theorem in this simple way.

Inserting Eq. (15) into the Stokes weak formulation, Eq. (11), and approximating the integrals on the left-hand side by ${\int }_{{\mathrm{\Omega }}^{k+\mathit{\theta }}}\cdot \mathrm{d}\mathrm{\Omega }\approx {\int }_{{\mathrm{\Omega }}^k}\cdot \mathrm{d}\mathrm{\Omega }$ yields the FSSA-stabilized weak formulation appropriate for glaciology: find $({\stackrel{\mathrm{̃}}{\mathit{u}}}^{k+\mathit{\theta }},{\stackrel{\mathrm{̃}}{p}}^{k+\mathit{\theta }})\in \mathcal{X}\times \mathcal{Q}$ such that

$$\begin{array}{}\text{(16)}& \begin{array}{rl}\left(\dot{\mathit{\epsilon }}\right(\mathit{v})& :\mathrm{2}\mathit{\eta }\dot{\mathit{\epsilon }}\left({\stackrel{\mathrm{̃}}{\mathit{u}}}^{k+\mathit{\theta }}\right){)}_{{\mathrm{\Omega }}^k}-(\mathrm{\nabla }\cdot \mathit{v},{\stackrel{\mathrm{̃}}{p}}^{k+\mathit{\theta }}{)}_{{\mathrm{\Omega }}^k}-(q,\mathrm{\nabla }\cdot {\stackrel{\mathrm{̃}}{\mathit{u}}}^{k+\mathit{\theta }}{)}_{{\mathrm{\Omega }}^k}\\ & +(\mathit{v},\mathit{\beta }{\stackrel{\mathrm{̃}}{\mathit{u}}}^{k+\mathit{\theta }}{)}_{{\mathrm{\Gamma }}_{\mathrm{b}}^k}+(\mathit{\rho }g\mathit{\theta }\mathrm{\Delta }t\mathit{v}\cdot \widehat{\mathit{z}},{\stackrel{\mathrm{̃}}{\mathit{u}}}^{k+\mathit{\theta }}\cdot \widehat{\mathit{n}}{)}_{{\mathrm{\Gamma }}_{\mathrm{s}}^k}\\ & =-(\mathit{v},\mathit{\rho }g\widehat{\mathit{z}}{)}_{{\mathrm{\Omega }}^k}-(\mathit{\rho }g\mathit{\theta }\mathrm{\Delta }t{a}_{\mathrm{s}}\mathit{v}\cdot \widehat{\mathit{z}},\widehat{\mathit{n}}\cdot \widehat{\mathit{z}}{)}_{{\mathrm{\Gamma }}_{\mathrm{s}}^k},\end{array}\end{array}$$

for all $(\mathit{v},q)\in \mathcal{X}\times \mathcal{Q}$. For the FSSA, letting θ=1, the velocity ${\stackrel{\mathrm{̃}}{\mathit{u}}}^{k+\mathrm{1}}$ is only an approximation of the solution u^k+1 obtained from the fully implicit scheme. The validity of the FSSA follows from the fact that the gravitational force is driving the ice flow.

To better understand the effect of the stabilization term, insight can be gained by applying the FSSA to the SIA approximation of the Stokes equations, for which it can be shown that the FSSA approximately coincides the evaluation of the pressure at the end of the time integration (see the Appendix in Löfgren et al., 2022).

3.4 Weak formulation of the free-surface equation

The time-discretized free-surface equation, Eq. (13), is discretized spatially using the FEM, which requires recasting it to its weak formulation: find $z_{\mathrm{s}}^{k+\mathrm{1}}\in \mathcal{V}$ such that

$$\begin{array}{}\text{(17)}& \begin{array}{rl}(v,z_{\mathrm{s}}^{k+\mathrm{1}}{)}_{{\mathrm{\Gamma }}_{\mathrm{s}}^{⟂}}& +\mathrm{\Delta }t{\left(v,({\stackrel{\mathrm{̃}}{u}}_x^{\mathrm{s}}{)}^{k+\mathit{\theta }}{\displaystyle \frac{\partial z_{\mathrm{s}}^{k+\mathrm{1}}}{\partial x}}\right)}_{{\mathrm{\Gamma }}_{\mathrm{s}}^{⟂}}\\ & +\mathrm{\Delta }t{\left(v,({\stackrel{\mathrm{̃}}{u}}_y^{\mathrm{s}}{)}^{k+\mathit{\theta }}{\displaystyle \frac{\partial z_{\mathrm{s}}^{k+\mathrm{1}}}{\partial y}}\right)}_{{\mathrm{\Gamma }}_{\mathrm{s}}^{⟂}}=(v,z_{\mathrm{s}}^k{)}_{{\mathrm{\Gamma }}_{\mathrm{s}}^{⟂}}\\ & +\mathrm{\Delta }t{\left(v,({\stackrel{\mathrm{̃}}{u}}_z^{\mathrm{s}}{)}^{k+\mathit{\theta }}\right)}_{{\mathrm{\Gamma }}_{\mathrm{s}}^{⟂}}\\ & +\mathrm{\Delta }t{\left(v,a_{\mathrm{s}}^{k+\mathrm{1}}\right)}_{{\mathrm{\Gamma }}_{\mathrm{s}}^{⟂}},\forall v\in \mathcal{V},\end{array}\end{array}$$

where ${\mathrm{\Gamma }}_{\mathrm{s}}^{⟂}\subset {\mathbb{R}}^{d-\mathrm{1}}$ is the projection of the free surface Γ_s onto the underlying plane (or line in two dimensions) with z=0. This advection-type equation is stabilized in Elmer/Ice by either the residual-free bubbles (RFB) method (Baiocchi et al., 1993) or streamline upwind Petrov–Galerkin (SUPG) stabilization (Franca and Frey, 1992). The stabilizing impact of the transport stabilization is investigated in this study.

4 Numerical experiments

4.1 Overview

In this section two experiments using varying bedrock slopes and sliding conditions are presented to demonstrate the applicability of the FSSA method to glacier modeling and to assess its stabilizing properties. In the first experiment, the method is applied to a 2D flow-line case, with an undulating bedrock generated using gradient noise (see Appendix A), superimposed on a sloping bedrock. The FSSA method is investigated with regards to accuracy and stability of different bedrock slopes, mesh resolutions, and upwinding schemes. The second experiment applies the FSSA method to the real-world glacier of Midtre Lovénbreen, Svalbard, and is also evaluated based on stability and accuracy.

4.2 Experiment 1: 2D “Perlin” glacier

4.2.1 Setup: advancing glacier

This experiment consists of a 2D glacier geometry with a sloping, undulating bedrock, where accumulation and sliding conditions are present. The bedrock is generated by superimposing three gradient-noise octaves (see Appendix A) on a parabola such that

$$\begin{array}{}\text{(18)}& \begin{array}{rl}{z}_{\mathrm{b}}\left(x\right)& ={\displaystyle \frac{\mathit{\alpha }}{L_x}}{\left(x-L_x\right)}^{\mathrm{2}}+C^{\mathrm{1}}{\text{octave}}_{\mathrm{1}}(\mathrm{\Delta }{x}^{\mathrm{1}},x)\\ & +C^{\mathrm{2}}{\text{octave}}_{\mathrm{2}}(\mathrm{\Delta }{x}^{\mathrm{2}},x)+C^{\mathrm{3}}{\text{octave}}_{\mathrm{3}}(\mathrm{\Delta }{x}^{\mathrm{3}},x),\end{array}\end{array}$$

where α is the average slope, L_x=8000 m is the horizontal extent of the domain, and the octaves represent noise of different frequencies. The coefficient of the parabola has been chosen such that $\frac{\mathrm{d}{z}_{\mathrm{b}}}{\mathrm{d}x}\left(\mathrm{0}\right)=\mathrm{2}\mathit{\alpha }$ and $\frac{\mathrm{d}{z}_{\mathrm{b}}}{\mathrm{d}x}\left(L_x\right)=\mathrm{0}$. The noise amplitudes are set to C¹=300 m, C²=500 m, and C³=600 m, and the respective octave frequencies are set to Δx¹=2000 m, Δx²=1000 m, and Δx³=500 m. The resulting bedrock is visible in Fig. 3.

The initial ice surface is a thin layer of ice:

$$\begin{array}{}\text{(19)}& z_{\mathrm{s}}(x,\mathrm{0})=z_{\mathrm{b}}\left(x\right)+\mathrm{10}\text{m}.\end{array}$$

To build up a glacier on the bedrock, a non-negative accumulation function that is linearly decaying (with the horizontal coordinate) with a maximum at x=0 is used:

$$\begin{array}{}\text{(20)}& a\left(x\right)=max\left(\mathrm{1}-{\displaystyle \frac{\mathrm{3}x}{L_x}},\mathrm{0}\right).\end{array}$$

The BCs imposed are impenetrability, Eq. (6), on the west- and east boundaries, Γ_W and Γ_E (see Fig. 1). On the surface Γ_s, the free-surface condition, Eq. (5), is applied. Lastly, on the bedrock Γ_b, the impenetrability BC, Eq. (6), is combined with the linear Weertman sliding law of Eq. (7), with a drag coefficient given by

$$\begin{array}{}\text{(21)}& \mathit{\beta }\left(x\right)={\mathit{\beta }}_{min}+{\displaystyle \frac{{\mathit{\beta }}_{max}-{\mathit{\beta }}_{min}}{\mathrm{1}+\exp \left(\frac{x-\mathit{\mu }}{\mathit{\sigma }}\right)}},\end{array}$$

where ${\mathit{\beta }}_{max}=\mathrm{1000}$ MPa yr m⁻¹, ${\mathit{\beta }}_{min}=\mathrm{0.01}$ MPa yr m⁻¹, σ=200 m, and μ=3000 m. This drag coefficient should be viewed as a transition from a no-slip condition when x≪μ to a free-slip condition when x≫μ, with the length of the transition zone controlled by σ.

Firstly, a study is performed to investigate how both stability and accuracy of the FSSA method are influenced by increasing bedrock slopes for an advancing glacier. Simulations are performed on three domains with different average bedrock slopes in Eq. (18): α=0.05 (≈2.9°) (gently sloping glacier), α=0.1 (≈5.7°) (moderately sloping glacier), and α=0.2 (≈11.5°) (steeply sloping glacier).

To estimate the error, a reference solution is obtained for all three cases by performing simulations using a short time-step size Δt=0.05 years and a fine mesh resolution with $(N_x,N_z)=(\mathrm{1000},\mathrm{10})$, where N_x and N_z are the number of layers in the horizontal and vertical directions, respectively. The reference simulations are performed until final times t=900, 700, and 500 years, for the respective cases. Subsequent simulations using coarser temporal resolutions are then started from an intermediate glacier surface obtained from the reference simulation, starting from times t=500, 400, and 300 years. The error is then estimated by comparing the ice thickness $H\left(x\right)=z_{\mathrm{s}}\left(x\right)-z_{\mathrm{b}}\left(x\right)$ of the coarse solution to the thickness of the reference simulation H_ref at the final times for the respective slopes. The ice thickness error ϵ_H is computed as

$$\begin{array}{}\text{(22)}& {\mathit{ϵ}}_H={\displaystyle \frac{\left|\right|H_{\mathrm{ref}}-H|{|}_{\mathrm{2}}}{\left|\right|H_{\mathrm{ref}}|{|}_{\mathrm{2}}}},\end{array}$$

where $\left|\right|\cdot |{|}_{\mathrm{2}}$ denotes the discrete L2-norm. The error is then measured for multiple coarse-temporal-resolution simulations, with and without the FSSA, for various time-step sizes Δt.

A stability study is conducted for estimating the LST for different mesh resolutions for the case of α=0.1. Starting from an intermediate glacier surface obtained at time t=300 years, 15 time steps are performed using a fixed Δt. In each time step, stability is examined by calculating the infinity norm of the velocity field: the solver is said to be unstable if $\left|\right|u_z|{|}_{\mathrm{\infty }}\ge \mathrm{100}$ m yr⁻¹. If the simulation is deemed stable, then Δt is incremented by 1 year and the simulation is restarted.

Following Gong et al. (2017) and Löfgren et al. (2022), all simulations use a constant rate factor A=100 yr⁻¹ MPa⁻³ in Eq. (3) and a constant ice density ρ=910 kg m⁻³. The elements used to solve the Stokes equations are the inf-sup stable Taylor–Hood elements (Taylor and Hood, 1973), and RFB is used to stabilize the transport problem.

4.2.2 Setup: retreating glacier

A study is also conducted to compare the stabilizing impact of the FSSA method to the standard upwinding schemes of RFB and SUPG, both of which are available in Elmer/Ice. In this case the methods are applied to a glacier subject to a negative net mass balance. Investigating such a case is of interest for two reasons: firstly, since the glacier approaches a steady state, it allows for performing very long-term simulations without the glacier reaching the end of the domain; secondly, a negative accumulation could potentially affect the stabilizing impact of the FSSA due to triggering the surface limiter imposing the minimum ice thickness, for which in this study the FSSA has not been adapted to take into consideration.

Stability is, as in the previous of the advancing-glacier study, evaluated based on the relative increase in the LST as compared to not using the FSSA. For this reason the same stability study as in previous experiment is performed for the three cases: no upwinding, RFB, and SUPG. The starting glacier surface is the final surface obtained at time t=700 years of the reference simulation. All simulations use a mesh resolution of $(N_x,N_z)=(\mathrm{1000},\mathrm{10})$.

To obtain a retreating glacier, melting is introduced into the accumulation function by simply letting

$$\begin{array}{}\text{(23)}& a\left(x\right)=\mathrm{1}-{\displaystyle \frac{\mathrm{3}x}{L_x}},\end{array}$$

which is essentially the same accumulation function as Eq. (20) but allows for negative values.

In order for the glacier to experience sliding throughout the simulation, the drag coefficient is modified so that slip occurs predominantly in the interior of the domain:

$$\begin{array}{}\text{(24)}& \mathit{\beta }\left(x\right)={\mathit{\beta }}_{min}+{\displaystyle \frac{{\mathit{\beta }}_{max}-{\mathit{\beta }}_{min}}{\mathrm{1}+\exp \left(\frac{x-\mathit{\mu }}{\mathit{\sigma }}\right)}}+{\displaystyle \frac{{\mathit{\beta }}_{max}-{\mathit{\beta }}_{min}}{\mathrm{1}+\exp \left(\frac{L_x-x-\mathit{\mu }}{\mathit{\sigma }}\right)}},\end{array}$$

where ${\mathit{\beta }}_{max}=\mathrm{1000}$ MPa yr m⁻¹, ${\mathit{\beta }}_{min}=\mathrm{0.001}$ MPa yr m⁻¹, σ=100 m, and μ=1500 m. Since this experiment is designed to study the effect of upwinding, the minimum drag coefficient is reduced by a factor of 10, compared to the previous advancing case, in order for the Stokes-coupled free-surface system to admit a more transport-like behavior, i.e., large horizontal velocities with a shear-to-slip ratio close to 0.

4.2.3 Results: advancing glacier

The reported estimated ice thickness error ϵ_H, as calculated by Eq. (22), are shown in Table 1 for different bedrock slopes α and FSSA stabilization parameters θ. The FSSA method (θ=1) is seen to be stable for time-step sizes Δt up to 25 years for α=0.05 and α=0.1, while for α=0.2 it was stable for all tested Δt values, up to 50 years. In the case of θ=0 the largest stable time-step size (LST) is between 5 and 10 years for α=0.05 and α=0.1 and between 10 and 25 years for α=0.1 and α=0.2. The FSSA method thus increases the LST by at least a factor of 2 for all cases and may even be as large as 5 times for α=0.05. Compared to ice-sheet simulations in, e.g., Löfgren et al. (2022), the time-step sizes are large even without stabilization.

In Table 2 the LST is reported for different mesh sizes for the intermediate sloping case of α=0.1. It is seen that without the FSSA that the LST is by all practical means mesh independent, which is in agreement with the mesh studies carried out by Löfgren et al. (2022). On the other hand, for the FSSA a slight mesh dependence of the LST ∼Δx^0.4 is seen. This means that the stabilizing effect of the FSSA decreases for higher mesh resolutions, where the relative increase in the LST reduces from 3 for the coarser mesh $(N_x,N_z)=(\mathrm{400},\mathrm{10})$ to 2.1 for the finer mesh $(N_x,N_z)=(\mathrm{1000},\mathrm{10})$.

The LST is larger for the steep-bedrock case, despite the velocity field also having a larger magnitude (see Fig. 3d–f). The reason for this might be related to the cases of low and moderately inclined bedrock having greater ice thickness, for which analytical expressions derived using zeroth-order approximation, e.g., SIA, have shown a strong inverse relation between the LST and ice thickness (Gong et al., 2017; Robinson et al., 2022). In Löfgren et al. (2022) it was shown that the characteristics of the instabilities for Stokes-coupled free-surface flow are related to the domain aspect ratio: thicker domains tend to give rise to long wavelength sloshing instabilities, while thin domains give rise to numerical oscillation of shorter wavelengths.

In the current experiment, thick domains are represented by α=0.05 and α=0.1, where velocities are low enough for a thicker ice to develop, and the thin domain is represented by α=0.2. Indeed from Fig. 4, which shows a time series of vertical velocity profiles for unstable time-step sizes, it is seen that for low and moderate bedrock slopes, Fig. 4a–c and d–f, respectively, the instabilities behave differently than for the case of steep bedrock, Fig. 4g–i. In the former cases the instabilities emerge as result of the vertical velocity profiles shifting in sign and growing in magnitude between time steps, resulting in the glacier surface sloshing around the stable reference surface (dashed black lines in Fig. 4). For the steep-bedrock case, while sloshing is possible to discern in the interior, the most prominent feature of the instabilities is the high-frequency numerical oscillations occurring at the glacier front, seen in Fig. 4i. Despite the different characteristic of the instabilities arising for the different cases, it is clear from Table 1 and Fig. 3 that the FSSA method mitigates instabilities in both cases, as was also concluded by Löfgren et al. (2022).

Comparing the errors in the FSSA method (θ=1) to no stabilization (θ=0) in Table 1, it is seen that the FSSA method generally yields a much more accurate solution, and as Δt increases, the discrepancy grows in favor of the FSSA. For Δt≥5 years the error for all slopes is almost twice as large when θ=0 compared to θ=1. For some cases, even using a time-step size twice as large with the FSSA compared to no FSSA, the increase in the error is much less than the expected 100 %.

Generally, the error increases linearly with the time-step size Δt, i.e., ϵ_H=𝒪(Δt), as is expected of the semi-implicit Euler time-stepping scheme in Eq. (13), as long as stability restrictions are satisfied. This observation holds regardless of slope α and stabilization parameter θ.

Figure 3 shows the glacier surfaces at the final times for different Δt and θ values. In all cases it is seen from the zoom-in plots that too large a time-step size gives too thin a glacier front, compared to the reference solution. Comparing Table 1 and Fig. 3, it is seen that a larger error corresponds to a thinner glacier front. Consequently, the FSSA method, which is generally more accurate in this experiment, yields a faster moving front for the same Δt. This is expected based on the fact that the FSSA method is a quasi-implicit time-stepping scheme, meaning it uses an estimate of the velocity from the next time step to update the glacier surface in Eq. (9). Since the glacier is growing in size, the velocity field is expected to increase in magnitude over the duration of the simulation such that $\left|\right|{\mathit{u}}^{k+\mathrm{1}}\left|\right|\ge \left|\right|{\mathit{u}}^k\left|\right|$ (k denotes the time step), meaning that the FSSA method gives a larger velocity coefficient in the free-surface equation and thus yields a faster moving front.

The error also increases with the bedrock slope α, despite the simulation times being shorter, which is expected given that larger slopes give a higher velocity coefficient in the free-surface equation, Eq. (9). From the large error seen for the steep-bedrock case of α=0.2, it is not even clear that stability considerations are the limiting factor for the time-step size, it might as well in practical applications come down to accuracy – depending on whether Eq. (22) is deemed a satisfactory error metric. From a practical point of view, this has the implication that using a time-step size close to the LST may not be a wise strategy when employing the FSSA method as the error in the final ice thickness is quite large for, e.g., the stable case of Δt=50 years, especially for the long-term simulations considered in this experiment. On the other hand, the smaller errors observed for α=0.05 and α=0.1 indicate that stability considerations are more important for the time-step size. Thus, for a given error tolerance, the FSSA method seems to offer the greatest potential for speedup for simulating glaciers that are on top of bedrock with a topography that is gently to moderately inclined. Regardless of the surface slope, it is obvious that the FSSA method not only gives a more stable solver but also increases its accuracy.

In summary, the main finding is that the FSSA method allows for using larger time-step sizes by increasing both accuracy and stability, for all slopes angles investigated. However, the method seems to offer the greatest benefit for gently to moderately sloping glaciers, where the time-step size seems to be mainly limited by stability considerations.

4.2.4 Results: retreating glacier

The LST for different numerical schemes and FSSA stabilization parameters θ is reported in Table 3. Compared with the LST in the advancing case (Table 2), which uses RFB, it is seen that the LST in this case is about 3 times smaller. This may be explained by the larger velocities of 90 m yr⁻¹ (see Fig. 5a) compared to 20 m yr⁻¹ (see Fig. 3e). Despite this, it is interesting that the larger velocities do not seem to have an impact on the stability of the FSSA, and in fact the LST is larger in this case (39 years compared to 27 years). However, it should be noted that the FSSA method yields smaller velocities as Δt increases (see Fig. 5), which for the case of Δt=39 years decreased to 45 m yr⁻¹. The fact that the LST is larger in this case demonstrates that the FSSA is also applicable to cases where the free surface is constrained by a minimum ice thickness.

Furthermore, it is seen that upwinding alone appears to have a negligible impact on the stability of the solver, despite surface velocities being dominated by sliding (see Fig. 5). However, combining the FSSA with upwinding increases the LST slightly for RFB (∼12 %) and substantially for SUPG (∼60 %), compared to the FSSA alone. The best choice stability-wise thus seems to be combining the FSSA and SUPG.

It was also found for the large time-step sizes Δt>25 years allowed for by the FSSA that non-zero surface velocities started arising in the deglaciated areas. However, these could effectively be mitigated by setting the accumulation part of the FSSA to 0 in deglaciated areas. This did not compromises stability nor accuracy as all simulations deemed stable were found to approach the same steady state.

Table 1 Relative error in the ice thickness, as defined by Eq. (22), for different time-step sizes Δt, the FSSA parameter θ, and bedrock slopes α. The error is calculated after final times of 900, 700, and 500 years for the respective slopes of α=0.05, 0.1, and 0.2. Entries marked with an X are unstable cases.

Table 2 Largest stable time-step size (LST) with and without FSSA stabilization for different mesh resolutions (N_x,N_z), where N_x and N_z denote the number of horizontal and vertical layers, respectively. The resolution for the reported LST is 1 year, meaning that it is between 12 and 13 years for, e.g., the case of θ=0 and $(N_x,N_z)=(\mathrm{400},\mathrm{10})$. The last column shows the relative increase in the LST.

Figure 3 (a–c) Glacier surfaces and (d–f) horizontal velocities of the reference simulations at final times T for different bedrock slopes α, FSSA stabilization parameters θ, and time-step sizes Δt. The solid gray lines in (a)–(c) are the reference surface z_s for the respective case, and the solid black lines are the bedrock z_b. The solid blue lines in (d)–(f) are the surface velocities $u_x^{\mathrm{s}}$, and the solid red lines are the bedrock slip velocities $u_x^{\mathrm{b}}$. The reference solution was obtained using a short time-step size Δt=0.05 years.

Figure 4 Vertical velocity profiles u_z from simulations using no FSSA (θ=0) and unstable time-step sizes of (a–c) Δt=10 years and (d–i) Δt=25 years shown at the indicated times t. The profiles are shown for the different bedrock slopes of (a–c) α=0.05, (d–f) α=0.1, and (g–i) α=0.2. The dashed black line in each figure is the glacier surface obtained from the stable reference simulation.

Table 3 Largest stable time-step size (LST) with and without FSSA stabilization for the numerical schemes: no upwinding, residual-free bubbles (RFB), and streamline upwind Petrov–Galerkin (SUPG). The resolution for the reported LST is 1 year, meaning that it is between 3 and 4 years for, e.g., the cases of θ=0. The last column shows the relative increase in the LST.

Figure 5 Horizontal velocity profiles (a) without the FSSA and (b) with the FSSA using stabilization parameter θ=1 and time-step size Δt=10 years. The solid blue line is the surface velocity $u_x^{\mathrm{s}}$, and the solid red line is the bedrock slip velocity $u_x^{\mathrm{b}}$.

4.3 Experiment 2: Midtre Lovénbreen

4.3.1 Setup

This experiment aims to demonstrate how the FSSA works for a real-world, three-dimensional (3D) glacier simulation. For this purpose, the valley glacier Midtre Lovénbreen, Svalbard (78.53° N, 12.04° E), is chosen as it is a thoroughly studied glacier that has been modeled using Elmer/Ice previously (Zwinger and Moore, 2009; Välisuo et al., 2017). The glacier is in this study classified as a moderately to steeply sloping glacier with an average surface slope α of 0.14 (8°). The initial geometry is shown in Fig. 6, where colors indicate ice thickness (Fig. 6a), surface mass balance (SMB; Fig. 6b), and velocity magnitude (Fig. 6c) and where gray contour lines represent surface heights in meters above sea level.

Figure 6 (a) Ice thickness, (b) surface mass balance (SMB) excluding the deglaciated area (gray), and (c) modeled surface velocity and surface elevation of Midtre Lovénbreen at the start of the simulation, in the year 1995. The ice flows from top to bottom. The computational domain is given by the 1962 extent of the glacier, and the white line indicates the outline of the glacier at the start of the simulation. The figures were created using the open-source visualization and analysis toolkit PyVista (Sullivan and Kaszynski, 2019).

The basal BC is the linear sliding law of Eq. (7) with the same drag coefficient as in Välisuo et al. (2017):

$$\begin{array}{}\text{(25)}& \mathit{\beta }(x,y,t)=\left\{\begin{array}{l}\mathrm{0.04}{\text{MPa\hspace{0.17em}yr\hspace{0.17em}m}}^{-\mathrm{1}}\text{if }H(x,y,t)\ge \mathrm{120}\text{m},\\ \mathrm{10}{\text{MPa\hspace{0.17em}yr\hspace{0.17em}m}}^{-\mathrm{1}}\text{if }H(x,y,t)<\mathrm{120}\text{m},\end{array}\right.\end{array}$$

where H is the ice thickness. This drag coefficient imposes high slip velocities at parts where the ice is thick (≥120 m) and essentially imposes a no-slip condition at the shallower parts (<120 m). Note that this also means that sliding velocities decrease as the glacier thins. The values for the drag coefficient were determined in Välisuo et al. (2017) by a manual inversion using observed surface velocities as input data. As in Välisuo et al. (2017), the bedrock elevation is given by a bedrock DEM (© Norwegian Polar Institute) created from ground-penetrating radar data (Rippin et al., 2003; Zwinger and Moore, 2009), the 1995 surface elevation DEM (© Norwegian Polar Institute) is based on digital photogrammetry from vertical aerial photographs, and the 2005 surface elevation DEM is a product derived from airborne lidar (light detection and ranging) data (James et al., 2006; Kohler et al., 2007).

Following Välisuo et al. (2017), the SMB used has been estimated from surface DEMs of the ice surface for the years 1995 and 2005. The SMB is obtained by solving the Stokes equations, Eqs. (1)–(2), on the domain defined by the 1995 surface DEM; from the time-discretized free-surface equation, Eq. (13), the SMB is given by

$$\begin{array}{}\text{(26)}& \begin{array}{rl}{a}_{\mathrm{s}}& ={\displaystyle \frac{z_{\mathrm{s}}^{\mathrm{2005}}-z_{\mathrm{s}}^{\mathrm{1995}}}{\mathrm{10}\text{ yr}}}+u_x^{\mathrm{1995}}{\displaystyle \frac{\partial z_{\mathrm{s}}^{\mathrm{2005}}}{\partial x}}\\ & +u_y^{\mathrm{1995}}{\displaystyle \frac{\partial z_{\mathrm{s}}^{\mathrm{2005}}}{\partial y}}-u_z^{\mathrm{1995}}.\end{array}\end{array}$$

The 3D mesh is generated by extruding a 2D footprint mesh into five layers (N_z=5) with a horizontal resolution of about 25 m. The footprint is large enough to cover the glacier at the size it was in 1962, which implies that a large deglaciated area is included in the domain.

For this experiment, linear equal-order elements are used in conjunction with Galerkin least-squares stabilization (Hughes et al., 1986) to circumvent the inf-sup stability condition. Furthermore, the free surface is constrained by a minimum ice thickness of $H_{min}=\mathrm{10}$ m. The transport problem is stabilized using SUPG.

The glacier evolution is simulated from the year 1995 to 2195 with and without the FSSA and compared and evaluated in regards to the LST and accuracy. To measure accuracy a reference solution is obtained by performing a simulation from the year 1995 to 2195 using a small time-step size of Δt=1 year and no FSSA (θ=0). The error is estimated by comparing the ice thickness in the final solutions to the reference solution, for various time-step sizes. Stability is for this experiment evaluated qualitatively, based on the presence of spurious shifts in the sign of the vertical velocity, i.e., sloshing, since it was found that checking the norm of the velocity alone did not accurately predict the presence of instability in this case.

To remedy convergence issues of the Newton solver (which appeared both with and without the FSSA) the derivative of the viscosity appearing in the Newton linearization is relaxed by a factor of $\mathrm{2}/\mathrm{3}$.

Figure 7 Midtre Lovénbreen in the year 2195. (a) Glacier outlines of the reference solution (white) using a short time-step size Δt=1 year, as well as outlines for a simulation using the larger Δt=40 years, with the FSSA (dark gray) and without the FSSA (orange). The thickness of the glacier as given by the reference solution is indicated with colors in panel (a). (b) Ice thickness along the black line in panel (a) for the reference solution (dashed black line in b) and simulations with and without the FSSA for Δt=20 and 40 years (solid lines in b). (c) Ice thickness errors as compared to the reference solution. Panel (a) was created using the open-source visualization and analysis toolkit PyVista (Sullivan and Kaszynski, 2019).

Figure 8 Vertical velocity profiles u_z at indicated times t on Midtre Lovénbreen for cases with (a–c) a stable Δt=5 years without the FSSA, (d–f) an unstable Δt=25 years without the FSSA, and (g–h) a stable Δt=25 years with the FSSA. The figure was created using the open-source visualization and analysis toolkit PyVista (Sullivan and Kaszynski, 2019).

4.3.2 Results

The final glacier surface after 200 years is shown in Fig. 7a, where the color denotes the ice thickness of the reference simulation. Compared with the thickness of the initial glacier in Fig. 6a, it is seen that the glacier surface has retreated a distance of about 1.5 km uphill. In addition the glacier has also experienced thinning from initially having a maximum thickness of about 180 m down to a maximum thickness of 100 m. The retreat of the glacier is expected given the considerably negative SMB, as can be seen from Fig. 6b. Note, however, that predicting the actual retreat of the glacier is not the objective of this study but rather demonstrates the stabilizing properties of the FSSA method for an SMB derived from experimental data. The true SMB is in reality likely to change substantially over the simulation period considered.

Figure 7a shows the glacier outlines for the stable reference simulation (white line), an FSSA-stabilized simulation with Δt=40 years (dark gray line), and an unstabilized simulation with Δt=40 years (orange line). It is seen that the unstabilized case deviates substantially from the reference simulation, due to instability, while the glacier outline from the FSSA simulation is to a large extent indistinguishable from the outline of the reference simulation. The largest time-step size tested for stability in each case is shown in Table 4, where it is seen that the FSSA simulation is stable for time-step sizes at least twice as large. From Fig. 7b it is seen, similarly to the 2D Perlin experiment, that instabilities arise in the deep parts of the domain, where the error is large – indicative of the sloshing-type instability encountered in the previous experiment. Indeed from Fig. 8d–f, instability is observed, in the sense of spurious shifts in the sign of the vertical velocity, for the case with Δt=25 years and θ=0 in the upper parts of the domain. These are effectively mitigated by reducing the time-step size (Fig. 8a–c) or by means of stabilization (Fig. 8g–i).

Furthermore, despite the unstable surface deviating substantially from the reference, the norm of the velocity never grew unboundedly in this case. Thus determining the presence of instability based on the norm of the velocity alone may not accurately predict the presence of instabilities. The cause of this may be the significantly negative SMB which causes glacier thinning and in a sense “stabilizes” the solver over time. However, as is evident from Fig. 7b the solution is still polluted by the initial instability. For this reason, the presence of the vertical velocity shifting in sign spuriously was taken as a criterion for instability in this case.

Figure 7c shows that for the stable time-step size Δt=20 years, the error using the FSSA is larger than without the FSSA at the deep parts of the glacier, but the accuracy near edges is higher so that the glacier area is more accurate (Fig. 7a and c). This is contrary to what was found in the previous 2D experiment, where the FSSA method was overall more accurate. The reason for this discrepancy could be explained by the fact that the estimated reference solution in this case is only refined temporally and therefore may not represent the analytical solution accurately enough. Figure 7b and c also show that the less accurate the solution is, the smaller the glacier retreat and larger ice thickness are compared to the reference solution. This should be contrasted with the first experiment which considered an advancing glacier, where the less accurate the solution is, the less the glacier had advanced.

In order to ensure computation times are not negatively affected by the FSSA method, an experiment was performed to measure the CPU time of the FSSA method compared to no stabilization. The computation times for different Δt values are shown in Table 5. It is seen that the FSSA method was faster for all cases. The FSSA method was about 10 % faster for the case with Δt=20 years, while the difference is only slightly in favor of the FSSA for Δt=10 years and Δt=5 years. The difference in computation times seems to be related in this case to the average number of nonlinear iterations needed for convergence. For example, for the case with Δt=20 years the FSSA method required about 10 % fewer iterations for convergence, which explains the 10 % difference in computation times. The increase in nonlinear iterations might be due to stability issues of the unstabilized solver as the relative difference becomes smaller with a decreasing Δt.

In summary, the FSSA method gave a more stable solution, increasing the LST by at least a factor of 2, without negatively impacting computation times. In regards to accuracy, the FSSA method yielded larger ice thickness errors in the interior, while the error was reduced at the glacier front. This has the implication that the FSSA method may be more suitable for predicting future glacier extent, while the non-FSSA method is more accurate for determining future glacier thinning.

Table 4 Largest stable time-step size of a 200-year-long simulation with and without FSSA stabilization for two different mesh resolutions Δx.

Table 5 CPU times for a 200-year-long simulation for different time-step sizes Δt, measured with and without FSSA stabilization. The simulations were run on a mesh with a horizontal resolution Δx∼25 m. Note that case with Δt=40 years and θ=0 is unstable.

5 General discussion and conclusions

The FSSA was implemented into Elmer/Ice and tested on simulations of synthetic glaciers as well as on Midtre Lovénbreen, Svalbard. The FSSA increased the largest stable time-step size (LST) by a factor of 2 for the simulation of Midtre Lovénbreen and up to a factor of 5 for the synthetic Perlin glacier with a low surface slope. Low glacier surface slopes were correlated with a shorter LST for the unstabilized method and were also the cases where the stabilization had the greatest effect. This may be due to the thicker ice that developed on the glaciers with low bedrock inclination, for which stability restrictions on lower order models have shown a strong inverse dependence (e.g., Gong et al., 2017; Robinson et al., 2022). For a case approaching a steady state, it was also found that despite the large time-step sizes allowed for by the FSSA, the same steady state was approached compared to a reference simulation using a short time-step size. Even without stabilization the LST was already quite large in many cases; for instance, for the advancing Perlin glacier with an intermediate slope, an LST of 13 years was observed, compared to 2 to 4 years for the ice-sheet simulations in, e.g., Löfgren et al. (2022). The larger LST observed here likely stems from a combination of lower flow speeds (∼20 m yr⁻¹ compared to ∼40 m yr⁻¹) and the considerably thinner ice (∼200 m compared to ∼2000 m).

It was also found that when compared to other means of stabilization, e.g., residual-free bubbles (RFB) and streamline upwind Petrov–Galerkin (SUPG), the FSSA was the only one when used alone that increased the LST. However, combining the FSSA and SUPG was the most stable choice. Furthermore, the mesh study revealed that the LST of the unstabilized case was mesh independent, in agreement with Löfgren et al. (2022), while the FSSA admitted a slight mesh dependence at high spatial resolutions. However, even for the finest mesh resolution investigated, the FSSA had twice the LST compared to no stabilization.

The FSSA mitigated instabilities and improved the accuracy overall, with the only exception being in the deep parts of Midtre Lovénbreen. As the computational cost is also low, the FSSA can be added as a security measure in simulations to prevent sloshing instability polluting the result, without negatively impacting accuracy or simulation times. For glacier simulations this may offer the greatest benefit of the FSSA, given the already large time-step sizes of the unstabilized method. Another potential usage are spin-up simulations, where a very large time-step size of, e.g., 50 years could be used. Climate data could then be incorporated on shorter timescales using semi-implicit time stepping of the surface.

In the previous study by Löfgren et al. (2022) the FSSA was tested in synthetic ice-sheet experiments, largely disregarding the effect of variable bedrock and sliding conditions and simplifying the contact problem near glacier fronts. This paper demonstrates that the FSSA method is applicable to more complex real-world simulations and the new implementation in Elmer/Ice makes the method accessible to a broad user base.

Some limitations of the current study are the generally low flow speeds (<100 m yr⁻¹) and that only linear sliding laws were considered. The FSSA method also remains to be adapted to higher-order time-stepping schemes, which could possibly yield better stability properties due to stronger coupling between the geometry and velocity, as has been demonstrated by Wirbel and Jarosch (2020) for a second-order Runge–Kutta scheme. Considering higher-order time stepping is an ongoing project for the authors.

Appendix A: Bedrock generation

This section presents an algorithm that is used to randomly generate bedrock topographies. It is based on a method that is common practice in computer graphics as a computationally inexpensive and flexible way of randomly generating visually appealing landscapes, clouds, textures, etc., known as gradient noise. The first application of gradient noise, the so-called Perlin noise, was developed by Perlin (1985) to model fire, water, and wrinkled surfaces and later adapted by Musgrave et al. (1989) for landscape generation. As the name suggests, gradient noise is based on generating a set of pseudo-random gradients at predefined vertices and then interpolating polynomials matching these gradients such that the resulting global interpolating function is smooth. In this study, the mesh vertices are assumed to be equally spaced with a spatial period Δx^f. The superscript f is used to denote the fact that the final noise function may consist of multiple spatial frequencies, so-called octaves – similar to a Fourier decomposition. See Fig. A1 for an example of a randomly generated bedrock consisting of three octaves.

Figure A1 Gradient noise consisting of three octaves (solid orange, blue, and yellow lines in b–d) superimposed on a parabola (dashed gray line in a) to generate a natural-looking bedrock topography (solid black line in a). (b–d) Spatial periods of the octaves are Δx¹=2000 m (solid orange line in b), Δx²=1000 m (solid blue line in c), and Δx³=500 m (solid yellow line in d). The black dots and solid gray lines in (b)–(d) denote the nodal values (zero in this case) and the matched pseudo-randomly generated tangents.

For the interpolation cubic Hermite polynomials are used:

$$\begin{array}{}\text{(A1)}& {\displaystyle }{p}^k\left(x\right){\displaystyle }=\sum _{i=\mathrm{0}}^{\mathrm{3}}{c}_i{x}^i,\end{array}$$

where p^k is the polynomial interpolation in cell k and c_i is coefficients determined by matching gradients and function values at the vertices. This gives the set of equations for each vertex x_i in cell k the following:

$$\begin{array}{}\text{(A2)}& {\displaystyle }{p}^k\left({\mathit{x}}_i\right)& {\displaystyle }=\mathrm{0},\text{(A3)}& {\displaystyle \frac{d{p}^k\left({\mathit{x}}_i\right)}{\mathrm{d}x}}& {\displaystyle }=f_{x_i},\end{array}$$

where $f_{x_i}$ is uniformly randomly generated gradients over the interval $[-\mathrm{1},\mathrm{1}]$, corresponding to slope angles between −45 and 45°. Since the number of vertices in each cell is two, this leads to a total of four equations, matching the total number of unknowns. Solving the resulting linear system then gives the octave characterized by the spatial period Δx^f of the noise function. The final noise function is obtained as a linear combination of all the octaves such that

$$\begin{array}{}\text{(A4)}& \text{noise}(x,y)=\sum _{f=\mathrm{1}}^n{A}^{\mathrm{f}}\text{octave}(\mathrm{\Delta }{x}^{\mathrm{f}},x,y),\end{array}$$

where A^f and Δx^f are the amplitude and spatial period of the fth octave, respectively, and n is the number of octaves.