Global mantle flow at ultra-high resolution: The competing influence of faulted plate margins, the strength of bending plates, and large-scale, nonlinear flow

A full understanding of the dynamics of plate motions requires numerical models with a realistic, nonlinear rheology and a mesh resolution sufficiently high to resolve large variations in viscosity over short length scales. We suspect that resolutions as fine as 1 km locally in global models of the whole mantle and lithosphere are necessary. We use the adaptive mesh mantle convection code Rhea to model convection in the mantle with plates in both regional and global domains. Rhea is a new generation parallel finite element mantle convection code designed to scale to hundreds of thousands of compute cores. It uses forest-of-octree-based adaptive meshes via the p4est library. With Rhea's adaptive capabilities we can create local resolution down to ~ 1 km around plate boundaries, while keeping the mesh at a much coarser resolution away from small features. The global models in this study have approximately 160 million elements, a reduction of ~ 2000x compared to a uniform mesh of the same high resolution. The unprecedented resolution in these global models allows us, for the first time, to resolve viscous dissipation in the bending plate as well as observe the trade-off between this process and the strength of slabs and the resistance of dipping thrust faults. Since plate velocities and 'plateness' are dynamic outcomes of numerical modeling, we must carefully incorporate both the full buoyancy field and the details of all plate boundaries at a fine scale. The global models were constructed with detailed maps of the age of the plates and a thermal model of the seismicity-defined slabs which grades into the more diffuse buoyancy resolved with tomography. In the regional models, the thermal model consists of plates following a halfspace cooling model, and slabs for which buoyancy is conserved at every depth. A composite formulation of Newtonian and non-Newtonian rheology along with yielding is implemented; plate boundaries are modeled as very narrow weak zones. Plate boundary parameters such as the width and strength of the weak zones are varied, as well as the characteristics of the yielding law. The effect of these variations on slab strength, plate velocities, the state of stress in the slab and viscosity in the hinge of the subducting plate is examined, and trade-offs determined. The results of the regional study are then used to fine-tune the instantaneous dynamic global mantle convection models predicting plate velocities. These are tested by assessing the plateness of the surface velocity field, and its misfit with measured surface velocities. The global mantle flow models allow us to address the cause of changes in plate motions and the distribution of energy dissipation within the convective system.