Viscosity jump in Earth’s mid-mantle

See allHide authors and affiliations

Science  11 Dec 2015:
Vol. 350, Issue 6266, pp. 1349-1352
DOI: 10.1126/science.aad1929

A mysterious mid-mantle slowdown

The viscosity of Earth's deep interior plays a key role in mediating plate tectonics. Rudolph et al. combined several geophysical data sets to model the viscosity of the mantle. Mantle viscosity abruptly increases below 1000 km. The increase could explain the stalling of subducting slabs and the deflections of hot upwelling plumes around this depth. Although the viscosity increase explains some recent unexpected observations, the origin of the jump itself remains a mystery.

Science, this issue p. 1349


The viscosity structure of Earth’s deep mantle affects the thermal evolution of Earth, the ascent of mantle plumes, settling of subducted oceanic lithosphere, and the mixing of compositional heterogeneities in the mantle. Based on a reanalysis of the long-wavelength nonhydrostatic geoid, we infer viscous layering of the mantle using a method that allows us to avoid a priori assumptions about its variation with depth. We detect an increase in viscosity at 800- to 1200-kilometers depth, far greater than the depth of the mineral phase transformations that define the mantle transition zone. The viscosity increase is coincident in depth with regions where seismic tomography has imaged slab stagnation, plume deflection, and changes in large-scale structure and offers a simple explanation of these phenomena.

The viscosity of Earth’s mantle controls the rate and pattern of mantle convection and, through it, the dynamics of our planet’s deep interior, including degassing of and heat transport from the interior, mixing of compositional heterogeneity, plume ascent and passive upwelling, and slab descent. The long-wavelength nonhydrostatic geoid is a key geophysical constraint on Earth’s internal viscosity structure. At the largest spatial scales (spherical harmonic degrees 2 to 7), the geoid is most sensitive to density structure and viscosity contrasts in the lower mantle. At smaller scales, the geoid becomes increasingly sensitive to upper mantle structure, which is primarily associated with subducting slabs. Because lateral viscosity variations have minor effects on the geoid at large spatial scales (1, 2)—though they may become more important at shorter length scales (3)—it is possible to infer deep mantle viscous layering from geoid observations. However, most studies of Earth’s mantle viscosity structure impose layer interfaces to be coincident with seismic velocity discontinuities. Thus, these studies may not resolve viscous layering whose origin is distinct from that of pressure-induced phase changes (e.g., at 410- and 660-km depth), or may miss phase transitions not clearly associated with seismic discontinuities.

We use the long-wavelength nonhydrostatic geoid to infer the mantle radial viscosity structure in a manner distinct from that of previous attempts in three key ways. First, we employ a transdimensional, hierarchical, Bayesian inversion procedure (4) that does not specify at the outset the number or location of interfaces in our layered viscosity structure. The Bayesian approach is very attractive for this inverse problem because it yields a posterior probability distribution that can be analyzed to quantify uncertainties of and trade-offs between model parameters (e.g., layer depth and viscosity contrast). Second, we explore various choices for the conversion between seismic velocity anomalies and density anomalies, including depth-dependent conversion factors based on thermodynamic principles, calculated using HeFESTo (5). Finally, we use a recent whole-mantle tomographic model, SEMUCB-WM1 (6), developed with waveform tomography using highly accurate wave propagation computations, to infer mantle density structure and a modern geoid model based on 10 years of GRACE satellite observations, combined with revised estimates of the hydrostatic flattening of Earth (7, 8).

A posterior probability density function for the radial profile of viscosity is shown in Fig. 1, where the mean (taken in log-space) viscosity at each depth is shown as a purple curve. In this particular inversion, we find evidence for relatively uniform viscosity throughout the upper mantle and transition zone. Below the mantle transition zone, there is a region of lower viscosity and an increase in viscosity between 670- and 1000-km depth. The preferred depth of this viscosity increase can be inferred from Fig. 1B and is centered about 1000 km.

Fig. 1 Properties of ensemble solution.

Viscosity inversion using depth-dependent Rρ,S from HeFESTo, lmax = 3, and assumption of uncorrelated errors yields radial viscosity profiles with a viscosity increase at 1000-km depth and a lower-viscosity channel between 670 and 1000 km. (A) A 2D histogram showing the posterior likelihood of viscosity and depth values. Horizontal dotted lines indicate depths of 670 and 1000 km. (B) A 2D histogram showing the posterior likelihood of layer interface depth and viscosity increase (>1 means viscosity increases with increasing depth). (C) Posterior likelihood of having a layer interface at each depth. (D) Distribution of residuals of solutions in ensemble solution. (E) Distribution of number of layers in models in the ensemble solution.

We carried out multiple inversions to explore the effects of (i) our treatment of data and model uncertainty, (ii) the degree of truncation of the spherical harmonic expansion of the geoid used to constrain our models, and (iii) the density scaling Embedded Image (Fig. 1). We consider features of the viscosity profiles to be robust if they are common among the separate inversions. We find that all solutions place the depth of viscosity increase considerably below 670-km depth, most often near 1000-km depth. This result appears to be independent of assumptions made, including maximum spherical harmonic degree Embedded Image, choice of depth-dependent or constant Embedded Image, or treatment of data and model covariance (7). Other features of the solutions are sensitive to these choices and, therefore, their robustness is proportional to the likelihood of the assumptions from which they result. Inversions with Embedded Image (dashed curves in Fig. 2) generally have a more pronounced peak in viscosity in the mid-mantle, underlain by a weaker region between 1500- and 2500-km depth and an increase in viscosity in the lowermost mantle. Several solutions, using depth-dependent Embedded Image or Embedded Image, feature a lower-viscosity layer between 670- and 1000-km depth. Some solutions include a high-viscosity “hill” in the mid-mantle between 1000- and 1500-km depth, separating upper and lower mantles of lower viscosity.

Fig. 2 Results from multiple inversions.

Mean radial profiles of viscosity obtained in eight inversions varying Rρ,S,lmax, and eliminating buoyancy contributions from the lowermost 1000 km of the mantle (denoted by a superscript “a”) all exhibit an increase in viscosity between 670- and 1000-km depth. Models with lmax = 7 are characterized by low viscosity in the mid–lower mantle.

Many early studies advocated for layered mantle convection with an interface at or somewhat below 670-km depth, and in particular Wen and Anderson (9) noted that the amplitude and pattern of the long-wavelength geoid and surface topography could be well reproduced using mantle flow models with an imposed barrier to flow about 250 km deeper than the 670-km seismic discontinuity. However, tomographic images of relict Farallon and Tethys slabs in the lower mantle suggest that the concept of layered mantle convection is at best incomplete, and we emphasize that our mantle flow calculations do not impose layered convection.

Our results favor viscosity structures in which the overall increase in viscosity is a factor of 10 to 150, in agreement with previous studies. All of our results favor the location (interface depth) of this viscosity increase lying below 670-km depth, and most models place this viscosity increase deeper still, in the vicinity of 1000-km depth. This result is particularly intriguing given the observation that most actively subducting slabs stagnate below the 670-km seismic discontinuity, at depths of 1000 km (10). For instance, both the GAP-P4 model (11) and SEMUCB-WM1 reveal slabs stagnating above the 670-km discontinuity in the Northern Honshu arc, but passing through the 670-km discontinuity and stagnating above 1000-km depth along the Tonga and Kermadec arcs. In at least one region, Central America, the slab appears to enter the lower mantle without stagnation. The mechanism responsible for this slab stagnation is unclear, as there is no velocity discontinuity at this depth in one-dimensional (1D) seismic models (12), nor a known phase transition.

Two mechanisms have been recently suggested for slab stagnation in the mid-mantle. First, King et al. (13) have suggested that the pyroxene to majoritic garnet phase transition in subducted slabs is kinetically hindered, and thus older, colder, slabs are more prone to stagnation. Marquardt and Miyagi (14), based on high-pressure deformation experiments of (Mg,Fe)O, argued that viscosity in the regions surrounding settling slabs in the shallow-most 900 km of the lower mantle may be about two orders of magnitude higher than previously expected, causing slabs to spread laterally and to settle very slowly through this region. Our results indicate that there may be a viscosity increase in the mid-mantle, and many of our inversions have viscosity contrasts at depths comparable to those suggested (14). However, we note that the observation of regional differences in slab behavior, and in particular the speculation that old, cold, slabs preferentially stagnate, cannot be explained using our 1D viscosity structure or by a viscosity contrast that would occur in the mantle surrounding all slabs, irrespective of age, without invoking additional mantle dynamic processes or subduction zone histories, such as the prevalence of trench rollback.

Previous inversions for layered viscosity structure with prescribed layer interfaces depths revealed some indication of an increase in viscosity at or around 1000-km depth. In particular, King and Masters (15) inverted for layered viscosity structure constrained by the geoid using a uniform velocity to density conversion factor, with velocity anomalies inferred from S-wave tomographic models, and found evidence for a viscosity increase of ~20 at 670-km depth and a second increase of ~5 at 1022-km depth. Forte and Peltier (16) also found, using a combination of a slab density model and lower-mantle tomographic model, that the agreement between modeled and observed geoid was better for a layered viscosity structure with an interface at 1200-km depth than at 670-km depth. Kido et al. (17) performed inversions for layered mantle viscosity structure (with prescribed layer depths) using a genetic algorithm and found evidence for a decrease in viscosity at 670-km depth and subsequent increase in viscosity at 1000-km depth. Our study is different in that we do not prescribe at the outset the number or locations of layer interfaces in our layered viscosity structure and as a result, we place the largest viscosity contrast in the model somewhat deeper than previous studies.

Many studies from the 1980s and 1990s employed layered structures with layering identical to that of the tomographic models then available (~11 layers), or layered structures with layers at the major seismic discontinuities. Subsequent models have introduced additional layers [for instance, 25 in (18)]. To justify such parameterizations, either additional observational constraints, such as rates of glacial isostatic adjustment, plate motions, or patterns of seismic anisotropy, or additional assumptions about the smoothness of the mantle viscosity structure, are required. Paulson et al. (19, 20) used geoid and relative sea-level data as constraints on a Monte-Carlo inversion for mantle viscosity structure with one, two, and three layers. One of the central conclusions was that the GRACE and relative sea-level data cannot be used to uniquely constrain a layered mantle viscosity structure with more than two layers. Two markedly different two-layer models were permitted by these inversions (with prescribed interface depth at 670 km), one having an upper mantle with viscosity around 5 × 1020 Pa-s and a lower mantle ~4.33 more viscous and the other having an upper mantle viscosity about an order of magnitude smaller and a viscosity contrast of ~1500, similar to what was found by Ricard et al. (21). Our results generally support the suggestion that the geoid alone cannot uniquely constrain the viscosity of more than a handful of layers. Indeed, many individual models in the posterior population for each of our inversions do have more than five layers (e.g., Fig. 1), but owing to trade-offs, the layer properties of these more complex structures cannot be uniquely constrained. The posterior distribution of solutions inherently captures these trade-offs between model parameters, and the precise viscosity structures of these inversions are largely dependent on assumptions in the inversion (7).

A viscosity contrast at 1000-km depth has important implications for the dynamics of convection in Earth’s mantle, including its thermal and chemical evolution. As ascending plumes encounter abrupt changes in viscosity (in numerical models), they can be laterally deflected and thinned. Similarly, downwellings in numerical simulations become elongated laterally and compressed vertically as they encounter viscosity increases. Deflection of upwellings is observed in some tomographic models. For instance, recent tomographic images obtained by full waveform tomography with sophisticated forward-modeling approaches reveal apparent deflection at 1000-km depth of the seismically slow structures both regionally beneath the Iceland hotspot (22) and globally (23). Indeed, examples of apparent deflected upwellings, such as the feature beneath the Macdonald hotspot in the South Pacific (Fig. 3), are globally not uncommon (23). In both studies (22, 23), the apparent radius of plumes also decreases from the lower to the upper mantle. The decrease in radius appears to be coincident with the deflection at 1000-km depth. Upwelling structures in numerical simulations of mantle convection with an imposed increase in viscosity at 1000-km depth show similar behavior (Fig. 3).

Fig. 3 Observed and modeled upwellings.

(A) Shear velocity anomaly isocontours delineate deflected downwellings at 1000-km depth (horizontal line) near McDonald hotspot in SEMUCB-WM1. (B) Dimensionless temperature (Tʹ) anomaly isocontours (and pseudocolor) show similar deflection and thinning of upwellings in a numerical geodynamic model with a viscosity increase at 1000-km depth. Cool and warm colors trace dimensionless temperature variations in (B) and denote seismically fast or slow regions in (A).

Other studies use the mantle radial correlation function (24) to analyze tomographic models and to compare tomographic and geodynamic models (24, 25). Radial correlation functions calculated for SEMUCB-WM1, as well as for the global P-wave tomographic model GAP-P4 (10) for spherical harmonic degrees 1 to 3 (Fig. 4, A and B), show a high degree of correlation throughout the lower mantle at depths greater than 1000 km and a rapid decrease in correlation at 1000-km depth. Nearly identical behavior is also present in the average of S-wave tomographic models SMEAN (25) (fig. S10). Other tomographic models show a change in radial correlation around this depth as well as a change in velocity heterogeneity, particularly at spherical harmonic degree 4 (25), and an independent test based on voxel tomography favors a vertical coherence minimum around 800-km depth, below the base of the transition zone (26).

Fig. 4 Radial correlation functions of tomographic and geodynamic models.

(A) Radial correlation functions for spherical harmonic degrees 1 to 3 from SEMUCB-WM1 and (B) GAP-P4 show an abrupt decorrelation of structure across 1000-km depth. Very similar radial correlation functions are seen in the temperature field from numerical mantle convection simulations with imposed plate motions, including a viscosity contrast at 1000-km depth (C), but not when the viscosity contrast is smaller and shallower, at 670-km depth (D).

Changes in the radial correlation function may be related to changes in viscosity. Numerical simulations of convection in spherical shell geometry show that endothermic phase changes (24) and depth-dependent viscosity can both cause corresponding changes in the radial correlation. We find that a viscosity increase at 1000 km (Fig. 4C) yields a radial correlation structure much more similar to that found in tomographic models (Fig. 4, A and B) than does a viscosity increase at 670 km (Fig. 4D). The rapid change in radial correlation at 1000-km depth in tomographic models thus suggests a contrast in viscosity, because no change in phase is known to occur at this depth. We emphasize that these models include simplified representations of mantle viscosity structure (fig. S7) and that a more gradual increase in viscosity may also be compatible with the observations. Other, more complex viscosity structures can also alter the behavior of upwellings and downwellings and consequently change the radial correlation structure. Convection simulations run with a “second asthenosphere,” a weak zone extending from 670- to 1000-km depth as suggested in some of our inversions (Fig. 1) as well as in inversions by Kido et al. (17), show a greater tendency toward layered convection (27), which promotes decorrelation.

The viscosity contrast at a 1000-km depth provides a physical mechanism for the observation that slabs and plumes stagnate or become deflected deeper than the transition zone in the absence of a pervasive compositional barrier or another endothermic phase change. It may also reconcile observations of changes in seismic structure (28) that led to a proposed hot abyssal layer (29), though this was originally placed at greater depths. Given the present state of understanding in mineral physics, no unique mechanism can be identified for this increase in viscosity, and our observation should motivate further experimental and computational studies. First principles calculations have indicated a continuous though gentle increase in the viscosity of bridgmanite due to greater vacancy diffusion starting at around 40 GPa (~1000 km) and continuing until the postperovskite phase transition (30). The increase in the strength of ferropericlase observed by Marquardt and Miyagi (14) is the first positive experimental evidence for a possible change in rheology at these depths. Whether this effect, which is localized in high–strain-rate regions (surrounding slabs), should be expected to contribute to the viscosity inferred on the basis of the very-long-wavelength components of the geoid, remains to be determined. The spin transition in ferropericlase occurs at much greater depths, and first-principles simulations suggest that the higher-pressure phase (low spin) should have increased diffusion and lower viscosity (31), with a viscosity minimum near 1500-km depth (32).

Two possible intriguing (though speculative) solutions remain. Changes in the relative abundance of ferric versus ferrous iron due to disproportionation (33) at these depths or gradually over a depth range might change the bonding strength in bridgmanite enough to markedly strengthen it. Perhaps of greater interest and of more pervasive dynamical consequence might be the gradual drying of the bridgmanite perovskite as the solubility of water in the structure decreases with pressure (34), becoming more viscous at 1000-km depth.

Supplementary Materials

Materials and Methods

Figs. S1 to S10

Table S1

References (3556)

References and Notes

  1. Methods can be found in the supplementary materials on Science Online.
  2. Acknowledgments: We thank T. Becker and C. O’Neill for developing and releasing the source code to HC and the Computational Infrastructure for Geodynamics ( for distributing software, Y. Ricard for enlightening discussion, and three anonymous reviewers. This project was initiated at the 2014 CIDER workshop at the Kavli Institute for Theoretical Physics, University of California, Santa Barbara. This work was supported by NSF grant EAR/1135452 and Natural Environment Research Council NE/K006061/1, as well as a Packard Science and Engineering Fellowship to V.L. All data are available in the manuscript and supplementary materials. The GRACE gravity model GGM05S can be obtained at:
View Abstract

Navigate This Article