# Tectonic predictions with mantle convection models

Nicolas Coltice, Grace E. Shephard

► To cite this version:

Nicolas Coltice, Grace E. Shephard. Tectonic predictions with mantle convection models. Geophysical Journal International, 2018, 213 (1), pp.16-29. 10.1093/gji/ggx531. hal-02327174

HAL Id: hal-02327174

https://univ-lyon1.hal.science/hal-02327174

Submitted on 13 Oct 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Tectonic predictions with mantle convection models

Nicolas Coltice<fnref n="1" /> and Grace E. Shephard<fnref n="2" />

Accepted 2017 December 11. Received 2017 November 27; in original form 2017 July 4

## SUMMARY

Over the past 15 yr, numerical models of convection in Earth's mantle have made a leap forward: they can now produce self-consistent plate-like behaviour at the surface together with deep mantle circulation. These digital tools provide a new window into the intimate connections between plate tectonics and mantle dynamics, and can therefore be used for tectonic predictions, in principle. This contribution explores this assumption. First, initial conditions at 30, 20, 10 and 0 Ma are generated by driving a convective flow with imposed plate velocities at the surface. We then compute instantaneous mantle flows in response to the guessed temperature fields without imposing any boundary conditions. Plate boundaries self-consistently emerge at correct locations with respect to reconstructions, except for small plates close to subduction zones. As already observed for other types of instantaneous flow calculations, the structure of the top boundary layer and upper-mantle slab is the dominant character that leads to accurate predictions of surface velocities. Perturbations of the rheological parameters have little impact on the resulting surface velocities. We then compute fully dynamic model evolution from 30 and 10 to 0 Ma, without imposing plate boundaries or plate velocities. Contrary to instantaneous calculations, errors in kinematic predictions are substantial, although the plate layout and kinematics in several areas remain consistent with the expectations for the Earth. For these calculations, varying the rheological parameters makes a difference for plate boundary evolution. Also, identified errors in initial conditions contribute to first-order kinematic errors. This experiment shows that the tectonic predictions of dynamic models over 10 My are highly sensitive to uncertainties of rheological parameters and initial temperature field in comparison to instantaneous flow calculations. Indeed, the initial conditions and the rheological parameters can be good enough for an accurate prediction of instantaneous flow, but not for a prediction after 10 My of evolution. Therefore, inverse methods (sequential or data assimilation methods) using short-term fully dynamic evolution that predict surface kinematics are promising tools for a better understanding of the state of the Earth's mantle.

**Key words:** Numerical modelling; Self-organization; Dynamics: convection currents, and mantle plumes; Dynamics of lithosphere and mantle; Kinematics of crustal and mantle deformation; Rheology: crust and lithosphere.

## 1 INTRODUCTION

In the theory of plate tectonics, the surface of the Earth is assumed to be divided into perfectly rigid plates, such that sufficient geological observations combined with geometric principles describe a coherent kinematic state. However, this revolutionary theory is not dynamic, hence it cannot be used to predict future and past states of the planet for which observations are too sparse or absent. Reconstructing past tectonics is therefore a difficult task (Gurnis et al. 2012), especially in areas where geological observations are lacking. For instance, 50 per cent of the world's present-day ocean floor is younger than 55 Ma, and a large fraction of the Pacific Ocean had disappeared prior to 60 Ma (Rowley 2008). Interpretation of

mantle seismic tomography can provide additional constraints, but the assumptions used still require testing (van Der Meer et al. 2010; Domeier et al. 2016). Unfortunately, even quantifying forces acting on plates today (Forsyth & Uyeda 1975) does not give access to how plate boundaries are generated and evolve. Analysing the plate velocity in tectonic reconstructions, for instance in terms of toroidal–poloidal partitioning brings more questions on the origins of plate velocity changes (Lithgow-Bertelloni et al. 1993).

As a consequence, dynamic models are needed to fill observational gaps. They can also handle diffuse deformation, extending the concept of plate tectonics beyond that of pure rigidity. These models consider that the plates and mantle constitute a single complex system (Bercovici 2003). Over the past 20 yr, numerical

<footnote n="1">¹Laboratoire de Géologie de Lyon, Université Claude Bernard Lyon 1, Ecole Normale Supérieure de Lyon, F-69100 CNRS, France.</footnote>

<footnote>E-mail: coltice@univ-lyon1.fr</footnote>

<footnote n="2">²Centre for Earth Evolution and Dynamics (CEED), Department of Geosciences, University of Oslo, 1028 Blindern, NO-0315 Oslo, Norway</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 1. Flow chart of the methodology used to generate the fully dynamic convection flows leading to the tectonic predictions.</figcaption></figure>

models of mantle convection have improved significantly through a better description of the rheology of the lithosphere (Trompert & Hansen 1998; Moresi & Solomatov 1998; Tackley 1998). The level of precision and sophistication is not at that of regional lithospheric models, but already allows for the localization of stress and strain in narrow regions surrounding stiff and coherent areas. The pseudo-plastic approximation produces plate-like behaviour self-consistently over a restricted range of parameters (van Heck & Tackley 2008; Foley & Becker 2009). Such models reveal the dynamic origin of some fundamental properties of plate tectonics on Earth at the present day, such as the size distribution of plates (Mallard et al. 2016) and the seafloor age versus area distribution (Coltice et al. 2012, 2013). However, their potential for tectonic predictions and reconstruction remains unexploited. Only Yoshida (2014) has explored the conditions required for Pangea breakup, with limited success. Indeed, uncertainties in the initial temperature field 200 My ago together with the intrinsic limit of predictability of mantle convection (Bello et al. 2014) restrict the possibility to realistically simulate the breakup of Pangea.

The following work presents tectonic predictions of instantaneous and dynamic evolution of 3-D spherical models of convection with plate-like behaviour. The goal is to explore the conditions of these models to reproduce plate boundaries and surface velocities of the Earth. Model errors and uncertainties on initial conditions play different roles whether instantaneous or dynamic predictions are considered.

## 2 METHOD

In this section, we detail how we generate the predictions of tectonic structures and kinematics (see also flow chart in Fig. 1). We use 3-D spherical models of mantle convection with plate-like behaviour, but at lower convective vigour than the mantle so it can be computationally tractable. First, we produce a guess of the thermal evolution of the mantle through imposing plate motions at the surface of the model. Then, we compute instantaneous and time-dependent flows starting from the guessed thermal states, without imposing any additional plate structure. Then, we analyse the deformation at the surface of the models in terms of plate boundaries and kinematics.

### 2.1 Physical and numerical model

We model the evolution of temperature, pressure and flow velocity in the Earth's mantle by an approximation of its dynamics. Numerical solutions of the equations of conservation of mass, momentum and energy, and advection of material properties are computed, together

with a pseudo-plastic rheology and a Boussinesq approximation for the equation of state. The physics of phase changes, compressibility, melting and deep dense chemical anomalies are neglected and the rheology is simplified. Such a model is already at the limit of current computational capabilities. Computing the guess of the thermal evolution, once parameters were fixed, took about two months on a supercomputer.

We use the code StagYY (Tackley 2008) to solve the set of equations in a 3-D spherical geometry over a Yin-Yang grid (Kageyama & Sato 2004). StagYY handles several orders of magnitude of viscosity contrasts between adjacent nodes (Tackley 2008) and has been benchmarked for pseudo-plasticity in 2-D (Tosi et al. 2015). The average resolution is 30 km, refined in the vertical direction close to boundary layers of up to 10 km, the lateral resolution being 35 km at the surface and 19 km at the core-mantle boundary. Improving the average resolution to 20 km produced consistent results in the dynamic predictions over 30 My of evolution. Viscosity increases with depth by a factor of 20 according to an activation volume. We impose a viscosity jump by a factor of 30 at 660 km, consistent with the viscosity structure of the Earth inferred from geoid anomalies (Ricard et al. 1993). An additional viscosity increase at around 1000 km depth has been proposed (Rudolph et al. 2015), but is not incorporated here. Uncertainties in the radial viscosity structure translate into errors in the modeling of deep mantle heterogeneity, especially in the sinking rate of slabs.

Viscosity is temperature-dependent:

$$ \eta(z, T) = \eta_0(z) \exp \left( \frac{E_a}{RT} \right), $$

with an activation energy $E_a$ of 142 kJ mol$^{-1}$. R is the gas constant and T the absolute dimensional temperature. Accounting for the full complexity of mantle rheology (King 2016) in such 3-D spherical models is a computational challenge, since extreme viscosity contrasts are difficult to resolve accurately.

The non-dimensional reference viscosity of 1 corresponds to a non-dimensional temperature of 0.64 at zero pressure. This value is chosen before the calculation is realized such as to correspond to the expected temperature at the base of the upper boundary layer. We set a cut-off for the maximum value of the non-dimensional viscosity at $10^4$ to limit viscosity variations. As a consequence, the viscosity contrast across the upper boundary layer is expected to be $10^4$, before the calculation is performed. After the calculation, the average value of the non-dimensional temperature at the base of the upper boundary layer is 0.75, that is, hotter than expected a priori. However, it is stable in the initial stage without imposed plate motions and in the stage with imposed plate motions (see the next subsection). Hence, the typical non-dimensional viscosity
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Viscosity (left) and temperature (right) profiles within the final snapshot of the convection model, in non-dimensional units. The viscosity profiles correspond to that generated by the geotherm on the right, and to the horizontally averaged viscosity of the last snapshot of the convection reconstruction (stiff slabs dominate the average).</figcaption></figure>

in the upper mantle (except in slabs) is around $10^{-1}$ as seen from Fig. 2.

We consider a stress dependence of the viscosity through a pseudo-plastic approximation in order to produce plate boundaries surrounding strong plate interiors (see, for instance Rolf *et al.* 2012). This choice leads to stiff slabs and one-sided subduction with imposed plate kinematics, as described by Bello *et al.* (2015). Viscosity also depends on the type of material, which is tracked with markers. We use three types of materials. Ambient mantle corresponds to the largest fraction of the spherical shell. Continental nuclei are 175 km thick, approximating continental shields (Fig. 3.) They are buoyant, with their buoyancy number being -0.4 ($200 \text{ kg m}^{-3}$ lighter than underlying mantle). They are 100 times more viscous than ambient mantle and their non-dimensional yield stress is 10 times larger than ambient mantle. The continental lithosphere that immediately surrounds the continent nuclei are 115 km thick and their buoyancy number is -0.3 ($150 \text{ kg m}^{-3}$ lighter than underlying mantle). They are 50 times more viscous than underlying mantle and they have 10 times larger yield stress. The Tibetan region of Eurasia, prior to collision, is similarly thick and buoyant as the surrounding belts. This specific continental block is modeled here by 50 times more viscous material, but 2.5 times larger yield stress than ambient mantle. The goal here is to parametrize efficient ductile deformation

<table><caption>**Table 1.** Non-dimensional and dimensional parameters of the reference convection model, also used to generate the Rayleigh number.</caption><thead><tr><th>Parameter</th><th>Non-dimensional value</th><th>Dimensional value</th></tr></thead><tbody><tr><td>Rayleigh number</td><td>10<sup>6</sup></td><td></td></tr><tr><td>Heat production rate</td><td>20</td><td>4.6 × 10<sup>-12</sup> W kg<sup>-1</sup></td></tr><tr><td>Surface temperature</td><td>0.12</td><td>255 K</td></tr><tr><td>Basal temperature</td><td>1.12</td><td>2390 K</td></tr><tr><td>Reference density</td><td>1</td><td>4400 kg m<sup>-3</sup></td></tr><tr><td>Thermal expansivity</td><td>1</td><td>4.5 × 10<sup>-5</sup> K<sup>-1</sup></td></tr><tr><td>Thermal diffusivity</td><td>1</td><td>10<sup>-6</sup> m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>Thermal conductivity</td><td>1</td><td>4 W m<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>Reference viscosity</td><td>1</td><td>10<sup>23</sup> Pa s</td></tr><tr><td>Viscosity jump factor at 660 km</td><td>30</td><td></td></tr><tr><td>Activation energy</td><td>8</td><td>142 kJ mol<sup>-1</sup></td></tr><tr><td>Yield stress at the surface</td><td>2 × 10<sup>4</sup></td><td>230 MPa</td></tr><tr><td>Yield stress depth derivative</td><td>2.5 × 10<sup>5</sup></td><td>1030 Pa m<sup>-1</sup></td></tr><tr><td>Continent nuclei viscosity factor</td><td>100</td><td></td></tr><tr><td>Continent nuclei yield stress</td><td>2 × 10<sup>5</sup></td><td>2300 MPa</td></tr><tr><td>Buoyancy number for continent nuclei</td><td>-0.4</td><td></td></tr><tr><td>Continent belts viscosity factor</td><td>50</td><td></td></tr><tr><td>Continent belts yield stress</td><td>2 × 10<sup>5</sup></td><td>2300 MPa</td></tr><tr><td>Buoyancy for continent belt</td><td>-0.3</td><td></td></tr><tr><td>Tibet viscosity factor</td><td>50</td><td></td></tr><tr><td>Tibet yield stress</td><td>5 × 10<sup>4</sup></td><td>590 MPa</td></tr><tr><td>Buoyancy number for Tibet</td><td>-0.3</td><td></td></tr><tr><td>Weak crust viscosity factor</td><td>0.1</td><td></td></tr><tr><td>Weak crust yield stress</td><td>2 × 10<sup>3</sup></td><td>23 MPa</td></tr><tr><td>Buoyancy number for weak crust</td><td>0.</td><td></td></tr><tr><td>Maximum viscosity cut-off</td><td>10<sup>4</sup></td><td>10<sup>27</sup> Pa s</td></tr></tbody></table>

during the collision (Zhang *et al.* 2004). The physical parameters of the model are listed in Table 1.

The solution is computed with an energy contribution from the core of 25 per cent of the total surface heat flux, the rest being internal heating. Both the surface and the bottom are isothermal, defining the temperature drop for the Rayleigh number $Ra$ of $10^6$, based on the reference viscosity defined above. The effective Rayleigh number based on averaged viscosity is here $5.9 \times 10^6$. The average surface velocity obtained with these physical parameters at statistical steady state, without imposing surface velocities, is 1.2 cm yr$^{-1}$ when scaled with a thermal diffusivity of $10^{-6} \text{ m}^2 \text{ s}^{-1}$. This is a factor of three lower than the Earth today. Unfortunately, computational cost limits the study to a lower $Ra$ than that which would produce Earth-like velocities. Since convective velocities are proportional to $Ra^{2/3}$, this factor of three suggests that increasing $Ra$ by a factor of five would generate appropriate Earth-like velocities with our approximation and keeping our dimensional value of

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Update of the shape of continents at 80 Ma. Shape of the continent boundaries and continent nuclei (in purple) at 80 Ma before (left) and after (right) the update of their shape.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Selected 3-D view state of the model corresponding to present day after Step 2 in the flow chart. Continental material is highlighted in yellow. South America is visible on the front side. The cold isotherm surface in blue (non-dimensional temperature 0.6) visualizes downwelling currents. The hot isotherm surface in red (non-dimensional temperature 0.9) shows plumes coming from the base of the model.</figcaption></figure>

thermal diffusivity. Another consequence of our low Rayleigh number is that convective structures are larger than for the Earth. The dimensional time is then scaled: dimensional velocities produced by the model are multiplied by three and the model time is divided by three, so that the values of velocities and time/age can directly be compared to the Earth for practical purposes.

## 2.2 Building guessed temperature fields with a convection reconstruction

The goal here is to build guessed temperature fields at 30, 20, 10 and 0 Ma using a numerical model of convection and plate reconstructions as information on the state of the mantle today and in the past. We use the methodology explained in more detail in Coltice et al. (2017) and illustrated in the flow chart (Fig. 1): (Step 1) we build a temperature field for the continent configuration at 200 Ma based on free convection with imposed and fixed continent configuration, (Step 2) we impose plate velocities as boundary conditions of the numerical model between 200 and 30, 20, 10 and 0 Ma in increments of 1 My, updating the continent shapes at 80 Ma to account for the moderate changes which happened in terms of continental growth and deformation (Fig. 3). We use the plate reconstructions of Seton et al. (2012), but since we performed the computations presented here, Müller et al. (2016) have published updates and improvements. Because convection in our model is less vigourous than on Earth,

the imposed velocities at present-day are scaled to be consistent with the convective vigour of our model (Bello et al. 2015): the rms value of imposed present-day velocities equals the rms surface velocity of the model without imposed kinematics. Imposing plate motion history generates artificial stresses at the surface, contrary to more realistic free-slip boundary conditions (Lowman 2011). A 3-D snapshot of the thermal state of the reconstruction at 0 Ma is depicted in Fig. 4.

In the following paragraphs, we compare the lateral temperature anomalies of the convection model at present day to seismic anomalies in tomographic models. Such a comparison is limited because seismic velocity is dependent on the local mineralogy and physical properties, rather than temperature alone, but is an appropriate first-order comparison for evaluating predicted mantle structure. Our model does not explicitly take into account for phase equilibria, melting and variable mantle chemistry. Therefore, regions where chemical anomalies in seismic velocities cannot be reproduced. Water and phases changes contribute substantially to seismic anomalies in the transition zone (Tauzin et al. 2013, for instance). Close to the core–mantle boundary regions, a combination of thermal and compositional effects results in broad regions of seismic velocity anomalies (Garnero et al. 2016). Considering these issues, we first compare the power spectrum of the tomographic model S40RTS (Ritsema et al. 2011) to that of the present-day velocity anomaly structure generated by our convection model. Because of the limitations explained above, we simply multiply the temperature
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Left: power spectrum of the seismic anomalies for the tomographic model S40RTS. Centre: power spectrum of the tomography-filtered velocity anomalies (proportional to temperature anomalies) in the convection model at present day. Right: correlation between the two fields. The amplitude of the power spectra is in logscaled.</figcaption></figure>

anomalies in our model by a factor of -0.4 per cent/100 K to obtain the *S* velocity anomaly field. We chose S40RTS because we could apply its resolution operator, hence taking into account the lower resolution of the tomographic model, the uneven distribution of earthquakes and seismic stations, and the parametrization of the tomographic inversion, as in Davies *et al.* (2012) and Koelemeijer *et al.* (2017). Since the resolution of the convection model is substantially finer than that of S40RTS (by more than a factor of 10), we refer to Coltice *et al.* (2017) for a discussion of structures of wavelength smaller than 1000 km, that is, harmonic degree >40. Both power spectra show strong degree two and strong degrees <10 in the upper mantle. The principal disagreements we interpret are within the deepest mantle and the transition zone, where degree 2 heterogeneities are very strong in S40RTS (see Fig. 5, left and centre). Indeed, the convection model does not involve deep chemical anomalies that are suspected to generate a strong seismic signature in the lower 1000 km of the mantle (Garnero *et al.* 2016, for a review). The convection model does not account for phase changes, mineralogical complexity (Nakagawa *et al.* 2012) and the water cycle (Richard *et al.* 2002), that would all otherwise produce seismic anomalies in the transition zone. Also, the stronger power in odd degrees in the deep mantle compared to S40RTS are expected from the filtering because of the tomographic models normal-mode data, giving strong power in even degrees, have a substantial weight. In the spectrum of the convection model, the temperature field displays a peak at long and intermediate wavelengths around 1500 km depth, which corresponds to the region where slabs start to fold and accumulate. This feature could change if we would take compressibility and phase transitions into account. The correlation between the filtered velocity anomalies computed from the convection model and S40RTS is high for the degree 2, and decreases with depth (see Fig. 5, right). The latter is expected, since the shallow mantle is more influenced than the deep mantle by the imposed plate motions, and because compressible effects, stronger as pressure increases, are not taken into account in our calculations.

We compare the location of slabs in the convection model to fast seismic anomalies in tomographic models. But tomographic models substantially differ: some are based on *S* waves, some on *P* waves; they use different 1-D reference model, seismic sources,

seismograms and picking of phases in seismograms; some use finite-frequency approximation and some ray theory only; they use different inversion domain decompositions, methods and parametrizations of the physics. Therefore, we use the vote map description of Shephard *et al.* (2017), for fast and slow seismic anomalies. The number of votes at a given location corresponds to the number of models in which a seismic velocity anomaly faster than the average of fast anomalies at a given depth is present. Shephard *et al.* (2017) described a method for fast seismic anomalies, which we extend to slow velocity anomalies. As a consequence, this tool provides the robust features of 14 tomographic models, seven for *P* waves (Montelli *et al.* 2006; Amaru 2007; Houser *et al.* 2008; Simmons *et al.* 2010, 2012; Burdick *et al.* 2012; Obayashi *et al.* 2013) and seven for *S* waves (Grand 2002; Montelli *et al.* 2006; Houser *et al.* 2008; Simmons *et al.* 2010; Ritsema *et al.* 2011; Auer *et al.* 2014; French & Romanowicz 2014).

Fig. 6 shows horizontal slices at depths of 500, 1500 and 2500 km. At 500 km, robust fast anomalies correspond to the cold sinking slabs in the convection model. Some robust cold anomalies beneath Africa do not correspond to strong cold features in the convection model. The slow robust anomalies which are not associated with plumes do not correspond to any features in the convection model. One possibility is that the slow features represent chemical heterogeneities. At 1500 km depth, the agreement between robust fast anomalies and cold slabs is weaker. For instance, below North America, the position of the Farallon slab in the model is ~1000 km west of that in the vote map. This is a common feature of such convection models, in which low-angle subduction is sometimes difficult to obtain (Bunge & Grand 2000). Another source of error can come from the radial viscosity distribution in our model, because it dictates how fast slabs sink in the lower mantle (Butterworth *et al.* 2014). At 2500 km depth, the disagreement is stronger. At this depth, the model lacks chemical heterogeneity, which is thought to be the source of the large slow velocity provinces, clearly seen on the corresponding vote map. The deepest structure in the convection model suffers the most from the approximation in initial conditions, hypothesis of incompressibility and from uncertainties of past subduction locations in plate reconstructions.
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Comparison between slices in the convection model and vote maps computed from 14 tomographic models. Left-hand column: maps at 500, 1500 and 2500 km depth of the non-dimensional temperature anomalies in the convection models. Central column: vote maps at the same depth for fast seismic anomalies in seven tomographic models of $V_s$ and seven tomographic models of $V_p$. Right-hand column: vote maps at the same depth for slow seismic anomalies in seven tomographic models of $V_s$ and seven tomographic models of $V_p$.</figcaption></figure>

Fig. 7 shows cross-sections for the Farallon, Tonga and Tethyan slabs. The Farallon slab is continuous in the convection model, but its dip angle seems to low compared to the vote map of fast anomalies. Therefore, the convection model predicts an erroneous cold structure below North America and East Atlantic in the lower mantle. The Tonga slab shows some similar patterns in both the convection model and vote maps of fast anomalies. However, the slab is broken into different pieces in the convection model, and sinks as isolated chunks. We attribute this artefact to the method of imposing plate motions. Imposing velocities at the surface of convection models violates the free-slip constraint, generating tangential stresses at the boundary (Nettelfield & Lowman 2007). These velocity gradients can breakup downwellings into several pieces at the trench, especially in intraoceanic domain because both sides of the subduction can yield (Bello et al. 2015). Below India, the location and geometry of the Tethyan slab in the convection model matches that expected from the vote map of fast anomalies. The slow seismic anomalies restricted to the transition zone do not correspond to hot anomalies in the convection model.

Overall, the computed temperature fields involve intrinsic errors. Convection structures are too thick (by a factor of two) because of the convective vigour being lower than that of the Earth. Also, the geometry of slabs is consistent with tomography models in the upper mantle but at the first order only, because of artificial breakoffs. The position of slabs is less accurate, relative to that of tomographic models, as the depth increases. The location of plumes in the numerical solution does not necessarily correspond to hotspots on Earth (see Fig. 8) because plumes emerge freely from the basal boundary

layer without a priori constraint. Finally, the deep mantle thermal structure retains a memory of the initial temperature field chosen at 200 Ma, which is uncertain. Errors therefore come from uncertainties and approximation of (1) the physics of the model, (2) initial conditions and (3) imposed plate kinematics. Therefore, we limit the prediction time frame to 30 My.

## 2.3 Instantaneous and dynamic predictions

We compute instantaneous flows in response to the guessed temperature fields provided by the convection reconstruction. We do not impose mechanically any pre-existing plate boundaries or surface velocities. Continents are the only pre-existing structures that exist in the models. In the relevant models, a 15 km weak crust at the surface of the ocean floor may also be incorporated. The weak crust is constantly created and disappears when it sinks into the mantle below 300 km depth. The viscosity and yield stress of the weak crust are 10 times lower than that of ambient mantle (see Table 1). It approximates hydrothermally altered rocks that are softer because of the presence of hydrated silicates like chlorite, amphibole and serpentine. The viscosity and the yield stress of this layer are set to 0.1 times the values of the ambient mantle. Such a layer is fundamental to the development of asymmetric subduction (Gerya et al. 2008; Crameri & Tackley 2014, 2015). It is here thicker than expected on Earth because the model has a lower Rayleigh number, hence thicker structures.

We also compute time-dependent convection evolution forward in time using guessed thermal states at 30 and 10 Ma as initial
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Comparison between cross-sections through the convection model and the votes computed from 14 tomographic models. Left-hand column: non-dimensional temperature anomalies in the convection model. Central column: votes for fast seismic anomalies in 14 tomographic models. Right-hand column: votes for slow seismic anomalies in 14 tomographic models. Top: cross-sections of the Farallon slab below South California. Middle: cross-sections of the Tonga slab. Bottom: cross-sections of the slab below the Himalayan collision.</figcaption></figure>

conditions. The system is chaotic: model and initial condition errors propagate in time (Bello et al. 2014, 2015). In test cases, Bocher et al. (2016) showed that the interval between corrections in a sequential data assimilation scheme (using surface velocities and seafloor age distribution as the data to match) has to be ≤15 My for accurate inversions of the convective temperature field. Therefore, we limit the prediction time frame to 30 My.

To study the role of the viscosity parameters, we compute numerical solutions for the instantaneous and dynamic models for (1) the same physical parameters as the convection reconstruction, (2) the same as the reference but with a lower yield stress (10⁴, i.e. 115 MPa in dimensional units) for the ambient mantle and (3) the same as the reference but with the weak crust.

To evaluate the quality of the predictions, the viscosity field just below the surface (5 km) is compared with the plate boundaries of the plate model used for the convection reconstruction (Seton et al. 2012). We also compare the kinematics emerging from the numerical model with that of the plate model, computing the mean-squared error on the velocity field:

$$ \text{MSE} = \frac{1}{N} \sum_{i=1}^{N} (\vec{V}(\mathbf{x}_i, T) - \vec{V}_{\text{plates}}(\mathbf{x}_i, T)) \cdot (\vec{V}(\mathbf{x}_i, T) - \vec{V}_{\text{plates}}(\mathbf{x}_i, T)), $$

where N is the number of nodes (414,144), $\vec{V}(\mathbf{x}_i, T)$ the predicted velocity vector at position $\mathbf{x}_i$ and age $T$, $\vec{V}_{\text{plates}}(\mathbf{x}_i, T)$ the velocity vector in the plate model (Seton et al. 2012). We note MSE, the tectonic mean-squared error which measures the mean-squared

difference between the average and plate velocities. Therefore, it is exactly the mean-squared plate velocity in the no-net rotation reference frame (the average velocity vector being the null vector).

# 3 RESULTS

## 3.1 Instantaneous predictions

We compute instantaneous flows in response to the reconstructed temperature fields at present day for the three parametrizations of the viscosity described above. Fig. 9 shows the surface viscosity fields and kinematics of the three solutions, compared to the plate tectonic reconstruction at present day. The three models show plate-like behaviour with 90 per cent of the deformation being concentrated in 11, 10 and 8 per cent of the surface for the low yield stress, reference and weak crust models, respectively. In the models, the network of very low (<10⁻¹) viscosity bands corresponds to the plate boundaries emerging from the model. In the three models, ridges located away from trenches match the plate reconstructions. But ridges in backarc basins do not emerge, or not at the right places. The location of trenches is also consistent with those of the Earth when subduction occurs below a continent. Intraoceanic trenches are less accurately predicted close to New Zealand, Japan and the Caribbean. The model with the weak crust produces the strongest viscosity contrast between plate interiors and boundaries. The model with the low yield stress produces a slightly more diffuse viscosity distribution, because yielding may occur over a broader area of high stresses. Overall, the layout of large plates self-consistently emerges when imposing this temperature field, as long as pseudo-plasticity is
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Temperature anomalies at 370 km depth and plume locations in the model at present day after Step 2 in the flow chart. Non-dimensional temperature anomalies respective to the laterally averaged temperature. The black triangles represent the location of plumes reaching the top boundary layer in the model. The white triangles represent hotspot locations from the GPlates 2.0 database (Whittaker *et al.* 2015). The model was not designed to produce plumes at the same location as on Earth. The elongate negative anomalies represent the location of subducted slabs. Small-scale convection below continents and ocean forms networks of alternating positive and negative anomalies (green and yellow-orange colours).</figcaption></figure>

introduced with the strong temperature dependence of the viscosity. The layout of small plates does not emerge here, whatever the viscosity parametrization.

The same figure shows the differences between the predicted and expected plate velocities of Seton *et al.* (2012). To the first order, the predicted velocity directions and magnitudes are consistent with the expected ones. As shown in Fig. 10, the lower value of MSE/MSE<sub>t</sub> is for the model with weak crust, being 0.32 (equivalent to the difference between plate velocities at 10 Ma and at present day), while it is 0.39 for the low yield stress model and 0.66 for the reference. MSE/MSE<sub>t</sub> for instantaneous flows produced with the

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Ratio of the mean-squared error (MSE) of the computed velocity field (relative to the expected values for the Earth at 0, 10, 20 and 30 Ma) to the mean-squared error of the velocity of the tectonic reconstruction MSE<sub>t</sub> (i.e. mean-squared surface velocity in the no net rotation reference frame). Open circles represent the instantaneous flow calculations, while filled squares represent the fully dynamic evolution.</figcaption></figure>

weak crust model modestly increases with the age of the convection reconstruction within the past 30 My. Some specific plates have systematically lower predicted velocities than expected: the Pacific, Nazca and Indian plates. The model with weak crust produces the highest velocities for these domains. The model with lower yield stress displays the stronger errors on velocity directions (15°) for the Pacific. However, the directions of the Nazca Plate are more accurate for this latter model than the others.

Fig. 8 shows the residual temperature at 370 km depth in the model together with the location of 21 plumes emerging from the reconstructed flow. These plumes emerge at locations that are not

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Viscosity field and kinematics of instantaneous flow models versus plate boundaries and kinematics on Earth, at present day. Top row: viscosity field at 10 km depth emerging from the instantaneous flow calculation, and the expected plate layout of the Earth, as indicated by plate boundaries in black and based on the reconstruction of Seton *et al.* (2012). The reference model is in the middle column. The model with a factor of two lower yield stress in on the left, and the model with a weak crust is on the right. Bottom row: for the same models, black arrows represent model velocities and green arrows represent the expected velocities, as derived from the plate reconstruction of Seton *et al.* (2012).</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Viscosity field and kinematics of dynamic models started at 10 Ma vesrus plate boundaries and kinematics on Earth. Top row: viscosity field at 10 km depth emerging from the calculation after 10 My of evolution, and the expected plate layout of the Earth. The model having a factor of two lower yield stress in on the left, and the model with a weak crust is on the right. Bottom row: for the same models, black arrows represent model velocities and green arrows represent the expected velocities at present day.</figcaption></figure>

imposed and therefore do not necessarily match those on Earth. However, they often correspond to regions of existing hotspots although the impact of deep chemical heterogeneities on plume onset is not taken into account. Indeed, the structure of downwellings already strongly constrains the onset locations of plumes (Davies & Davies 2009). The errors in the predicted plate boundaries and velocities do not correlate with the presence of plumes nearby or in terms of the numbers of plumes beneath a plate.

## 3.2 Dynamic flow predictions

We compute a dynamic model evolution starting from the convection reconstruction at 10 Ma. From 10 to 0 Ma, the flow is self-organized and we do not impose any plate boundaries or tectonic constraints. After 10 My of evolution, Fig. 11 shows the present-day viscosity field at the surface and the predicted kinematics for the low yield stress model and the model with weak crust. Both models show ridges at the expected locations except in backarc basins. The major discrepancy comes from the North Atlantic ridge, which is no longer a ridge after 10 My of evolution, but rather a shear band localizing incipient convergence (Fig. 11). The model with a weak crust still displays the ridges surrounding the Bauer Plate close to the East Pacific Rise, while they should stop spreading. The Chile Ridge is progressively fading out in both models. Trenches are located at, or close to the expected locations. Backarc basins develop in the western Pacific, but with differences in plate boundary locations relative to the Earth. The plate boundaries in these regions differ from one model to the other, the weak crust model displaying sharper bands of low viscosity and smaller plates.

The kinematics of both models show similar errors in terms of velocity direction and amplitude for most plates. The direction of the Pacific is off by <20° for both models, but the model with weak crust predicts faster velocities, which are more consistent with the observations. The velocities of Africa and Antarctica are larger than expected for the Earth, especially for the weak crust model.

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Erroneous prediction of slab breakoff under Cascadia. The residual non-dimensional temperature field at 0 Ma for the dynamic evolution started at 10 Ma. The slab has broken off from the surface, which is not expected on Earth as based on seismic models.</figcaption></figure>

Predicted kinematics for North America is the major issue of both models. The direction is more than 90° off, leading to a closing of the North Atlantic ocean basin. It comes from the breakoff of the slab as seen in the cross-section (Fig. 12). It profoundly modifies the kinematics beyond the region whatever the rheological parameters. The value of MSE/MSE~t~ at the final time is more than four times the initial value (1.2 and 1.87, respectively) for the weak crust model and the low yield stress model, respectively.
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Viscosity field and kinematics of the dynamic models started at 30 Ma versus plate boundaries and kinematics on Earth at the corresponding time steps. Left-hand column: viscosity field at 10 km depth emerging from dynamic evolution of the model with weak crust, and the expected plate layout of the Earth over the time evolution based on plate boundaries from the model of Seton *et al.* (2012, black lines). Right-hand column: black arrows represent model velocities and green arrows represent the expected velocities over the time evolution.</figcaption></figure>

We compute a longer dynamic evolution for the weak crust model, which has the lower MSE/MSE, for both the instantaneous and 10 My evolution tests. The numerical solution corresponds to a free evolution started from the initial condition set by the convection reconstruction at 30 Ma, as depicted in Fig. 13. Over this time, the predictions of the locations of several plate boundaries degrade quickly. Only the South Atlantic ridge and the South Indian ridges remain precise, moving in the appropriate directions. The Galapagos ridge initiates as expected but further south of the location on Earth. The India-Eurasia collision continues, thanks to the low resistance of the Tibet block, and subduction on the West Pacific operates as well as under South America. However, subduction under North America quickly stops, because of the early breakoff (between 30 and 20 Ma) of the slab as for the 10 My dynamic evolution. Again, the North Atlantic starts to be in compression after the breakoff, shutting down the ridge system. Also, the subduction system north and east of Australia quickly retreats until it reaches the ocean-

continent boundary, instead of remaining at a similar position in the expected plate layout. As for the preceding calculations, backarc basins are generated with rapidly evolving ridge systems in connection with the moving trench. However, the small plate pattern does not match the expected one on Earth.

The predicted kinematics show a progressive 20° change of direction of the Pacific Plate towards the south, while it is expected to remain constant on Earth. The direction of the Australian Plate also changes direction progressively to reach a 30° offset towards the east, leading to the opening of a ridge system south of Southeast Asia. These changes of directions correlate with the retreat of the trench in the South-East Pacific described above, modifying the force balance on the Pacific and Australian plates that are converging. As for the 10 My evolution, the North American motion is quickly inconsistent with Earth evolution, before changing back again at the end to produce kinematics more consistent with the expectations. However, the relative motion between North America
and Eurasia still corresponds to a slowly converging boundary in-
stead of a slowly diverging one. The MSE/MSE, in Fig. 10 quickly
grows as for the 10 My model, and stabilizes at about twice its initial
value, and four times the value of the instantaneous flow calculation
at 0 Ma. The change of direction of the Pacific and Australian plates,
as well as the incorrect kinematics of North America, produce the
early peak of errors because of inaccurate trench evolution (fast re-
treat in the South East of the Pacific and slab breakoff under North
America).

4 DISCUSSION

In this study, we compute first a reconstruction of convection in
the mantle consistent with the physics and approximations used
for the subsequent instantaneous and dynamic predictions. Most of
the limitations are caused by computational power that is not yet
sufficient to reach more realistic parametrizations of the physics.
From reconstructed thermal fields, we compute instantaneous flows
where plate boundaries and surface kinematics are not prescribed.

The plate layouts emerging from these flows are consistent with
the ones expected for the Earth, except close to subduction zones
where the plate fragmentation does not produce the observed plate
boundaries. A substantial decrease of the yield stress or a weak crust
at the surface of the ocean floor has a minor impact on the resulting
plate configuration. The predicted kinematics follow the same con-
clusions for the instantaneous models: velocities have directions and
magnitudes close to what is expected on Earth. Discrepancies are
again related to selected subduction regions: the Pacific and Nazca
plates are slower in the prediction that expected, while they are of
the correct magnitude elsewhere. Introducing a weak crust speeds
up these plates, by reducing the coupling between the sinking and
upper plates. The direction of the Nazca Plate can slightly vary with
rheological parameters, but by an angle <30°. These results are
confirmed for instantaneous calculations at 30, 20, 10 and 0 Ma.
Therefore, surface kinematics and plate boundary emergence are
first-order outcomes of the temperature field in these models. The
rheological parameters are second order. Extreme perturbations of
the rheological parameters used to build the guessed temperature
fields would certainly change this result, but would be inconsistent
with the approach we develop, which aims at keeping consistent
physics for both guessing initial conditions and realizing predic-
tions.

A clear observation is that plumes have no influence on the instan-
taneous kinematics and plate boundaries here. They neither produce
effortless plate boundaries nor alter surface kinematics. The vis-
cosity contrast (6 orders of magnitude here) is so large between the
surface and hot plumes that in most cases they easily spread below
the cold boundary layer, slightly changing their thermal structure
without modifying the force balance as proposed by Monnereau
et al. (1993).

Stadler et al. (2010) and Alisic et al. (2012) worked on models
comparable to the ones presented here since they also incorporated
strong slabs and large lateral viscosity variations. They proposed
similar conclusions: the direction and magnitude of plate veloci-
ties remain consistent varying the rheological parameters, except
for the Nazca Plate and for small plates. These models belong
to a larger class of models, which differ from those presented
in this paper because (1) rigid plates or plate boundaries are im-
posed while they self-consistently emerge in this paper, and (2) the
guessed temperature field at present day derives from conversion
of seismic anomalies or imposed location of slabs in the interior
of the mantle whereas they are outputs of the models here. Within

this class of geodynamic models (i.e. imposed mantle initial con-
ditions and/or plate kinematics), substantial differences in rheo-
logical parametrizations produce successful kinematic predictions.
Ghosh & Holt (2012) predict accurate plate motions from a guess of
the temperature field derived from seismology, taking into account
lateral viscosity variations in the lithosphere and asthenosphere
only. Ricard et al. (1989), Becker & O'Connell (2001) and Conrad
& Lithgow-Bertelloni (2002) also predict accurate plate motions
without lateral variations of viscosity, and with different types of
guessed density inside the Earth's mantle [these types of density
models correlating with each other—see Becker & Boschi (2002)].
Becker & O'Connell (2001) showed that plate motions are mostly
sensitive to the structure of the lithosphere and upper-mantle slabs.
Taking into account the contribution of lower mantle slabs slightly
improves the predictions (Becker & O'Connell 2001; Conrad &
Lithgow-Bertelloni 2002; Alisic et al. 2012). Since all these mod-
els have a diversity of rheological parameters for slabs and the
lithosphere, the results agree with the observation made here that
rheology is second order for the instantaneous predictions of surface
velocities.

The results from the instantaneous predictions contrast with the
dynamical evolution started from guesses of past temperature fields.
The models started at 10 and 30 Ma display discrepancies in slab
evolution that quickly arise within the first 10 My. The trench east
of Australia retreats faster than expected. Considering the presence
of continental Zealandia instead of pure oceanic floor (Mortimer
et al. 2017) would certainly impede the retreat. The subduction
under North America breaks off, whereas it is expected to persist to
the present day on Earth. It is certainly artificially generated by the
errors in the reconstructed temperature field because of the recurrent
chopping off of slabs by imposing plate velocities at the surface. The
breakoff of the Farallon slab, and the low angle of the sinking slab
contract the forces that drag North America westwards. Therefore,
the North Atlantic Ridge starts to localize incipient convergence.
This change of force balance in the East Pacific, combined with
the strong subduction in the west are responsible for the westward
motion of the Pacific Plate instead of being north-westward.

The fast growth of errors comes from feedbacks between errors in
the initial temperature field, which are stronger in the lower mantle
than the upper mantle, and errors of parametrization of the physics.
Unfortunately, the initial temperature field contains errors coming
from (1) errors in the initial condition at 200 Ma (Step 1 of the chart
flow in Fig. 1), (2) errors in physical parameters used for Step 2
(Fig. 1) since, for instance, slab sinking rate depends on the radial
viscosity structure and (3) uncertainties in plate reconstructions. As
yet, we do not have a way in which to correct all these issues, which
all point the deep mantle as the major source of errors.

The lower mantle is also the region where our parametriza-
tion of convection fails the most. Indeed, we neglect compress-
ibility, that is, the decrease of thermal expansivity with pressure
(Chopelas & Boehler 1992). When taken into account, it slows down
slabs, which are consequentially more stagnant (Tosi et al. 2013).
Another limitation of our models is that deep chemical heterogene-
ity is not incorporated. Furthermore, the top of the lower mantle
is also the location of phase transitions. Depending on the density
change and Clapeyron slope of the transitions, mostly at 660 km
depth, sinking slabs can stagnate and lie for some time at a phase
boundary (Christensen & Yuen 1984; Tackley et al. 1993).

Compared to the instantaneous models, dynamic calculations
demonstrate stronger discriminating power for sources of errors
in kinematic predictions. Therefore, they have rich potential for
inversions of rheology and guessed temperature fields, even over
short timescales. Indeed, the initial conditions and the rheological parameters can be good enough for an accurate prediction of instantaneous flow, but not for a prediction after 10 My of evolution. We suggest here that using inversions of dynamical evolution using surface velocities as data constraints rather than inputs should lead to improved rheologies and resulting mantle flow. Methods like sequential data assimilation (Bocher et al. 2016, 2017) and adjoint-based inversions (Li et al. 2017) are under development for that very purpose.

Nonetheless, the dynamical framework we used has strong limitations. The physics is approximated since compressibility is not taken into account, and the rheology is empirical instead of being defined by properties at the mineralogical scale. The vigour of convection is lower than that of Earth, therefore convective structures are probably about twice larger than expected for our planet. Increasing the convective vigour could also increase the time dependence and the chaotic nature of the flow. Most of these limitations are caused by the computational cost of the time-dependent calculations. Parallelization in time could be a solution (Samuel 2012), however, it is then difficult to simultaneously test a variety of initial conditions at 200 Ma and parametrizations of the physics. With all these simplifications, the models presented here already generate tectonics consistent at first order with what is expected, even for the dynamic evolution.

## 5 CONCLUSIONS

We compare the tectonic predictions (kinematics and plate boundary locations) of 3-D spherical convection models with plate-like behaviour with tectonic reconstructions for the Earth. We show that calculation of instantaneous flows generate plate boundaries and kinematics consistent with what is expected for present day and in the past, except for small plates close to subduction zones. Perturbing the rheological parameters does not significantly modify the results although a weaker coupling between subducting plates and continents improves the predictions. Lithosphere structure and upper-mantle slabs overcome rheological approximations and errors in the temperature field of the lower mantle. Plumes and small-scale convection have imperceptible effects on the plate layout and kinematics. The models evolving freely over several tens of million years show a rapid growth of errors. In the models presented here, errors in the guessed past states interact with errors on rheological parameters. These calculations show that short-term (10-30 My) dynamical evolution models are more suitable experiments than instantaneous flow calculations for the inversion of the temperature field and rheological parameters. Such methods based on adjoint codes (Li et al. 2017) and bayesian approaches (Bocher et al. 2016, 2017) are under development.

## ACKNOWLEDGEMENTS

We thank two anonymous reviewers for fruitful comments, questions and suggestions. We thank Paula Koelemeijer for discussions on the tomographic models, and for providing the tools to compute the filtered velocity anomalies in the convection model. We thank Thorsten Becker for encouragements and discussions, and Mélanie Gérault and Claire Mallard for their help on a former version of the manuscript. NC is funded by European Research Council within the framework of the SP2-Ideas Programme ERC-2013-CoG, under ERC grant agreement 617588. GES is funded by VISTA – a basic research program in collaboration between The Norwegian Academy of Science and Letters, and Statoil (project 6268, DEFMOD). GES acknowledges support from the Research Council of Norway through its Centers of Excellence funding scheme, project number 223272. Calculations were performed at P2CHPD Lyon.

## REFERENCES

* Alisic, L., Gurnis, M., Stadler, G., Burstedde, C. & Ghattas, O., 2012. Multiscale dynamics and rheology of mantle flow with plates, *J. geophys. Res.*, **117**, doi:10.1029/2012JB009234.

* Amaru, M., 2007. Global travel time tomography with 3-D reference models, *PhD thesis*, Utrecht University.

* Auer, L., Boschi, L., Becker, T., Nissen-Meyer, T. & Giardini, D., 2014. Savani: a variable resolution whole-mantle model of anisotropic shear velocity variations based on multiple data sets, *J. geophys. Res.*, **119**(4), 3006–3034.

* Becker, T.W. & Boschi, L., 2002. A comparison of tomographic and geodynamic mantle models, *Geochem. Geophys. Geosys.*, **3**(1), doi:10.1029/2001GC000168.

* Becker, T.W. & O'Connell, R.J., 2001. Predicting plate velocities with geodynamic models, *Geochem., Geophys., Geosys.*, **2**(12), doi:10.1029/2001GC000171.

* Bello, L., Coltice, N., Rolf, T. & Tackley, P.J., 2014. On the predictability limit of convection models of the earth's mantle, *Geochem. Geophys. Geosyst.*, **15**, 2319–2328.

* Bello, L., Coltice, N., Tackley, P.J., Müller, R.D. & Cannon, J., 2015. Assessing the role of slab rheology in coupled plate-mantle convection models, *Earth planet. Sci. Lett.*, **430**, 191–201.

* Bercovici, D., 2003. The generation of plate tectonics from mantle convection, *Earth planet. Sci. Lett.*, **205**(3), 107–121.

* Bocher, M., Coltice, N., Fournier, A. & Tackley, P., 2016. A sequential data assimilation approach for the joint reconstruction of mantle convection and surface tectonics, *Geophys. J. Int.*, **204**(1), 200–214.

* Bocher, M., Fournier, A. & Coltice, N., 2017. Ensemble kalman filter for the reconstruction of the earth's mantle circulation, *Nonlin. Process. Geophys.*, in press.

* Bunge, H.-P. & Grand, S.P., 2000. Mesozoic plate-motion history below the northeast pacific ocean from seismic images of the subducted farallon slab, *Nature*, **405**(6784), 337–340.

* Burdick, S. *et al.*, 2012. Model update march 2011: upper mantle heterogeneity beneath north america from traveltime tomography with global and usarray transportable array data, *Seismol. Res. Lett.*, **83**, 23–28.

* Butterworth, N., Talsma, A., Müller, R., Seton, M., Bunge, H.-P., Schuberth, B., Shephard, G. & Heine, C., 2014. Geological, tomographic, kinematic and geodynamic constraints on the dynamics of sinking slabs, *J. Geodyn.*, **73**, 1–13.

* Chopelas, A. & Boehler, R., 1992. Thermal expansivity in the lower mantle, *Geophys. Res. Lett.*, **19**(19), 1983–1986.

* Christensen, U.R. & Yuen, D.A., 1984. The interaction of a subducting lithospheric slab with a chemical or phase boundary, *J. geophys. Res.*, **89**, 4389–4402.

* Coltice, N., Rolf, T., Tackley, P.J. & Labrosse, S., 2012. Dynamic causes of the relation between area and age of the ocean floor, *Science*, **336**, 335–338.

* Coltice, N., Seton, M., Rolf, T., Müller, R. & Tackley, P., 2013. Convergence of tectonic reconstructions and mantle convection models for significant fluctuations in seafloor spreading, *Earth planet. Sci. Lett.*, **383**, 92–100.

* Coltice, N., Larrouturou, G., Debayle, E. & Garnero, E.J., 2017. Interactions of scales of convection in the earth's mantle, *Tectonophysics*, in press.

* Conrad, C.P. & Lithgow-Bertelloni, C., 2002. How mantle slabs drive plate tectonics, *Science*, **298**, 207–209.

* Crameri, F. & Tackley, P.J., 2014. Spontaneous development of arcuate single-sided subduction in global 3-d mantle convection models with a free surface, *J. geophys. Res.*, **119**(7), 5921–5942.

* Crameri, F. & Tackley, P.J., 2015. Parameters controlling dynamically self-consistent plate tectonics and single-sided subduction in global models of mantle convection, *J. geophys. Res.*, **120**(5), 3680–3706.
Davies, D.R. & Davies, J.H., 2009. Thermally-driven mantle plumes rec-
conde multiple hot-spot observations, *Earth planet. Sci. Lett.*, **278**(1),
50–54.

Davies, D.R., Goes, S., Davies, J.H., Schuberth, B., Bunge, H.-P. & Ritsema, J., 2012. Reconciling dynamic and seismic models of earth's lower mantle: the dominant role of thermal heterogeneity, *Earth planet. Sci. Lett.*, **353**, 253–269.

Domeier, M., Doubrovine, P.V., Torsvik, T.H., Spakman, W. & Bull, A.L., 2016. Global correlation of lower mantle structure and past subduction, *Geophys. Res. Lett.*, **43**, 4945–4953.

Foley, B.J. & Becker, T.W., 2009. Generation of plate-like behavior and mantle heterogeneity from a spherical, viscoplastic convection model, *Geochem. Geophys. Geosyst.*, **10**(8), doi:10.1029/2009GC002378.

Forsyth, D.W. & Uyeda, S., 1975. On the relative importance of the driving forces of plate motion, *Geophys. J. R. astr. Soc.*, **43**, 163–200.

French, S. & Romanowicz, B., 2014. Whole-mantle radially anisotropic shear velocity structure from spectral-element waveform tomography, *Geophys. J. Int.*, **199**(3), 1303–1327.

Garnero, E.J., McNamara, A.K. & Shim, S.-H., 2016. Continent-sized anomalous zones with low seismic velocity at the base of earth's mantle, *Nature Geosci.*, **9**, 481–489.

Gerya, T.V., Connolly, J.A. & Yuen, D.A., 2008. Why is terrestrial subduction one-sided?, *Geology*, **36**(1), 43–46.

Ghosh, A. & Holt, W.E., 2012. Plate motions and stresses from global dynamic models, *Science*, **335**, 839–843.

Grand, S.P., 2002. Mantle shear-wave tomography and the fate of subducted slabs, *Phil. Trans. R. Soc. Lond., A: Math. Phys. Eng. Sci.*, **360**(1800), 2475–2491.

Gurnis, M. et al., 2012. Plate tectonic reconstructions with continuously closing plates, *Comput. Geosci.*, **38**(1), 35–42.

Houser, C., Masters, G., Shearer, P. & Laske, G., 2008. Shear and compres-
sional velocity models of the mantle from cluster analysis of long-period
waveforms, *Geophys. J. Int.*, **174**(1), 195–212.

Kageyama, A. & Sato, T., 2004. 'yin-yang grid': an overset grid in spherical geometry, *Geochem. Geophys. Geosyst.*, **5**, Q09005.

King, S.D., 2016. Reconciling laboratory and observational models of man-
tle rheology in geodynamic modelling, *J. Geodyn.*, **100**, 33–50.

Koelemeijer, P., Schuberth, B., Davies, D., Deuss, A. & Ristema, J., 2017. Constraints on the presence of post-perovskite in earth's lowermost mantle from tomographic-geodynamic model comparisons, *Earth planet. Sci. Lett.*, in press.

Li, D., Gurnis, M. & Stadler, G., 2017. Towards adjoint-based inversion of time-dependent mantle convection with nonlinear viscosity, *Geophys. J. Int.*, **209**(1), 86–105.

Lithgow-Bertelloni, C., Richards, M.A., Ricard, Y., O'Connell, R.J. & Engebretson, D.C., 1993. Toroidal-poloidal partitioning of plate motions since 120 Ma, *Geophys. Res. Lett.*, **20**, 375–378.

Lowman, J.P., 2011. Mantle convection models featuring plate tectonic be-
havior: an overview of methods and progress, *Tectonophysics*, **510**(1),
1–16.

Mallard, C., Coltice, N., Seton, M., Müller, R.D. & Tackley, P.J., 2016. Subduction controls the distribution and fragmentation of earth's tectonic plates, *Nature*, **535**, 140–143.

Monnereau, M., Rabinowicz, M. & Arquis, E., 1993. Mechanical erosion and reheating of the lithosphere: a numerical model for hotspot swells, *J. geophys. Res.*, **98**(B1), 809–823.

Montelli, R., Nolet, G., Dahlen, F. & Masters, G., 2006. A catalogue of deep mantle plumes: new results from finite-frequency tomography, *Geochem. Geophys. Geosyst.*, **7**(11), doi:10.1126/science.1092485.

Moresi, L. & Solomatov, V., 1998. Mantle convection with a brittle litho-
sphere: thoughts on the global tectonic styles of the earth and venus,
*Geophys. J. Int.*, **133**, 669–682.

Mortimer, N. et al., 2017. Zealandia: Earth's hidden continent, *GSA Today*, **27**(3), doi:10.1130/GSATG321A.1.

Müller, R.D. et al., 2016. Ocean basin evolution and global-scale plate reorganization events since pangea breakup, *Annu. Rev. Earth planet. Sci.*, **44**, 107–138.

Nakagawa, T., Tackley, P.J., Deschamps, F. & Connolly, J.A., 2012. Radial 1-d seismic structures in the deep mantle in mantle convection simula-tions with self-consistently calculated mineralogy, *Geochem. Geophys. Geosyst.*, **13**(11), doi:10.1029/2012GC004325.

Nettelfield, D. & Lowman, J.P., 2007. The influence of plate-like surface motion on upwelling dynamics in numerical mantle convection models, *Phys. Earth planet. Inter.*, **161**(3), 184–201.

Obayashi, M., Yoshimitsu, J., Nolet, G., Fukao, Y., Shiobara, H., Sugioka, H., Miyamachi, H. & Gao, Y., 2013. Finite frequency whole mantle p wave tomography: improvement of subducted slab images, *Geophys. Res. Lett.*, **40**, 5652–5657.

Ricard, Y., Vigny, C. & Froidevaux, C., 1989. Mantle heterogeneities, geoid, and plate motion: a Monte Carlo inversion, *J. geophys. Res.*, **94**(B10), 13 739–13 754.

Ricard, Y., Richards, M., Lithgow-Bertelloni, C. & Le Stunff, Y., 1993. A geodynamic model of mantle density heterogeneity, *J. geophys. Res.*, **98**, 21 895–21 909.

Richard, G., Monnereau, M. & Ingrin, J., 2002. Is the transition zone an empty water reservoir? inferences from numerical model of mantle dynamics, *Earth planet. Sci. Lett.*, **205**(1), 37–51.

Ritsema, J., Deuss, A., van Heijst, H. & Woodhouse, J., 2011. S40rts: a degree-40 shear-velocity model for the mantle from new rayleigh wave dispersion, teleseismic traveltime and normal-mode splitting function measurements, *Geophys. J. Int.*, **184**, 1223–1236.

Rolf, T., Coltice, N. & Tackley, P., 2012. Linking continental drift, plate tectonics and the thermal state of the earth's mantle, *Earth planet. Sci. Lett.*, **351**, 134–146.

Rowley, D.B., 2008. Extrapolating oceanic age distributions: lessons from the pacific region, *J. Geol.*, **116**, 587–598.

Rudolph, M.L., Lekić, V. & Lithgow-Bertelloni, C., 2015. Viscosity jump in earth's mid-mantle, *Science*, **350**(6266), 1349–1352.

Samuel, H., 2012. Time domain parallelization for computational geodynamics, *Geochem. Geophys. Geosyst.*, **13**(1), doi:10.1029/2011GC003905.

Seton, M. et al., 2012. Global continental and ocean basin reconstructions since 200ma, *Earth-Sci. Rev.*, **113**, 212–270.

Shephard, G., Matthews, K., Hosseini, K. & Domeier, M., 2017. On the consistency of seismically imaged lower mantle slabs, *Sci. Rep.*, **7**, 1–17.

Simmons, N.A., Forte, A.M., Boschi, L. & Grand, S.P., 2010. Gypsum: a joint tomographic model of mantle density and seismic wave speeds, *J. geophys. Res.*, **115**, doi:10.1029/2010JB007631.

Simmons, N.A., Myers, S.C., Johannesson, G. & Matzel, E., 2012. Llnl-g3dv3: global p wave tomography model for improved regional and teleseismic travel time prediction, *J. geophys. Res.*, **117**, doi:10.1029/2012JB009525.

Stadler, G., Gurnis, M., Burstedde, C., Wilcox, L.C., Alisic, L. & Ghattas, O., 2010. The dynamics of plate tectonics and mantle flow: from local to global scales, *Science*, **329**, 1033–1038.

Tackley, P.J., 1998. Self-consistent generation of tectonics plates in three-dimensional mantle convection, *Earth planet. Sci. Lett.*, **157**(1), 9–22.

Tackley, P.J., 2008. Modelling compressible mantle convection with large viscosity contrasts in a three-dimensional spherical shell using the yin-yang grid, *Phys. Earth planet. Inter.*, **171**, 7–18.

Tackley, P.J., Stevenson, D.J., Glatzmaier, G.A. & Schubert, G., 1993. Effects of an endothermic phase transition at 670 km depth in a spherical model of convection in the Earth's mantle, *Nature*, **361**, 699–704.

Tauzin, B., van Der Hilst, R.D., Wittlinger, G. & Ricard, Y., 2013. Multi-
ple transition zone seismic discontinuities and low velocity layers below
western united states, *J. geophys. Res.*, **118**(5), 2307–2322.

Tosi, N., Yuen, D.A., de Koker, N. & Wentzcovitch, R.M., 2013. Mantle dynamics with pressure-and temperature-dependent thermal expansivity and conductivity, *Phys. Earth planet. Inter.*, **217**, 48–58.

Tosi, N. et al., 2015. A community benchmark for viscoplastic thermal convection in a 2-D square box, *Geochem. Geophys. Geosyst.*, **16**, 2175–2196.
Trompert, R. & Hansen, U., 1998. Mantle convection simulations with rheologies that generate plate-like behaviour, *Nature*, **395**, 686–689.

van Der Meer, D.G., Spakman, W., van Hinsbergen, D.J., Amaru, M.L. & Torsvik, T.H., 2010. Towards absolute plate motions constrained by lower-mantle slab remnants, *Nature Geosci.*, **3**, 36–40.

van Heck, H. & Tackley, P., 2008. Planforms of self-consistently generated plates in 3d spherical geometry, *Geophys. Res. Lett.*, **35**(19), doi:10.1029/2008GL035190.

Whittaker, J., Afonso, J., Masterton, S., Müller, R., Wessel, P., Williams, S. & Seton, M., 2015. Long-term interaction between mid-ocean ridges and mantle plumes, *Nature Geosci.*, **8**(6), 479–483.

Yoshida, M., 2014. Effects of various lithospheric yield stresses and different mantle-heating modes on the breakup of the pangea supercontinent, *Geophys. Res. Lett.*, **41**, 3060–3067.

Zhang, P.-Z. et al., 2004. Continuous deformation of the tibetan plateau from global positioning system data, *Geology*, **32**, 809–812.