Plate Motions and Stresses from Global Dynamic Models

See allHide authors and affiliations

Science  17 Feb 2012:
Vol. 335, Issue 6070, pp. 838-843
DOI: 10.1126/science.1214209


Delineating the driving forces behind plate motions is important for understanding the processes that have shaped Earth throughout its history. However, the accurate prediction of plate motions, boundary-zone deformation, rigidity, and stresses remains a difficult frontier in numerical modeling. We present a global dynamic model that produces a good fit to such parameters by accounting for lateral viscosity variations in the top 200 kilometers of Earth, together with forces associated with topography and lithosphere structure, as well as coupling with mantle flow. The relative importance of shallow structure versus deeper mantle flow varies over Earth’s surface. Our model reveals where mantle flow contributes toward driving or resisting plate motions. Furthermore, subducted slabs need not act as strong stress guides to satisfy global observations of plate motions and stress.

Predicting plate motions correctly, along with stresses within the plates, has been a challenge for global dynamic models. Accurate predictions are vitally important for understanding the forces responsible for the movement of plates, mountain building, rifting of continents, and strain accumulation released in earthquakes. Previous studies have investigated these driving forces by either predicting stresses in the plates alone (1, 2) or plate motions alone (35). Other studies have taken the important step of predicting both plate motions and stresses in a single model (68). However, in addition to predicting plate motions, a successful global dynamic model must also explain plate rigidity and plate boundary-zone deformation, as well as intraplate stress patterns. Furthermore, the presence of lateral viscosity variations within the top 200 km of Earth influences the coupling between lithosphere and mantle convection. A systematic investigation of this influence is needed to improve our understanding of the driving mechanisms for plate tectonics.

We used global dynamic models to investigate the influence of lateral viscosity variations in the lithosphere and asthenosphere on both surface motions and stresses within the plates and plate boundary zones. Our models include incorporation of the effects of topography and lithosphere structure and a lithosphere coupled with whole-mantle convection, driven by density buoyancies within the mantle. Our modeling reveals the lateral viscosity variations that are necessary for matching observations. The results further emphasize the relative contributions of (i) topography and lithosphere structure and (ii) coupling with whole-mantle convection, both of which vary over Earth’s surface.

We solved the three-dimensional (3D) force balance equations after depth-integrating them from a surface of variable elevation to a common depth reference level (100 km below sea level) to obtain deviatoric stresses, strain rates, and horizontal velocities within the top 100 km of the planet (9). The body forces in these equations were derived from two sources: (i) topography and lithosphere density structure and (ii) density-driven convection within the mantle constrained by tomography and history of subduction. Benchmarking tests have demonstrated that despite the simplification used in this method, we are able to recover the horizontal components of stress, strain rate, and velocity in the upper 100 km of a full 3D whole-mantle convection model with better than 99% accuracy (10). We tested different radial and lateral viscosity variations in the lithosphere and asthenosphere, where the lateral variations were assigned based on positions of cratons and weak plate boundary zones (Fig. 1A). A relatively narrow range of viscosity models gave acceptable fits to the observations. Viscosity models that simultaneously gave a good fit to both plate motions and stresses required a stiff lithosphere (1023 Pa · s) with stiffer (1024 Pa · s) cratons (white regions in Fig. 1A) and weaker plate boundary zones (1020 to 1022 Pa · s) in the top 100 km and a moderately strong asthenosphere (300-km thickness, 1020 Pa · s). The successful models had keels beneath the cratons with viscosities less than 1023 Pa · s between depths of 100 to 200 km.

Fig. 1

(A) Absolute viscosity model (top 100 km) that provided a best fit to our observations. (B) Kinematic NNR model from (11) (blue arrows), along with predicted velocities from our global dynamic model (red arrows) in an NNR frame. The dynamic model includes contributions from both coupling with whole-mantle convection and lithosphere structure and topography. (C) Same as in (B), except the predicted velocities (red arrows) are from mantle tractions only. (D) Average poles of rotation of major tectonic plates (yellow stars) predicted by the dynamic model on top of individually inferred poles from relatively undeformed patches on the respective plates (blue dots). The NNR MORVEL poles from (12) are shown as black dots within their respective 95% confidence level error ellipses (in red). PAC, Pacific; NAM, North America; SAM, South America; ARB, Arabia; NUB, Nubia; NAZ, Nazca; EUR, Eurasia; and AUS, Australia.

The velocity field predicted by our best-fit dynamic model in a no-net-rotation (NNR) frame shows a remarkably good fit to the NNR plate motion model defined by Global Positioning System (GPS) (Fig. 1B) (11). The root mean square misfit of the velocity field from our complete dynamic model (mantle flow–associated tractions plus lithosphere structure and topography) compared at 63,000 spaced points (1° by 1°) with the kinematic NNR model is ~1 cm/year. The relative contribution of motions associated with coupling with whole-mantle flow versus topography and lithosphere structure can be understood by inspection of Fig. 1C, which is based on the contribution from mantle circulation tractions only. The relative driving mechanisms of topography and lithosphere structure versus coupling with mantle flow varies from plate to plate. The India and Nazca plates have a dominant influence from coupling with mantle flow, whereas other plates and regions approach parity in the relative contribution, with mantle-flow tractions dominating slightly. It is obvious, however, that the contribution from coupling with mantle circulation alone fails to predict surface motions.

We calculated the poles of rotation (table S1) for the angular velocities of the major tectonic plates that were predicted by the dynamic model and compared them with the latest NNR kinematic model, MORVEL (Fig. 1D) (12). The velocity of any given patch on the surface of our dynamic model was parameterized by an angular velocity possessing a pole position. The small scatter in the pole positions for these patches (blue dots) shows that the plates are behaving almost rigidly at the stress levels output by the dynamic model and for an effective viscosity of the plates of 1 × 1023 Pa · s. A comparison of the poles of rotation from the NNR MORVEL model (red 95% confidence ellipse) with the average pole positions (yellow stars) shows that the predictions for the North and the South American plates are nearly perfect (Fig. 1D and table S1). The predicted poles of the Pacific, Europe, Nubia, Nazca, and Arabia plates also lie fairly close to the respective poles of the kinematic MORVEL model.

Comparison of relative surface motions provided by our best-fit dynamic model in selected frames of reference with GPS observations (11) shows that the dynamic model is predicting motions in both the plates and plate boundary zones (Fig. 2, A to C). The poles of rotation for the relative angular velocities predicted by the dynamic model (red stars) are, in most cases, close to angular velocities (blue stars) from the latest kinematic plate-motion estimates (13). The general agreement of the predicted velocities from the dynamic model to the GPS vectors demonstrates that the model is predicting the correct deformation tensor field within plate boundary zones, including diffuse plate boundary zones, which is a difficult problem for global dynamic models.

Fig. 2

(A) Model velocity vectors (red arrows) from our global dynamic model plotted along with GPS vectors (11) (blue arrows) over North America in a Pacific fixed reference frame. A zoom-in view of the western U.S. region is shown in the inset map. Poles predicted by the dynamic model (red star) and the MORVEL plate-motion model (13) (blue star) are shown for PAC-NAM relative motion. (B) Same as in (A), but with Nubia fixed. Poles are for ARB-EUR relative motion. (C) Same as in (A), but with India fixed. Poles are for AUS-India relative motion.

Earth’s lithospheric stress field gives an indication of the driving forces that cause continental deformation and form mountain ranges and plateaus (1, 14, 15). We compared the orientation and style of our predicted deviatoric stresses with the World Stress Map (WSM) data (16) in the intraplate areas. Our predicted most compressive principal stresses (Fig. 3B) are in good agreement with the WSM horizontal most compressive stress (SHmax) directions and style (Fig. 3A). We show the stress results in three important continental deformation zones: western North America, the India-Asia collision zone, and the central Mediterranean. In western North America, deviatoric stresses from the dynamic model predict the opening of the Basin and Range Province, strike-slip along the San Andreas system, compression within the Juan de Fuca trench, and north-south compression over the Cascadia forearc (Fig. 4A). In the central Mediterranean and eastern Turkey regions, the modeled stresses (Fig. 4B) are compatible with findings of the known deformation field (17). The Hellenic arc displays trench-perpendicular compression, whereas clear strike-slip deformation is predicted along the North Anatolian fault. Improvement in prediction for these continental regions is possible through incorporation of the influence of smaller-scale convection (18). The predicted deviatoric stresses in Tibet show a predominantly strike-slip style of deformation (also mixed with normal fault-style deviatoric stress) and a rotation of SHmax within Tibet around the Eastern Himalayan Syntaxis region (Fig. 4C), similar to what is observed there. The contribution from topography and lithosphere structure plays an important role for these areas of continental deformation. We also compared our predicted stress orientation and style with strain-rate data from the Global Strain Rate Map (GSRM) (19) in the deforming areas by computing a correlation coefficient (20) between the predicted stress tensors and the GSRM strain-rate tensors. The comparison shows a good match in the mid-oceanic ridges, continental Africa and the Indo-Australia oceanic plate boundary zone, as well as in the Andes (fig. S1).

Fig. 3

(A) SHmax directions from the World Stress Map averaged within 1° by 1° areas. Red indicates a normal fault regime, blue indicates a thrust regime, and green denotes a strike-slip regime. Only SHmax directions where the regime was known have been used. (B) Most compressive horizontal principal deviatoric stress axes from our best-fitting dynamic model. The colors indicate the strain environment predicted by the deviatoric stresses in the dynamic model. Red, blue, green: same as in (A).

Fig. 4

(A) Deviatoric stress prediction in western North America plotted on top of topography. Red denotes tensional stresses, whereas black denotes compressive stresses. (B) Same as in (A), but in the central Mediterranean. (C) Same as in (A), but in the India-Asia collision zone. (D) Average horizontal tractions at the reference-level depth of 100 km below sea level, along with 95% confidence error ellipses, derived from nine mantle-flow models that fit the observations. The nine models possessed the narrow range of acceptable lateral viscosity variations described in the text and in (9). (E) Predicted second invariant of strain rates from our best-fitting dynamic model.

The traction field (21) that contributes to our best-fit dynamic model is long wavelength (Fig. 4D) and shows convergent flow in areas of downgoing slabs and divergent flow in places such as central Africa and the Pacific. Comparison of these tractions with surface velocities (Fig. 1B) indicates whether the tractions are driving or resisting. If mantle flow is leading plate motion, tractions are driving; if mantle flow is trailing the plate, then tractions are resistive. Tractions are driving in areas like the Nazca plate, eastern North America, the North Atlantic and Western Europe, northern and eastern Siberia, northern Africa, the Indian and Australian plates, and the Pacific plate in regions approaching the subduction zones. In these places, tractions act in similar directions to plate motions in an NNR frame, and thus mantle flow is leading the plate motion. On the other hand, in western North America, the northern part of South America, and southern Africa, tractions are resistive, as they are in a direction opposite to the NNR surface velocity. This is an important conclusion of our study that addresses the hugely controversial issue of whether mantle tractions are driving or resistive (14, 22, 23).

We also calibrated the absolute deviatoric stress magnitudes in the lithosphere verified through benchmarking using whole-mantle convection models (10). The average stress levels (second invariant of deviatoric stresses) within the lithosphere in our best-fit dynamic models are between 20 and 80 MPa (fig. S2). The higher stresses occur within the plates, whereas the plate boundary zones have lower values. At those stress levels, and given the predicted stress and strain-rate tensor fields, we obtain near-plate rigidity (Fig. 4E) and a close prediction to surface motions. The rigidity of the plates is evident from the low strain rates predicted in the intraplate areas, 1 to 4 × 10−9 per year (white to red areas in Fig. 4E). A goal of this study is to investigate the relative contribution of mantle flow versus lithosphere structure and topography. Although this ratio of the relative contribution varies depending on location, an average of 70% of the magnitude of lithosphere deviatoric stresses is associated with coupling with mantle flow, and the remaining 30% is associated with lithosphere structure and topography.

The issue of whether slabs remain strong (8) or weaken (2427) as they subduct is still unresolved (28). Our convection model is solely density-driven with Newtonian viscosity; no nonlinear rheology or stiff slabs have been considered, and yet our model predicts global plate motions, as well as motions within most of the world’s diffuse plate boundary zones. Our results cannot rule out the need for stiff slabs; we can only infer that they are not a necessary condition for predicting plate motions and plate boundary-zone deformation.

Supporting Online Material


SOM Text

Figs. S1 and S2

Table S1

References (2937)

Computer Codes

References and Notes

  1. Methods are available as supporting material on Science Online.
  2. By tractions we mean –τrθ and –τrϕ, which are the effective force per unit area imposed on the base of the lithosphere by mantle flow.
  3. Acknowledgments: This study was supported by NSF grant EAR-0911300. Maps were prepared using Generic Mapping Tools version 4.5.7. We thank D. Argus and C. Kreemer for sharing their kinematic models and GPS data for comparison with the dynamic models, L. Wen for help with using his code, and three anonymous reviewers for their comments that substantially improved the manuscript. We are grateful to UNAVCO, Incorporated Research Institutions for Seismology, NSF-EarthScope, International GNSS Service, NASA, and countless researchers for facilitation and availability of geodetic and seismic data. The codes are available as a zip file as part of the supporting online material.

Stay Connected to Science

Navigate This Article