  Numerical models for predicting future ice-mass loss of the Antarctic and
Greenland ice sheet requires 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 since
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, but linear, 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 nonlinearity of the problem, 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 presented
is deemed accurate, efficient and straightforward to implement into existing
ice-sheet solvers.
