A new free-surface stabilization algorithm for geodynamical modelling: Theory and numerical tests

Miguel Andres-Martinez, Jason Morgan, Marta Perez-Gussinye, Lars Ruepke

Research output: Contribution to journalArticlepeer-review

129 Downloads (Pure)


The surface of the solid Earth is effectively stress free in its subareal portions, and hydrostatic beneath the oceans. Unfortunately, this type of boundary condition is difficult to treat computationally, and for computational convenience, numerical models have often used simpler approximations that do not involve a normal stress-loaded, shear-stress free top surface that is free to move. Viscous flow models with a computational free surface typically confront stability problems when the time step is bigger than the viscous relaxation time. The small time step required for stability (< 2 Kyr) makes this type of model computationally intensive, so there remains a need to develop strategies that mitigate the stability problem by making larger (at least 10 ∼Kyr) time steps stable and accurate. Here we present a new free-surface stabilisation algorithm for finite element codes which solves the stability problem by adding to the Stokes formulation an intrinsic penalization term equivalent to a portion of the future load at the surface nodes. Our algorithm is straightforward to implement and can be used with both Eulerian or Lagrangian grids. It includes α and β parameters to respectively control both the vertical and the horizontal slope-dependent penalization terms, and uses Uzawa-like iterations to solve the resulting system at a cost comparable to a non-stress free surface formulation. Four tests were carried out in order to study the accuracy and the stability of the algorithm: 1) a decaying first-order sinusoidal topography test, 2) a decaying high-order sinusoidal topography test, 3) a Rayleigh-Taylor instability test, and 4) a steep-slope test. For these tests, we investigate which α and β parameters give the best results in terms of both accuracy and stability. We also compare the accuracy and the stability of our algorithm with a similar implicit approach recently developed by Kaus et al. (2010). We find that our algorithm is slightly more accurate and stable for steep slopes, and also conclude that, for longer time steps, the optimal α controlling factor for both approaches is ∼2/3, instead of the 1/2 Crank-Nicolson parameter inferred from a linearized accuracy analysis. This more-implicit value coincides with the velocity factor for a Galerkin time discretization applied to our penalization term using linear shape functions in time.
Original languageEnglish
Pages (from-to)41-51
Number of pages11
JournalPhysics of the Earth and Planetary Interiors
Early online date20 Jul 2015
Publication statusPublished - Sept 2015


  • Free-surface stabilization
  • Geodynamic modelling

Cite this