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 (

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

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

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

One way of stabilizing the problem is to use a fully implicit time-stepping scheme, which has been demonstrated by

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

These issues were addressed by

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

The rest of the paper is structured as follows: in Sect.

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

Cross section of a generic glacier domain

In order to specify the appropriate boundary condition (BC), the glacier boundary

The explanation of each BC is as follows: Eq. (

The time evolution of a glacier (or an ice sheet) is determined by its surface position

A first-order time-stepping approach for solving the Stokes equations coupled to the free-surface equation is shown in Fig.

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

Examples of first-order

The Stokes equations, Eqs. (

The nonlinear nature of Eq. (

The weak formulation of the Stokes-coupled free-surface equations, Eqs. (

The explicit nature, when

Inserting Eq. (

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

The time-discretized free-surface equation, Eq. (

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

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

The initial ice surface is a thin layer of ice:

The BCs imposed are impenetrability, Eq. (

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. (

To estimate the error, a reference solution is obtained for all three cases by performing simulations using a short time-step size

A stability study is conducted for estimating the LST for different mesh resolutions for the case of

Following

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

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

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:

The reported estimated ice thickness error

In Table

The LST is larger for the steep-bedrock case, despite the velocity field also having a larger magnitude (see Fig.

In the current experiment, thick domains are represented by

Comparing the errors in the FSSA method

Generally, the error increases linearly with the time-step size

Figure

The error also increases with the bedrock slope

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.

The LST for different numerical schemes and FSSA stabilization parameters

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.

It was also found for the large time-step sizes

Relative error in the ice thickness, as defined by Eq. (

Largest stable time-step size (LST) with and without FSSA stabilization for different mesh resolutions

Vertical velocity profiles

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

Horizontal velocity profiles

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 (

The basal BC is the linear sliding law of Eq. (

Following

The 3D mesh is generated by extruding a 2D footprint mesh into five layers (

For this experiment, linear equal-order elements are used in conjunction with Galerkin least-squares stabilization

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

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

Midtre Lovénbreen in the year 2195.

Vertical velocity profiles

The final glacier surface after 200 years is shown in Fig.

Figure

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.

Figure

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

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.

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

CPU times for a 200-year-long simulation for different time-step sizes

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

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

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

Some limitations of the current study are the generally low flow speeds (

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

Gradient noise consisting of three octaves (solid orange, blue, and yellow lines in

For the interpolation cubic Hermite polynomials are used:

Elmer version

AL: conceptualization, project administration, investigation, methodology, software, validation, visualization, writing (original draft preparation). CH: software, investigation, writing (review and editing). JA: conceptualization, project administration, supervision, funding acquisition, writing (original draft preparation). PR and TZ: project administration, methodology, investigation, validation, software, resources, funding acquisition, writing (review and editing).

At least one of the (co-)authors is a member of the editorial board of

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

Josefin Ahlkrona received funding from the Swedish Research Council (grant no. 2021-04001), which funded the work of Josefin Ahlkrona, Thomas Zwinger, and André Löfgren. Thomas Zwinger's contribution was within the framework of the Academy of Finland COLD consortium (grant no. 322978). Peter Råback received funding from the European High Performance Computing Joint Undertaking (EuroHPC JU) and Spain, Italy, Iceland, Germany, Norway, France, Finland, and Croatia (grant agreement no. 101093038). Josefin Ahlkrona and Christian Helanow received funding from the Swedish e-Science Research Centre (SeRC). Computations on the CSC – IT Center for Science Ltd. platforms Mahti and cPouta were supported by the HPC-Europa3 program, part of the European Union's Horizon 2020 research and innovation program (grant agreement no. 730897). The publication of this article was funded by Swedish Research Council, Forte, Formas, and Vinnova.

This paper was edited by Johannes J. Fürst and reviewed by two anonymous referees.