Research ArticleGEOPHYSICS

What drives tectonic plates?

See allHide authors and affiliations

Science Advances  30 Oct 2019:
Vol. 5, no. 10, eaax4295
DOI: 10.1126/sciadv.aax4295


Does Earth’s mantle drive plates, or do plates drive mantle flow? This long-standing question may be ill posed, however, as both the lithosphere and mantle belong to a single self-organizing system. Alternatively, this question is better recast as follows: Does the dynamic balance between plates and mantle change over long-term tectonic reorganizations, and at what spatial wavelengths are those processes operating? A hurdle in answering this question is in designing dynamic models of mantle convection with realistic tectonic behavior evolving over supercontinent cycles. By devising these models, we find that slabs pull plates at rapid rates and tear continents apart, with keels of continents only slowing down their drift when they are not attached to a subducting plate. Our models show that the tectonic tessellation varies at a higher degree than mantle flow, which partly unlocks the conceptualization of plate tectonics and mantle convection as a unique, self-consistent system.


How does the surface of Earth deform? In the framework of plate tectonics, the surface is represented by rigid plates separated along discrete plate boundary faults (1, 2), placing the focus on the forces acting at the edge of plates (ridge push, slab pull, and shear along transform faults) and at their base (active or passive mantle drag) (3). In the past 40 years, scores of dynamic models have been deployed to evaluate these forces at present day, from fits between plate velocities and seafloor age (4) to global force balances solving for mantle flow (58). Results point to a prevalence of slab pull force over mantle drag at the base of plates, which suggests that tectonic plates drive mantle flow (3). Albeit extremely powerful, this could be a misleading representation, as a tectonic plate is a conceptualization itself. Mantle rocks together and oceanic and continental crust belong to a single convective system, in which the lithosphere is a thermal boundary layer. Isolating tectonic plates from the surrounding bulk Earth becomes hazardous when accounting for the complexity of Earth’s deformation and rock rheology (9, 10), as (i) diffuse deformation exists in most tectonically active regions of the world (11) and closing plate boundaries remains challenging (12); (ii) mechanical properties are continuously depth dependent and laterally variable, making the lithosphere-asthenosphere boundary difficult to define; and (iii) the present-day or recent tectonic pattern can be misleading, with the overall force balance potentially having changed over time, typically during a supercontinental cycle.

Here, we present global three-dimensional spherical mantle convection models that lift the constraint of perfectly rigid plates while capturing tectonics at a regional scale. These models, which account for plate-like behavior and involve aggregation and dispersal of continents, reveal collective interactions of dynamic structures at the global scale throughout a supercontinental cycle. On this basis, we propose a new paradigm to go beyond the recurrent opposition between mantle convection and plate tectonics of either exclusively top-down (slab-driven) or bottom-up (mantle upwellings and plumes) driving forces. We instead evaluate the mechanisms that make the interior actively drag the surface and vice versa.


How the interior drags the surface

In the framework of continuum mechanics, forces are expressed locally. The buoyancy force is radial, driving the flow upward or downward, with viscous stresses expressing resistance or drag. Dynamic pressure gradients drive the horizontal flow from pressure highs to pressure lows in two ways. First, poloidal flow arises from the downward suction of the near surface by negatively buoyant material (or the upward push from positively buoyant material), with gradients of gravitational potential energy, like in the lithosphere, also accounting for this. Second, toroidal flow is the faster horizontal flow in low-viscosity regions surrounding high-viscosity regions like cold subducting slabs (and, as a corollary, is nonexistent in the absence of lateral viscosity variations). In the convective system, plates represent a fragmented thermal boundary layer where heat diffuses. Within this layer, lateral flow from horizontal dynamic pressure gradients dominates, with trenches representing archetypical dynamic pressure lows and the thickened continental lithosphere representing pressure highs. Interior flow imposes the motion at the surface when surface velocity is lower than the velocity below the boundary layer. In practice, the sublithospheric mantle only drags the surface in areas where lateral pressure gradients are sufficiently small in the boundary layer to allow viscous stresses to dominate. Below the boundary layer, important lateral pressure gradients driving horizontal flow (see fig. S1) are essentially linked to (i) deep slabs that are disconnected from the top boundary layer and generate negative pressure regions above them, (ii) retreating and advancing slabs that are associated with pressure gradients crossing them within the asthenosphere moving toward the pressure low (13) or toroidal flow circumventing viscous slabs from pressure high to pressure low, and (iii) low-viscosity mantle plumes, which spread below the boundary layer and generate pressure highs, pushing the mantle above them during the course of their ascent. Overall, lateral pressure gradients within the asthenosphere correspond to slab suction (i and ii) (6) and plume push (iii) (14). This original descriptive mechanical framework is a convenient approach to revisit surface- or interior- driven motion while remaining immune from the restrictive view that slab pull versus mantle convection considerations invariably impose.

To evaluate the strength of internal drag, we quantify the areal fraction FD of the surface, which is dragged by the interior for a series of spherical convection models. The goal is to evaluate the impact of the following set of parameters: heating mode, rheology, and the presence of continents. Heating purely from within generates a unique boundary layer at the top of the system, whereas basal heating produces active upwellings in the form of plumes that have the potential to be an additional force at the surface. Rheological contrasts facilitate decoupling between surface and interior flows. Continental roots have been proposed to act as keels, enhancing viscous coupling between the lithosphere and the mantle (15). First, we define a drag coefficient D as the normalized difference between the surface velocity and the sublithospheric velocity (240 km depth) in the direction of surface motion (equation in the Supplementary Materials). D is not a measure of shear stress per se, but negative values strictly correspond to surface plates that are being dragged by the interior, and positive values correspond to surface-driven drag. FD is the areal fraction of negative D values. To evaluate how the heating mode, rheology, and presence of continents modulate FD, we compute 10 models (Fig. 1 and fig. S2), at a lower convective vigor than that of Earth, with a Rayleigh number of 106 and an average resolution of 45 km (see the Supplementary Materials and fig. S2 for more details on the calculations). We consider a radial viscosity jump by a factor of 30 at 660 km, which is identified as the most important structure from inversions of mantle viscosity profiles (16). Although alternative viscosity profiles could have been chosen (17), the goal here is to adopt a simple model including the presence of a viscosity jump that stalls downwellings.

Fig. 1 Time-averaged areal fraction FD of the surface being dragged by the mantle for a variety of 3D spherical convection models.

The first two models have no lateral viscosity variations, whereas the other models have plate-like behavior. See fig. S2 for details.

Radially viscous models (i.e., with no lateral viscosity variations) always yield time-averaged FD that are close to 0, whether heated purely from within or below. This is consistent with semianalytical models of convection, wherein surface velocities are larger than at depth (18). In models that disregard the temperature dependence of the viscosity, upwellings and plumes have the same viscosity as the boundary layer, imposing a pressure high within the boundary layer when they approach the surface. In that case, the surface drags the interior. When accounting for plate-like behavior, we use a modified Arrhenius law approximation for the temperature dependence of viscosityη(T)exp (A(TTs)/ΔT+1A2)(1)where A is a constant that determines the viscosity contrast, akin to the activation energy, T is the temperature, Ts is the surface temperature, and ΔT is the temperature difference between the bottom and top. This formulation produces smaller viscosity gradients than the classic Arrhenius law that we use thereafter (19), making it easier to solve numerically. A is 30 for the reference model, leading to a maximum viscosity contrast of six orders of magnitude (see table S1). The reference model without continents is similar to the reference model of Mallard et al. (20), with intermediate yield stress producing a plate size distribution similar to that of Earth. Introducing self-consistent plate-like generation allows >10% of the surface to be dragged by the interior flow (reference model). The difference of this result relative to radially viscous models arises from two causes. First, the temperature dependence of viscosity results in strong lithospheric plates, and hence prevents hot upwellings and plumes from reaching the very surface. Instead, they spread below the stiff boundary layer, often at rates greater than that of plate motion. Second, subducting slabs excite large-scale convecting cells that interact with multiple plate fragments. Plate fragments connected to slabs drag the mantle, and conversely, fragments that are not attached to any slab are dragged by the interior flow, with fragments pulled by connected slabs moving faster. Increasing basal heating from 4 to 43% leads to higher FD, as it boosts the first cause. Increasing the nondimensional yield stress by 30% generates larger fragments (20), hence larger FD as a result of the second cause (21). Yet, the impact is limited in both cases by decreasing the value of A by a factor of 2 or removing the viscosity jump at 660 km, resulting in mantle drag occurring at less than 25% of the surface area.

In this formulation, continents play the most important role. We introduce them in the reference model as buoyant viscous rafts, each composed of a 200-km-thick interior surrounded by a 125-km-thick and slightly weaker outer rim. In this model, viscosities and yield stresses are 100 and 10 times higher in continental interiors and rims than in the mantle, respectively. In the numerical solution comprising a unique continent covering 10% of the surface (comparable to the size of Eurasia), FD for both total and continental area is comparable with the reference model without continents. However, when multiple continents with present-day Earth geometries, covering ∼40% of the total area, are accounted for, FD increases to 35% on average, with even higher FD values over continental areas (Fig. 1). The major distinction between the two numerical solutions comes from the collision and rifting that only occurs when multiple continents are present. When continental collision occurs above a downwelling that breaks off, slab suction is triggered (22), which increases pressure gradients in the asthenosphere and ultimately drags continents one against the other. The presence of a warm upwelling further strengthens pressure gradients below the surface. A present-day analog of this scenario could be represented by the India-Asia collision (23). In addition, multiple continents also favor the existence of several large plates (almost equal in number to the number of continents themselves).

How mantle drag evolves over supercontinent cycles

To assess how the surface and interior reciprocally drag each other, we now compute the numerical solution for a more accurate and realistic, yet challenging, model building upon Arnould et al. (24). All parameters, resulting to a Rayleigh number of 107, are listed in table S1. Rheology is set with a classic Arrhenius viscosityη(T)exp (Ea+PVaRTEaRT0)(2)where Ea is the activation energy, Va is the activation volume, R is the gas constant, and P is the lithostatic pressure. T0 and η0 are the reference temperature and viscosity for the Rayleigh number. A factor of 30 viscosity jump imposed at the depth of 660 km, and the value of the yield stress at the surface is 61 MPa. This rheological law, combined with the systematic formation of a 14-km weak top layer in oceanic areas, generates asymmetric to one-sided subduction zones (with only very large and stable subduction systems occasionally appearing two-sided), as proposed previously (25). Continental rafts are introduced as above. Two antipodal 500-km-thick provinces, chemically dense and 10 times more viscous than the surrounding mantle, mimic the Large Low shear velocity Provinces (LLSVPs) atop the core at present day (26). The initial continental configuration is designed to resemble Pangea, located above the Atlantic LLSVP (see fig. S3) (27). The resolution is 23 km on average (about 50,000,000 volumetric cells) and is refined to a radial resolution of <10 km in boundary layers, with markers allowed to move. We obtained a 1500–million year (Ma) duration dynamic evolution of this system (after nine full months and 390,200 time steps on a parallel supercomputer). For computational reasons, other parameter combinations could not be tested (such as phase changes and alternative viscosity profiles).

Our numerical solution gives values for surface velocities (~4 cm year−1), heat flow (50 TW), and hypsometry (−20 km to +11 km), all within the expected range for Earth. The characteristics of plumes are also consistent with observations (28), ranging from 12 to 25, and remaining stable through time (with some persisting for more than 100 Ma). Plume temperature and buoyancy fluxes are 100 to 200 K in excess of the surrounding asthenospheric temperature and 102 to 103 kg s−1, respectively. In all cases but one, plumes are anchored at the edges or the top of the dense material lying at the base. The deep piles move very slowly and could be considered stable over periods of 300 Ma but substantially change over time scales of 500 Ma (see movies S1 to S3). Over the entire run, the average temperature decreases smoothly by about 100 K, and continental area is progressively reduced to half of its initial value due to the ablation of continental material at convergent boundaries (subduction and collision). This corresponds to a flux of 0.15 km2 year−1, consistent with estimates of present-day fluxes of recycling at both erosional and accretionary margins (29). As there is no continent production in our model to compensate for the loss, we stopped further calculations after 1500 Ma of integration. The computed thermal heterogeneity (fig. S4) is large scale, with peaks of power in the shallowest and lowermost mantle, as in tomographic power spectra (30). The radial correlation (fig. S4) shows a minimum value around 800 km (30). Overall, the power spectrum of the temperature field reveals not only that higher degrees are excited in the lithosphere as opposed to within the asthenosphere but also that plate tectonics behavior including at high harmonic degrees appears to act independently of the underlying planform of asthenospheric flow. This is expressed by a high degree of plate fragmentation that does not mirror the lower degrees of mantle flow, both in our models (e.g., movie S1) and on Earth (for instance, the African rift and China sea opening during closure of the Tethys).

As a result of the level of complexity of the model, we will not endeavor to describe individual regions over such a long temporal evolution. The model reproduces the most salient features of Earth tectonics, including subduction, back-arc extension (Fig. 2) continental drift through supercontinent assembly via introversion and extroversion modes and dispersal. Two types of rifting environments followed by seafloor spreading are observed. The first rift type involves the fast and abrupt breakup of a continent operating over distances of >1000 km (see initial stage of continent dispersal in movie S2). Once continents are fully separated, spreading rates increase up to ∼10 cm year−1 for 10 to 30 Ma. The second type of rift is a zipper mode, with propagation of the tip through the weaker regions of the continent (Fig. 3A and second stage of continental dispersal in movie S2). Opening and propagation of this type of rifting can be either slow (∼1 cm year−1) or fast (>3 cm year−1), with failed episodes of zipper rifting also existing.

Seafloor spreading operates by developing ridge-like structures (Figs. 4 and 5). In spectacular rifting events, these ridge-like structures propagate in the newly opened oceanic domains (seen as the most linear bathymetric highs in Fig. 3A and movie S2). Pure transform motion is rare but does exist in a transient form (Fig. 3B). However, transtensive/transpressive shear zones are numerous (see examples in Fig. 4), with a minimal component of extension/compression. Trenches lie either within oceanic domains or at continent-ocean boundaries (see bathymetric lows in Figs. 4 and 5). They initiate within a variety of settings, favoring areas containing buoyancy contrasts over short distances, such as continent-ocean boundaries, areas close to oceanic ridges (see movie S3) where rising plumes may trigger new subduction zones (31), or at ancient transform shear zones where instabilities can emerge (32). Intraoceanic trenches are geometrically arcuate and mostly retreat. Trench advance does occur but only where several trenches interact, resulting in the termination of short-lived basins. Back-arc basins expand in overriding plates, typically in intraoceanic settings, while others arise at continent-ocean boundaries (Fig. 2 and movies S1 and S2), as instabilities of elongated subducting systems. Small back-arc basins may develop within larger ones, closing domains between aggregating continents (see movie S1). Most slabs stagnate below the 660-km viscosity jump (33) before penetrating deeper (as in Fig. 2). Slab edges are often wrapped over themselves in the upper mantle domain while remaining connected to other plate boundaries at the surface across transtensive or transpressive domains (movie S3). The Sumatra-Banda subduction zone at present day corresponds to this type of setting (34). The replication of Earth-like features by our model validates its predictive ability and allows us to derive general behavior laws from its predictions.

Fig. 2 Surface and cross-section view across a collision (1) and a back-arc system (2), 350 Ma after the start of the model.

Topography is shown for continents only, while oceanic domains are depicted by the depth of the 1670 K isotherm, which represents the base of the thermal boundary layer. Label 3 points to an incipient subduction zone at a continent-ocean boundary. The subduction zone of back-arc system 2 stagnates below 660 km. There, the subducting plate moving at about 3 cm year−1 drags the interior, while, conversely, the overriding plate is dragged by underlying mantle flow that moves at about 8 cm year−1. In the area of collision 1, mantle velocities in the vicinity of the detaching slab are faster on the right-hand side than on the left-hand side, showing persisting mantle-driven convergence.

Fig. 3 Two examples of major features of plate tectonics in the convection model.

(A) Rift propagation within a supercontinent. Depth in oceanic domains refers to the 1670 K isotherm, which represents the base of the cold boundary layer. Short-wavelength features depict small-scale convection. Spreading rate is 10 cm year−1. Strain rate is displayed for continents only. (B) Regions of the divergence field with absolute values higher than 10−16 s−1 are superimposed on the vorticity field. Purple represents right-lateral strike-slip motion and clockwise spin. Green represents the opposite. The absence of divergence in areas of localized right-lateral motion are akin to transform faults (although more diffuse than actual transform boundaries on Earth). The occurrence of spinning plates shows the importance of toroidal motion at the surface. Continents are in yellow.

Fig. 4 Surface maps of model snapshots.

Left: Total topography (calculated without water loading with an average depth set to −2300 m). Red arrows represent surface velocities. The red dotted lines show the location of the transcurrent component during supercontinent formation. Right: Drag coefficient D. Yellow dots indicate plume locations. The blue lines show the contour of −500 K temperature anomalies at 240 km depth, highlighting the position of shallow slabs. For the lower snapshots, the black dashed polygons identify two plates, defined using the strain rate map. D is variable within the interior of each.

Fig. 5 Surface maps of model snapshots.

Left: Total topography (calculated without water loading with an average depth set to −2300 m). Red arrows represent surface velocities. Right: Drag coefficient D. Yellow dots indicate plume locations. The blue lines show the contour of −500 K temperature anomalies at 240 km depth. For the lower snapshots, the black dashed polygon identifies a plate, defined using the strain rate map. D is variable within its interior. The red dashed polygon shows the opening of a new ocean with very positive values of D.

Fig. 6 Evolution of FD, mean surface velocity (bulk and continents only), and number of hot plumes.

Large-scale tectonic events (collisions, rifting, and assembled supercontinent) are color-coded above. Time 0 is the start of the calculation with freely moving continents.

Overall, FD can be conveniently extracted as a single metric, which enables the comparison of FD during various large-scale tectonic events (Fig. 6). It is consistent with the results obtained from simpler models (values of 25% on average and 35% for continents only; Fig. 1). FD is minimal following rifting events and generally increases after continental collisions. Most of the signal comes from variations within continental areas, where FD reaches values of less than 10% after full breakup and as high as 70% for the majority of the time before supercontinent aggregation at 710 Ma. Rifting facilitates continental drift and plate reorganizations, with initial phases of oceanic spreading corresponding to high continental (surface) root mean square (rms) velocity and lower FD. Conversely, collisions and aggregation prevent rapid drift, hinder the underlying mantle flow, and increase FD. These alternating behaviors characterize the dynamics of plates together with mantle convection throughout the Wilson cycle. As shown in Figs. 5 and 6, D varies laterally below a given plate, shifting from interior drag to interior resistance. The plates highlighted in Figs. 4 and 5 are unambiguous examples of this. This behavior can also quickly flip in time, switching in less than 50 Ma from drag to resistance as observed between 300 and 350 Ma during collision.

Evaluating this in greater detail, the numerical solution starts with a 300-Ma highly dynamic burn-in period, owing to the release of the initial geometrical setup constraints on continents and LLSVPs. During this stage, continents rapidly disperse and aggregate in successive sequences (Figs. 4 and 6), with introversion-mode continental assembly occurring between 330 and 410 Ma (Fig. 4). In addition, FD exceeds 0.8, consistent with reorganization of the flow following the birth of new subduction zones and the death of earlier ones. At around 350 Ma, an introversion-mode collision initiates closure of an ocean with a strong transcurrent component (Fig. 4). The collision propagates along the edges of the continents, creating a protracted period of intercontinental dragging lasting up to 100 Ma. The change of drag patterns across the collision stage correlates with the transition from rapid continental motion at velocities up to 10 cm year−1 to a slowed drift velocity of around 1 cm year−1. This transition corresponds to a major shift for the continent, from being on a subducting plate shearing the mantle to being dragged by flow in the asthenosphere powered by the detaching slab. From 430 Ma, a long period of FD increase begins (Fig. 6), remaining stable thereafter for more than 100 Ma at approximately 0.4 (0.6 for continents only). This episode corresponds to lower rates of global surface fragmentation and the recurrent development of back-arc basins at the edges of the central continent from a continent-ocean subduction zone. This slow phase of aggregation is predominantly due to the extroversion-mode consumption of the internal ocean (movie S1). During this ocean closure, new basins with small-scale subduction systems are formed and evolve on time scales of 10 Ma (movie S1). Continents are located on overriding plates and are dragged against each other during trench retreat and oceanic closure. When continents sit on upper plates, like between 350 and 710 Ma, continental motion is dominantly the result of mantle drag.

A supercontinent eventually aggregates at 710 Ma (Fig. 5). Newly formed subduction zones progressively merge at the periphery of continents 220 Ma after assembly. Once merged trenches constitute elongated rings, rifting starts at 930 Ma (Fig. 5), splitting the supercontinent through repeated pulses into four pieces (movie S2) that drift toward each other until the end of the calculation. Prior to the onset of rifting, the number of plumes reaching the surface quickly increases from less than 20 to almost 30 (Fig. 6), before dropping down to a total of 16 when rifting ends. Although plumes are more frequent at 930 Ma, FD reaches minimal values. This observation might seem paradoxical but shows that instead of pushing the boundary layer, the plumes thermally weaken it, with some surging into the rift and newly formed ridge. The length of trenches and ridges increases synchronously (Fig. 6) throughout this period, with subsequent aggregation underway but not complete by 1500 Ma. Over the entire evolution, the topography of continents is, as expected, depressed in the vicinity of subduction zones (35) and is elevated during supercontinent configuration, with high elevations forming in convergent zones and persisting over time. Because we do not model erosion, the topography only disappears owing to dynamic pull or gravitational collapse.


Despite several limitations such as strictly one-sided subduction zones, compressibility, and strain-dependent rheology, our model consistently shows that mantle convection operates with pressure gradients located within the thermal boundary layers: the surface seemingly dragging the interior to first order. However, our results also point to 20 to 40% of the surface instead being dragged by the interior. In continental areas, fluctuations of FD are more extreme. No systematic behaviors emerge for areas where the mantle drags the surface, except at the very location of plumes within stable areas. They spread below the surface often in the direction of motion but at a faster rate as a result of their low viscosity. Lithospheric keels underneath continents are thought to increase the coupling with the interior (15). Here, we find that subducting lithosphere drags the mantle irrespective of the presence of continental keels and, hence, that continental plates are able to reach high velocities, as shown in Figs. 4 and 5 or as observed for Australia today. In this case, the viscous keels of continents do little to limit continental drift. The velocity is dependent on the overall buoyancy of the slab, supporting the importance of the length of trenches rather than the depth of the slab (which does not vary as extensively as trench length), as proposed by (3). When continents are located on overriding plates, they are often dragged by the mantle and drift slower than the average global velocity (36). In this scenario, continental keels play a role, suggesting that defining a speed limit for continents (37) could be valid only when they are dragged.

Over a supercontinent cycle, the surface drag of the interior peaks during young oceanic spreading. In our model, continental-scale rifting is initiated when elongated slabs (along-strike) activate large-scale mantle flow that tears the lithosphere (Fig. 5 and movie S2 and S3). Rifting can start with the propagation of a mid-ocean ridge through the continent (movie S2). Plumes weaken the boundary layer, which localizes the fragmentation path. Following rifting, continental blocks quickly disperse during several tens of million years (38). In our models, the dynamic force from the plume is negligible; the mantle often resists surface displacement (positive D) on both sides of the rift (Fig. 5). This is primarily due to the fact that the buoyancy flux of plumes rarely exceeds 104 kg s−1, while that of a slab is typically 105 kg s−1. Also, because plume viscosities are up to seven orders of magnitude lower than that of slabs, they spread below the lithosphere instead of reaching the very surface. Slabs break off in the aftermath of continental collisions, and mantle drag may dominate until new breakup commences. Our results differ from Zhang et al. (39), as our temperature dependence of viscosity is more than three orders of magnitude larger, which is closer to what is expected for Earth’s mantle. As a result of this, plumes in our models quickly spread out, curbing the associated overpressure, while in models with smaller viscosity contrasts, plumes are strong enough to push against the boundary layer. Yoshida and Santosh (40) used models similar to ours, with the primary differences from our calculation being long enough to capture one supercontinent cycle, a finer grid resolution in our models, and our inclusion of basal thermochemical piles.

Global mantle convection models with plate-like behavior allowed us to investigate the sources of lateral motions at the surface of Earth within a framework in which the mantle and lithosphere constitute a single self-organized system. Slabs connected to the surface generate the strongest pressure gradients within the lithosphere. These gradients trigger the fastest velocities regardless of the presence of continental blocks embedded within the pulled plate. The large-scale mantle flow slowly drags the surface only when the pressure gradient is sufficiently low within the boundary layer. Plumes spread below the lithosphere, playing a rheological role in weakening the asthenosphere rather than pushing plates. Following a collision, the mantle drags the surface, with the flow being excited by the sinking detached slab and supported by warm upwellings.

Our models, particularly the one reproducing salient features of Earth, show that the structure of the surface to asthenospheric drag pattern varies at a higher degree than revealed by the tessellation of plates. Continents reinforce that behavior, making the surface more sensitive to large-scale mantle flow for plates that are not connected to slabs. Keels limit drift only when the mantle is dragging them and not when they are embedded within a subducting plate. On Earth, as in our numerical predictions of an Earth-like planet, the large-scale flow and regional patterns of deformation are often unrelated; for instance, rifting or back-arc opening may occur within large continental blocks themselves undergoing collision. While such geological observations have restricted the conceptualization of plate tectonics and mantle convection as a unique self-consistent system for decades, models of mantle convection with plate-like tectonics show that this deadlock can now be resolved.


We solved the conservation equations of mass, momentum, energy, and composition for mantle convection (41), accounting for variable viscosity on a yin-yang grid (42) with the code StagYY (43). We worked with the Boussinesq approximation, neglecting the dynamic effects of compressibility. The set of coupled equations the code solves isv=0 (3)(η(v+(v)T))p=Ra(α(z)α0T+Bi Ci)er (4)tT+vT=2T+H (5)tCi+vCi=0 (6)where T is the temperature, v is the velocity, p is the dynamic pressure, η is the viscosity, α(z) is the thermal expansion at depth z and α0 the thermal expansion at the surface, Ci is the composition for the material i (being mantle, continental lithosphere, dense basal mantle, and weak surface layer), H is the internal heating rate, and er is the unit vector in the radial direction.

The convective vigor of the system is measured by the Rayleigh number RaRa=ρ0gα0ΔTh3κη0 (7)where ρ0 is the reference density, g is the gravitational acceleration. ΔT is the temperature drop across the shell of the thickness h, κ is the thermal diffusivity, and η0 is the reference viscosity.

Bi is the buoyancy ratio of material i that corresponds to the density ratio (ρi − ρ0)/(ρ0αΔT), where ρi is the density of material i given for the temperature at the surface.

In terms of boundary conditions, the surface and the base are both free slip. Temperatures at the surface and at the base are fixed over the calculation.

To generate plate-like behavior, we used a pseudoplastic approximation: When stresses reach the yield stress limit σY, viscosity drops as σY/(2ϵ˙II), where ϵ˙II is the second invariant of the strain rate tensor. A well-known drawback is that subduction is often two sided. A weak surface layer and a free surface (which we do not model here) generate asymmetric and one-sided subduction (44). The nondimensional value of σY in our reference model at low Ra is 1.5 × 104 and is comparable with the model and parameters of (20) (H being higher here), displaying Earth-like size distribution of plates. We tracked compositional fields corresponding to continental material, weak top layer, and dense basal material with markers advected in the flow (45).

We selected nondimensional parameters for the models at low convective vigor following the studies of Rolf et al. (46) and Mallard et al. (20). Within their proposed range of yield stress values, and allowing for little increase with depth, convection models produce plate-like behavior and plate size distributions with both large and small plates. For the model at high convective vigor, we built on the work of (24) in spherical annulus geometry, who determined a range of Ra and rheological parameters that generate plate-like behavior, and both heat flow and topography that are similar to Earth. In our high Ra model, thermal expansion decreases with depth by a factor of 3, as suggested by mineral physics (47).

The initial conditions of the calculations correspond to statistical steady states, meaning that we first computed a solution running a model for more than 5000 Ma with fixed markers (continents, weak layer, and dense basal material), progressively increasing resolution to reach a state where temperature, vrms, and heat flow fluctuate around their average values, respectively. These solutions are our initial conditions. 3D snapshots of the 11 models are shown in fig. S1. Profiles of temperature and viscosity for the high Ra model are shown in fig. S5.

The drag coefficient D is computed asD=vsurfvdvsurfvsurfvrms (8)where vsurf is the local surface velocity vector and vd is the horizontal velocity vector at depth d. vrms is the typical velocity scale of tectonic velocity, being the time-averaged surface rms velocity (4 cm year−1). As a consequence, D does not measure a stress per se but a velocity gradient factor. Very negative values of D mean that the interior moves substantially faster than the surface.

We detected plume conduits using a threshold of 160 km year−1 on the radial heat advection field (the product of temperature and vertical velocity) at a depth of 350 km and using a spatial connectivity scheme as described by (48).


Supplementary material for this article is available at

Fig. S1. Close-up of relative dynamic pressure field from the cross-section shown in Fig. 2.

Fig. S2. 3D snapshots of the 10 low Ra models and the model at Earth-like convective vigor (high Rayleigh number).

Fig. S3. Initial distribution of continental material (yellow), the stiffer roots being thicker, and dense basal piles (red) in the model with high convective vigor.

Fig. S4. Properties of the mantle convection model for the snapshot of Fig. 2 and fig. S1.

Fig. S5. Temperature and viscosity profiles in the high Ra calculation.

Movie S1. Movie showing the closing of an ocean and supercontinent assembly.

Movie S2. Movie showing breakup sequences of a supercontinent.

Movie S3. Movie showing subduction initiations.

Table S1. Physical parameters of the reference model (low Rayleigh number) and the high Rayleigh number model.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: Calculations were performed on the AUGURY supercomputer at P2CHPD Lyon. We thank F. Dubuffet for help on the supercomputer. We thank M. Tetley and A. Holt who helped us improve the English writing. Funding: The support for this research was provided by the European Union’s Horizon 2020 research and innovation program under ERC grant agreement no. 617588. Author contributions: All authors contributed to the development and analysis of the models, preparation of the figures, and writing of the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: The code StagYY is the property of P. J. Tackley and ETH-Zürich and is available for collaborative studies from P. J. Tackley ( N.C. provides the minor modifications of the code StagYY, which are specific to this study upon request. Data used to generate maps and time series are available at Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article