In prior work we found that precise approximation of the continuity constraint is crucial for accurate propagation of tracer data when advected through a background incompressible velocity field (Sime et al., 2021). Here we extend this investigation to compressible flows using the anelastic liquid approximation (ALA) and address four related issues: 1. exact conservation of tracer discretized fields through a background compressible velocity; 2. exact mass conservation; 3. addition and removal of tracers without affecting (exact) conservation to preserve a consistent number of tracers per cell; and 4. the diffusion of tracer data, for example, as induced by thermal or chemical effects. In this process we also present an abstract formulation of the interior penalty hybrid discontinuous Galerkin (HDG) finite element formulation for diffusion problems, and apply it to the advection-diffusion and compressible Stokes systems. Finally we present numerical experiments exhibiting the HDG compressible Stokes momentum formulation’s superconvergent compressibility approximation and reproduce community numerical benchmarks from the literature for the ALA.

Tracer methods are widespread in computational geodynamics for modeling the advection of chemical data. However, they present certain numerical challenges, especially when used over long periods of simulation time. We address two of these in this work: the necessity for mass conservation of chemical composition fields and the need for the velocity field to be pointwise divergence free to avoid gaps in tracer coverage. We do this by implementing the hybrid discontinuous Galerkin (HDG) finite element method combined with a mass-conserving constrained projection of the tracer data. To demonstrate the efficacy of this system we compare it to other common finite element formulations of the Stokes system and projections of the chemical composition. We provide a reference of the numerical properties and error convergence rates which should be observed by using these various discretization schemes. This serves as a tool for verification of existing or new implementations. We summarize these data in a reproduction of a published Rayleigh–Taylor instability benchmark, demonstrating the importance of careful choices of appropriate and compatible discretization methods for all aspects of geodynamics simulations.

England and May (2021) present a model comparison for the forearc thermal structure in subduction zones that employs a finite element modeling technique and a new variation on approximate equations. We have some comments on various claims made in this published paper.