Increasing stable time step sizes in ice sheet modelling

In order to understand the rate at which an ice sheet is losing mass one has to consider its dynamics. Ice is a very slow moving, highly viscous, non-newtonian fluid and as such is most accurately described by the full Stokes equation. Time dependence is taken into account by coupling the Stokes equation to the so called free surface equation, which describes how the free surface boundary of the ice sheet is advected due to the Stokes velocity field.

A problem with this system is that it is numerically quite unstable and has a very strict time step constraint, where very small time steps are needed in order to have a stable solver. This constitutes a severe limitation for making long term predictions as the expensive nonlinear Stokes equation has to be solved in each time step.

By adding an additional term to the weak form of the Stokes equation we achieve stability for time steps 10-20 times larger than without stabilization. This stabilization technique is straightforward to implement into existing code and does not result in significantly larger computation times or memory usage.