A stabilization algorithm for geodynamic numerical simulations with a free surface

Geodynamic processes are in many cases coupled to surface processes such as erosion of mountains or deposition of sediments in basins. Numerical models that study the effect of such surface processes on the dynamics of the lithosphere use either a self-consistent free surface approach or approximate it by adding a weak, low density “air” layer on top of the model domain. A main problem with such setups, however, is that they are numerically unstable, which is caused by the fact that the density difference that drives geodynamic processes (∼30–100 kg/m3) is much smaller than the density difference between rocks and air (∼3000 kg/m3). In order to retain isostatic balance, one should thus employ a time step that is at least 30 times smaller than in models with a free slip upper boundary condition, which makes such models computationally very expensive.Here we describe a new free surface stabilization approach (FSSA) that largely overcomes this time step restriction. In this work, we apply the stabilization technique to consistent free surface models utilizing a finite element discretization. The approach is based on the observation that most geodynamical simulations perform time-dependent simulations by solving the static Stokes equations and then updating the material or temperature fields in a separate step. If, however, the time-dependency of these quantities is taken into account during the discretization of the momentum equations, additional surface traction terms appear in the weak form of the FE formulation. These terms, which depend on the employed time step and velocity, essentially correct for the change of forces in an element due to advection or distortion of every element. By time-discretizing them, they act as a stabilizing term that largely increases the time steps that can be employed (by a factor 20–100 depending on the time stepping criteria that is used).We illustrate the accuracy and power of this method with examples of a free surface Rayleigh–Taylor instability and with a free subduction experiment.