Search

Earth’s surface topography is a direct physical expression of our planet’s dynamics. Most is isostatic, controlled by thickness and density variations within the crust and lithosphere, but a substantial proportion arises from forces exerted by underlying mantle convection. This dynamic topography directly connects the evolution of surface environments to Earth’s deep interior, but predictions from mantle flow simulations are often inconsistent with inferences from the geological record, with little consensus about its spatial pattern, wavelength and amplitude. Here, we demonstrate that previous comparisons between predictive models and observational constraints have been biased by subjective choices. Using measurements of residual topography beneath the oceans, and a hierarchical Bayesian approach to performing spherical harmonic analyses, we generate a robust estimate of Earth’s oceanic residual topography power spectrum. This indicates water-loaded power of 0.5 ± 0.35 km2 and peak amplitudes of up to ~0.8 ± 0.1 km at long wavelengths (~104 km), decreasing by roughly one order of magnitude at shorter wavelengths (~103 km). We show that geodynamical simulations can be reconciled with observational constraints only if they incorporate lithospheric structure and its impact on mantle flow. This demonstrates that both deep (long-wavelength) and shallow (shorter-wavelength) processes are crucial, and implies that dynamic topography is intimately connected to the structure and evolution of Earth’s lithosphere. The structure of the lithosphere and its impact on mantle flow significantly influence the impact of Earth’s interior dynamics on surface topography, suggest statistical analyses of Earth’s topography.

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 (S...

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 (FE) 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 FE 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.

Mineral grain size in the mantle affects fluid migration by controlling mantle permeability; the smaller the grain size, the less permeable the mantle is. Mantle shear viscosity also affects fluid migration by controlling compaction pressure; high mantle shear viscosity can act as a barrier to fluid flow. Here, we investigate for the first time their combined effects on fluid migration in the mantle wedge of subduction zones over ranges of subduction parameters and patterns of fluid influx using a 2-D numerical fluid migration model. Our results show that fluids introduced into the mantle wedge beneath the forearc are first dragged down-dip by the mantle flow due to small grain size (< 1 mm) and high mantle shear viscosity that develop along the base of the mantle wedge. Increasing grain size with depth allows upward fluid migration out of the high shear viscosity layer at sub-arc depths. Fluid introduced into the mantle wedge at post-arc depths migrate upwards due to relatively large grain size in the deep mantle wedge, forming secondary fluid pathways behind the arc. Fluids that reach the shallow part of the mantle wedge spread trench-wards due to the combined effect of high mantle shear viscosity and advection by the inflowing mantle and eventually pond at 55–65 km depths. These results show that grain size and mantle shear viscosity together play an important role in focusing fluids beneath the arc.

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.

At mid-ocean ridges, oceanic crust is emplaced in a narrow neovolcanic region on the seafloor, whereas basaltic melt that forms this oceanic crust is generated in a wide region beneath as suggested by a few geophysical surveys. The combined observations suggest that melt generated in a wide region at depths has to be transported horizontally to a small region at the surface. We present results from a suite of two-phase models applied to the mid-ocean ridges, varying half-spreading rate and intrinsic mantle permeability using new openly available models, with the goal of understanding melt focusing beneath mid-ocean ridges and its relevance to the lithosphere-asthenosphere boundary (LAB). Three distinct melt focusing mechanisms are recognized in these models: 1) melting pressure focusing, 2) decompaction layers and 3) ridge suction, of which the first two play dominant roles in focusing melt. All three of these mechanisms exist in the fundamental two phase flow formulation but the manifestation depends largely on the choice of rheological model. The models also show that regardless of spreading rates, the amount of melt and melt transport patterns are sensitive to changes in intrinsic permeability, K0. In these models, the LAB is delineated by the melt-rich decompaction layers, which are essentially defined by the temperature dependent rheological and freezing boundaries. Geophysical observations place the LAB at a steeper incline as compared to the gentler profile suggested by most of our models. The models suggest that one way to reconcile this discrepancy is to have stronger melting pressure focusing mechanism as it is the only mechanism in these models that can focus melt before reaching the typical model thermal LAB. The apparent lack of observable decompaction layers in the geophysical observations hints at the possibility that melting pressure focusing could be significant. These models help improve our understanding of melt focusing beneath mid-ocean ridges and could provide new constraints for mantle rheology and permeability.

The paleomagnetic record is central to our understanding of the history of the Earth. The orientation and intensity of magnetic minerals preserved in ancient rocks indicate the geodynamo has been alive since at least the Archean and possibly the Hadean. A paleomagnetic signature of the initial solidification of the inner core, arguably the singular most important event in core history, however, has remained elusive. In pursuit of this signature we investigate the assumption that the field is a geocentric axial dipole (GAD) over long time scales. We study a suite of numerical dynamo simulations from a paleomagnetic perspective to explore how long the field should be time-averaged to obtain stable paleomagnetic pole directions and intensities. We find that running averages over 20 − 40 kyr are needed to obtain stable paleomagnetic poles with α95

It has been proposed that the spatial variations recorded in the geochemistry of hotspot lavas, such as the bilateral asymmetry recorded at Hawaii, can be directly mapped as the heterogeneous structure and composition of their deep-mantle source. This would imply that source-region heterogeneities are transported into, and preserved within, a plume conduit, as the plume rises from the deep-mantle to Earth's surface. Previous laboratory and numerical studies, which neglect density and rheological variations between different chemical components, support this view. However, in this paper, we demonstrate that this interpretation cannot be extended to distinct chemical domains that differ from surrounding mantle in their density and viscosity. By numerically simulating thermo-chemical mantle plumes across a broad parameter space, in 2-D and 3-D, we identify two conduit structures: (i) bilaterally asymmetric conduits, which occur exclusively for cases where the chemical effect on buoyancy is negligible, in which the spatial distribution of deep-mantle heterogeneities is preserved during plume ascent; and (ii) concentric conduits, which occur for all other cases, with dense material preferentially sampled within the conduit's centre. In the latter regime, the spatial distribution of geochemical domains in the lowermost mantle is not preserved during plume ascent. Our results imply that the heterogeneous structure and composition of Earth's lowermost mantle can only be mapped from geochemical observations at Earth's surface if chemical heterogeneity is a passive component of lowermost mantle dynamics (i.e. its effect on density is outweighed by, or is secondary to, the effect of temperature). The implications of our results for: (i) why oceanic crust should be the prevalent component of ocean island basalts; and (ii) how we interpret the geochemical evolution of Earth's deep-mantle are also discussed.

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 $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.

Arc volcanoes tend to occur at locations where the slab is at approximately 100 km depth but most models of fluid production from the downgoing slab suggest fluids are released over a wide range of depths. Reconciling the models with the observations suggests that focusing of slab-produced fluids is necessary if flux-melting is a primary mechanism for the production of arc magmas. This paper investigates one possible mechanism for inducing focusing of fluid flow toward the sub-arc mantle. Through a series of simplified models we explore the role of compaction pressure gradients in modifying fluid flow. These gradients are produced by variations in fluid flux interacting with the permeability and viscosity structure of the solid mantle. When these gradients are neglected, high-permeability systems are dominated by buoyancy and fluid flow is primarily vertical. However, when included, compaction pressure terms have three principal effects: (i) enhancement of upslope flow within high-permeability layers in the slab produced by local dehydration reactions, (ii) deflection of fluids along the sloping rheologically strengthening region in the upper thermal boundary layer, and (iii) production of non-linear porosity waves that locally modulate the flow of fluids and can allow significant transient accumulation of fluids. We demonstrate significant localization of fluid flux toward the sub-arc region due to the permeability and solid viscosity structure. We also estimate the amount of melting expected among the different models and show that models with compaction pressure could produce ∼10% flux melting, whereas distributed fluid flow produces ≲1% in most cases.

Small-scale experiments of volcanic ash particle settling in water have demonstrated that ash particles can either settle slowly and individually, or rapidly and collectively as a gravitationally unstable ash-laden plume. This has important implications for the emplacement of tephra deposits on the seabed. Numerical modelling has the potential to extend the results of laboratory experiments to larger scales and explore the conditions under which plumes may form and persist, but many existing models are computationally restricted by the fixed mesh approaches that they employ. In contrast, this paper presents a new multiphase flow model that uses an adaptive unstructured mesh approach. As a simulation progresses, the mesh is optimized to focus numerical resolution in areas important to the dynamics and decrease it where it is not needed, thereby potentially reducing computational requirements. Model verification is performed using the method of manufactured solutions, which shows the correct solution convergence rates. Model validation and application considers 2-D simulations of plume formation in a water tank which replicate published laboratory experiments. The numerically predicted settling velocities for both individual particles and plumes, as well as instability behaviour, agree well with experimental data and observations. Plume settling is clearly hindered by the presence of a salinity gradient, and its influence must therefore be taken into account when considering particles in bodies of saline water. Furthermore, individual particles settle in the laminar flow regime while plume settling is shown (by plume Reynolds numbers greater than unity) to be in the turbulent flow regime, which has a significant impact on entrainment and settling rates. Mesh adaptivity maintains solution accuracy while providing a substantial reduction in computational requirements when compared to the same simulation performed using a fixed mesh, highlighting the benefits of an adaptive unstructured mesh approach.

Identifying the dominant controls on Earth’s surface topography is of critical importance to understanding both the short- and long-term evolution of geological processes and past- and present-day dynamics of Earth’s coupled mantle–lithosphere system. The ability to simulate a stress free — or a so-called ‘free surface’ — boundary condition is required to examine such processes via numerical models. However, at present, geodynamical models incorporating a free surface are limited, as most underlying free surface algorithms place severe restrictions on the computational timestep. Consequently, the simulations are often intractable. In this study, we introduce a new approach for incorporating a free surface within geodynamical models: an algorithm, in which free surface elevation is treated as an independent variable and is solved for in conjunction with the momentum and continuity equation, using implicit time integration. We demonstrate that the method is straightforward to implement in existing models and, using a series of analytical and benchmark comparisons, we show that it does not suffer from the timestep constraints of previous schemes. Furthermore, the scheme can be made second order accurate in time, at no additional cost. The method therefore dramatically improves the computational efficiency of geodynamical simulations including a free surface, whilst maintaining solution accuracy.

International audience; Transition zone slab deformation influences Earth's 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 subduct-ing 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 dictates 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 (330), our regimes reproduce the diverse range of seismically imaged slab morphologies.

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.