Search

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.

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 II the finite element techniques that can be used to predict thermal structure are discussed in an introductory fashion along with their verification and validation. This is followed by a comparison of the predicted thermal structure to geophysical, geochemical, and petrological observations.

Computational models of mantle convection must accurately represent curved boundaries and the associated boundary conditions of a 3-D spherical shell, bounded by Earth's surface and the core–mantle boundary. This is also true for comparable models in a simplified 2-D cylindrical geometry. It is of fundamental importance that the codes underlying these models are carefully verified prior to their application in a geodynamical context, for which comparisons against analytical solutions are an indispensable tool. However, analytical solutions for the Stokes equations in these geometries, based upon simple source terms that adhere to physically realistic boundary conditions, are often complex and difficult to derive. In this paper, we present the analytical solutions for a smooth polynomial source and a delta-function forcing, in combination with free-slip and zero-slip boundary conditions, for both 2-D cylindrical- and 3-D spherical-shell domains. We study the convergence of the Taylor–Hood (P2–P1) discretisation with respect to these solutions, within the finite element computational modelling framework Fluidity, and discuss an issue of suboptimal convergence in the presence of discontinuities. To facilitate the verification of numerical codes across the wider community, we provide a Python package, Assess, that evaluates the analytical solutions at arbitrary points of the domain.

The plate interface undergoes two transitions between seismogenic depths and subarc depths. A brittle-ductile transition at 20–50 km depth is followed by a transition to full viscous coupling to the overlying mantle wedge at ∼80 km depth. We review evidence for both transitions, focusing on heat-flow and seismic-attenuation constraints on the deeper transition. The intervening ductile shear zone likely weakens considerably as temperature increases, such that its rheology exerts a stronger control on subduction-zone thermal structure than does frictional shear heating. We evaluate its role through analytic approximations and two-dimensional finite-element models for both idealized subduction geometries and those resembling real subduction zones. We show that a temperature-buffering process exists in the shear zone that results in temperatures being tightly controlled by the rheological strength of that shear zone’s material for a wide range of shear-heating behaviors of the shallower brittle region. Higher temperatures result in weaker shear zones and hence less heat generation, so temperatures stop increasing and shear zones stop weakening. The net result for many rheologies are temperatures limited to ≤350–420 °C along the plate interface below the cold forearc of most subduction zones until the hot coupled mantle is approached. Very young incoming plates are the exception. This rheological buffering desensitizes subduction-zone thermal structure to many parameters and may help explain the global constancy of the 80 km coupling limit. We recalculate water fluxes to the forearc wedge and deep mantle and find that shear heating has little effect on global water circulation.

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

Subduction is a key driver of plate tectonics on Earth1. A range of observations indicate that the termination of subduction leads to diverse and unexplained tectonic and geological activity, including anomalous magmatism, exhumation and topographic subsidence, followed by rapid uplift. However, the mechanism driving this complex surface response remains enigmatic. A prime example of recent subduction termination can be found in northern Borneo (Malaysia), where subduction ceased in the late Miocene2 and was followed by a puzzling tectonic response3,4,5,6,7,8 that cannot be reconciled with our current understanding of post-subduction tectonics. Here, we use new passive-seismic data to image, in unprecedented detail, a sub-vertical lithospheric drip that developed as a Rayleigh-Taylor gravitational instability9 from the root of a volcanic arc. We use thermo-mechanical simulations to reconcile these images with time-dependent dynamical processes within the crust and underlying mantle following subduction termination. Our model predictions illustrate how significant extension from a downwelling lithospheric drip can thin the crust in an adjacent orogenic belt, facilitating lower crustal melting and possible exhumation of subcontinental material. Our study provides a new paradigm for core-complex formation in other areas of recent subduction termination.

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

Melt transport beneath the lithosphere is elusive. With a distinct viscosity and density from the surrounding mantle, magmatic melt moves on a different time scale as the surrounding mantle. To resolve the temporal scale necessary to accurately capture melt transport in the mantle, the model simulations become numerically expensive quickly. Recent computational advances make possible two-phase numerical explorations to understand magma transport in the mantle. We review results from a suite of two-phase models applied to the mid-ocean ridges, where we varied 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). Here, we highlight the importance of viscosities for the melt focusing mechanisms. In addition, magmatic porosity waves that are a natural consequence of these two-phase flow formulations. We show that these waves could explain long-period temporal variations in the seafloor bathymetry at the Southeast Indian Ridge.

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.

Tsunami generated by submarine slides are arguably an under-considered risk in comparison to earthquake-generated tsunami. Numerical simulations of submarine slide-generated waves can be used to identify the important factors in determining wave characteristics. Here we use Fluidity, an open source finite element code, to simulate waves generated by deformable submarine slides. Fluidity uses flexible unstructured meshes combined with adaptivity which alters the mesh topology and resolution based on the simulation state, focussing or reducing resolution, when and where it is required. Fluidity also allows a number of different numerical approaches to be taken to simulate submarine slide deformation, free-surface representation, and wave generation within the same numerical framework. In this work we use a multi-material approach, considering either two materials (slide and water with a free surface) or three materials (slide, water and air), as well as a sediment model (sediment, water and free surface) approach. In all cases the slide is treated as a viscous fluid. Our results are shown to be consistent with laboratory experiments using a deformable submarine slide, and demonstrate good agreement when compared with other numerical models. The three different approaches for simulating submarine slide dynamics and tsunami wave generation produce similar waveforms and slide deformation geometries. However, each has its own merits depending on the application. Mesh adaptivity is shown to be able to reduce the computational cost without compromising the accuracy of results.

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.

International audience; The migration pathways of hydrous fluids in the mantle wedge are influenced by the compaction of the porous mantle matrix, which depends on the matrix permeability, fluid viscosity, and fluid density. Experimental studies show that when fluids are interconnected, the permeability depends on mineral grain size and porosity, the latter of which depends on the amount of fluids introduced into the system (fluid influx). Here, we investigate the role of fluid influx, fluid viscosity, and fluid density in controlling fluid migration in the mantle wedge, using a 2-D numerical model accounting for the effects of grain-size variation and matrix compaction. Our models predict that fluid influx and fluid viscosity are key controls on fluid pathways, while fluid density plays a secondary role. Temperature dependence of fluid viscosity promotes downdip drag of fluids at the base of the forearc mantle toward the subarc region. High fluid influx at postarc depths promotes updip flow near the base of the mantle wedge, guiding the fluids arcward. The model that is applied to northern Cascadia predicts upward fluid migration focused beneath the arc but cannot explain high electrical conductivity observed slightly west of the upward fluid migration. We estimate the amount of hydrous melt that can be produced in the mantle wedge using calculated fluid distributions. Up to a few percent partial melting is predicted in a relatively small region in the core part of the subarc mantle wedge in most subduction settings, including northern Cascadia, and beneath the backarc in old-slab subduction zones.

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.

We consider multiphysics applications from algorithmic and architectural perspectives, where “algorithmic” includes both mathematical analysis and computational complexity, and “architectural” includes both software and hardware environments. Many diverse multiphysics applications can be reduced, en route to their computational simulation, to a common algebraic coupling paradigm. Mathematical analysis of multiphysics coupling in this form is not always practical for realistic applications, but model problems representative of applications discussed herein can provide insight. A variety of software frameworks for multiphysics applications have been constructed and refined within disciplinary communities and executed on leading-edge computer systems. We examine several of these, expose some commonalities among them, and attempt to extrapolate best practices to future systems. From our study, we summarize challenges and forecast opportunities.

The emergence of geographically and geochemically distinct double volcanic chains on the Pacific plate coincides with a recent azimuthal change in the motion of the plate. Volcanic hotspots in regions that are not associated with plate tectonics are widely thought to be surface expressions of thermochemical upwellings from the deep Earth. But the emergence of the double Loa and Kea tracks of volcanism along the Hawaiian–Emperor hotspot island chain, and the systematic geochemical differences between them, remain enigmatic. Timothy Jones et al. show that the emergence of these and other double volcanic tracks on the Pacific plate coincides with a recent change in Pacific plate motion. They propose that 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, producing the geographically and geochemically distinct volcanic chains. 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 chain1, Loa and Kea, and the systematic geochemical differences between them2,3 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 peridotite4, dominates the low-pressure melt region beneath Loa-track volcanism, yielding the systematic geochemical differences observed between Loa- and Kea-type lavas3,5,6,7,8,9. 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.

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.

Verification is a process crucially important for the final users of a computational model: code is useless if its results cannot be relied upon. Typically, verification is seen as a discrete event, performed once and for all after development is complete. However, this does not reflect the reality that many geoscientific codes undergo continuous development of the mathematical model, discretisation and software implementation. Therefore, we advocate that in such cases verification must be continuous and happen in parallel with development: the desirability of their automation follows immediately. This paper discusses a framework for automated continuous verification of wide applicability to any kind of numerical simulation. It also documents a range of test cases to show the possibilities of the framework.

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.

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.

The interfaces by which users specify the scenarios to be simulated by scientific computer models are frequently primitive, under-documented and ad-hoc text files which make using the model in question difficult and error-prone and significantly increase the development cost of the model. In this paper, we present a model-independent system, Spud, which formalises the specification of model input formats in terms of formal grammars. This is combined with an automated graphical user interface which guides users to create valid model inputs based on the grammar provided, and a generic options reading module, libspud, which minimises the development cost of adding model options. Together, this provides a user friendly, well documented, self validating user interface which is applicable to a wide range of scientific models and which minimises the developer input required to maintain and extend the model interface.