Increasing stable time-step sizes of the free-surface problem arising in ice-sheet simulations

Numerical models for predicting future ice mass loss of the Antarctic and Greenland ice sheets require accurately representing their dynamics. Unfortunately, ice-sheet models suffer from a very strict time-step size constraint, which for higher-order models constitutes a severe bottleneck; in each time step a nonlinear and computationally demanding system of equations has to be solved. In this study, stable time-step sizes are increased for a full-Stokes model by implementing a so-called free-surface stabilization algorithm (FSSA). Previously this stabilization has been used successfully in mantle-convection simulations where a similar viscous-flow problem is solved. By numerical investigation it is demonstrated that instabilities on the very thin domains required for ice-sheet modeling behave differently than on the equal-aspect-ratio domains the stabilization has previously been used on. Despite this, and despite the different material properties of ice, it is shown that it is possible to adapt FSSA to work on idealized ice-sheet domains and increase stable time-step sizes by at least one order of magnitude. The FSSA method presented is deemed accurate, efficient and straightforward to implement into existing ice-sheet solvers.