# Tectonic predictions with mantle convection models

Tectonic predictions with mantle convection models 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. 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 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. Figure 1. View largeDownload slide Flow chart of the methodology used to generate the fully dynamic convection flows leading to the tectonic predictions. Figure 1. View largeDownload slide Flow chart of the methodology used to generate the fully dynamic convection flows leading to the tectonic predictions. 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:   \begin{equation*} \eta (z,T)=\eta _0(z)\, \exp {\left(\frac{E_a}{RT}\right)} , \end{equation*} with an activation energy Ea 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 104 to limit viscosity variations. As a consequence, the viscosity contrast across the upper boundary layer is expected to be 104, 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 in the upper mantle (except in slabs) is around 10−1 as seen from Fig. 2. Figure 2. View largeDownload slide Viscosity (left) and temperature (right) profiles within the final snapshot of the convection model, in non-dimensional units. The viscosity profiles corresponds 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). Figure 2. View largeDownload slide Viscosity (left) and temperature (right) profiles within the final snapshot of the convection model, in non-dimensional units. The viscosity profiles corresponds 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). 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 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 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 during the collision (Zhang et al. 2004). The physical parameters of the model are listed in Table 1. Figure 3. View largeDownload slide 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. Figure 3. View largeDownload slide 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. Table 1. Non-dimensional and dimensional parameters of the reference convection model, also used to generate the Rayleigh number. Parameter  Non-dimensional value  Dimensional value  Rayleigh number  106    Heat production rate  20  4.6 × 10−12 W kg−1  Surface temperature  0.12  255 K  Basal temperature  1.12  2390 K  Reference density  1  4400 kg m−3  Thermal expansivity  1  4.5 × 10−5 K−1  Thermal diffusivity  1  10−6 m2 s−1  Thermal conductivity  1  4 W m−1 K−1  Reference viscosity  1  1023 Pa s  Viscosity jump factor at 660 km  30    Activation energy  8  142 kJ mol−1  Yield stress at the surface  2 × 104  230 MPa  Yield stress depth derivative  2.5×105  1030 Pa m−1  Continent nuclei viscosity factor  100    Continent nuclei yield stress  2 × 105  2300 MPa  Buoyancy number for continent nuclei  −0.4    Continent belts viscosity factor  50    Continent belts yield stress  2 × 105  2300 MPa  Buoyancy for continent belt  −0.3    Tibet viscosity factor  50    Tibet yield stress  5 × 104  590 MPa  Buoyancy number for Tibet  −0.3    Weak crust viscosity factor  0.1    Weak crust yield stress  2 × 103  23 MPa  Buoyancy number for weak crust  0.    Maximum viscosity cut-off  104  1027 Pa s  Parameter  Non-dimensional value  Dimensional value  Rayleigh number  106    Heat production rate  20  4.6 × 10−12 W kg−1  Surface temperature  0.12  255 K  Basal temperature  1.12  2390 K  Reference density  1  4400 kg m−3  Thermal expansivity  1  4.5 × 10−5 K−1  Thermal diffusivity  1  10−6 m2 s−1  Thermal conductivity  1  4 W m−1 K−1  Reference viscosity  1  1023 Pa s  Viscosity jump factor at 660 km  30    Activation energy  8  142 kJ mol−1  Yield stress at the surface  2 × 104  230 MPa  Yield stress depth derivative  2.5×105  1030 Pa m−1  Continent nuclei viscosity factor  100    Continent nuclei yield stress  2 × 105  2300 MPa  Buoyancy number for continent nuclei  −0.4    Continent belts viscosity factor  50    Continent belts yield stress  2 × 105  2300 MPa  Buoyancy for continent belt  −0.3    Tibet viscosity factor  50    Tibet yield stress  5 × 104  590 MPa  Buoyancy number for Tibet  −0.3    Weak crust viscosity factor  0.1    Weak crust yield stress  2 × 103  23 MPa  Buoyancy number for weak crust  0.    Maximum viscosity cut-off  104  1027 Pa s  View Large 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 106, based on the reference viscosity defined above. The effective Rayleigh number based on averaged viscosity is here 5.9 × 106. 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 m2 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 Ra2/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 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. Figure 4. View largeDownload slide 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. Figure 4. View largeDownload slide 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. 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 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. Figure 5. View largeDownload slide 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. Figure 5. View largeDownload slide 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. 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 6. View largeDownload slide 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 Vs and seven tomographic models of Vp. Right-hand column: vote maps at the same depth for slow seismic anomalies in seven tomographic models of Vs and seven tomographic models of Vp. Figure 6. View largeDownload slide 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 Vs and seven tomographic models of Vp. Right-hand column: vote maps at the same depth for slow seismic anomalies in seven tomographic models of Vs and seven tomographic models of Vp. 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. Figure 7. View largeDownload slide 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. Figure 7. View largeDownload slide 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. 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. Figure 8. View largeDownload slide 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). Figure 8. View largeDownload slide 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). 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 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 (104, 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:   \begin{eqnarray*} {\rm MSE}&=&\frac{1}{N}\sum _{i=1}^{N}\left(\vec{V}(\mathbf {x_i},T)-\vec{V}_{{\rm plates}}(\mathbf {x_i},T)\right) \nonumber\\ &&\cdot \left(\vec{V}(\mathbf {x_i},T)-\vec{V}_{{\rm plates}}(\mathbf {x_i},T)\right), \end{eqnarray*} where N is the number of nodes (414,144), $$\vec{V}(\mathbf {x_i},T)$$ the predicted velocity vector at position xi and age T, $$\vec{V}_{{\rm plates}}(\mathbf {x_i},T)$$ the velocity vector in the plate model (Seton et al. 2012). We note MSEt 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−1) 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 introduced with the strong temperature dependence of the viscosity. The layout of small plates does not emerge here, whatever the viscosity parametrization. Figure 9. View largeDownload slide 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). Figure 9. View largeDownload slide 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). 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/MSEt 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/MSEt for instantaneous flows produced with the 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. Figure 10. View largeDownload slide 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 MSEt ( 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. Figure 10. View largeDownload slide 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 MSEt ( 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. 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 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. Figure 11. View largeDownload slide 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. Figure 11. View largeDownload slide 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. 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. 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/MSEt 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 12. View largeDownload slide 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. Figure 12. View largeDownload slide 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. We compute a longer dynamic evolution for the weak crust model, which has the lower MSE/MSEt 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. Figure 13. View largeDownload slide 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. Figure 13. View largeDownload slide 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. 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 instead of a slowly diverging one. The MSE/MSEt 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 retreat 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 conclusions 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 predictions. A clear observation is that plumes have no influence on the instantaneous kinematics and plate boundaries here. They neither produce erroneous plate boundaries nor alter surface kinematics. The viscosity 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 velocities 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 imposed 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 conditions and/or plate kinematics), substantial differences in rheological 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 models 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 parametrization of convection fails the most. Indeed, we neglect compressibility, 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 heterogeneity 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. Multi-scale 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. Google Scholar CrossRef Search ADS   Becker T.W., Boschi L., 2002. A comparison of tomographic and geodynamic mantle models, Geochem. Geophys. Geosyst. , 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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   Bercovici D., 2003. The generation of plate tectonics from mantle convection, Earth planet. Sci. Lett. , 205( 3), 107– 121. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   Bocher M., Fournier A., Coltice N., 2017. Ensemble kalman filter for the reconstruction of the earth's mantle circulation, Nonlin. Process. Geophys. , in press. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS PubMed  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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   Chopelas A., Boehler R., 1992. Thermal expansivity in the lower mantle, Geophys. Res. Lett. , 19( 19), 1983– 1986. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS PubMed  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. Google Scholar CrossRef Search ADS   Coltice N., Larrouturou G., Debayle E., Garnero E.J., 2017. Interactions of scales of convection in the earth's mantle, Tectonophysics , in press. Google Scholar CrossRef Search ADS   Conrad C.P., Lithgow-Bertelloni C., 2002. How mantle slabs drive plate tectonics, Science , 298, 207– 209. Google Scholar CrossRef Search ADS PubMed  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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   Davies D.R., Davies J.H., 2009. Thermally-driven mantle plumes reconcile multiple hot-spot observations, Earth planet. Sci. Lett. , 278( 1), 50– 54. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   French S., Romanowicz B., 2014. Whole-mantle radially anisotropic shear velocity structure from spectral-element waveform tomography, Geophys. J. Int. , 199( 3), 1303– 1327. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   Gerya T.V., Connolly J.A., Yuen D.A., 2008. Why is terrestrial subduction one-sided?, Geology , 36( 1), 43– 46. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   Gurnis M.et al.  , 2012. Plate tectonic reconstructions with continuously closing plates, Comput. Geosci. , 38( 1), 35– 42. Google Scholar CrossRef Search ADS   Houser C., Masters G., Shearer P., Laske G., 2008. Shear and compressional velocity models of the mantle from cluster analysis of long-period waveforms, Geophys. J. Int. , 174( 1), 195– 212. Google Scholar CrossRef Search ADS   Kageyama A., Sato T., 2004. ‘yin-yang grid’: an overset grid in spherical geometry, Geochem. Geophys. Geosyst. , 5, Q09005. Google Scholar CrossRef Search ADS   King S.D., 2016. Reconciling laboratory and observational models of mantle rheology in geodynamic modelling, J. Geodyn. , 100, 33– 50. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   Lowman J.P., 2011. Mantle convection models featuring plate tectonic behavior: an overview of methods and progress, Tectonophysics , 510( 1), 1– 16. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS PubMed  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. Google Scholar CrossRef Search ADS   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 lithosphere: thoughts on the global tectonic styles of the earth and venus, Geophys. J. Int. , 133, 669– 682. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   Nakagawa T., Tackley P.J., Deschamps F., Connolly J.A., 2012. Radial 1-d seismic structures in the deep mantle in mantle convection simulations 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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   Rowley D.B., 2008. Extrapolating oceanic age distributions: lessons from the pacific region, J. Geol. , 116, 587– 598. Google Scholar CrossRef Search ADS   Rudolph M.L., Lekić V., Lithgow-Bertelloni C., 2015. Viscosity jump in earth’s mid-mantle, Science , 350( 6266), 1349– 1352. Google Scholar CrossRef Search ADS PubMed  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. Google Scholar CrossRef Search ADS   Shephard G., Matthews K., Hosseini K., Domeier M., 2017. On the consistency of seismically imaged lower mantle slabs, Sci. Rep. , 7, 1– 17. Google Scholar CrossRef Search ADS PubMed  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. Google Scholar CrossRef Search ADS PubMed  Tackley P.J., 1998. Self-consistent generation of tectonics plates in three-dimensional mantle convection, Earth planet. Sci. Lett. , 157( 1), 9– 22. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   Tauzin B., van Der Hilst R.D., Wittlinger G., Ricard Y., 2013. Multiple transition zone seismic discontinuities and low velocity layers below western united states, J. geophys. Res. , 118( 5), 2307– 2322. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   Tosi N.et al.  , 2015. A community benchmark for viscoplastic thermal convection in a 2-D square box, Geochem. Geophys. Geosyst. , 16, 2175– 2196. Google Scholar CrossRef Search ADS   Trompert R., Hansen U., 1998. Mantle convection simulations with rheologies that generate plate-like behaviour, Nature , 395, 686– 689. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   Zhang P.-Z.et al.  , 2004. Continuous deformation of the tibetan plateau from global positioning system data, Geology , 32, 809– 812. Google Scholar CrossRef Search ADS   © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# Tectonic predictions with mantle convection models

, Volume 213 (1) – Apr 1, 2018
14 pages

/lp/ou_press/tectonic-predictions-with-mantle-convection-models-r9bdpHYfbo
Publisher
The Royal Astronomical Society
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx531
Publisher site
See Article on Publisher Site

### Abstract

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. 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 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. Figure 1. View largeDownload slide Flow chart of the methodology used to generate the fully dynamic convection flows leading to the tectonic predictions. Figure 1. View largeDownload slide Flow chart of the methodology used to generate the fully dynamic convection flows leading to the tectonic predictions. 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:   \begin{equation*} \eta (z,T)=\eta _0(z)\, \exp {\left(\frac{E_a}{RT}\right)} , \end{equation*} with an activation energy Ea 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 104 to limit viscosity variations. As a consequence, the viscosity contrast across the upper boundary layer is expected to be 104, 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 in the upper mantle (except in slabs) is around 10−1 as seen from Fig. 2. Figure 2. View largeDownload slide Viscosity (left) and temperature (right) profiles within the final snapshot of the convection model, in non-dimensional units. The viscosity profiles corresponds 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). Figure 2. View largeDownload slide Viscosity (left) and temperature (right) profiles within the final snapshot of the convection model, in non-dimensional units. The viscosity profiles corresponds 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). 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 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 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 during the collision (Zhang et al. 2004). The physical parameters of the model are listed in Table 1. Figure 3. View largeDownload slide 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. Figure 3. View largeDownload slide 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. Table 1. Non-dimensional and dimensional parameters of the reference convection model, also used to generate the Rayleigh number. Parameter  Non-dimensional value  Dimensional value  Rayleigh number  106    Heat production rate  20  4.6 × 10−12 W kg−1  Surface temperature  0.12  255 K  Basal temperature  1.12  2390 K  Reference density  1  4400 kg m−3  Thermal expansivity  1  4.5 × 10−5 K−1  Thermal diffusivity  1  10−6 m2 s−1  Thermal conductivity  1  4 W m−1 K−1  Reference viscosity  1  1023 Pa s  Viscosity jump factor at 660 km  30    Activation energy  8  142 kJ mol−1  Yield stress at the surface  2 × 104  230 MPa  Yield stress depth derivative  2.5×105  1030 Pa m−1  Continent nuclei viscosity factor  100    Continent nuclei yield stress  2 × 105  2300 MPa  Buoyancy number for continent nuclei  −0.4    Continent belts viscosity factor  50    Continent belts yield stress  2 × 105  2300 MPa  Buoyancy for continent belt  −0.3    Tibet viscosity factor  50    Tibet yield stress  5 × 104  590 MPa  Buoyancy number for Tibet  −0.3    Weak crust viscosity factor  0.1    Weak crust yield stress  2 × 103  23 MPa  Buoyancy number for weak crust  0.    Maximum viscosity cut-off  104  1027 Pa s  Parameter  Non-dimensional value  Dimensional value  Rayleigh number  106    Heat production rate  20  4.6 × 10−12 W kg−1  Surface temperature  0.12  255 K  Basal temperature  1.12  2390 K  Reference density  1  4400 kg m−3  Thermal expansivity  1  4.5 × 10−5 K−1  Thermal diffusivity  1  10−6 m2 s−1  Thermal conductivity  1  4 W m−1 K−1  Reference viscosity  1  1023 Pa s  Viscosity jump factor at 660 km  30    Activation energy  8  142 kJ mol−1  Yield stress at the surface  2 × 104  230 MPa  Yield stress depth derivative  2.5×105  1030 Pa m−1  Continent nuclei viscosity factor  100    Continent nuclei yield stress  2 × 105  2300 MPa  Buoyancy number for continent nuclei  −0.4    Continent belts viscosity factor  50    Continent belts yield stress  2 × 105  2300 MPa  Buoyancy for continent belt  −0.3    Tibet viscosity factor  50    Tibet yield stress  5 × 104  590 MPa  Buoyancy number for Tibet  −0.3    Weak crust viscosity factor  0.1    Weak crust yield stress  2 × 103  23 MPa  Buoyancy number for weak crust  0.    Maximum viscosity cut-off  104  1027 Pa s  View Large 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 106, based on the reference viscosity defined above. The effective Rayleigh number based on averaged viscosity is here 5.9 × 106. 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 m2 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 Ra2/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 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. Figure 4. View largeDownload slide 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. Figure 4. View largeDownload slide 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. 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 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. Figure 5. View largeDownload slide 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. Figure 5. View largeDownload slide 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. 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 6. View largeDownload slide 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 Vs and seven tomographic models of Vp. Right-hand column: vote maps at the same depth for slow seismic anomalies in seven tomographic models of Vs and seven tomographic models of Vp. Figure 6. View largeDownload slide 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 Vs and seven tomographic models of Vp. Right-hand column: vote maps at the same depth for slow seismic anomalies in seven tomographic models of Vs and seven tomographic models of Vp. 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. Figure 7. View largeDownload slide 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. Figure 7. View largeDownload slide 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. 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. Figure 8. View largeDownload slide 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). Figure 8. View largeDownload slide 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). 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 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 (104, 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:   \begin{eqnarray*} {\rm MSE}&=&\frac{1}{N}\sum _{i=1}^{N}\left(\vec{V}(\mathbf {x_i},T)-\vec{V}_{{\rm plates}}(\mathbf {x_i},T)\right) \nonumber\\ &&\cdot \left(\vec{V}(\mathbf {x_i},T)-\vec{V}_{{\rm plates}}(\mathbf {x_i},T)\right), \end{eqnarray*} where N is the number of nodes (414,144), $$\vec{V}(\mathbf {x_i},T)$$ the predicted velocity vector at position xi and age T, $$\vec{V}_{{\rm plates}}(\mathbf {x_i},T)$$ the velocity vector in the plate model (Seton et al. 2012). We note MSEt 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−1) 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 introduced with the strong temperature dependence of the viscosity. The layout of small plates does not emerge here, whatever the viscosity parametrization. Figure 9. View largeDownload slide 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). Figure 9. View largeDownload slide 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). 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/MSEt 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/MSEt for instantaneous flows produced with the 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. Figure 10. View largeDownload slide 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 MSEt ( 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. Figure 10. View largeDownload slide 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 MSEt ( 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. 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 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. Figure 11. View largeDownload slide 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. Figure 11. View largeDownload slide 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. 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. 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/MSEt 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 12. View largeDownload slide 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. Figure 12. View largeDownload slide 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. We compute a longer dynamic evolution for the weak crust model, which has the lower MSE/MSEt 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. Figure 13. View largeDownload slide 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. Figure 13. View largeDownload slide 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. 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 instead of a slowly diverging one. The MSE/MSEt 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 retreat 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 conclusions 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 predictions. A clear observation is that plumes have no influence on the instantaneous kinematics and plate boundaries here. They neither produce erroneous plate boundaries nor alter surface kinematics. The viscosity 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 velocities 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 imposed 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 conditions and/or plate kinematics), substantial differences in rheological 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 models 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 parametrization of convection fails the most. Indeed, we neglect compressibility, 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 heterogeneity 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. Multi-scale 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. Google Scholar CrossRef Search ADS   Becker T.W., Boschi L., 2002. A comparison of tomographic and geodynamic mantle models, Geochem. Geophys. Geosyst. , 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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   Bercovici D., 2003. The generation of plate tectonics from mantle convection, Earth planet. Sci. Lett. , 205( 3), 107– 121. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   Bocher M., Fournier A., Coltice N., 2017. Ensemble kalman filter for the reconstruction of the earth's mantle circulation, Nonlin. Process. Geophys. , in press. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS PubMed  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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   Chopelas A., Boehler R., 1992. Thermal expansivity in the lower mantle, Geophys. Res. Lett. , 19( 19), 1983– 1986. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS PubMed  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. Google Scholar CrossRef Search ADS   Coltice N., Larrouturou G., Debayle E., Garnero E.J., 2017. Interactions of scales of convection in the earth's mantle, Tectonophysics , in press. Google Scholar CrossRef Search ADS   Conrad C.P., Lithgow-Bertelloni C., 2002. How mantle slabs drive plate tectonics, Science , 298, 207– 209. Google Scholar CrossRef Search ADS PubMed  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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   Davies D.R., Davies J.H., 2009. Thermally-driven mantle plumes reconcile multiple hot-spot observations, Earth planet. Sci. Lett. , 278( 1), 50– 54. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   French S., Romanowicz B., 2014. Whole-mantle radially anisotropic shear velocity structure from spectral-element waveform tomography, Geophys. J. Int. , 199( 3), 1303– 1327. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   Gerya T.V., Connolly J.A., Yuen D.A., 2008. Why is terrestrial subduction one-sided?, Geology , 36( 1), 43– 46. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   Gurnis M.et al.  , 2012. Plate tectonic reconstructions with continuously closing plates, Comput. Geosci. , 38( 1), 35– 42. Google Scholar CrossRef Search ADS   Houser C., Masters G., Shearer P., Laske G., 2008. Shear and compressional velocity models of the mantle from cluster analysis of long-period waveforms, Geophys. J. Int. , 174( 1), 195– 212. Google Scholar CrossRef Search ADS   Kageyama A., Sato T., 2004. ‘yin-yang grid’: an overset grid in spherical geometry, Geochem. Geophys. Geosyst. , 5, Q09005. Google Scholar CrossRef Search ADS   King S.D., 2016. Reconciling laboratory and observational models of mantle rheology in geodynamic modelling, J. Geodyn. , 100, 33– 50. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   Lowman J.P., 2011. Mantle convection models featuring plate tectonic behavior: an overview of methods and progress, Tectonophysics , 510( 1), 1– 16. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS PubMed  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. Google Scholar CrossRef Search ADS   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 lithosphere: thoughts on the global tectonic styles of the earth and venus, Geophys. J. Int. , 133, 669– 682. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   Nakagawa T., Tackley P.J., Deschamps F., Connolly J.A., 2012. Radial 1-d seismic structures in the deep mantle in mantle convection simulations 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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   Rowley D.B., 2008. Extrapolating oceanic age distributions: lessons from the pacific region, J. Geol. , 116, 587– 598. Google Scholar CrossRef Search ADS   Rudolph M.L., Lekić V., Lithgow-Bertelloni C., 2015. Viscosity jump in earth’s mid-mantle, Science , 350( 6266), 1349– 1352. Google Scholar CrossRef Search ADS PubMed  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. Google Scholar CrossRef Search ADS   Shephard G., Matthews K., Hosseini K., Domeier M., 2017. On the consistency of seismically imaged lower mantle slabs, Sci. Rep. , 7, 1– 17. Google Scholar CrossRef Search ADS PubMed  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. Google Scholar CrossRef Search ADS PubMed  Tackley P.J., 1998. Self-consistent generation of tectonics plates in three-dimensional mantle convection, Earth planet. Sci. Lett. , 157( 1), 9– 22. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   Tauzin B., van Der Hilst R.D., Wittlinger G., Ricard Y., 2013. Multiple transition zone seismic discontinuities and low velocity layers below western united states, J. geophys. Res. , 118( 5), 2307– 2322. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   Tosi N.et al.  , 2015. A community benchmark for viscoplastic thermal convection in a 2-D square box, Geochem. Geophys. Geosyst. , 16, 2175– 2196. Google Scholar CrossRef Search ADS   Trompert R., Hansen U., 1998. Mantle convection simulations with rheologies that generate plate-like behaviour, Nature , 395, 686– 689. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   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. Google Scholar CrossRef Search ADS   Zhang P.-Z.et al.  , 2004. Continuous deformation of the tibetan plateau from global positioning system data, Geology , 32, 809– 812. Google Scholar CrossRef Search ADS   © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Apr 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations