On the speed and numerical stability of ice-dynamics approximations

The Stokes solution to ice dynamics is computationally expensive, and in many cases unnecessary. Many approximations have been developed that reduce the complexity of the problem and thus reduce computational cost. Most approximations can generally be tuned to give reasonable solutions to ice-dynamics problems, depending on the domain and scale being simulated. However, the inherent numerical stability of time-stepping with different solvers has not been studied in detail. Here we investigate how different approximations lead to limits on the maximum timestep in mass conservation calculations for both idealized and realistic geometries. The ice-sheet models Yelmo and CISM are used to compare the following approximations: the shallow-ice approximation (SIA), the shallow-shelf approximation (SSA), the SIA+SSA approximation (Hybrid) and two variants of the L1L2 solver, namely one that reduces to SIA in the case of no-sliding (dubbed L1L2-SIA here) and the so-called depth-integrated viscosity approximation (DIVA). We find that these approaches vary significantly with respect to numerical stability. The extreme dependence on the local surface gradient of the SIA-based approximations (SIA, Hybrid, L1L2-SIA) leads to an amplified local velocity response and greater potential for instability, especially as grid resolution increases. In contrast, the SSA and DIVA approximations allow for longer time steps, because numerical oscillations in ice thickness are damped with increasing resolution. Given its high fidelity to the Stokes solution and its favorable stability properties, we demonstrate the strong case for using the DIVA approximation in many contexts.