On the Application of Compatible Finite Elements for Divergence-Free and Mesh-Independent Viscoelastic-Plastic Rheology in Geodynamics Simulations

The movement of the Earth's crust and mantle in geodynamics is typically modeled as the flow of a viscous fluid governed by the Stokes equations. Incorporating plasticity into material rheology often results in mesh-dependent behavior, which poses challenges for accurate numerical simulations. Several approaches have been proposed to mitigate mesh dependence and develop solvers that decouple solution errors from viscosity and mesh size.

Traditionally, Taylor-Hood (TH) and Crouzeix-Raviart (CR) elements of order 2 are used for geodynamics simulations. In this study, we examine the numerical solution of variable-viscosity Stokes flow with plasticity and Drucker-Prager type yielding using Scott-Vogelius (SV) compatible finite elements in combination with pseudo-Jacobian and augmented Lagrange methods. The Scott-Vogelius element is unique among finite elements for the mixed formulation of Stokes flow, as it has an associated De Rham complex. This theoretically ensures a divergence-free velocity field. We investigate the degree of decoupling between velocity errors, pressure errors, and viscosity-induced errors in a viscoelasto-plastic case study.

Our results show that Taylor-Hood elements (CG2 × CG1 for velocity and pressure) fail to provide accurate solutions in such cases. While the low-order CR elements perform better, the higher-order SV elements (CG4 × DG3) yield the best results. 

We conclude that due to the inherent mesh-dependent behavior and viscosity dependent errors in TH elements, CR or SV elements should be preferred for geodynamics simulations.