Tectonic predictions with mantle convection models

varying the rheological parameters makes a difference for plate boundary evolution. Also, identiﬁed errors in initial conditions contribute to ﬁrst-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 ﬁeld in comparison to instantaneous ﬂow calculations. Indeed, the initial conditions and the rheological parameters can be good enough for an accurate prediction of instantaneous ﬂow, 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.


I N T RO D U C T I O N
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 , 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 .
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. , 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.

M E T H O D
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 timedependent 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.

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 ). 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: with an activation energy E a of 142 kJ mol −1 . R is the gas constant and T the absolute dimensional temperature. Accounting for the full complexity of mantle rheology (King 2016) in such 3-D spherical models is a computational challenge, since extreme viscosity contrasts are difficult to resolve accurately. The non-dimensional reference viscosity of 1 corresponds to a non-dimensional temperature of 0.64 at zero pressure. This value is chosen before the calculation is realized such as to correspond to the expected temperature at the base of the upper boundary layer. We set a cut-off for the maximum value of the non-dimensional viscosity at 10 4 to limit viscosity variations. As a consequence, the viscosity contrast across the upper boundary layer is expected to be 10 4 , before the calculation is performed. After the calculation, the average value of the non-dimensional temperature at the base of the upper boundary layer is 0.75, that is, hotter than expected a priori. However, it is stable in the initial stage without imposed plate motions and in the stage with imposed plate motions (see the next subsection). Hence, the typical non-dimensional viscosity 18 N. Coltice and G.E. Shephard  in the upper mantle (except in slabs) is around 10 −1 as seen from Fig. 2. We consider a stress dependence of the viscosity through a pseudo-plastic approximation in order to produce plate boundaries surrounding strong plate interiors (see, for instance Rolf et al. 2012). This choice leads to stiff slabs and one-sided subduction with imposed plate kinematics, as described by Bello et al. (2015). Viscosity also depends on the type of material, which is tracked with markers. We use three types of materials. Ambient mantle corresponds to the largest fraction of the spherical shell. Continental nuclei are 175 km thick, approximating continental shields (Fig. 3.) They are buoyant, with their buoyancy number being −0.4 (200 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  (Zhang et al. 2004). The physical parameters of the model are listed in Table 1. The solution is computed with an energy contribution from the core of 25 per cent of the total surface heat flux, the rest being internal heating. Both the surface and the bottom are isothermal, defining the temperature drop for the Rayleigh number Ra of 10 6 , based on the reference viscosity defined above. The effective Rayleigh number based on averaged viscosity is here 5.9 × 10 6 . The average surface velocity obtained with these physical parameters at statistical steady state, without imposing surface velocities, is 1.2 cm yr −1 when scaled with a thermal diffusivity of 10 −6 m 2 s −1 . This is a factor of three lower than the Earth today. Unfortunately, computational cost limits the study to a lower Ra than that which would produce Earth-like velocities. Since convective velocities are proportional to Ra 2/3 , this factor of three suggests that increasing Ra by a factor of five would generate appropriate Earth-like velocities with our approximation and keeping our dimensional value of  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. 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.

Building guessed temperature fields with a convection reconstruction
The goal here is to build guessed temperature fields at 30, 20, 10 and 0 Ma using a numerical model of convection and plate reconstructions as information on the state of the mantle today and in the past. We use the methodology explained in more detail in Coltice et al. (2017) and illustrated in the flow chart ( Fig. 1): (Step 1) we build a temperature field for the continent configuration at 200 Ma based on free convection with imposed and fixed continent configuration, (Step 2) we impose plate velocities as boundary conditions of the numerical model between 200 and 30, 20, 10 and 0 Ma in increments of 1 My, updating the continent shapes at 80 Ma to account for the moderate changes which happened in terms of continental growth and deformation (Fig. 3). We use the plate reconstructions of Seton et al. (2012), but since we performed the computations presented here, Müller et al. (2016) have published updates and improvements. Because convection in our model is less vigourous than on Earth, the imposed velocities at present-day are scaled to be consistent with the convective vigour of our model (Bello et al. 2015): the rms value of imposed present-day velocities equals the rms surface velocity of the model without imposed kinematics. Imposing plate motion history generates artificial stresses at the surface, contrary to more realistic free-slip boundary conditions (Lowman 2011). A 3-D snapshot of the thermal state of the reconstruction at 0 Ma is depicted in Fig. 4.
In the following paragraphs, we compare the lateral temperature anomalies of the convection model at present day to seismic anomalies in tomographic models. Such a comparison is limited because seismic velocity is dependent on the local mineralogy and physical properties, rather than temperature alone, but is an appropriate first-order comparison for evaluating predicted mantle structure. Our model does not explicitly take into account for phase equilibria, melting and variable mantle chemistry. Therefore, regions where chemical anomalies in seismic velocities cannot be reproduced. Water and phases changes contribute substantially to seismic anomalies in the transition zone (Tauzin et al. 2013,for instance). Close to the core-mantle boundary regions, a combination of thermal and compositional effects results in broad regions of seismic velocity anomalies (Garnero et al. 2016). Considering these issues, we first compare the power spectrum of the tomographic model S40RTS (Ritsema et al. 2011) to that of the present-day velocity anomaly structure generated by our convection model. Because of the limitations explained above, we simply multiply the temperature Figure 5. Left: power spectrum of the seismic anomalies for the tomographic model S40RTS. Centre: power spectrum of the tomography-filtered velocity anomalies (proportional to temperature anomalies) in the convection model at present day. Right: correlation between the two fields. The amplitude of the power spectra is in logscaled. anomalies in our model by a factor of −0.4 per cent/100 K to obtain the S velocity anomaly field. We chose S40RTS because we could apply its resolution operator, hence taking into account the lower resolution of the tomographic model, the uneven distribution of earthquakes and seismic stations, and the parametrization of the tomographic inversion, as in Davies et al. (2012) and Koelemeijer et al. (2017). Since the resolution of the convection model is substantially finer than that of S40RTS (by more than a factor of 10), we refer to Coltice et al. (2017) for a discussion of structures of wavelength smaller than 1000 km, that is, harmonic degree >40. Both power spectra show strong degree two and strong degrees <10 in the upper mantle. The principal disagreements we interpret are within the deepest mantle and the transition zone, where degree 2 heterogeneities are very strong in S40RTS (see Fig. 5, left and centre). Indeed, the convection model does not involve deep chemical anomalies that are suspected to generate a strong seismic signature in the lower 1000 km of the mantle (Garnero et al. 2016,for a review). The convection model does not account for phase changes, mineralogical complexity (Nakagawa et al. 2012) and the water cycle (Richard et al. 2002), that would all otherwise produce seismic anomalies in the transition zone. Also, the stronger power in odd degrees in the deep mantle compared to S40RTS are expected from the filtering because of the tomographic models normal-mode data, giving strong power in even degrees, have a substantial weight. In the spectrum of the convection model, the temperature field displays a peak at long and intermediate wavelengths around 1500 km depth, which corresponds to the region where slabs start to fold and accumulate. This feature could change if we would take compressibility and phase transitions into account. The correlation between the filtered velocity anomalies computed from the convection model and S40RTS is high for the degree 2, and decreases with depth (see Fig. 5, right). The latter is expected, since the shallow mantle is more influenced than the deep mantle by the imposed plate motions, and because compressible effects, stronger as pressure increases, are not taken into account in our calculations.
We compare the location of slabs in the convection model to fast seismic anomalies in tomographic models. But tomographic models substantially differ: some are based on S waves, some on P waves; they use different 1-D reference model, seismic sources, seismograms and picking of phases in seismograms; some use finitefrequency 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. 2010Simmons et al. , 2012Burdick 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. Comparison between slices in the convection model and vote maps computed from 14 tomographic models. Left-hand column: maps at 500, 1500 and 2500 km depth of the non-dimensional temperature anomalies in the convection models. Central column: vote maps at the same depth for fast seismic anomalies in seven tomographic models of V s and seven tomographic models of V p . Right-hand column: vote maps at the same depth for slow seismic anomalies in seven tomographic models of V s and seven tomographic models of V p . Fig. 7 shows cross-sections for the Farallon, Tonga and Tethyan slabs. The Farallon slab is continuous in the convection model, but its dip angle seems to low compared to the vote map of fast anomalies. Therefore, the convection model predicts an erroneous cold structure below North America and East Atlantic in the lower mantle. The Tonga slab shows some similar patterns in both the convection model and vote maps of fast anomalies. However, the slab is broken into different pieces in the convection model, and sinks as isolated chunks. We attribute this artefact to the method of imposing plate motions. Imposing velocities at the surface of convection models violates the free-slip constraint, generating tangential stresses at the boundary (Nettelfield & Lowman 2007). These velocity gradients can breakup downwellings into several pieces at the trench, especially in intraoceanic domain because both sides of the subduction can yield (Bello et al. 2015). Below India, the location and geometry of the Tethyan slab in the convection model matches that expected from the vote map of fast anomalies. The slow seismic anomalies restricted to the transition zone do not correspond to hot anomalies in the convection model.
Overall, the computed temperature fields involve intrinsic errors. Convection structures are too thick (by a factor of two) because of the convective vigour being lower than that of the Earth. Also, the geometry of slabs is consistent with tomography models in the upper mantle but at the first order only, because of artificial breakoffs. The position of slabs is less accurate, relative to that of tomographic models, as the depth increases. The location of plumes in the numerical solution does not necessarily correspond to hotspots on Earth (see Fig. 8) because plumes emerge freely from the basal boundary layer without a priori constraint. Finally, the deep mantle thermal structure retains a memory of the initial temperature field chosen at 200 Ma, which is uncertain. Errors therefore come from uncertainties and approximation of (1) the physics of the model, (2) initial conditions and (3) imposed plate kinematics. Therefore, we limit the prediction time frame to 30 My.

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 22 N. Coltice and G.E. Shephard  conditions. The system is chaotic: model and initial condition errors propagate in time (Bello et al. 2014(Bello et al. , 2015. In test cases, Bocher et al. (2016) showed that the interval between corrections in a sequential data assimilation scheme (using surface velocities and seafloor age distribution as the data to match) has to be ≤15 My for accurate inversions of the convective temperature field. Therefore, we limit the prediction time frame to 30 My.
To study the role of the viscosity parameters, we compute numerical solutions for the instantaneous and dynamic models for (1) the same physical parameters as the convection reconstruction, (2) the same as the reference but with a lower yield stress (10 4 , 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 meansquared error on the velocity field: where N is the number of nodes (414,144), V (x i , T ) the predicted velocity vector at position x i and age T, V plates (x i , T ) the velocity vector in the plate model (Seton et al. 2012). We note MSE t 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).

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. 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 weak crust model modestly increases with the age of the convection reconstruction within the past 30 My. Some specific plates have systematically lower predicted velocities than expected: the Pacific, Nazca and Indian plates. The model with weak crust produces the highest velocities for these domains. The model with lower yield stress displays the stronger errors on velocity directions (15 • ) for the Pacific. However, the directions of the Nazca Plate are more accurate for this latter model than the others. Fig. 8 shows the residual temperature at 370 km depth in the model together with the location of 21 plumes emerging from the reconstructed flow. These plumes emerge at locations that are not Figure 9. Viscosity field and kinematics of instantaneous flow models versus plate boundaries and kinematics on Earth, at present day. Top row: viscosity field at 10 km depth emerging from the instantaneous flow calculation, and the expected plate layout of the Earth, as indicated by plate boundaries in black and based on the reconstruction of Seton et al. (2012). The reference model is in the middle column. The model with a factor of two lower yield stress in on the left, and the model with a weak crust is on the right. Bottom row: for the same models, black arrows represent model velocities and green arrows represent the expected velocities, as derived from the plate reconstruction of Seton et al. (2012).

24
N. Coltice and G.E. Shephard Figure 11. Viscosity field and kinematics of dynamic models started at 10 Ma vesrus plate boundaries and kinematics on Earth. Top row: viscosity field at 10 km depth emerging from the calculation after 10 My of evolution, and the expected plate layout of the Earth. The model having a factor of two lower yield stress in on the left, and the model with a weak crust is on the right. Bottom row: for the same models, black arrows represent model velocities and green arrows represent the expected velocities at present day.
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.

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 selforganized and we do not impose any plate boundaries or tectonic constraints. After 10 My of evolution, Fig. 11 shows the present-day viscosity field at the surface and the predicted kinematics for the low yield stress model and the model with weak crust. Both models show ridges at the expected locations except in backarc basins. The major discrepancy comes from the North Atlantic ridge, which is no longer a ridge after 10 My of evolution, but rather a shear band localizing incipient convergence (Fig. 11). The model with a weak crust still displays the ridges surrounding the Bauer Plate close to the East Pacific Rise, while they should stop spreading. The Chile Ridge is progressively fading out in both models. Trenches are located at, or close to the expected locations. Backarc basins develop in the western Pacific, but with differences in plate boundary locations relative to the Earth. The plate boundaries in these regions differ from one model to the other, the weak crust model displaying sharper bands of low viscosity and smaller plates.
The kinematics of both models show similar errors in terms of velocity direction and amplitude for most plates. The direction of the Pacific is off by <20 • for both models, but the model with weak crust predicts faster velocities, which are more consistent with the observations. The velocities of Africa and Antarctica are larger than expected for the Earth, especially for the weak crust model. Predicted kinematics for North America is the major issue of both models. The direction is more than 90 • off, leading to a closing of the North Atlantic ocean basin. It comes from the breakoff of the slab as seen in the cross-section (Fig. 12). It profoundly modifies the kinematics beyond the region whatever the rheological parameters. The value of MSE/MSE t at the final time is more than four times the initial value (1.2 and 1.87, respectively) for the weak crust model and the low yield stress model, respectively. Figure 13. Viscosity field and kinematics of the dynamic models started at 30 Ma versus plate boundaries and kinematics on Earth at the corresponding time steps. Left-hand column: viscosity field at 10 km depth emerging from dynamic evolution of the model with weak crust, and the expected plate layout of the Earth over the time evolution based on plate boundaries from the model of Seton et al. (2012 , black lines). Right-hand column: black arrows represent model velocities and green arrows represent the expected velocities over the time evolution.
We compute a longer dynamic evolution for the weak crust model, which has the lower MSE/MSE t for both the instantaneous and 10 My evolution tests. The numerical solution corresponds to a free evolution started from the initial condition set by the convection reconstruction at 30 Ma, as depicted in Fig. 13. Over this time, the predictions of the locations of several plate boundaries degrade quickly. Only the South Atlantic ridge and the South Indian ridges remain precise, moving in the appropriate directions. The Galapagos ridge initiates as expected but further south of the location on Earth. The India-Eurasia collision continues, thanks to the low resistance of the Tibet block, and subduction on the West Pacific operates as well as under South America. However, subduction under North America quickly stops, because of the early breakoff (between 30 and 20 Ma) of the slab as for the 10 My dynamic evolution. Again, the North Atlantic starts to be in compression after the breakoff, shutting down the ridge system. Also, the subduction system north and east of Australia quickly retreats until it reaches the ocean-continent boundary, instead of remaining at a similar position in the expected plate layout. As for the preceding calculations, backarc basins are generated with rapidly evolving ridge systems in connection with the moving trench. However, the small plate pattern does not match the expected one on Earth.
The predicted kinematics show a progressive 20 • change of direction of the Pacific Plate towards the south, while it is expected to remain constant on Earth. The direction of the Australian Plate also changes direction progressively to reach a 30 • offset towards the east, leading to the opening of a ridge system south of Southeast Asia. These changes of directions correlate with the retreat of the trench in the South-East Pacific described above, modifying the force balance on the Pacific and Australian plates that are converging. As for the 10 My evolution, the North American motion is quickly inconsistent with Earth evolution, before changing back again at the end to produce kinematics more consistent with the expectations. However, the relative motion between North America and Eurasia still corresponds to a slowly converging boundary instead of a slowly diverging one. The MSE/MSE t 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).

D I S C U S S I O N
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 Tectonic predictions with convection models 27 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(Bocher et al. , 2017 and adjointbased 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.

C O N C L U S I O N S
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(Bocher et al. , 2017 are under development.

A C K N O W L E D G E M E N T S
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 VISTAa basic research program in collaboration between The Norwegian Academy of Science and Letters, and Statoil (project 6268, DEFMOD). GES acknowledges support from the Research Coun-cil of Norway through its Centers of Excellence funding scheme, project number 223272. Calculations were performed at P2CHPD Lyon.