Simulations Of Mantle Convection Research Articles

3D spherical models of Martian mantle convection constrained by melting history

While most of Tharsis rise was in place by end of the Noachian period, at least one volcano on Tharsis swell (Arsia Mons) has been active within the last 2 Ma. This places an important constraint on mantle convection and on the thermal evolution of Mars. The existence of recent volcanism on Mars implies that adiabatic decompression melting and, hence, upwelling convective flow in the mantle remains important on Mars at present. The thermal history on Mars can be constrained by the history of melt production, specifically generating sufficient melt in the first billion years of the planets history to produce Tharsis rise as well as present day melt to explain recent volcanism. In this work, mantle convection simulations were performed using finite element code CitcomS in a 3D sphere starting from a uniformly hot mantle and integrating forward in time for the age of the solar system. We implement constant and decaying radioactive heat sources; and vary the partitioning of heat sources between the crust and mantle, and consider decreasing core–mantle boundary temperature and latent heat of melting. The constant heat source calculations produce sufficient melt to create Tharsis early in Martian history and continue to produce significant melt to the present. Calculations with decaying radioactive heat sources generate excessive melt in the past, except when all the radiogenic elements are in the crust, and none produce melt after 2 Gyr. Producing a degree-1 or degree-2 structure may not be pivotal to explain the Tharsis rise: we present multi-plume models where not every plume produces melt. The Rayleigh number controls the timing of the first peak of volcanism while late-stage volcanism is controlled more by internal mantle heating. Decreasing the Rayleigh number increases the lithosphere thickness (i.e., depth), and increasing lithosphere thickness increases the mean mantle temperature. Increasing pressure reduces melt production while increasing temperature increases melt production; hence predicting melt production from convection parameters is not straightforward. Generating enough melt in the mantle to create Tharsis early on and also to explain recent volcanism may require other mechanisms such as small-scale convection or lowering the thermal conductivity of the crust.

The role of harzburgite layers in the morphology of subducting plates and the behavior of oceanic crustal layers

AbstractNumerical simulations of mantle convection with semidynamic plate subduction in 3‐D geometry is performed to investigate the role of harzburgite layers in the morphology of subducting plates and the behavior of oceanic crustal layers. The results show that a “buckled” crustal layer is observed under the “heel” of the stagnant slab that begins to penetrate into the lower mantle, regardless of the magnitude of the viscosity of the harzburgite layer, ∆ηhrz, when the factor of viscosity increase at the boundary of the upper and lower mantle, ∆ηlm, is larger than 60–100. As ∆ηhrz increases, the curvature of buckling is larger. When ∆ηlm and ∆ηhrz are larger, the volumes of crustal and harzburgite materials trapped in the mantle transition zone (MTZ) are also larger, although almost all of the materials penetrate into the lower mantle. These materials are trapped in the MTZ for over tens of millions of years.

Mantle convection modeling of the supercontinent cycle: Introversion, extroversion, or a combination?

The periodic assembly and dispersal of continental fragments, referred to as the supercontinent cycle, bear close relation to the evolution of mantle convection and plate tectonics. Supercontinent formation involves complex processes of “introversion” (closure of interior oceans), “extroversion” (closure of exterior oceans), or a combination of these processes in uniting dispersed continental fragments. Recent developments in numerical modeling and advancements in computation techniques enable us to simulate Earth's mantle convection with drifting continents under realistic convection vigor and rheology in Earth-like geometry (i.e., 3D spherical-shell). We report a numerical simulation of 3D mantle convection, incorporating drifting deformable continents, to evaluate supercontinent processes in a realistic mantle convection regime. Our results show that supercontinents are assembled by a combination of introversion and extroversion processes. Small-scale thermal heterogeneity dominates deep mantle convection during the supercontinent cycle, although large-scale upwelling plumes intermittently originate under the drifting continents and/or the supercontinent.

Mantle dynamics in super-Earths: Post-perovskite rheology and self-regulation of viscosity

The discovery of extra-solar “super-Earth” planets with sizes up to twice that of Earth has prompted interest in their possible lithosphere and mantle dynamics and evolution. Simple scalings suggest that super-Earths are more likely than an equivalent Earth-sized planet to be undergoing plate tectonics. Generally, viscosity and thermal conductivity increase with pressure while thermal expansivity decreases, resulting in lower convective vigour in the deep mantle, which, if extralopated to the largest super-Earths might, according to conventional thinking, result in no convection in their deep mantles due to the very low effective Rayleigh number. Here we evaluate this. First, as the mantle of a super-Earth is made mostly of post-perovskite we here extend the density functional theory (DFT) calculations of post-perovskite activation enthalpy of to a pressure of 1TPa, for both slowest diffusion (upper-bound rheology) and fastest diffusion (lower-bound rheology) directions. Along a 1600K adiabat the upper-bound rheology would lead to a post-perovskite layer of a very high (∼1030Pas) but relatively uniform viscosity, whereas the lower-bound rheology leads to a post-perovskite viscosity increase of ∼7 orders of magnitude with depth; in both cases the deep mantle viscosity would be too high for convection. Second, we use these DFT-calculated values in statistically steady-state numerical simulations of mantle convection and lithosphere dynamics of planets with up to ten Earth masses. The models assume a compressible mantle including depth-dependence of material properties and plastic yielding induced plate-like lithospheric behaviour. Results confirm the likelihood of plate tectonics for planets with Earth-like surface conditions (temperature and water) and show a self-regulation of deep mantle temperature. The deep mantle is not adiabatic; instead feedback between internal heating, temperature and viscosity regulates the temperature such that the viscosity has the value needed to facilitate convective loss of the radiogenic heat, which results in a very hot perovskite layer for the upper-bound rheology, a super-adiabatic perovskite layer for the lower-bound rheology, and an azimuthally-averaged viscosity of no more than 1026Pas. Convection in large super-Earths is characterised by large upwellings (even with zero basal heating) and small, time-dependent downwellings, which for large super-Earths merge into broad downwellings. In the context of planetary evolution, if, as is likely, a super-Earth was extremely hot/molten after its formation, it is thus likely that even after billions of years its deep interior is still extremely hot and possibly substantially molten with a “super basal magma ocean” – a larger version of the proposal of Labrosse et al. (Labrosse, S., Hernlund, J.W., Coltice, N. [2007]. Nature 450, 866–869), although this depends on presently unknown melt–solid density contrast and solidus.

Three-dimensional visualization of numerically simulated, present-day global mantle flow

The behavior of mantle convection plays a primary role in the earth’s dynamic evolution. Numerical simulation of mantle convection is one of the most powerful tools to resolve the thermal and mechanical evolution in the earth’s mantle. The developments in numerical techniques and increases in computational power over the last 10–15 years have facilitated great progress in simulation studies of mantle convection. In addition, with increases in computational resolution, one can expect greater focus on the use of three-dimensional (3D) visualization of mantle convection in gaining a deeper understanding of mantle dynamics. This paper presents simulations of the present-day, instantaneous global mantle flow based on geophysical observations such as seismic tomography models and 3D visualizations of present-day mantle convection. Numerical results are summarized in terms of two aspects of geodynamics: (1) the mixing properties and topological structure of the mantle interior, and (2) the presence of strong mantle upwelling beneath the South Pacific enhanced by the low-viscosity D″ layer.

Mantle dynamics with pressure- and temperature-dependent thermal expansivity and conductivity

In numerical simulations of mantle convection it is commonly assumed that the coefficients of thermal expansion α and thermal conduction k are either constant or pressure-dependent. Pressure changes are generally computed using parametrizations that rely on extrapolations of low-pressure data for a single upper-mantle phase. Here we collect data for both the pressure and temperature dependence of α from a database of first-principles calculations, and of k from recent experimental studies. We use these data-sets to construct analytical parametrizations of α and k for the major upper- and lower-mantle phases that can be easily incorporated into exisiting convection codes. We then analyze the impact of such parametrizations on Earth’s mantle dynamics by employing two-dimensional numerical models of thermal convection. When α is the only variable parameter, both its temperature and pressure dependence enhance hot plumes and tend to inhibit the descent of cold downwellings. Taking into account a variable k leads to a strong increase of the bulk mantle temperature, which reduces the buoyancy available to amplify bottom boundary layer instabilities and causes mantle flow to be driven primarily by the instability of cold plates whose surface velocity also tends to rise. When both parameters are considered together, we observe an increased propensity to local layering which favors slab stagnation in the transition zone and subsequent thickening in the lower mantle. Furthermore, the values of k near the core-mantle boundary ultimately control the effect of this physical property on convection, which stresses the importance of determining the thermal conductivity of the post-perovskite phase.

Mantle temperature under drifting deformable continents during the supercontinent cycle

AbstractThe thermal heterogeneity of the Earth's mantle under the drifting continents during a supercontinent cycle is a controversial issue in earth science. Here, a series of numerical simulations of mantle convection are performed in 3‐D spherical‐shell geometry, incorporating drifting deformable continents and self‐consistent plate tectonics, to evaluate the subcontinental mantle temperature during a supercontinent cycle. Results show that the laterally averaged temperature anomaly of the subcontinental mantle remains within several tens of degrees (±50°C) throughout the simulation time. Even after the formation of the supercontinent and the development of subcontinental plumes due to the subduction of the oceanic plates, the laterally averaged temperature anomaly of the deep mantle under the continent is within +10°C. This implies that there is no substantial temperature difference between the subcontinental and suboceanic mantles during a supercontinent cycle. The temperature anomaly immediately beneath the supercontinent is generally positive owing to the thermal insulation effect and the active upwelling plumes from the core‐mantle boundary. In the present simulation, the formation of a supercontinent causes the laterally averaged subcontinental temperature to increase by a maximum of 50°C, which would produce sufficient tensional force to break up the supercontinent. The regular periodicity of the supercontinent cycles observed in previous 2‐D and 3‐D simulation models with rigid nondeformable continents and without self‐consistent plate tectonics is not confirmed.

Large-scale adaptive mantle convection simulation

SUMMARY A new generation, parallel adaptive-mesh mantle convection code, Rhea, is described and benchmarked. Rhea targets large-scale mantle convection simulations on parallel computers, andthushasbeendevelopedwithastrongfocusoncomputationalefficiencyandparallelscalabilityofbothmeshhandlingandnumericalsolvers.Rheabuildsmantleconvectionsolversona collection ofparallel octree-based adaptivefiniteelement librariesthatsupport new distributed data structures and parallel algorithms for dynamic coarsening, refinement, rebalancing and repartitioning of the mesh. In this study we demonstrate scalability to 122 880 compute cores and verify correctness of the implementation. We present the numerical approximation and convergence properties using 3-D benchmark problems and other tests for variable-viscosity Stokes flow and thermal convection.

Radial 1‐D seismic structures in the deep mantle in mantle convection simulations with self‐consistently calculated mineralogy

Numerical thermo‐chemical mantle convection simulations in a spherical annulus geometry with self‐consistently calculated mineralogy and mineral physics are used to predict detailed deep mantle seismic structures, particularly local radial profiles of shear wave velocity (Vs) and bulk sound velocity (Vb). The predicted structures are compared to seismological observations and used to guide the interpretation of seismic observations and to test the model. The mantle composition is described as a mixture of MORB (Mid‐Oceanic‐Ridge‐Basalt) and harzburgitic end‐members in the Na2O‐CaO‐FeO‐MgO‐Al2O3‐SiO2system. To assess the influence of chemical variability, four different sets of end‐member compositions are evaluated. Results confirm that the post‐perovskite (pPv) phase causes anti‐correlated S wave and bulk sound velocities in the deep mantle, due to pPv being fast in Vs but slow in bulk sound velocity. Local 1‐D seismic profiles display great lateral variability, and often have multiple discontinuities in the deep mantle due to MORB layers in folded slabs, with a positive Vs anomaly and negative bulk sound anomaly, or the perovskite‐pPv phase transition. The pPv transition is not visible inside piles of segregated MORB because of the high temperature and small velocity contrast of pPv in MORB. Piles of segregated MORB are seismically slow in both Vs and bulk sound despite being expected to be fast in Vs, because they are hotter than the surrounding material. Anelasticity has a significant influence on profiles of Vs only in the lower thermal boundary layer, which corresponds to below 2600 to 2800 km depth depending on region, where temperatures are higher than the extrapolated adiabat. These results indicate the importance of using a joint geodynamical‐mineralogical approach to predict and aid in the interpretation of deep mantle seismic structure, because interpretations based on seismology and mineral physics alone may be misleading and do not capture the strong lateral variability in 1‐D structure obtained here: for example, multiple reflections arising from folded slabs and the precise balance between thermal and compositional influences on seismic structure.

Impact of grain size on the convection of terrestrial planets

This article presents a set of simulations of mantle convection, using a new model of grain size‐dependent rheology. In the present paper, it is shown that this rheology behaves in many ways as a visco‐plastic rheology. I use a model of grain size evolution which has been calibrated on experimental data in a previous paper. In this physical model, the grain size is directly related to the stress state, following a temperature‐dependent piezometric law. The rheology used here allows both diffusion and dislocation creep, depending on the grain size. At low stress, the grain size is high and forces the rheology to be dislocation dominated. For sufficiently high stresses, the equilibrium regime reached by the grains is located in the diffusion creep. In this case, the viscosity is linked to a stress‐dependent grain size, which actually makes the rheology more non‐Newtonian than it is in dislocation creep. This experimentally calibrated model allowed me to perform a set of numerical experiments of convection in which the rheology may be diffusion or dislocation creep dominated, depending on the state of the stress tensor. Then, The stress exponent varies from 3 to 5 because of grain size, which has a large impact on the temperature dependence of the viscosity. The present paper shows the impact of this new model on the convection regimes of terrestrial planets. In particular, for a wide range of parameters, I observe the episodic regime which is thought to govern the dynamics of Venus. This process of episodic resurfacing was obtained in previous simulations using visco‐plastic rheologies is a tight range of parameters. I obtain it here without using an ad hoc plasticity law, only using a viscous rheology based on laboratory measurements. In these simulations, I show that the cooling rate of the terrestrial planets may be largely modified by the consideration of a grain size‐dependent rheology.

Plume’s buoyancy and heat fluxes from the deep mantle estimated by an instantaneous mantle flow simulation based on the S40RTS global seismic tomography model

It is still an open question as to how much heat is transported from the deep mantle to the upper mantle by mantle upwelling plumes, which would impose a strong constraint on models of the thermal evolution of the earth. Here I perform numerical computations of instantaneous mantle flow based on a recent highly resolved global seismic tomography model (S40RTS), apply new simple fluid dynamics theories to the plume’s radius and velocity, considering a Poiseuille flow assumption and a power-law relationship between the boundary layer thickness and Rayleigh number, and estimate the plume’s buoyancy and heat fluxes from the deep lower mantle under varying plume viscosity. The results show that for some major mantle upwelling plumes with localized strong ascent velocity under the South Pacific and Africa, the buoyancy fluxes of each plume beneath the ringwoodite to perovskite+magnesiowüstite (“660-km”) phase decomposition boundary are comparable to those inferred from observed hotspot swell volumes on the earth, i.e., on the order of 1Mgs−1, when the plume viscosity is 1019–1020Pas. This result, together with previous numerical simulations of mantle convection and the gentle Clausius–Clapeyron slope for the 660-km phase decomposition derived from recent high-pressure measurements under dehydrated/hydrated conditions in the mantle transition zone, implies that mantle upwelling plumes in the lower mantle penetrate the 660-km phase decomposition boundary without significant loss in thermal buoyancy because of the weak thermal barrier at the 660-km boundary. The total plume heat flux under the South Pacific is estimated to be about 1TW beneath the 660-km boundary, which is significantly smaller than the core–mantle boundary heat flux. Previously published scaling laws for the plume’s radius and velocity based on a plume spacing theory, which explains well plume dynamics in three-dimensional time-dependent mantle convection, suggest that these plume fluxes depend critically on plume viscosity. This result implies that the predicted plume viscosity may not be as low as expected, i.e., on the order of 1020Pas, considering the weak thermal barrier for upwelling plumes at the 660-km boundary, and the plume heat fluxes in the deep lower mantle under the South Pacific may be comparable to the lower limit for estimated core-mantle boundary heat flux inferred from a number of independent geodynamic studies.

Depth-dependent viscosity and mantle stress amplification: implications for the role of the asthenosphere in maintaining plate tectonics

Boundary layer theory is used to derive scaling relationships for plate stresses in a mantle convection system with a low-viscosity asthenosphere. The theory assumes a plate tectonic like mode of mantle convection with flow driven by an active upper boundary layer. The theory predicts that the confinement of horizontal mantle flow within a low-viscosity, sublithospheric channel can lead to an increase in plate stress compared to the case lacking a channel (even if the absolute viscosity of the sublithosphere mantle does not change between the two cases). The theory further predicts increasing shear stress with decreasing low-viscosity channel thickness. If the thickness of tectonic plates is determined dominantly by a dehydrated chemical lithosphere, then the plate normal stress is predicted to also increase with decreasing channel thickness. We use 3-D spherical shell simulations of mantle convection with temperature-, depth- and stress dependent rheology to test scaling trends. The simulations and theoretical scalings demonstrate that a low-viscosity layer (asthenosphere) can amplify convective stresses. If the level of convective stress plays a role in maintaining and/or reactivating plate boundaries, this suggests that a relatively thin low viscosity layer may help to maintain plate tectonics. The numerical simulations support this suggestion as they show that an increase in the thickness of a low viscosity channel can cause the system to transition from an active-lid mode of convection to a stagnant lid state. Collectively, the simulations and theoretical scalings lead to the conclusion that the role of the asthenosphere in maintaining plate tectonics does not come principally from a basal lubrication effect, associated with a low absolute asthenosphere viscosity, but, instead, from a mantle flow channelization effect, associated with a high viscosity contrast from the asthenosphere to the mantle below.

High accuracy mantle convection simulation through modern numerical methods

The first author was supported by the Graduate School in Mathematics and Computing (FMB) at Uppsala University, Sweden. The second and third authors are supported in part through the Computational Infrastructure in Geodynamics initiative (CIG), through the National Science Foundation under Award No. EAR-0949446 and The University of California-Davis. This publication is based in part on work supported by Award No. KUS-C1-016-04, made by King Abdullah University of Science and Technology (KAUST). The third author is also supported in part by an Alfred P. Sloan Research Fellowship.Some computations for this paper were performed on the 'Ranger' and 'Lonestar' clusters at the Texas Advanced Computing Center (TACC), and the 'Brazos' and 'Hurr' clusters at the Institute for Applied Mathematics and Computational Science (IAMCS) at Texas A&M University. Ranger was funded by NSF award OCI-0622780, and we used an allocation obtained under NSF award TG-MCA04N026. The authors acknowledge the Texas A&M Supercomputing Facility for providing computing resources on 'Lonestar' useful in conducting the research reported in this paper. Part of Brazos was supported by NSF award DMS-0922866. Hurr is supported by Award No. KUS-C1-016-04 made by King Abdullah University of Science and Technology (KAUST). Some computations were performed on resources provided by SNIC through Uppsala Multidisciplinary Center for Advanced Computational Science (UPPMAX) under project p2010002.

Benchmarking FEniCS for mantle convection simulations

This paper evaluates the usability of the FEniCS Project for mantle convection simulations by numerical comparison to three established benchmarks. The benchmark problems all concern convection processes in an incompressible fluid induced by temperature or composition variations, and cover three cases: (i) steady-state convection with depth- and temperature-dependent viscosity, (ii) time-dependent convection with constant viscosity and internal heating, and (iii) a Rayleigh–Taylor instability. These problems are modeled by the Stokes equations for the fluid and advection–diffusion equations for the temperature and composition. The FEniCS Project provides a novel platform for the automated solution of differential equations by finite element methods. In particular, it offers a significant flexibility with regard to modeling and numerical discretization choices; we have here used a discontinuous Galerkin method for the numerical solution of the advection–diffusion equations. Our numerical results are in agreement with the benchmarks, and demonstrate the applicability of both the discontinuous Galerkin method and FEniCS for such applications.

Hysteresis in mantle convection: Plate tectonics systems

We use simulations of mantle convection with surface yielding to show that multiple tectonic regimes are possible for equivalent system parameter values. Models with the same lithospheric strength parameters and the same vigor of convection can display different modes of tectonics. Within the region of multiple solutions, the evolutionary pathway of the system is the dominant factor that determines the tectonic mode (e.g., whether mantle convection operates in a plate tectonic like mode). The extent of the multiple regimes window is found to increase with the temperature‐dependent viscosity contrast across the system. The implication for models that seek to predict the tectonic regimes of planets is that the temporal evolution of the planet needs to be taken into account. A further implication is that modeling studies can lead to different conclusions regarding the tectonic state of a planet, extra‐solar planets in particular, despite the final model parameter values remaining equivalent.

The conditions for plate tectonics on super-Earths: Inferences from convection models with damage

Numerical simulations of mantle convection with a damage-grainsize feedback are used to develop scaling laws for predicting conditions at which super-Earths would have plate tectonics. In particular, the numerical simulations are used to determine how large a viscosity ratio between pristine lithosphere and mantle (μl/μm) can be offset by damage to allow mobile (plate-like) convection. Regime diagrams of μl/μm versus the damage number (D) show that the transition from stagnant lid to mobile convection occurs for higher μl/μm as D increases; a similar trend occurs for increasing Rayleigh number. We hypothesize a new criterion for the onset of plate tectonics on terrestrial planets: that damage must reduce the viscosity of shear zones in the lithosphere to a critical value equivalent to the underlying mantle viscosity; a scaling law based on this hypothesis reproduces the numerical results. For the Earth, damage is efficient in the lithosphere and provides a viable mechanism for the operation of plate tectonics. We scale our theory to super-Earths and map out the transition between plate-like and stagnant-lid convection with a “planetary plate-tectonic phase” diagram in planet size-surface temperature space. Both size and surface conditions are found to be important, with plate tectonics being favored for larger, cooler planets. This gives a natural explanation for Earth, Venus, and Mars, and implies that plate tectonics on exoplanets should correlate with size, incident solar radiation, and atmospheric composition.

Influence of magmatism on mantle cooling, surface heat flow and Urey ratio

Two-dimensional thermo-chemical mantle convection simulations are used to investigate the influence of melting-inducted differentiation on the thermal evolution of Earth's mantle, focussing in particular on matching the present-day surface heat flow and the ‘Urey ratio’. The influence of internal heating rate, initial mantle temperature and partitioning of heat-producing elements into basaltic crust are studied. High initial mantle temperatures, which are expected following Earth's accretion, cause major differences in early mantle thermo-chemical structures, but by the present-day surface heat flux and internal structures are indistinguishable from cases with a low initial temperature. Assuming three different values of mantle heat production that vary by more than a factor of two results in small differences in present-day heat flow, as does assuming different partitioning ratios of heat-producing elements into crust. Indeed, all of the cases presented here, regardless of exact parameters, have approximately Earth's present-day heat flow, with substantial fractions coming from the core and from mantle cooling. As a consequence of the model present-day surface heat flow varying only slightly with parameters, the Urey ratio (the ratio of total heat production to the total surface heat flow) is highly dependent on the amount of internal heat production, and due to the large uncertainty in this, the Urey ratio is considered to be a much poorer constraint on thermal evolution than the heat flow. The range of present-day Urey ratio observed in simulations here is about 0.3 to 0.5, which is consistent with observational and geochemical constraints (Jaupart et al., 2007). Magmatic heat transport contributes an upper bound of 9% to Earth's present-day heat loss but a much higher fraction at earlier times—often more than convective heat loss—so neglecting this causes an overestimation of the Urey ratio. Magmatic heat transport also plays an important role in mantle cooling. Considering these points, it is important to include magmatic effects when attempting to understand the thermal evolution of the Earth.

Influence of the Ringwoodite-Perovskite transition on mantle convection in spherical geometry as a function of Clapeyron slope and Rayleigh number

Abstract. We investigate the influence on mantle convection of the negative Clapeyron slope ringwoodite to perovskite and ferro-periclase mantle phase transition, which is correlated with the seismic discontinuity at 660 km depth. In particular, we focus on understanding the influence of the magnitude of the Clapeyron slope (as measured by the Phase Buoyancy parameter, P) and the vigour of convection (as measured by the Rayleigh number, Ra) on mantle convection. We have undertaken 76 simulations of isoviscous mantle convection in spherical geometry, varying Ra and P. Three domains of behaviour were found: layered convection for high Ra and more negative P, whole mantle convection for low Ra and less negative P, and transitional behaviour in an intervening domain. The boundary between the layered and transitional domain was fit by a curve P = α Raβ where α = −1.05, and β = −0.1, and the fit for the boundary between the transitional and whole mantle convection domain was α = −4.8, and β = −0.25. These two curves converge at Ra ≈ 2.5 × 104 (well below Earth mantle vigour) and P ≈ −0.38. Extrapolating to high Ra, which is likely earlier in Earth history, this work suggests a large transitional domain. It is therefore likely that convection in the Archean would have been influenced by this phase change, with Earth being at least in the transitional domain, if not the layered domain.

A note on continents and the Earth’s Urey ratio

The Urey ratio (mantle heat production divided by heat loss) is a key constraint for thermal history models. Recent Urey ratio estimates are in the range of 0.21–0.49. It has been well noted in the literature that classic thermal history models predict a larger value. The extent of the discrepancy is considered significant enough to be deemed a paradox. The ‘Urey ratio paradox’ has motivated a number of studies and some considerable debates. Classic thermal history models, and the majority of those that have followed, do not explicitly consider the effects of continents. For consistency, the present day observations they seek to match should be adjusted to account for continents. We use mantle convection simulations with continents to argue that the adjusted ratio shifts to 0.33–0.76. Classic thermal history models predict a Urey ratio that, although at the high end, is within this range. This suggests that the ‘paradox’ has been exaggerated. It also points toward a first-order role for continents in the Earth’s thermal evolution

Two-dimensional mantle convection simulations using an internal state variable model: the role of a history dependent rheology on mantle convection

SUMMARY We apply the Bammann inelastic internal state variable model (BIISV) to a mantle convection code TERRA2D to investigate the influence of a history dependent solid mechanics model on mantle convection. We compare and contrast the general purpose BIISV model to the commonly used power-law model. We implemented the BIISV model using a radial return algorithm and tested it against previously published mantle convection simulation results for verification. Model constants for the BIISV are used based on experimental stress–strain behaviour found in the literature. After implementation we give illustrative simulation examples were the BIISV produces hardened areas on the cold thermal boundary layer that the power-law model cannot produce. The hardened boundary layers divert material downward giving a plausible reason for the current subduction zones that are present on the Earth.