TY - JOUR
T1 - A new free-surface stabilization algorithm for geodynamical modelling
T2 - Theory and numerical tests
AU - Andres-Martinez, Miguel
AU - Morgan, Jason
AU - Perez-Gussinye, Marta
AU - Ruepke, Lars
PY - 2015/9
Y1 - 2015/9
N2 - 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.
AB - 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.
KW - Free-surface stabilization
KW - Geodynamic modelling
UR - http://www.sciencedirect.com/science/article/pii/S0031920115000990
U2 - 10.1016/j.pepi.2015.07.003
DO - 10.1016/j.pepi.2015.07.003
M3 - Article
SN - 0031-9201
VL - 246
SP - 41
EP - 51
JO - Physics of the Earth and Planetary Interiors
JF - Physics of the Earth and Planetary Interiors
ER -