Research Article

The Dynamics of Plate Tectonics and Mantle Flow: From Local to Global Scales

See allHide authors and affiliations

Science  27 Aug 2010:
Vol. 329, Issue 5995, pp. 1033-1038
DOI: 10.1126/science.1191223


Plate tectonics is regulated by driving and resisting forces concentrated at plate boundaries, but observationally constrained high-resolution models of global mantle flow remain a computational challenge. We capitalized on advances in adaptive mesh refinement algorithms on parallel computers to simulate global mantle flow by incorporating plate motions, with individual plate margins resolved down to a scale of 1 kilometer. Back-arc extension and slab rollback are emergent consequences of slab descent in the upper mantle. Cold thermal anomalies within the lower mantle couple into oceanic plates through narrow high-viscosity slabs, altering the velocity of oceanic plates. Viscous dissipation within the bending lithosphere at trenches amounts to ~5 to 20% of the total dissipation through the entire lithosphere and mantle.

Mantle convection and associated plate tectonics are principal controls on the thermal and geological evolution of the Earth. These processes are central to our understanding of the origin and evolution of tectonic deformation, the evolution of the thermal and compositional states of the mantle, and ultimately the evolution of Earth as a whole. Plate creation and motion largely govern the loss of heat from the solid earth (1), and the strength of plates may control the energy dissipation and hence heat loss over geological time (2). However, despite the central importance of plate dynamics, there are fundamental uncertainties on the forces resisting and driving plate motions.

Although there is consensus that the 1- to 10-cm-per-year motion of plates is driven largely by the thermal buoyancy within subducted slabs (3) and perturbed by upper-mantle solid-solid phase transitions (4) and cooling of oceanic lithosphere from ridge to trench, the importance of the aseismic extension of slabs within the lower mantle (5) remains unresolved. The strength of subducted slabs probably regulates the velocity of plate tectonics. The vast majority of available negative buoyancy driving plates is within the transition zone and lower mantle, and if slabs are strong, then this force can be coupled directly into the edges of oceanic plates at trenches (6). However, if the oceanic lithosphere is strong during initial subduction as it bends below the trench, then the dissipation within the narrow high-viscosity slab could limit plate velocity (7). Although the importance of plate margin and slab strength has been studied in two- and three-dimensional Cartesian models aimed at understanding the physics of subduction (4, 811) and in limited regional models that assimilate observed structure (12, 13), the incorporation of realistic rheologies into models with narrow slabs and plate boundaries has remained an elusive goal of global geodynamics. Whether slabs are weak or strong remains unresolved (4).

With the incorporation of strong slabs and realistic treatment of plate margins, the ability to observationally constrain models would increase substantially. Observations constrain the deformation of slabs and will prove useful in global models: Examples include the strain rate and state of stress within slabs from deep focus earthquakes (14) and the kinematics of slab rollback in subduction zones with present-day back-arc extension (15). Large fractions of Earth’s surface (~15% globally) do not follow a rigid plate tectonic model but undergo deformation close to trenches and farther from plate margins (16). Some oceanic plates are deforming diffusively within their interiors, especially the Indo-Australian plates (17). The rich array of geodetic, topographic, gravitational, and seismic observations from local to regional scales constrains these deformations and could validate global dynamics if we are able to capture the commensurate scales in models with realistic rheologies.

Arguably, the biggest limitation on current progress is not observational, but computational: Solution through models that incorporate realistic rheologies and local geological structure has historically been prohibitive because of limitations in numerical methods and computational resources. Taken as a whole, the current generation of models poorly exploits the observational constraints on present-day deformation. For example, models of plate motion do not use observations of deformation at plate margins and interiors because numerical simulation of global mantle convection down to the scale of faulted plate boundaries has been intractable because of the wide range of time and length scales involved. Using adaptive mesh refinement (AMR) on highly parallel computers has allowed us to incorporate realistic rheologies and fine-scale observations in global mantle convection simulations and in turn reach fundamental conclusions about the forces driving the plates and the energy dissipation throughout the solid earth by intimately linking these models to observations.

Parallel adaptive mesh refinement. Capturing the large viscosity variations occurring at plate boundaries requires a mesh with about 1-km local resolution. A globally uniform mesh with 1-km mesh size would require ~1012 mesh elements, beyond both the capacity of contemporary supercomputers and the reach of numerical solution methods. With AMR, we achieve 1-km resolution near plate boundaries while using a coarser 5-km resolution within thermal boundary layers (including the oceanic lithosphere) and 15- to 50-km resolution for the rest of the mantle, saving a factor of over 103 as compared with a uniform mesh. The resulting reduction in problem size to a few hundred million elements is critical to making the simulations tractable on petascale supercomputers.

However, scaling AMR to thousands of processors is a challenge (18). Adaptively refined meshes entail irregular and dynamically changing topological mesh relations, whereas solution on parallel computers makes it necessary to store just a small part of the mesh on each processor. These mesh partitions must be changed after each refinement and coarsening step to ensure that the computational load on individual processors is balanced. We have developed scalable algorithms for these mesh operations as well as numerical solution of the mantle flow equations in our AMR finite element framework (1921). Our parallel AMR algorithms are based on a forest of adaptive octrees, in which multiple warped cubes are joined to represent general geometries. The spherical shell used to represent the mantle is composed of 24 cubes (Fig. 1). Each cube is adaptively subdivided by using an octree data structure, which allows for fast algorithms to manage the mesh adaptivity and to construct the mesh connectivity information required in numerical simulations (22). These algorithms, which have scaled to over 200,000 processor cores, are applicable to a broad spectrum of multiscale scientific and engineering problems that require high resolution in localized (possibly dynamically evolving) regions, such as near fronts, discontinuities, material interfaces, reentrant corners, and boundary and interior layers.

Fig. 1

(A) Splitting of Earth’s mantle into 24 warped cubes. The effective viscosity field is shown; the narrow low-viscosity zones corresponding to plate boundaries are seen as red lines on Earth’s surface. (B) Zoom into the hinge zone of the Australian plate [as indicated by the box in (C)] showing the adaptively refined mesh with a finest resolution of about 1 km. (C) Cross section showing the refinement that occurs both around plate boundaries and dynamically in response to the nonlinear viscosity, with plastic failure in the region from the New Hebrides to Tonga in the SW Pacific. Plates are labeled Australian, AUS; New Hebrides, NH; Tonga, TO; and Pacific, PAC.

Rheology and assimilated constraints. We solved for present-day, instantaneous mantle flow as a balance of mass and momentum (Stokes flow), in which flow velocities, strain rate, and state of stress are determined while a distribution of temperature (density) throughout the solid earth is assumed. Diffusion- and dislocation-creep and plastic-yielding govern lithosphere and mantle viscosity. Experimental constraints provide us with the form of the constitutive relation (22) and the ranges for uncertain parameters, over which we then vary. With these relations, dislocation-creep dominates in the upper mantle, whereas diffusion-creep is present throughout the whole mantle but dominates deformation in the lower mantle. Plastic-yielding is incorporated and defines a maximum strength of a material.

We incorporated plate margins as narrow zones with reduced viscosity through a spatially dependent prefactor to the constitutive relation, not unlike the approach used by others (4). Ridge and transform boundaries are treated as shallow vertical planes with locations from published compilations (23). For subduction zones, we used a global synthesis of seismic focal mechanisms for interplate thrust events, constraining a surface that continuously changes shape while following the strike of the plate boundary (22). The depth of these zones can extend below the depth of seismic coupling into the mantle wedge because water release from subducting slabs may result in a dipping low-viscosity zone above the slab within the mantle wedge (24).

The global models are driven entirely by temperature variations with three components, one based on the thermal age of the surface, a second based on seismicity-defined slabs, and a third on seismic tomography (22), which is similar to the input used in recent, lower-resolution models (25, 26). We created upper-mantle structure with seismicity and not with tomography so as to maximize the sharpness of features in the upper mantle. Globally, slabs are too poorly resolved in current tomography models because mantle wedges and slabs often blur together. Resolving the slabs as high-viscosity stress guides is critical for correctly resolving plate motions (8).

Plate motions and plateness. Even when the rheology is laboratory-based and mantle structure is seismically constrained, rigid plates are not guaranteed with nonlinear and plastic rheologies. Consequently, besides the traditional use of a single Euler pole for each plate we assessed how closely surface velocities matched that of rigid blocks—that is, the plateness (22). However, to fit plate motions we found it necessary to substantially reduce the yield stress (to ~100 MPa) from experimental values, which is not surprising given prior results (8). A nominal case with only buoyancy associated with upper-mantle slabs and plate-cooling had plate motions generally matching those in a no-net rotation (NNR) frame of reference (27). However, we found large variations in surface viscosity and strain rate on small and large scales, especially near subducting plate boundaries with strain rates varying regionally from ~10−14 to ~10−17 s−1 (Fig. 2A). Bands of high strain rate bound the subducting plate boundary adjacent to a band of low strain rate on the overriding plate. Although globally small (Fig. 2A), these variations are much longer in wavelength than the surface expression of faults between plates. In models in which these dipping zones are comparable in dimension with plate thickness, the subducting plate can flux through the weak zone, which is a common artifact of low-resolution models. Instead, the lithosphere slides by the fault so that viscosity and strain rate of the slab are determined dynamically. Consequently, the nature of coupling between slabs and plates is qualitatively different from previous models. These larger-scale variations in strain rate and viscosity are generated by bending, with the width of the deformation region being dependent on the nonlinear rheology and degree of coupling between subducting slab and deeper mantle (Figs. 3D and 4D).

Fig. 2

Strain rate, plate velocities, and plateness for three cases centered at 180°W. (A, B and D) Case 1, with only plate cooling and upper mantle slabs. (C) Case 2, identical to case 1 except for lower-mantle lateral structure. (E and F) Case 4, similar to case 2, except that n = 3.5. (A) Second invariant of strain rate. (B), (C), and (F) Plate motions in a NNR from (27) as green arrows and predicted velocities as black arrows; actual plate margins are shown as red, gray, and blue symbols. (D) and (E) Plateness for PAC shown in two ways: vector difference between computed velocity and velocity from best-fitting Euler pole, P2 (22), as a raster field with color palette shown to the right of (D); and individually inferred Euler poles within spherical caps (radius 20°) with magnitude of rotation (ω ) denoted with color of pole [palette shown to the right of (E)]. The Nuvel1-NNR pole position is shown as a red triangle and best-fitting pole for all computed velocities within PAC as a black square.

Fig. 3

Velocity vectors and surface strain rate for case 1 in (A) and case 6 in (B) that show the effect of increasing viscosity in the weak zone within the Columbia, Peru, and Chile trenches. White arrows are from Nuvel1-NNR (27), and black arrows are from predictions. Position of cross sections in (C) and (D) is indicated with red line. (C) Down-dip compression (blue) and down-dip tension (red) from CMT solutions (22). (D) Predicted compressional axes for case 6. Plates are labeled Nazca, NAZ, and South America, SAM.

Fig. 4

(A) Map of surface strain rate with plate velocities for plates in a NNR frame of reference for the New Hebrides-Tonga region. Predictions are in black and white arrows from Nuvel1-NNR (27), with micro-plate motions from (23). Plates are labeled Australian, AUS; New Hebrides, NH; Tonga, TO; Kermadec, KE; and Pacific, PAC. Cross-sections in (B), (C), and (D) denoted by red line in (A). (B) Predicted velocity and viscosity. (C) Down-dip compression (blue) and down-dip tension (red) from CMT solutions (22). (D) Predicted compressional axes for case 6.

Plate velocities for the nominal case (Fig. 2B) better match observed plate motions as compared with those of linear models (5, 28). The major oceanic plates move faster than the overriding plates so that convergence is asymmetric, unlike symmetric convergence in linear models (5, 28). The complex rheology of slabs allows upper-mantle slabs to pull oceanic plates toward trenches. Although the end result is similar to parametrizations used in linear models (29), the pull results from the self-consistent balance between buoyancy and the nonlinear resisting forces in the bending plate and, hence, gives greater insight into the forces driving and resisting plate motion.

When strong slabs embed into the lower mantle, plate motions can change substantially, as evident with the addition of temperature variations in the transition zone and lower mantle inferred from seismic tomography. This addition causes overall mantle flow velocities and most oceanic plate velocities to decrease (Fig. 2, B and C) (22), which is a behavior opposite to that expected from a linear system (5). In addition, the range in surface strain rate decreases while the plate margins become narrower but higher in viscosity. A combination of slabs penetrating the transition zone and temperature-dependent viscosity within the lower mantle (22) are likely causes for these results. A cross section through the Chilean subducting boundary (Fig. 3D) shows a gap in the continuity of high viscosity between slab and lower mantle. When temperature variations in the transition zone and lower mantle are present, the high-viscosity slab embeds into the high-viscosity lower mantle so that the much higher viscosities at depth resist surface tectonic motion (22). The presence of lower-mantle structure leads to a breakdown in asymmetric convergence at trenches (Fig. 2C). The inclusion of high-viscosity narrow slabs causes oceanic plates to either speed up or slow down through a much greater coupling with the high-viscosity lower mantle and occurs despite the nonlinearity of slab rheology. This suggests a close balance between plate motions and deep mantle processes.

Conducting searches over the ranges of uncertain parameters in the constitutive relation [activation energy in upper and lower mantle, grain-size, yield stress, and the nonlinear exponent n (22)], we were able to achieve rigid plate motions (Fig. 2D). Nearly perfect plate motions with little intraplate strain for the Pacific plate (PAC) can be found, with only the region south of New Zealand experiencing substantial strain. The rigidity of PAC is a sensitive function of the nonlinear rheology; if the nonlinearity is too strong, the plateness of PAC degrades to unacceptable levels. For example, if we combine n = 3.5 and larger activation volumes with lower mantle buoyancy, substantial strain occurs within PAC, resulting in deformation extending >103 km from convergent boundaries and widely separated Euler poles for the northern and southern parts of PAC (Fig. 2E). Even though values closer to n = 3.5 are more consistent with experimental work (30), we find that intraplate strain increased rapidly for n > 3 with platenesses apparently inconsistent with observations. Plateness also improves with smaller viscosities in the upper mantle, which is consistent with generic models (8). The Indo-Australian plate separates into two zones with a broad area of deformation (Fig. 2A), which is consistent with observations (17); however, the range of platenesses was less sensitive to the nonlinear viscosity than it was for PAC.

Slab deformation and interplate coupling. The nominal model substantially over-predicts the velocity of the Nazca (NAZ) plate toward the east and South America (SAM) toward the west (Fig. 2B). Linear models also predict unrealistic symmetric convergence toward the Peru-Chile trench (5, 28), which is inconsistent with plate kinematics in an NNR frame of reference that shows SAM to be approximately stationary. The addition of the lower-mantle structure as described above degrades the fit as both plates slow down, creating symmetric convergence toward the trench (Fig. 2C). Cross sections through the Nazca slab beneath SAM show that the state of stress is inconsistent with deep focus earthquakes, which indicate that the deepest events are strongly in down-dip compression, whereas events in 200 to 300 km in depth are in tension (Fig. 3C). In the lithosphere, although the subducting plate is in strong tension near the trench, SAM is only in mild compression with the axis of principal compression approximately orientated north-south. Great earthquakes along the Chile trench suggest that this is amongst the most coupled plate boundaries (31). By increasing the viscosity of the weak zone between NAZ and SAM (22), we simultaneously change the broad-scale plate kinematics and the state of stress in upper-mantle slabs and the overriding plate. Compared with the nominal case, with only plate cooling and slabs, increasing the interplate coupling slows NAZ toward the east, stops the westward motion of SAM, and reorients compression within the Andes to be strongly east-west (Fig. 3), but otherwise global plate motions remain unchanged.

Besides fitting deformation of ocean-continent boundaries, we find models consistent with deformation, rapid trench rollback, and back-arc extension that characterize some ocean-ocean subduction zones. Predicting such deformation has remained elusive for global and regional models. Although most generic models of trench rollback remove the overriding plate so that hot mantle extends to the surface enveloping the subducting slab (10, 11), our approach includes the overriding plate. We find that the New Hebrides slab with the attendant New Hebrides (NH) plate nearly always rolls back (without penetration into the transition zone, the slab descends vertically into the mantle, pulling the coupled Australian plate). The slab is always in a state of down-dip tension, as observed seismically.

The better known rapid rollback of Tonga-Kermadec with its distinctive down-dip compression within the Benioff zone emerges in some models and may be a sensitive function of the rheological structure within the overriding plate. During the analysis of early simulations, we found that two data sets used to create input conditions—namely, the position of plate boundaries (23) and plate age (32)—were inconsistent with the ridges not overlapping and that the zone of youngest lithosphere within the Lau Basin and Hauvre Trough was less extensive than thought (22). We found that the Kermadec (KE) and Tonga (TO) micro-plates became locked to the Australian (AUS), and the Tonga trench did not roll back even for different n, yield stresses, and upper- and lower-mantle viscosities. However, by aligning the ridges and decreasing the age of the overriding plate to be consistent with regional geology (22), KE and TO became decoupled from AUS, allowing the TO-KE trench to rapidly roll back, as indicated by rapid eastward motion of the two micro-plates (Fig. 4A). In detail, KE pivots to the east around an Euler pole near the North Island of New Zealand, as observed.

In addition, the state of stress for models with rapid trench rollback also fits the state of stress from seismic focal mechanisms. The modeled slabs are in down-dip compression in the transition zone, strongly in compression between the transition zone and upper mantle along the top of the slab, but in tension within the hinge zone, all of which is consistent with seismic focal mechanisms (Fig. 4, C to D). Viscosity near the hinge zone is reduced by ~10 times over ~100 km (Fig. 4D) and even more so along the axis of the trench and is consistent with the extreme local reduction in plate rigidity in the Kermadec trench and outer trench wall (33). With the low viscosities in the wedge and large pressure gradients near slabs, flow velocities within the mantle wedge are ~5 times higher than that of the slab, which is consistent with a recent regional study (13).

Energy dissipation in global models. The location of energy dissipation within the solid earth remains unresolved; some analyses place nearly all dissipation within the hinge zone (7, 34), whereas others do not (35, 36). However, previous models made a priori assumptions about plate motions and hinge zone viscosity. By using seismic focal mechanisms and trench rollback, the state of stress and strain rate can be observationally constrained independently of plate motions.

We find high dissipation (22) within slabs, the mantle wedge, and slab regions within the lower mantle (Fig. 5). Within slabs, the highest values of dissipation occur in the hinge zones, but dissipation is high through the entire slab; another local maximum occurs within the transition zone, which is consistent with a peak in energy release from seismicity (37). The low-viscosity zone between plates is also a zone of high dissipation. Surprisingly, even though the viscosity is probably very low (12), the mantle wedge is a zone of high dissipation (Fig. 5C) because strain rates are high (Fig. 4B) (13).

Fig. 5

Pointwise dissipation for models in cross section through small-circle at 17°S for case 6 (without lower-mantle structure) in (A) and through Case 8 (with lower-mantle structure) in (B). The Tonga area is denoted with the letter T, and the Northern Chile area with letters Ch. (C) Zoom-in on the Tonga area for case 6. (D) Same zoom-in for case 8.

We found a smaller proportion of the dissipation within the hinge zones as compared with the entire mantle and lithosphere. Even when assuming no buoyancy within the lower mantle, the dissipation within the hinge zone is generally less than 20% of the total viscous dissipation within the mantle (22). The lower mantle could be a source of substantial dissipation (Fig. 5) (22), reaching as high as 70% of the total, depending on the viscosity of the lower mantle and the lateral variations driven by temperature-dependent viscosity. If that is the case, then the dissipation in the hinge zone could be <5%. By considering the whole mantle and lithosphere, in which the hinge zone is one of several components in a tightly coupled system, these estimates are smaller than previous estimates (7, 34) but consistent with recent energetic arguments (38).

Supporting Online Material

Materials and Methods

SOM Text

Figs. S1 to S6

Tables S1 to S5


References and Notes

  1. Materials and methods are available as supporting material on Science Online.
  2. We thank W. Leng, T. Becker, and J.-P. Avouac for helpful comments on the manuscript. This work was partially supported by NSF’s PetaApps program (OCI-0749334 and OCI-0748898), NSF Earth Sciences (EAR-0426271 and EAR-0810303), U.S. Department of Energy Office of Science’s SciDAC program (DE-FC02-06ER25782 and DE-SC0002710), NSF grants CCF-0427985 and DMS-072474, and the Caltech Tectonics Observatory (by the Gordon and Betty Moore Foundation). Computing resources on Texas Advanced Computing Center’s Ranger and Spur systems were provided through the NSF TeraGrid under grant TG-MCA04N026.
View Abstract

Stay Connected to Science

Navigate This Article