Density-contrast instabilities in finite element modelling of slow viscous flow

Finite element methods are often used to model Earth processes involving slow viscous or viscoelastic flow. Inertial terms of the Navier-Stokes equations are neglected in very slow flows, so timestep size is not limited by the Courant instability. However, where there is advection of density contrasts in a gravitational field, over-advection can lead to numerically induced flow oscillations. We derive analytic results for the maximum stable timestep size in two cases: a free surface over a fluid of uniform density, and a free surface kept level by sedimentation/erosion, but with a density gradient in the underlying medium. Using parameters appropriate to the Earth's crust we show that the density-contrast instability occurs for timesteps larger than 3000 years for the constant-density case. For a fluid with a density gradient of 10 kg/m³ per km the solution is stable for timesteps up to about 200,000 years if full erosion/sedimentation is implemented.