Mantle convection modeling with viscoelastic/brittle lithosphere: Numerical methodology and plate tectonic modeling

The earth's tectonic plates are strong, viscoelastic shells which make up the outermost part of a thermally convecting, predominantly viscous layer. Brittle failure of the lithosphere occurs when stresses are high. In order to build a realistic simulation of the planet's evolution, the complete viscoelastic/brittle convection system needs to be considered. A particle-in-cell finite element method is demonstrated which can simulate very large deformation viscoelasticity with a strain-dependent yield stress. This is applied to a plate-deformation problem. Numerical accuracy is demonstrated relative to analytic benchmarks, and the characteristics of the method are discussed.

Longevity and stability of cratonic lithosphere: Insights from numerical simulations of coupled mantle convection and continental tectonics

[1] The physical conditions required to provide for the tectonic stability of cratonic crust and for the relative longevity of deep cratonic lithosphere within a dynamic, convecting mantle are explored through a suite of numerical simulations. The simulations allow chemically distinct continents to reside within the upper thermal boundary layer of a thermally convecting mantle layer. A rheologic formulation, which models both brittle and ductile behavior, is incorporated to allow for plate-like behavior and the associated subduction of oceanic lithosphere. Several mechanisms that may stabilize cratons are considered. The two most often invoked mechanisms, chemical buoyancy and/or high viscosity of cratonic root material, are found to be relatively ineffective if cratons come into contact with subduction zones. High root viscosity can provide for stability and longevity but only within a thick root limit in which the thickness of chemically distinct, high-viscosity cratonic lithosphere exceeds the thickness of old oceanic lithosphere by at least a factor of 2. This end-member implies a very thick mechanical lithosphere for cratons. A high brittle yield stress for cratonic lithosphere as a whole, relative to oceanic lithosphere, is found to be an effective and robust means for providing stability and lithospheric longevity. This mode does not require exceedingly deep strength within cratons. A high yield stress for only the crustal or mantle component of the cratonic lithosphere is found to be less effective as detachment zones can then form at the crust-mantle interface which decreases the longevity potential of cratonic roots. The degree of yield stress variations between cratonic and oceanic lithosphere required for stability and longevity can be decreased if cratons are bordered by continental lithosphere that has a relatively low yield stress, i.e., mobile belts. Simulations that combine all the mechanisms can lead to crustal stability and deep root longevity for model cratons over several mantle overturn times, but the dominant stabilizing factor remains a relatively high brittle yield stress for cratonic lithosphere.

Elasticity, yielding and episodicity in simple models of mantle convection

We explore the implications of refinements in the mechanical description of planetary constituents on the convection modes predicted by finite-element simulations. The refinements consist in the inclusion of incremental elasticity, plasticity (yielding) and multiple simultaneous creep mechanisms in addition to the usual visco-plastic models employed in the context of unified plate-mantle models. The main emphasis of this paper rests on the constitutive and computational formulation of the model. We apply a consistent incremental formulation of the non-linear governing equations avoiding the computationally expensive iterations that are otherwise necessary to handle the onset of plastic yield. In connection with episodic convection simulations, we point out the strong dependency of the results on the choice of the initial temperature distribution. Our results also indicate that the inclusion of elasticity in the constitutive relationships lowers the mechanical energy associated with subduction events.

Anisotropy model for mantle convection

The paper presents a new theory for modeling flow in anisotropic, viscous rock. This theory has originally been developed for the simulation of large deformation processes including folding and kinking in multi-layered visco-elastic rock. The orientation of slip planes in the context of crystallographic slip is determined by the normal vector, the so-called director of these surfaces. The model is applied to simulate anisotropic natural mantle convection. We compare the evolution of the director and approximately steady states of isotropic and anisotropic convection. The isotropic case has a simple steady state solution, whereas the orthotropic convection model produces a continuously evolving patterning in tile core of the convection cell which makes only a near-steady condition possible, in which the thermal boundary layer appears to be well aligned with the flow and hence as observed in seismic tomomgraphy strong anistropic.

Anisotropic convection model for the earth's mantle

The paper presents a theory for modeling flow in anisotropic, viscous rock. This theory has originally been developed for the simulation of large deformation processes including the folding and kinking of multi-layered visco-elastic rock (Muhlhaus et al. [1,2]). The orientation of slip planes in the context of crystallographic slip is determined by the normal vector - the director - of these surfaces. The model is applied to simulate anisotropic mantle convection. We compare the evolution of flow patterns, Nusselt number and director orientations for isotropic and anisotropic rheologies. In the simulations we utilize two different finite element methodologies: The Lagrangian Integration Point Method Moresi et al [8] and an Eulerian formulation, which we implemented into the finite element based pde solver Fastflo (www.cmis.csiro.au/Fastflo/). The reason for utilizing two different finite element codes was firstly to study the influence of an anisotropic power law rheology which currently is not implemented into the Lagrangian Integration point scheme [8] and secondly to study the numerical performance of Eulerian (Fastflo)- and Lagrangian integration schemes [8]. It turned out that whereas in the Lagrangian method the Nusselt number vs time plot reached only a quasi steady state where the Nusselt number oscillates around a steady state value the Eulerian scheme reaches exact steady states and produces a high degree of alignment (director orientation locally orthogonal to velocity vector almost everywhere in the computational domain). In the simulations emergent anisotropy was strongest in terms of modulus contrast in the up and down-welling plumes. Mechanisms for anisotropic material behavior in the mantle dynamics context are discussed by Christensen [3]. The dominant mineral phases in the mantle generally do not exhibit strong elastic anisotropy but they still may be oriented by the convective flow. Thus viscous anisotropy (the main focus of this paper) may or may not correlate with elastic or seismic anisotropy.

Stirring in 3-d spherical models of convection in the Earth's mantle

On a global scale basalts from mid-ocean ridges are strikingly more homogeneous than basalts from intraplate volcanism. The observed geochemical heterogeneity argues strongly for the existence of distinct reservoirs in the Earth's mantle. It is an unresolved problem of Geodynamics as to how these findings can be reconciled with large-scale convection. We review observational constraints, and investigate stirring properties of numerical models of mantle convection. Conditions in the early Earth may have supported layered convection with rapid stirring in the upper layers. Material that has been altered near the surface is transported downwards by small-scale convection. Thereby a layer of homogeneous depleted material develops above pristine mantle. As the mantle cools over Earth history, the effects leading to layering become reduced and models show the large-scale convection favoured for the Earth today. Laterally averaged, the upper mantle below the lithosphere is least affected by material that has experienced near-surface differentiation. The geochemical signature obtained during the previous episode of small-scale convection may be preserved there for the longest time. Additionally, stirring is less effective in the high viscosity layer of the central lower mantle [1, 2], supporting the survival of medium-scale heterogeneities there. These models are the first, using 3-d spherical geometry and mostly Earth-like parameters, to address the suggested change of convective style. Although the models are still far from reproducing our planet, we find that proposal might be helpful towards reconciling geochemical and geophysical constraints.

Thermal effects in the evolution of initially layered mantle material

A simplified model for anisotropic mantle convection based on a novel class of rheologies, originally developed for folding instabilities in multilayered rock (MUHLHAUS et al., 2002), is extended ¨ through the introduction of a thermal anisotropy dependent on the local layering. To examine the effect of the thermal anisotropy on the evolution of mantle material, a parallel implementation of this model was undertaken using the Escript modelling toolkit and the Finley finite-element computational kernel (DAVIES et al., 2004). For the cases studied, there appears too little if any effect. For comparative purposes, the effects of anisotropic shear viscosity and the introduced thermal anisotropy are also presented. These results contribute to the characterization of viscous anisotropic mantle convection subject to variation in thermal conductivities and shear viscosities.

Core formation and mantle differentiation on Mars

Geochemical investigation of Martian meteorites (SNC meteorites) yields important constraints on the chemical and geodynamical evolution of Mars. These samples may not be representative of the whole of Mars; however, they provide constraints on the early differentiation processes on Mars. The bulk composition of Martian samples implies the presence of a metallic core that formed concurrently as the planet accreted. The strong depletion of highly siderophile elements in the Martian mantle is only possible if Mars had a large scale magma ocean early in its history allowing efficient separation of a metallic melt from molten silicate. The solidification of the magma ocean created chemical heterogeneities whose ancient origin is manifested in the heterogeneous 142Nd and 182W abundances observed in different meteorite groups derived from Mars. The isotope anomalies measured in SNC meteorites imply major chemical fractionation within the Martian mantle during the life time of the short-lived isotopes 146Sm and 182Hf. The Hf-W data are consistent with very rapid accretion of Mars within a few million years or, alternatively, a more protracted accretion history involving several large impacts and incomplete metal-silicate equilibration during core formation. In contrast to Earth early-formed chemical heterogeneities are still preserved on Mars, albeit slightly modified by mixing processes. The preservation of such ancient chemical differences is only possible if Mars did not undergo efficient whole mantle convection or vigorous plate tectonic style processes after the first few tens of millions of years of its history.

An anisotropic viscous representation of Mohr-Coulomb failure for use is modelling coupled mantle-continent dynamics

In mantle convection models it has become common to make use of a modified (pressure sensitive, Boussinesq) von Mises yield criterion to limit the maximum stress the lithosphere can support. This approach allows the viscous, cool thermal boundary layer to deform in a relatively plate-like mode even in a fully Eulerian representation. In large-scale models with embedded continental crust where the mobile boundary layer represents the oceanic lithosphere, the von Mises yield criterion for the oceans ensures that the continents experience a realistic broad-scale stress regime. In detailed models of crustal deformation it is, however, more appropriate to choose a Mohr-Coulomb yield criterion based upon the idea that frictional slip occurs on whichever one of many randomly oriented planes happens to be favorably oriented with respect to the stress field. As coupled crust/mantle models become more sophisticated it is important to be able to use whichever failure model is appropriate to a given part of the system. We have therefore developed a way to represent Mohr-Coulomb failure within a code which is suited to mantle convection problems coupled to large-scale crustal deformation. Our approach uses an orthotropic viscous rheology (a different viscosity for pure shear to that for simple shear) to define a prefered plane for slip to occur given the local stress field. The simple-shear viscosity and the deformation can then be iterated to ensure that the yield criterion is always satisfied. We again assume the Boussinesq approximation - neglecting any effect of dilatancy on the stress field. An additional criterion is required to ensure that deformation occurs along the plane aligned with maximum shear strain-rate rather than the perpendicular plane which is formally equivalent in any symmetric formulation. It is also important to allow strain-weakening of the material. The material should remember both the accumulated failure history and the direction of failure. We have included this capacity in a Lagrangian-Integration-point finite element code and will show a number of examples of extension and compression of a crustal block with a Mohr-Coulomb failure criterion, and comparisons between mantle convection models using the von Mises versus the Mohr-Coulomb yield criteria. The formulation itself is general and applies to 2D and 3D problems, although it is somewhat more complicated to identify the slip plane in 3D.

Anomalous compressibility of ferropericlase throughout the iron spin cross-over

The thermoelastic properties of ferropericlase Mg(1-x)Fe(x)O (x = 0.1875) throughout the iron high-to-low spin cross-over have been investigated by first principles at Earth`s lower mantle conditions. This cross-over has important consequences for elasticity such as an anomalous bulk modulus (K(S)) reduction. At room temperature the anomaly is somewhat sharp in pressure but broadens with increasing temperature. Along a typical geotherm it occurs across most of the lower mantle with a more significant K(S) reduction at approximate to 1,400-1,600 km depth. This anomaly might also cause a reduction in the effective activation energy for diffusion creep and lead to a viscosity minimum in the mid-lower mantle, in apparent agreement with results from inversion of data related with mantle convection and postglacial rebound.