C. R. Wilson

Interaction of subducted slabs with the mantle transition‐zone: A regime diagram from 2‐D thermo‐mechanical models with a mobile trench and an overriding plate

Transition zone slab deformation influences Earths thermal, chemical and tectonic evolution. However, the mechanisms responsible for the wide-range of imaged slab morphologies remain debated. Here, we use 2-D thermo-mechanical models with a mobile trench, an overriding plate, a temperature- and stress-dependent rheology, and a 10, 30 or 100-fold increase in lower mantle viscosity, to investigate the effect of initial subducting- and overriding-plate ages on slab transition-zone interaction. Four subduction styles emerge: (i) a “vertical folding” mode, with a quasi-stationary trench, near-vertical subduction and buckling/folding at depth (VF); (ii) slabs that induce mild trench retreat, which are flattened/“horizontally deflected” and stagnate at the upper-lower mantle interface (HD); (iii) inclined slabs, which result from rapid sinking and strong trench retreat (ISR); (iv) a two-stage mode, displaying backward-bent and subsequently inclined slabs, with late trench retreat (BIR). Transitions from regime (i) to (iii) occur with increasing subducting-plate age (i.e. buoyancy and strength). Regime (iv) develops for old (strong) subducting and overriding plates. We find that the interplay between trench motion and slab deformation at depth dictate the subduction style, both being controlled by slab strength, which is consistent with predictions from previous compositional subduction models. However, due to feedbacks between deformation, sinking rate, temperature and slab strength, the subducting-plate buoyancy, overriding-plate strength and upper-lower mantle viscosity jump are also important controls in thermo-mechanical subduction. For intermediate upper-lower mantle viscosity jumps (×30), our regimes reproduce the diverse range of seismically imaged slab morphologies.

Fluidity: A fully unstructured anisotropic adaptive mesh computational modeling framework for geodynamics

We present a new computational modeling framework, Fluidity, for application to a range of two- and three-dimensional geodynamic problems, with the focus here on mantle convection. The approach centers upon a finite element discretization on unstructured simplex meshes, which represent complex geometries in a straightforward manner. Throughout a simulation, the mesh is dynamically adapted to optimize the representation of evolving solution structures. The adaptive algorithm makes use of anisotropic measures of solution complexity, to vary resolution and allow long, thin elements to align with features such as boundary layers. The modeling framework presented differs from the majority of current mantle convection codes, which are typically based upon fixed structured grids. This necessitates a thorough and detailed validation, which is a focus of this paper. Benchmark comparisons are undertaken with a range of two- and three-dimensional, isoviscous and variable viscosity cases. In addition, model predictions are compared to experimental results. Such comparisons highlight not only the robustness and accuracy of Fluidity but also the advantages of anisotropic adaptive unstructured meshes, significantly reducing computational requirements when compared to a fixed mesh simulation.

Anisotropic mesh adaptivity for multi-scale ocean modelling

Research into the use of unstructured mesh methods in oceanography has been growing steadily over the past decade. The advantages of this approach for domain representation and non-uniform resolution are clear. However, a number of issues remain, in particular those related to the computational cost of models produced using unstructured mesh methods compared with their structured mesh counterparts. Mesh adaptivity represents an important means to improve the competitiveness of unstructured mesh models, where high resolution is only used when and where necessary. In this paper, an optimization-based approach to mesh adaptivity is described where emphasis is placed on capturing anisotropic solution characteristics. Comparisons are made between the results obtained with uniform isotropic resolution, isotropic adaptive resolution and fully anisotropic adaptive resolution.

A community benchmark for viscoplastic thermal convection in a 2‐D square box

Numerical simulations of thermal convection in the Earth’s mantle often employ a pseudoplastic rheology in order to mimic the plate-like behavior of the lithosphere. Yet the benchmark tests available in the literature are largely based on simple linear rheologies in which the viscosity is either assumed to be constant or weakly dependent on temperature. Here we present a suite of simple tests based on nonlinear rheologies featuring temperature, pressure, and strain rate-dependent viscosity. Eleven different codes based on the finite volume, finite element, or spectral methods have been used to run five benchmark cases leading to stagnant lid, mobile lid, and periodic convection in a 2-D square box. For two of these cases, we also show resolution tests from all contributing codes. In addition, we present a bifurcation analysis, describing the transition from a mobile lid regime to a periodic regime, and from a periodic regime to a stagnant lid regime, as a function of the yield stress. At a resolution of around 100 cells or elements in both vertical and horizontal directions, all codes reproduce the required diagnostic quantities with a discrepancy of at most

The Mantle Wedge's Transient 3-D Flow Regime and Thermal Structure

3% in the presence of both linear and nonlinear rheologies. Furthermore, they consistently predict the critical value of the yield stress at which the transition between different regimes occurs. As the most recent mantle convection codes can handle a number of different geometries within a single solution framework, this benchmark will also prove useful when validating viscoplastic thermal convection simula- tions in such geometries.

The concurrent emergence and causes of double volcanic hotspot tracks on the Pacific plate

Arc volcanism, volatile cycling, mineralization, and continental crust formation are likely regulated by the mantle wedges flow regime and thermal structure. Wedge flow is often assumed to follow a regular corner-flow pattern. However, studies that incorporate a hydrated rheology and thermal buoyancy predict internal small-scale-convection (SSC). Here, we systematically explore mantle-wedge dynamics in 3-D simulations. We find that longitudinal “Richter-rolls” of SSC (with trench-perpendicular axes) commonly occur if wedge hydration reduces viscosities to ≲1·1019 Pa s, although transient transverse rolls (with trench-parallel axes) can dominate at viscosities of ∼5·1018−1·1019 Pa s. Rolls below the arc and back arc differ. Subarc rolls have similar trench-parallel and trench-perpendicular dimensions of 100–150 km and evolve on a 1–5 Myr time-scale. Subback-arc instabilities, on the other hand, coalesce into elongated sheets, usually with a preferential trench-perpendicular alignment, display a wavelength of 150–400 km and vary on a 5–10 Myr time scale. The modulating influence of subback-arc ridges on the subarc system increases with stronger wedge hydration, higher subduction velocity, and thicker upper plates. We find that trench-parallel averages of wedge velocities and temperature are consistent with those predicted in 2-D models. However, lithospheric thinning through SSC is somewhat enhanced in 3-D, thus expanding hydrous melting regions and shifting dehydration boundaries. Subarc Richter-rolls generate time-dependent trench-parallel temperature variations of up to ∼150 K, which exceed the transient 50–100 K variations predicted in 2-D and may contribute to arc-volcano spacing and the variable seismic velocity structures imaged beneath some arcs.

TerraFERMA: The Transparent Finite Element Rapid Model Assembler for multiphysics problems in Earth sciences

Mantle plumes are buoyant upwellings of hot rock that transport heat from Earth’s core to its surface, generating anomalous regions of volcanism that are not directly associated with plate tectonic processes. The best-studied example is the Hawaiian–Emperor chain, but the emergence of two sub-parallel volcanic tracks along this chain, Loa and Kea, and the systematic geochemical differences between them have remained unexplained. Here we argue that the emergence of these tracks coincides with the appearance of other double volcanic tracks on the Pacific plate and a recent azimuthal change in the motion of the plate. We propose a three-part model that explains the evolution of Hawaiian double-track volcanism: first, mantle flow beneath the rapidly moving Pacific plate strongly tilts the Hawaiian plume and leads to lateral separation between high- and low-pressure melt source regions; second, the recent azimuthal change in Pacific plate motion exposes high- and low-pressure melt products as geographically distinct volcanoes, explaining the simultaneous emergence of double-track volcanism across the Pacific; and finally, secondary pyroxenite, which is formed as eclogite melt reacts with peridotite, dominates the low-pressure melt region beneath Loa-track volcanism, yielding the systematic geochemical differences observed between Loa- and Kea-type lavas. Our results imply that the formation of double-track volcanism is transitory and can be used to identify and place temporal bounds on plate-motion changes.

Reconciling mantle wedge thermal structure with arc lava thermobarometric determinations in oceanic subduction zones

We introduce and describe a new software infrastructure TerraFERMA, the\emph{Transparent Finite Element Rapid Model Assembler}, for the rapid and reproducible description and solution of coupled multiphysics problems. The design of TerraFERMA is driven by two computational needs in Earth sciences. The first is the need for increased flexibility in both problem description and solution strategies for coupled problems where small changes in model assumptions can lead to dramatic changes in physical behavior. The second is the need for software and models that are more transparent so that results can be verified, reproduced, and modified in a manner such that the best ideas in computation and Earth science can be more easily shared and reused. TerraFERMA leverages three advanced open-source libraries for scientific computation that provide high-level problem description (FEniCS), composable solvers for coupled multiphysics problems (PETSc), and an options handling system (SPuD) that allows the hierarchical management of all model options. TerraFERMA integrates these libraries into an interface that organizes the scientific and computational choices required in a model into a single options file from which a custom compiled application is generated and run. Because all models share the same infrastructure, models become more reusable and reproducible, while still permitting the individual researcher considerable latitude in model construction. TerraFERMA solves partial differential equations using the finite element method. It is particularly well suited for nonlinear problems with complex coupling between components. TerraFERMA is open-source and available at http://terraferma.github.io, which includes links to documentation and example input files. This article is protected by copyright. All rights reserved.

Compressible Flows on Adaptive and Unstrucured Meshes with FLUIDITY

Subduction zone mantle wedge temperatures impact plate interaction, melt generation, and chemical recycling. However, it has been challenging to reconcile geophysical and geochemical constraints on wedge thermal structure. Here we chemically determine the equilibration pressures and temperatures of primitive arc lavas from worldwide intra-oceanic subduction zones and compare them to kinematically driven thermal wedge models. We find that equilibration pressures are typically located in the lithosphere, starting just below the Moho, and spanning a wide depth range of ∼25 km. Equilibration temperatures are high for these depths, averaging ∼1300°C. We test for correlations with subduction parameters and find that equilibration pressures correlate with upper plate age, indicating overriding lithosphere thickness plays a role in magma equilibration. We suggest that most, if not all, thermobarometric pressure and temperature conditions reflect magmatic re-equilibration at a mechanical boundary, rather than reflecting the conditions of major melt generation. The magma re-equilibration conditions are difficult to reconcile, to a first order, with any of the conditions predicted by our dynamic models, with the exception of subduction zones with very young, thin upper plates. For most zones, a mechanism for substantially thinning the overriding plate is required. Most likely thinning is localised below the arc, as kinematic thinning above the wedge corner would lead to a hot forearc, incompatible with forearc surface heat flow and seismic properties. Localised sub-arc thermal erosion is consistent with seismic imaging and exhumed arc structures. Furthermore, such thermal erosion can serve as a weakness zone and affect subsequent plate evolution. This article is protected by copyright. All rights reserved.

Fluid flow in subduction zones: The role of solid rheology and compaction pressure

Fluidity is an open source, general purpose, multi‐phase CFD code capable of solving numerically the Navier‐Stokes and accompanying field equations on arbitrary unstructured finite element meshes in one, two and three dimensions. It uses a moving finite element/control volume method which allows arbitrary movement of the mesh in time dependent problems. It has a wide range of finite element/control volume element choices including mixed formulations. Here, continuous Galerkin (CG) and control volumes (CV) solutions of the compressible Navier‐Stokes (N‐S) equations are presented for the stratified tests cases of the rising thermal bubble and inertia gravity waves. Results show good agreement with previously published literature and novel result presented here is the ability to dynamically adapt the mesh to increase resolution in the region of interest, thus reducing the number of degrees of freedom in the problem without decreasing the accuracy of the result. Finally, results from the case of a fully three...