Weak form Shallow Ice Approximation models with an improved time step restriction

TL;DR

This work tackles the limited time-step size in Shallow Ice Approximation (SIA) simulations of ice sheets when advancing the free surface. It shows that rewriting SIA in weak form and adding Free Surface Stabilization Algorithm (FSSA) terms yields a linear-in-$Δx$ dt restriction, rather than the conventional quadratic scaling, and extends the approach to a weak-form linear Stokes formulation with the SIA viscosity. Through theoretical cost analyses and extensive 2D numerical experiments on slab-on-slope, idealized ice sheets, and a Greenland cross-section, the study demonstrates that weak-form formulations—especially W-SIAStokes and W-SIAStokes-FSSA—achieve superior stability and efficiency, often matching or surpassing full Stokes accuracy at a fraction of the cost. The results provide practical guidance on parameter choices (notably the FSSA parameter $θ$) and indicate strong potential for applying these formulations in ice-sheet spin-up and long-time, isothermal simulations, with future work extending to non-isothermal physics and 3D geometries.

Abstract

The Shallow Ice Approximation (SIA) model on strong form is commonly used for inferring the flow dynamics of grounded ice sheets. The solution to the SIA model is a closed-form expression for the velocity field. When that velocity field is used to advance the ice surface in time, the time steps have to take small values due to quadratic scaling in terms of the horizontal mesh size. In this paper we write the SIA model on weak form, and add in the Free Surface Stabilization Algorithm (FSSA) terms. We find numerically that the time step restriction scaling is improved from quadratic to linear, but only for large horizontal mesh sizes. We then extend the weak form by adding the initially neglected normal stress terms. This allows for a linear time step restriction across the whole range of the horizontal mesh sizes, leading to an improved efficiency. Theoretical analysis demonstrates that the inclusion of FSSA stabilization terms transitions the explicit time stepping treatment of second derivative surface terms to an implicit approach. Moreover, a computational cost analysis, combined with numerical results on stability and accuracy, advocates for preferring the SIA models written on weak form over the standard SIA model.

Figures (11)

• Figure 1: Propagation of the surface elevation function in time when computed as a solution to the nonlinear Stokes problem (reference), where the simulation time step is $\Delta t = 0.1$ years. Horizontal mesh size and vertical mesh size are $\Delta x = 250$ meters and $\Delta y = 90$ meters respectively.
• Figure 2: Surface elevations at simulation time $t=6$ years for different SIA problem formulations and the reference formulation (nonlinear Stokes problem). Horizontal and vertical mesh sizes are $\Delta x = 250$ meters and $\Delta y = 90$ meters respectively. The largest feasible time step $\Delta t_*$ for the given $\Delta x$ and $\Delta y$ is chosen for each formulation: $\Delta t_* = 6$ years, $\Delta t_* = 12$ years, $\Delta t_* = 1.8$ years, $\Delta t_* = 0.04$ years, $\Delta t_* = 0.008$ years, $\Delta t = 0.1$ years for (W-SIAStokes-FSSA), (W-SIA-FSSA), (W-SIAStokes), (W-SIA), (SIA), (Reference) respectively.
• Figure 3: Left plot: scaling of the largest feasible time step $\Delta t_*$ as a function of the horizontal mesh size $\Delta x$ when the vertical mesh size $\Delta y = 90$ meters is fixed. Right plot: the model error as a function of the wall clock runtime when the nonlinear Stokes solution is taken as a reference. In all cases the final time is fixed at $t=20$ years. Mesh sizes $\Delta x = 250$ meters and $\Delta y = 90$ meters are also fixed. The time step is refined starting at the formulation largest feasible time step $\Delta t = \Delta t_*, \Delta t_*/2, \Delta t_*/4,...$ The time step for the reference solution is fixed at $\Delta t = 0.1$ years.
• Figure 4: Largest feasible time step $\Delta t_*$ as a function of the horizontal mesh size $\Delta x$ for the W-SIAStokes-FSSA formulation, when the FSSA parameter $\theta$ varies. Vertical mesh size is fixed at $\Delta y = 83.3$ meters (left plot), and at $\Delta y = 41.67$ meters.
• Figure 5: Left plot: scaling of the largest feasible time step $\Delta t_*$ as a function of the horizontal mesh size $\Delta x$ when the vertical mesh size $\Delta y = 90$ meters is fixed. Right plot: the model error as a function of the wall clock runtime when the nonlinear Stokes solution is taken as a reference. In all cases the final time is fixed at $t=20$ years. Mesh sizes $\Delta x = 250$ meters and $\Delta y = 90$ meters are also fixed.