LAUSR

Abstract Geodynamic simulations increasingly rely on models with a true free surface to investigate questions of dynamic topography, tectonic deformation, gravity perturbations, and global mantle convection. However, implementations of free surface boundary conditions have proven challenging from a standpoint of accuracy, robustness, and stability. In particular, time integration of a free surface tends to suffer from a numerical instability that manifests as sloshing surface motions, also known as the “drunken sailor” instability. This instability severely limits stable timestep sizes to those much smaller than can be used in geodynamic simulations without a free surface. Several schemes have been proposed in the literature to deal with these instabilities. Here we analyze the problem of creeping viscous flow with a free surface and discuss the origin of these instabilities. We demonstrate their cause and how existing stabilization schemes work to damp them out. We also propose a new scheme for removing instabilities from free surface calculations. It does not require modifications to the system matrix, nor additional variables, but is instead an explicit scheme based on nonstandard finite differences. It relies on a single stabilization parameter which may be identified with the smallest relaxation timescale of the free surface. Finally, we present numerical results to show the effectiveness of the new approach and discuss the free surface implementation in the open source, community based mantle convection software ASPECT .