Mantle Convection Models Research Articles

Closure of the Mongol–Okhotsk Ocean: Insights from seismic tomography and numerical modelling

Abstract The existence of the Palaeozoic and Mesozoic Mongol–Okhotsk Ocean is evident from the Mongol–Okhotsk suture, which stretches from central Mongolia to the Sea of Okhotsk. A lack of sufficient palaeomagnetic data and an otherwise diffuse suture with an abrupt termination to the west has led to difficulties in reconstructing the history, geometry and closure of this ocean. Both the timing and style of the ocean's closure are unclear and have led to several alternative reconstructions. Closure timing ranges between the Late Jurassic (∼155 Ma) and beginning of the Early Cretaceous (∼120 Ma), and the proposed kinematics include contemporaneous subduction along two opposite margins, subduction along only one margin or with a component of left-lateral shear. In the present study, numerical models of mantle convection are coupled with global plate reconstructions to investigate ambiguities regarding the closure of the Mongol–Okhotsk Ocean. In order to decipher the tectonic history of this enigmatic region, two end-member scenarios of subduction location – either along the present-day northern or the southern margins of the Mongol–Okhotsk Ocean – are imposed as kinematic surface boundary conditions for the past 230 Myrs. Through a comparison to seismic tomography, the results indicate a preferred subduction history along the Siberian margin (relative northern margin) of the Mongol–Okhotsk Ocean. At present-day, the slab remnant is predicted to be located farther west than previously proposed. Furthermore, we find that the subducting slabs in this region generate a hot, dense pile at the same location and with a similar shape as the Perm Anomaly.

Global-scale modelling of melting and isotopic evolution of Earth's mantle: melting modules for TERRA

Abstract. Many outstanding problems in solid-Earth science relate to the geodynamical explanation of geochemical observations. Currently, extensive geochemical databases of surface observations exist, but satisfying explanations of underlying mantle processes are lacking. One way to address these problems is through numerical modelling of mantle convection while tracking chemical information throughout the convective mantle. We have implemented a new way to track both bulk compositions and concentrations of trace elements in a finite-element mantle convection code. Our approach is to track bulk compositions and trace element abundances via particles. One value on each particle represents bulk composition and can be interpreted as the basalt component. In our model, chemical fractionation of bulk composition and trace elements happens at self-consistent, evolving melting zones. Melting is defined via a composition-dependent solidus, such that the amount of melt generated depends on pressure, temperature and bulk composition of each particle. A novel aspect is that we do not move particles that undergo melting; instead we transfer the chemical information carried by the particle to other particles. Molten material is instantaneously transported to the surface layer, thereby increasing the basalt component carried by the particles close to the surface and decreasing the basalt component in the residue. The model is set to explore a number of radiogenic isotopic systems, but as an example here the trace elements we choose to follow are the Pb isotopes and their radioactive parents. For these calculations we will show (1) the evolution of the distribution of bulk compositions over time, showing the buildup of oceanic crust (via melting-induced chemical separation in bulk composition), i.e. a basalt-rich layer at the surface, and the transportation of these chemical heterogeneities through the deep mantle; (2) the amount of melt generated over time; (3) the evolution of the concentrations and abundances of different isotopes of the trace elements (U, Th, K and Pb) throughout the mantle; and (4) a comparison to a semi-analytical theory relating observed arrays of correlated Pb isotope compositions to melting age distributions (Rudge, 2006).

Formation of Australian continental margin highlands driven by plate–mantle interaction

Passive margin highlands occur on most continents on Earth and play a critical role in the cycle of weathering, erosion, and atmospheric circulation. Yet, in contrast to the well-developed understanding of collisional mountain belts, such as the Alps and Himalayas, the origin of less elevated (1–2 km) passive margin highlands is still unknown. The eastern Australian highlands are a prime example of these plateaus, but compared to others they have a well-documented episodic uplift history spanning 120 million years. We use a series of mantle convection models to show that the time-dependent interaction of plate motion with mantle downwellings and upwellings accounts for the broad pattern of margin uplift phases. Initial dynamic uplift of 400–600 m from 120–80 Ma was driven by the eastward motion of eastern Australia's margin away from the sinking eastern Gondwana slab, followed by tectonic quiescence to about 60 Ma in the south (Snowy Mountains). Renewed uplift of ∼700 m in the Snowy Mountains is propelled by the gradual motion of the margin over the edge of the large Pacific mantle upwelling. In contrast the northernmost portion of the highlands records continuous uplift from 120 Ma to present-day totalling about 800 m. The northern highlands experienced a continuous history of dynamic uplift, first due to the end of subduction to the east of Australia, then due to moving over a large passive mantle upwelling. In contrast, the southern highlands started interacting with the edge of the large Pacific mantle upwelling ∼40–50 million years later, resulting in a two-phase uplift history. Our results are in agreement with published uplift models derived from river profiles and the Cretaceous sediment influx into the Ceduna sub-basin offshore southeast Australia, reflecting the fundamental link between dynamic uplift, fluvial erosion and depositional pulses in basins distal to passive margin highlands.

Mantle‐induced subsidence and compression in SE Asia since the early Miocene

AbstractRift basins developed extensively across Sundaland, the continental core of Southeast Asia, since the Eocene. Beginning in the early Miocene, basins in southern Sundaland experienced widespread synchronous compression (inversion) and marine inundation, despite a large drop in long‐term global sea level. The mechanism for this large‐scale synchronous regional sea level rise, basin inversion, and subsidence is not well understood and contrary to expectations from traditional basin models and eustatic sea level trends. We present geodynamic models of mantle convection with both deformable and rigid plate reconstructions to investigate this enigma. Models suggest that a slab stagnates within the transition zone beneath Southeast Asia before the Miocene. The stagnant slab penetrated through the 660 km mantle discontinuity during the early Miocene and formed a slab avalanche event, due to continuous subduction and accumulation of negatively buoyant slabs. This avalanche may have induced large‐scale marine inundation, regional compression, and basin inversion across southern Sundaland. We argue mantle convection induced large‐scale basin compression, in contrast to conventional plate margin‐induced compression; this suggests mantle convection may exert a much stronger control on surface processes than previously recognized.

Constraining mantle viscosity structure for a thermochemical mantle using the geoid observation

AbstractLong‐wavelength geoid anomalies provide important constraints on mantle dynamics and viscosity structure. Previous studies have successfully reproduced the observed geoid using seismically inferred buoyancy in whole‐mantle convection models. However, it has been suggested that large low shear velocity provinces (LLSVPs) underneath Pacific and Africa in the lower mantle are chemically distinct and are likely denser than the ambient mantle. We formulate instantaneous flow models based on seismic tomographic models to compute the geoid and constrain mantle viscosity by assuming both thermochemical and whole‐mantle convection. Geoid modeling for the thermochemical model is performed by considering the compensation effect of dense thermochemical piles and removing buoyancy structure of the compensation layer in the lower mantle. Thermochemical models well reproduce the observed geoid, thus reconciling the geoid with the interpretation of LLSVPs as dense thermochemical piles. The viscosity structure inverted for thermochemical models is nearly identical to that of whole‐mantle models. In the preferred model, the lower mantle viscosity is ∼10 times higher than the upper mantle viscosity that is ∼10 times higher than the transition zone viscosity. The weak transition zone is consistent with the proposed high water content there. The geoid in thermochemical mantle models is sensitive to seismic structure at midmantle depths, suggesting a need to improve seismic imaging resolution there. The geoid modeling constrains the vertical extent of dense and stable chemical piles to be within ∼500 km above CMB. Our results have implications for mineral physics, seismic tomographic studies, and mantle convection modeling.

The compressible adjoint equations in geodynamics: derivation and numerical assessment

The adjoint method is a powerful means to obtain gradient information in a mantle convection model relative to past flow structure. While the adjoint equations in geodynamics have been derived for the conservation equations of mantle flow in their incompressible form, the applicability of this approximation to Earth is limited, because density increases by almost a factor of two (Dziewonski and Anderson, Phys Earth Planet Inter 25:297–356, 1981) from the surface to the Core Mantle Boundary. Here we introduce the compressible adjoint equations for the conservation equations in the anelastic-liquid approximation. Our derivation applies an operator formulation in Hilbert spaces, to connect to recent work in seismology (Fichtner et al., Phys Earth Planet Inter 157(1–2):86–104, 2006b) and geodynamics (Horbach et al., Int J Geomath 5(2):163–194, 2014), where the approach was used to derive the adjoint equations for the wave equation and incompressible mantle flow. We present numerical tests of the newly derived equations based on twin experiments, focusing on three simulations. A first, termed Compressible, assumes the compressible forward and adjoint equations, and represents the consistent means of including compressibility effects. A second, termed Mixed, applies the compressible forward equation, but ignores compressibility effects in the adjoint equations, where the incompressible equations are used instead. A third simulation, termed Incompressible, neglects compressibility effects entirely in the forward and adjoint equations relative to the reference twin. The compressible and mixed formulations successfully restore earlier mantle flow structure, while the incompressible formulation yields noticeable artifacts. Our results suggest the use of a compressible formulation, when applying the adjoint method to seismically derived mantle heterogeneity structure.

Pipelined, Flexible Krylov Subspace Methods

We present variants of the Conjugate Gradient (CG), Conjugate Residual (CR), and Generalized Minimal Residual (GMRES) methods which are both pipelined and flexible. These allow computation of inner products and norms to be overlapped with operator and nonlinear or nondeterministic preconditioner application.The methods are hence aimed at hiding network latencies and synchronizations which can become computational bottlenecks in Krylov methods on extreme-scale systems or in the strong-scaling limit. The new variants are not arithmetically equivalent to their base flexible Krylov methods, but are chosen to be similarly performant in a realistic use case, the application of strong nonlinear preconditioners to large problems which require many Krylov iterations. We provide scalable implementations of our methods as contributions to the PETSc package and demonstrate their effectiveness with practical examples derived from models of mantle convection and lithospheric dynamics with heterogeneous viscosity structure. These represent challenging problems where multiscale nonlinear preconditioners are required for the current state-of-the-art algorithms, and are hence amenable to acceleration with our new techniques. Large-scale tests are performed in the strong-scaling regime on a contemporary leadership supercomputer, where speedups approaching, and even exceeding $2\times$ can be observed. We conclude by analyzing our new methods with a performance model targeted at future exascale machines.

The mantle convection model with non-Newtonian rheology and phase transitions: The flow structure and stress fields

The mantle convection model with phase transitions, non-Newtonian viscosity, and internal heat sources is calculated for two-dimensional (2D) Cartesian geometry. The temperature dependence of viscosity is described by the Arrhenius law with a viscosity step of 50 at the boundary between the upper and lower mantle. The viscosity in the model ranges within 4.5 orders of magnitude. The use of the non-Newtonian rheology enabled us to model the processes of softening in the zone of bending and subduction of the oceanic plates. The yield stress in the model is assumed to be 50 MPa. Based on the obtained model, the structure of the mantle flows and the spatial fields of the stresses σ xz and σ xx in the Earth’s mantle are studied. The model demonstrates a stepwise migration of the subduction zones and reveals the sharp changes in the stress fields depending on the stage of the slab detachment. In contrast to the previous model (Bobrov and Baranov, 2014), the self-consistent appearance of the rigid moving lithospheric plates on the surface is observed. Here, the intense flows in the upper mantle cause the drift and bending of the top segments of the slabs and the displacement of the plumes. It is established that when the upwelling plume intersects the boundary between the lower and upper mantle, it assumes a characteristic two-level structure: in the upper mantle, the ascending jet of the mantle material gets thinner, whereas its velocity increases. This effect is caused by the jump in the viscosity at the boundary and is enhanced by the effect of the endothermic phase boundary which impedes the penetration of the plume material from the lower mantle to the upper mantle. The values and distribution of the shear stresses σ xz and superlithostatic horizontal stresses σ xx are calculated. In the model area of the subducting slabs the stresses are 60–80 MPa, which is by about an order of magnitude higher than in the other mantle regions. The character of the stress fields in the transition region of the phase boundaries and viscosity step by the plumes and slabs is analyzed. It is established that the viscosity step and endothermic phase boundary at a depth of 660 km induce heterogeneities in the stress fields at the upper/lower mantle boundary. With the assumed model parameters, the exothermic phase transition at 410 km barely affects the stress fields. The slab regions manifest themselves in the stress fields much stronger than the plume regions. This numerically demonstrates that it is the slabs, not the plumes that are the main drivers of the convection. The plumes partly drive the convection and are partly passively involved into the convection stirred by the sinking slabs.

A sequential data assimilation approach for the joint reconstruction of mantle convection and surface tectonics

With the progress of mantle convection modelling over the last decade, it now becomes possible to solve for the dynamics of the interior flow and the surface tectonics to first order. We show here that tectonic data (like surface kinematics and seafloor age distribution) and mantle convection models with plate-like behaviour can in principle be combined to reconstruct mantle convection. We present a sequential data assimilation method, based on suboptimal schemes derived from the Kalman filter, where surface velocities and seafloor age maps are not used as boundary conditions for the flow, but as data to assimilate. Two stages (a forecast followed by an analysis) are repeated sequentially to take into account data observed at different times. Whenever observations are available, an analysis infers the most probable state of the mantle at this time, considering a prior guess (supplied by the forecast) and the new observations at hand, using the classical best linear unbiased estimate. Between two observation times, the evolution of the mantle is governed by the forward model of mantle convection. This method is applied to synthetic 2-D spherical annulus mantle cases to evaluate its efficiency. We compare the reference evolutions to the estimations obtained by data assimilation. Two parameters control the behaviour of the scheme: the time between two analyses, and the amplitude of noise in the synthetic observations. Our technique proves to be efficient in retrieving temperature field evolutions provided the time between two analyses is 10 Myr. If the amplitude of the a priori error on the observations is large (30 per cent), our method provides a better estimate of surface tectonics than the observations, taking advantage of the information within the physics of convection.

Seismic evidence of on-going sublithosphere upper mantle convection for intra-plate volcanism in Northeast China

A 3-D crustal and upper mantle S-wave velocity model of NE China is constructed by inversion of phase velocity dispersion curves at 6–140 s periods from ambient noise tomography and two-plane surface wave tomography. The seismic data used in this study are collected from 120 China Earthquake Administration (CEA) permanent stations and 127 portable stations of NECESSArray. We observe strong low S-wave velocity beneath the Changbaishan volcano in the upper mantle to at least 200-km depth, which is interpreted as a mantle upwelling beneath the Changbaishan volcano that is consistent with the body wave tomographic image. The Songliao Basin is dominated by a high velocity extending to at least 200-km depth. Built upon the observed velocity anomalies, we propose a sub-lithosphere mantle convection model for NE China in which the upwelling of upper mantle materials from the mantle transition zone to the Changbaishan volcano could induce a local sub-lithosphere convection in the upper mantle and the strong high velocity of the upper mantle beneath the Songliao Basin corresponds to the downwelling limb of this convection cell. The downwelling beneath the Songliao Basin could also induce secondary local convection in the asthenosphere to the west, leading to local asthenospheric upwelling beneath the Abaga and Halaha volcanoes in the Xing'an–Mongolia Orogenic Belt.

On retrodictions of global mantle flow with assimilated surface velocities

Modeling past states of Earth's mantle and relating them to geologic observations such as continental-scale uplift and subsidence is an effective method for testing mantle convection models. However, mantle convection is chaotic and two identical mantle models initialized with slightly different temperature fields diverge exponentially in time until they become uncorrelated, thus limiting retrodictions (i.e., reconstructions of past states of Earth's mantle obtained using present information) to the recent past. We show with 3-D spherical mantle convection models that retrodictions of mantle flow can be extended significantly if knowledge of the surface velocity field is available. Assimilating surface velocities produces in some cases negative Lyapunov times (i.e., e-folding times), implying that even a severely perturbed initial condition may evolve toward the reference state. A history of the surface velocity field for Earth can be obtained from past plate motion reconstructions for time periods of a mantle overturn, suggesting that mantle flow can be reconstructed over comparable times.

Couette and Poiseuille flows in a low viscosity asthenosphere: Effects of internal heating rate, Rayleigh number, and plate representation

Mantle convection models with a low viscosity asthenosphere and high viscosity surface plates have been shown to produce very large aspect ratio convection cells like those inferred to exist in Earth’s mantle and to exhibit two asthenospheric flow regimes. When the surface plate is highly mobile, the plate velocity exceeds the flow velocities in the asthenosphere and the plate drives a Couette-type flow in the asthenospheric channel. For sluggish plates, the flow velocities in the asthenosphere exceed the plate velocity and the asthenospheric flow is more Poiseuille-like. It has been shown that under certain circumstances, flows become increasingly Couette-like as the aspect ratio of the plate is increased in numerical simulations. These models also show an increase in the average surface heat flux with aspect ratio which is counterintuitive, as one would expect that large aspect ratio models would result in older and colder oceanic lithosphere. Previous investigations have used single internal heating rates and Rayleigh numbers and a plate formulation that did not preclude significant deformation within the plate. In this paper, we investigate the conditions necessary for Couette and Poiseuille asthenospheric flows and for surface heat flux to increase with plate aspect ratio by varying the internal heating rate, the Rayleigh number and the representation of surface plates in 2D mantle convection models Plates are represented as a high viscosity layer with (1) a free-slip top surface boundary condition and (2) a force-balance boundary condition that imposes a constant surface velocity within the plate. We find that for models with a free-slip surface boundary condition, the internal heating rate and Rayleigh number do not strongly affect the dominance of Couette or Poiseuille flows in the asthenosphere but the increase in surface heat flux with model aspect ratio in the Poiseuille asthenospheric flow regime increases with internal heating rate. For models using the force-balance representation of surface plates, the flow regime in the asthenosphere is found to be almost independent of the plate aspect ratio and spatially averaged surface heat flux decreases with increasing aspect ratio. The dependence of heat flux on plate aspect ratio is an important consideration in Earth’s thermal history while the preponderance of Poiseuille or Couette flows in the asthenosphere may be important for understanding the mobility of surface plates and asthenospheric seismic anisotropy.

Formation of ridges in a stable lithosphere in mantle convection models with a viscoplastic rheology.

Numerical simulations of mantle convection with a viscoplastic rheology usually display mobile, episodic or stagnant lid regimes. In this study, we report a new convective regime in which a ridge can form without destabilizing the surrounding lithosphere or forming subduction zones. Using simulations in 2‐D spherical annulus geometry, we show that a depth‐dependent yield stress is sufficient to reach this ridge only regime. This regime occurs when the friction coefficient is close to the critical value between mobile lid and stagnant lid regimes. Maps of convective regime as a function of the parameters friction coefficients and depth dependence of viscosity are provided for both basal heating and mixed heating situations. The ridge only regime appears for both pure basal heating and mixed heating mode. For basal heating, this regime can occur for all vertical viscosity contrasts, while for mixed heating, a highly viscous deep mantle is required.

The long‐wavelength geoid from three‐dimensional spherical models of thermal and thermochemical mantle convection

AbstractThe Earth's long‐wavelength geoid anomalies have long been used to constrain the dynamics and viscosity structure of the mantle in an isochemical, whole mantle convection model. However, there is strong evidence that the seismically observed large low shear velocity provinces (LLSVPs) in the lower mantle underneath the Pacific and Africa are chemically distinct and likely denser than the ambient mantle. In this study, we have formulated dynamically self‐consistent 3‐D spherical mantle convection models to investigate how chemically distinct and dense piles above the core‐mantle boundary may influence the geoid. Our dynamic models with realistic mantle viscosity structure produce dominantly spherical harmonic degree‐2 convection, similar to that of the present‐day Earth. The models produce two broad geoid and topography highs over two major thermochemical piles in the lower mantle, consistent with the positive geoid anomalies over the Pacific and African LLSVPs. Our geoid analysis showed that the bottom layer with dense chemical piles contributes negatively to the total geoid, while the layer immediately above the chemical piles contributes positively to the geoid, canceling out the effect of the piles. Thus, the bottom part of the mantle, as a compensation layer, has zero net contribution to the total geoid, and the thickness of the compensation layer is ~1000 km or 2 to 3 times as thick as the chemical piles. Our results help constrain and interpret the large‐scale thermochemical structure of the mantle using surface observations of the geoid and topography, as well as seismic models of the mantle.

Iron-spin transition controls structure and stability of LLSVPs in the lower mantle

Seismic tomography models have revealed that there exist two large low shear velocity provinces (LLSVPs) beneath Africa and the Pacific Ocean in the Earth's lower mantle. Waveform modeling results suggest both LLSVPs have steep sharp side boundaries, which imply that they probably are compositionally heterogeneous from the ambient mantle. When applying the surface plate motion history in the last a few hundred million years (Ma) as the driving mechanism, numerical modeling has successfully reproduced the geographical distribution of the two LLSVPs in thermochemical mantle convection. However, two prominent seismic features of the LLSVPs, the steep side boundaries and the high elevation, can hardly be obtained in previous geodynamic models. Here, we include in our mantle convection model the effects of iron-spin transition of ferropericlase which substantially change physical properties of the mantle. Our results show that iron-spin transition plays a dominant role in controlling the structure and stability of LLSVPs. Large chemical blocks with steep side boundaries and high elevations up to ∼1200 km above core–mantle boundary (CMB) emerge in our models with the volume content of ferropericlase ∼20%. Such blocks cause a shear wave velocity decrease of ∼3.5% which is consistent with the seismic observations. Our results also show that these LLSVPs are transient structures in the lower mantle, which can typically last for a few hundred million years before destroyed by large-scale mantle motion.

Parameters controlling dynamically self‐consistent plate tectonics and single‐sided subduction in global models of mantle convection

AbstractRecent advances in numerical modeling allow global models of mantle convection to more realistically reproduce the behavior at convergent plate boundaries; in particular, the inclusion of a free surface at the outer boundary has been shown to facilitate self‐consistent development of single‐sided subduction. This allows for a more extensive study of subduction in the context of global mantle convection, as opposed to commonly used regional models. Our first study already indicated important differences between mantle convection with single‐sided subduction and mantle convection with double‐sided subduction. Here we further investigate the effect of various physical parameters and complexities on inducing Earth‐like plate tectonics and its evolution in time. Results reinforce the previous finding that using a free surface instead of a free‐slip outer boundary dramatically changes subduction style, with free surface cases displaying many episodes of single‐sided subduction, which leads to more realistic slab dip, stress state, trench retreat rate, and slab‐induced mantle flow. Longevity of single‐sided subduction is promoted by a layer of hydrated crust with a low yield strength to lubricate the subduction channel, a low‐viscosity asthenosphere, and a high strength of the slab (determined by a combination of high‐diffusion creep viscosity and intermediate friction coefficient), although its effective viscosity is in the observationally constrained range in the bending region. The time evolution displays interesting events including subduction polarity reversals, subduction shut‐off, and slab break‐off.

On the temporal evolution of long‐wavelength mantle structure of the Earth since the early Paleozoic

AbstractThe seismic structure of the Earth's lower mantle is characterized by a dominantly degree‐2 pattern with the African and Pacific large low shear velocity provinces (i.e., LLSVP) that are separated by circum‐Pacific seismically fast anomalies. It is important to understand the origin of such a degree‐2 mantle structure and its temporal evolution. In this study, we investigated the effects of plate motion history and mantle viscosity on the temporal evolution of the lower mantle structure since the early Paleozoic by formulating 3‐D spherical shell models of thermochemical convection. For convection models with realistic mantle viscosity and no initial structure, it takes about ∼50 Myr to develop dominantly degree‐2 lower mantle structure using the published plate motion models for the last either 120 Ma or 250 Ma. However, it takes longer time to develop the mantle structure for more viscous mantle. While the circum‐Pangea subduction in plate motion history models promotes the formation of degree‐2 mantle structure, the published pre‐Pangea plate motions before 330 Ma produce relatively cold lower mantle in the African hemisphere and significant degree‐1 structure in the early Pangea (∼300 Ma) or later times, even if the lower mantle has an initially degree‐2 structure and a viscosity as high as 1023 Pas. This suggests that the African LLSVP may not be stationary since the early Paleozoic. With the published plate motion models and lower mantle viscosity of 1022 Pas, our mantle convection models suggest that the present‐day degree‐2 mantle structure may have largely been formed by ∼200 Ma.

The effects of internal heating and large scale climate variations on tectonic bi-stability in terrestrial planets

We use 3D mantle convection and planetary tectonics models to explore the links between tectonic regimes and the level of internal heating within the mantle of a planet (a proxy for thermal age), planetary surface temperature, and lithosphere strength. At both high and low values of internal heating, for moderate to high lithospheric yield strength, hot and cold stagnant-lid (single plate planet) states prevail. For intermediate values of internal heating, multiple stable tectonic states can exist. In these regions of parameter space, the specific evolutionary path of the system has a dominant role in determining its tectonic state. For low to moderate lithospheric yield strength, mobile-lid behavior (a plate tectonic-like mode of convection) is attainable for high degrees of internal heating (i.e., early in a planet's thermal evolution). However, this state is sensitive to climate driven changes in surface temperatures. Relatively small increases in surface temperature can be sufficient to usher in a transition from a mobile- to a stagnant-lid regime. Once a stagnant-lid mode is initiated, a return to mobile-lid is not attainable by a reduction of surface temperatures alone. For lower levels of internal heating, the tectonic regime becomes less sensitive to surface temperature changes. Collectively our results indicate that terrestrial planets can alternate between multiple tectonic states over giga-year timescales. Within parameter space regions that allow for bi-stable behavior, any model-based prediction as to the current mode of tectonics is inherently non-unique in the absence of constraints on the geologic and climatic histories of a planet.

Tectonic speed limits from plate kinematic reconstructions

The motion of plates and continents on the planet's surface are a manifestation of long-term mantle convection and plate tectonics. Present-day plate velocities provide a snapshot of this ongoing process, and have been used to infer controlling factors on the speeds of plates and continents. However, present-day velocities do not capture plate behaviour over geologically representative periods of time. To address this shortcoming, we use a plate tectonic reconstruction approach to extract time-dependent plate velocities and geometries from which root mean square (RMS) velocities are computed, resulting in a median RMS plate speed of ∼4 cm/yr over 200 Myr. Linking tectonothermal ages of continental lithosphere to the RMS plate velocity analysis, we find that the increasing portions of plate area composed of continental and/or cratonic lithosphere significantly reduces plate speeds. Plates with any cratonic portion have a median RMS velocity of ∼5.8 cm/yr, while plates with more than 25% of cratonic area have a median RMS speed of ∼2.8 cm/yr. The fastest plates (∼8.5 cm/yr RMS speed) have little continental fraction and tend to be bounded by subduction zones, while the slowest plates (∼2.6–2.8 cm/yr RMS speed) have large continental fractions and usually have little to no subducting part of plate perimeter. More generally, oceanic plates tend to move 2–3 times faster than continental plates, consistent with predictions of numerical models of mantle convection. The slower motion of continental plates is compatible with deep keels impinging on asthenospheric flow and increasing shear traction, thus anchoring the plate in the more viscous mantle transition zone. We also find that short-lived (up to ∼10 Myr) rapid accelerations of Africa (∼100 and 65 Ma), North America (∼100 and 55 Ma) and India (∼130,80 and 65 Ma) appear to be correlated with plume head arrivals as recorded by large igneous province (LIPs) emplacement. By evaluating factors influencing plate speeds over the Mesozoic and Cenozoic, our temporal analysis reveals simple principles that can guide the construction and evaluation of absolute plate motion models for times before the Cretaceous in the absence of hotspot tracks and seafloor spreading histories. Based on the post-Pangea plate motions, one principle that can be applied to pre-Pangea times is that plates with less than ∼50% continental area can reach RMS velocities of ∼20 cm/yr, while plates with more than 50% continental fraction do not exceed RMS velocities of ∼10 cm/yr. Similarly, plates with large portions of continental or cratonic area with RMS velocities exceeding ∼15 cm/yr for more than ∼10 Myr should be considered as potential artefacts requiring further justification of plate driving forces in such scenarios.

Pangea breakup and northward drift of the Indian subcontinent reproduced by a numerical model of mantle convection.

Since around 200 Ma, the most notable event in the process of the breakup of Pangea has been the high speed (up to 20 cm yr−1) of the northward drift of the Indian subcontinent. Our numerical simulations of 3-D spherical mantle convection approximately reproduced the process of continental drift from the breakup of Pangea at 200 Ma to the present-day continental distribution. These simulations revealed that a major factor in the northward drift of the Indian subcontinent was the large-scale cold mantle downwelling that developed spontaneously in the North Tethys Ocean, attributed to the overall shape of Pangea. The strong lateral mantle flow caused by the high-temperature anomaly beneath Pangea, due to the thermal insulation effect, enhanced the acceleration of the Indian subcontinent during the early stage of the Pangea breakup. The large-scale hot upwelling plumes from the lower mantle, initially located under Africa, might have contributed to the formation of the large-scale cold mantle downwelling in the North Tethys Ocean.