On the choice of finite element for applications in geodynamics

Abstract. Geodynamical simulations over the past decades have widely been built on quadrilateral and hexahedral finite elements. For the discretisation of the key Stokes equation describing slow, viscous flow, most codes use either the unstable Q1 × P0 element, a stabilised version of the equal-order Q1 × Q1 element, or more recently the stable Taylor-Hood element with continuous (Q2 × Q1) or discontinuous (Q2 × P−1) pressure. However, it is not clear which of these choices is actually the best at accurately simulating typical geodynamic situations. Herein, we are providing for the first time a systematic comparison of all of these elements. 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 stabilised Q1 × Q1 element has great difficulty producing accurate solutions for buoyancy-driven flows – the dominant forcing for mantle convection flow – and that the Q1 × P0 element is too unstable and inaccurate in practice. As a consequence, we believe that the Q2 × Q1 and Q2 × 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.