The thermal structure of subduction zones is fundamental to our understanding of physical and chemical processes that occur at active convergent plate margins. These include magma generation and related arc volcanism, shallow and deep seismicity, and metamorphic reactions that can release fluids. Computational models can predict the thermal structure to great numerical precision when models are fully described but this does not guarantee accuracy or applicability. In a pair of companion papers the construction of thermal subduction zone models, their use in subduction zone studies, and their link to geophysical and geochemical observations is explored. In part I the motivation to understand the thermal structure is presented based on experimental and observational studies. This is followed by a description of a selection of thermal models for the Japanese subduction zones.

A document by Peter van Keken. Click on the document to view its contents.

A document by Peter van Keken. Click on the document to view its contents.

In prior work we found that precise approximation of the continuity constraint is crucial for accurate propagation of tracer data when advected through a background incompressible velocity field (Sime et al., 2021). Here we extend this investigation to compressible flows using the anelastic liquid approximation (ALA) and address four related issues: 1. exact conservation of tracer discretized fields through a background compressible velocity; 2. exact mass conservation; 3. addition and removal of tracers without affecting (exact) conservation to preserve a consistent number of tracers per cell; and 4. the diffusion of tracer data, for example, as induced by thermal or chemical effects. In this process we also present an abstract formulation of the interior penalty hybrid discontinuous Galerkin (HDG) finite element formulation for diffusion problems, and apply it to the advection-diffusion and compressible Stokes systems. Finally we present numerical experiments exhibiting the HDG compressible Stokes momentum formulation’s superconvergent compressibility approximation and reproduce community numerical benchmarks from the literature for the ALA.

The evolution of mantle composition can be viewed as process of destruction whereby the initial chemical state is overprinted and reworked with time. Analyses of ocean island basalts reveals that some portion of the mantle has survived this process, retaining a chemically ‘primitive’ signature. A question that remains is how this primitive signature has survived four and half billion years of vigorous convection. We hypothesize that some of Earth’s primitive mantle is buried within a slab graveyard at the core-mantle boundary. We explore this possibility using high‐resolution finite element models of mantle convection, in which oceanic lithosphere is produced at zones of plate spreading and subducted at zones of plate convergence. Upon subduction, dense oceanic crust sinks to the base of the mantle and gradually accumulates to form broad and robust thermochemical piles. Sinking oceanic crust entrains the surrounding mantle whose composition is predominantly primitive early in the model’s evolution. As a result, thermochemical piles are initially supplied with relatively high concentrations of primitive material –summing up to ~30% their total mass. The dense oceanic crust that dominates the piles resists efficient mixing and preserves the primitive material that it is intermingled with. The significance of this process is shown to be proportional the rate of mantle processing through time and the excess density of oceanic crust at mantle pressures and temperatures. Unlike existing theories for the survival of Earth’s primitive mantle, this one does not require the early Earth to have anomalously high density or large scale viscosity contrasts.

Tracer methods are widespread in computational geodynamics for modeling the advection of chemical data. However, they present certain numerical challenges, especially when used over long periods of simulation time. We address two of these in this work: the necessity for mass conservation of chemical composition fields and the need for the velocity field to be pointwise divergence free to avoid gaps in tracer coverage. We do this by implementing the hybrid discontinuous Galerkin (HDG) finite element method combined with a mass-conserving constrained projection of the tracer data. To demonstrate the efficacy of this system we compare it to other common finite element formulations of the Stokes system and projections of the chemical composition. We provide a reference of the numerical properties and error convergence rates which should be observed by using these various discretization schemes. This serves as a tool for verification of existing or new implementations. We summarize these data in a reproduction of a published Rayleigh–Taylor instability benchmark, demonstrating the importance of careful choices of appropriate and compatible discretization methods for all aspects of geodynamics simulations.

Seismic wave amplitudes have tremendous sensitivity to subduction structure; however, they are affected by attenuation, scattering and focusing, and have therefore been sparsely used compared with traveltimes. We measure and model teleseismic body wave amplitudes recorded at a dense broadband array in the Washington Cascades (iMUSH). These data show anomalous amplitude variations with complex azimuthal dependence at the low frequency of 0.05 Hz, accompanied by significant multipathing. We demonstrate using spectral-element numerical simulations that focusing of the teleseismic wavefield by the Juan de Fuca slab is responsible for some of the amplitude anomalies. The focusing effects can contaminate the apparent differential attenuation measurements and produce at least 20% of the inferred attenuation signal. The focusing results in complex azimuthal patterns that produce different phase and amplitude variations than does intrinsic attenuation, which should allow separation of elastic (focusing) and anelastic effects. Our results indicate that the amplitudes are sensitive to the subducting slab geometry and subduction structure, and can be used to refine seismic images. Ubiquitous focusing effects are observed along the arc, suggesting a continuous Juan de Fuca slab from Canada to northern California.

England and May (2021) present a model comparison for the forearc thermal structure in subduction zones that employs a finite element modeling technique and a new variation on approximate equations. We have some comments on various claims made in this published paper.