A comparison of numerical surface topography calculations in geodynamic modelling: an evaluation of the ‘sticky air’ method

TLDR

Surface topography in geodynamic models is typically computed using free‑surface approximations, with the “sticky air” method emerging as a popular free‑surface proxy. This study theoretically derives the physical criteria that determine when the sticky‑air approximation accurately represents a true free surface. The authors analyze two benchmark scenarios—cosine‑perturbation isostatic relaxation and plume‑induced topography—comparing six finite‑difference and finite‑element codes that compute topography via normal‑stress, body‑fitting, and Lagrangian tracking, and examine marker‑distribution effects. The sticky‑air method remains valid when the dimensionless ratio (η_st/η_ch)/(h_st/L)^3 is sufficiently small, ensuring that the low‑viscosity air layer does not distort the free‑surface response.

Abstract

Calculating surface topography in geodynamic models is a common numerical problem. Besides other approaches, the so-called ‘sticky air’ approach has gained interest as a free-surface proxy at the top boundary. The often used free slip condition is thereby vertically extended by introducing a low density, low viscosity fluid layer. This allows the air/crust interface to behave in a similar manner to a true free surface. We present here a theoretical analysis that provides the physical conditions under which the sticky air approach is a valid approximation of a true free surface. Two cases are evaluated that characterize the evolution of topography on different timescales: (1) isostatic relaxation of a cosine perturbation and (2) topography changes above a rising plume. We quantitatively compare topographies calculated by six different numerical codes (using finite difference and finite element techniques) using three different topography calculation methods: (i) direct calculation of topography from normal stress, (ii) body-fitting methods allowing for meshing the topography and (iii) Lagrangian tracking of the topography on an Eulerian grid. It is found that the sticky air approach works well as long as the term (ηst/ηch)/(hst/L)3 is sufficiently small, where ηst and hst are the viscosity and thickness of the sticky air layer, and ηch and L are the characteristic viscosity and length scale of the model, respectively. Spurious lateral fluctuations of topography, as observed in some marker-based sticky air approaches, may effectively be damped by an anisotropic distribution of markers with a higher number of markers per element in the vertical than in the horizontal direction.