On the choice of finite element for applications in geodynamics

Geodynamical simulations over the past decades have widely been built on quadrilateral and hexahedral finite elements. For the discretization of the key Stokes equation describing slow, viscous flow, most codes use either the unstable Q(1) x P-0 element, a stabilized version of the equalorder Q(1) x Q(1) element, or more recently the stable Taylor-Hood element with continuous (Q(2) x Q(1)) or discontinuous (Q(2) x P-1) pressure. However, it is not clear which of these choices is actually the best at accurately simulating "typical" geodynamic situations. Herein, we provide a systematic comparison of all of these elements for the first time. We use a series of benchmarks that illuminate different aspects of the features we consider typical of mantle convection and geodynamical simulations. We will show in particular that the stabilized Q(1) x Q(1) element has great difficulty producing accurate solutions for buoyancy-driven flows - the dominant forcing for mantle convection flow - and that the Q(1) x P-0 element is too unstable and inaccurate in practice. As a consequence, we believe that the Q(2) x Q(1) and Q(2) x P-1 elements provide the most robust and reliable choice for geodynamical simulations, despite the greater complexity in their implementation and the substantially higher computational cost when solving linear systems.