| 12 Nov 2025
Grounding-line dynamics in a Stokes ice-flow model (Elmer/Ice v9.0): Improved numerical stability allows larger time steps
Abstract. The efficient and accurate simulation of grounding line dynamics in marine ice-sheet models remains a challenge, largely due to restrictive time-step limitations. The restrictive time step size of ice-sheet simulations (∼ ∆t = 0.01−0.5 years) has led to the routine use of approximate models that compromise physical complexity compared to full Stokes models. To address the time-step restriction and enhance the applicability of full Stokes simulations, we implement a numerical stabilisation scheme at the ice-ocean interface, namely the Free-Surface Stabilisation Algorithm (FSSA). The FSSA acts by predicting the surface elevation at the next time step, resulting in a reduction in surface oscillations and an increase in the largest numerically stable time step. When applied to the ice-ocean interface, FSSA acts in combination with the sea spring numerical stabilisation scheme, allowing larger time steps to be taken.
In order to test the capabilities of the FSSA when applied to the ice-ocean interface, we perform the benchmark simulation of Experiment 3a from the Marine Ice Sheet Model Intercomparison Project (MISMIP). These simulations demonstrate the ability of the model to capture grounding line migration on both prograde slopes (oceanward sloping) and retrograde slopes (inland sloping). We find a time step size of ∆t = 10 years to be numerically stable and accurate in the MISMIP experiment, which is more than an order of magnitude larger than the small time steps traditionally used. In comparison, a time-step size of ∆t = 50 years can maintain numerical stability, but is not capable of capturing the full range of grounding-line motion in the MISMIP experiments. We further demonstrate the applicability of the FSSA to a 3D marine terminating model domain, finding that a time-step size of ∆t = 10 years is numerically stable. The increase in the largest numerically stable time step by greater than an order of magnitude in marine-terminating Stokes ice-flow problems through the inclusion of FSSA broadens the applicability of Stokes models, which have otherwise been deemed too computationally expensive for large-scale applications.
The preprint under consideration presents the free-surface stabilization scheme (FSSA), a stabilization scheme for the free-surface Stokes equation. The free surface Stokes equation is a time-dependent equation that is generally solved with explicit time-stepping schemes, due to the difficulty in solving the free-surface evolution and the Stokes equations simultaneoulsy. A consequence of using explicit schemes is that time steps have to be small. The FSSA has been shown to enable higher time steps without reducing the numerical accuracy signficantly. In this preprint, the authors study the impact of using the FSSA in the context of marine ice sheets, where the free surface detaches from the bedrock at the grounding line.
The preprint is well written and contains two well-documented numerical tests that demonstrate the power of the FSSA in the context of marine ice sheets. My main criticism of this paper is that it is very similar to other existing preprints and publications, including
- Josefin Ahlkrona, A. Clara J. Henry, and André Löfgren (2025), A Fully Implicit Second Order Method for Viscous Free Surface Stokes Flow – Application to Glacier Simulations (under consideration for publication in GMD),
- Tilda Westling Dolling, A. Clara J. Henry, and Josefin Ahlkrona (2025), A numerical stabilization scheme for the shallow shelf approximation (arXiv:2510.02943).
Since the need for stabilization schemes like this one is urgent in the glaciological community, I recommend this article for publication once the following comments are addressed.
General comments
- As indicated above, this paper is closely related to a series of papers on the FSSA (and its depth-integrated variants) written by some of the authors. I think the introduction should contain a more exhaustive description of what exactly is novel about this publication (its application to grounding line dynamics) and in what way it is related to existing publications.
- Section 4. For the MISMIP experiment, you only solve the unstabilized system ($\theta = 0$) for one timestep. Does the scheme work for higher timesteps if $\theta = 0$. If so, it would be important to include these results so that the reader understands the effects of not including the FSSA. If the scheme breaks down for higher timesteps when $\theta = 0$, please indicate this.
- The numbering of figures throughout the text seems to be incorrect. See for example reference to "Figure 4.2" in Line 281 and "Fig. 5.1" in Line 300.
Minor comments
- Line 3. "The restrictive time step size of ice-sheet simulations (∼ ∆t = 0.01−0.5 years) has led to the routine use of approximate models that compromise physical complexity compared to full Stokes models." In my opinion, approximations of the Stokes equations are used in glaciology because the Stokes equations, as opposed to e.g. depth-integrated models, are very expensive to solve. Please rewrite along the lines of ""The restrictive time step size of ice-sheet simulations (∼ ∆t = 0.01−0.5 years) is one of the reasons why... "
- Figure 1. This figure is not referred to in the text. Please refer to to it when introducing the notation for the domain $\Omega$ and its boundary.
- Line 81. Clarify that, despite the Stokes equations appearing to be stationary in time, time-dependence enters the problem because the domain Ω changes in time.
- Lines 105 and 110. Writing $v \in \mathbb{R}^3$ and $q \in \mathbb{R}$ is very confusing. Please rewrite as e.g. $v: \Omega \to \mathbb{R}^3$ and $q:\Omega \to \RR$.
- Figure 2. Please refer to this figure in the text.
- Line 137. "However, due to the elliptical nature of the Stokes equations, the coupled system of the Stokes and the free-surface equations cannot be fully implicit in the velocity." I do not understand this sentence. Are you not solving this coupled system implicity in the paper "A Fully Implicit Second Order Method for Viscous Free Surface Stokes Flow - Application to Glacier Simulations"? Please clarify and explain relationship to this paper.
- Line 185. What does the velocity with a tilde, $\tilde{u}$, refer to? It seems to be like this variable has not been defined.
- Line 221. Write "Table 1" as opposed to "Table (1)". Same for "Table (2)" in the next line.
- Line 245. Why are Picard iterations used? Is a Newton solver not implemented in Elmer/Ice? The p-Stokes equations are the optimality conditions for differentiable convex optimization problem and therefore work very well, in general, when Newton's method is applied.
- Line 308. "The isothermal ice fluidity, A, and the friction coefficient, C, were tuned..." Please indicate how they were tuned. Did you solve an inverse problem to find the corresponding scalar values?
- Equation (28). Is this choice of boundary condition for the ice divide common? Please reference other works where this is used. What value is chosen for $C_n$? I can imagine that, as $C_n$ becomes very large, an impenetrability boundary condition is effectively enforced (this is equivalent to penalizing u.n at the boundary). Is this true? If so, please indicate this in the text to facilitate understanding for the reader.
- Line 343. "Figs 5.1 and 5.1". Please correct.
- Line 363. "a time-step size of ∆ = 50 years". Please correct.
- Line 382. "We hypothesis...". Please write verb correctly.