Inferring Earth’s discontinuous chemical layering from the 660-kilometer boundary topography

See allHide authors and affiliations

Science  15 Feb 2019:
Vol. 363, Issue 6428, pp. 736-740
DOI: 10.1126/science.aav0822

Inferring blocked mantle convection

The boundaries between rocks with different physical properties in Earth's interior come from either a change in crystal structure or a change in chemical composition. Wu et al. examined the roughness of the boundary between Earth's upper and lower mantle, thought to form from a change in mineral structure (see the Perspective by Houser). To their surprise, in some locations, the boundary has small-scale roughness that requires some chemical difference above and below the boundary. This observation provides evidence of partially blocked mantle circulation that leads to some chemical differences between the upper and lower mantle.

Science, this issue p. 736; see also p. 696


Topography, or depth variation, of certain interfaces in the solid Earth can provide important insights into the dynamics of our planet interior. Although the intermediate- and long-range topographic variation of the 660-kilometer boundary between Earth’s upper and lower mantle is well studied, small-scale measurements are far more challenging. We found a surprising amount of topography at short length scale along the 660-kilometer boundary in certain regions using scattered P'P' seismic waves. Our observations required chemical layering in regions with high short-scale roughness. By contrast, we did not see such small-scale topography along the 410-kilometer boundary in the upper mantle. Our findings support the concept of partially blocked or imperfect circulation between the upper and lower mantle.

The globally observed 660-km seismic discontinuity defines the top of the lower mantle and is commonly understood to involve the phase transition of the mineral ringwoodite to bridgmanite and ferropericlase. The detailed nature of this interface provides constraints on the chemical and dynamic properties of the whole mantle. Several lines of evidence support the boundary being due to the phase transition alone. If this is the case, the natural conclusion is that the whole mantle is convecting on geologic time scales (1, 2), despite the mineralogical differences between the upper and lower mantle. However, not all observations support this picture of the discontinuity. Other geochemical and mineralogical lines of evidence suggest a chemical interface, which requires some sort of dominantly layered or compartmentalized convection in order to maintain chemical differences between the upper and lower mantle (36). Seismic waves can be used to measure many features of the discontinuity related to the physical properties at the boundary, including sharpness, density, elasticity contrast, and topographic variations (711).

Topographic variation provides essential clues for understanding the nature of the 660-km discontinuity. Topography is the result of dynamic processes and the heterogeneous distribution of density. The free surface and the core mantle boundary of Earth feature topography from scales of a thousand kilometers to a few kilometers (1214), leading to the expectation that the 660-km discontinuity might also be rough at many scales. The scale of roughness on a boundary provides insight into the dynamic processes responsible for the topographic variations.

Our current picture of the topography of the 660-km discontinuity comes from precursors of reflected body waves (such as PP, SS, and P'P') (7, 8, 10) or converted phases such as P660s (15, 16) and S660P (17, 18). These methods reveal the large-scale (~1000 km) topography from thermal variations (Fig. 1A) of up to tens of kilometers (7, 10). Intermediate-scale (~100 km) topography (Fig. 1A) has been mapped with receiver functions (15, 16) or converted phases (17, 18) because they have smaller Fresnel zones. Most of the intermediate- and large-scale topography results are interpreted as lateral temperature variations (19, 20), but some studies revealed that the 660-km topography may be associated with more complex mechanisms than just the phase transition (2123).

Fig. 1 PSD of the 660-km interface and ray path of P'•660•P'.

(A) PSD of 2D free-surface topography (blue dashed line) and 660-km discontinuity topography (red dashed lines) as a function of wave number k (km1). The dashed blue line is P = C2Dk−3 with C2D = 0.3 m (33), and the red dashed lines represent C2D = 10, 100, and 1000 m, respectively. The large, intermediate, and small lateral-scale ranges of the 660-km discontinuity topography are labeled as “1,” “2,” and “3” respectively. Intermediate-scale topography has not been thoroughly sampled because of the limited distribution of seismic stations. (B) Ray path of P'•660•P'. We set a fictitious source (red star) and receiver (blue triangle) on the equator. In contrast to most routinely reported seismic phases, which travel in a great circle plane (the equator in this figure), asymmetrical scattering permits out-of-plane scattering waves. (C) Cartoon ray paths of P'•660•P' (black lines) and P'SurfP' (gray lines) scattered from topography at the relevant interface. The black and gray dotted lines show ray paths of P'P' and P'660P', respectively, waves undergoing symmetrical reflections.

In order to determine the small-scale (~10 km) topographic variations of the 660-km and the 410-km seismic discontinuities, we used the scattering of short period waves. This strategy has been successful for the upper crust, where P'SurfP' waves were generated by asymmetric (out of plane) back-scattering of small-scale free-surface topography and/or heterogeneities in the upper crust (24). We used the P'•d•P' phase (25), in which a rough interface is at a depth d below the surface of Earth (Fig. 1, B and C), and the double amplifications of PKP near its caustic distance can substantially enhance weak signals.

We chose seismic waveforms at small angular epicentral distances (roughly, 0° to 40°), at which asymmetrical scattering waves P'•d•P' are readily observable (2426). Because the 660-km discontinuity has a much lower impedance contrast than the free surface (27), we carefully selected seismic stations with low noise levels and strong earthquakes (table S1) to enhance the chance of observing clear P'•660•P'. We also used high-quality small-aperture seismic arrays such as IL (Eielson), YK (Yellow Knife), and AS (Alice Springs) to enhance signal-to-noise ratio (SNR) and perform beam-forming analysis. These three arrays have been used extensively to detect various types of weak signals because they have superb performance (25, 26, 28, 29).

The raw seismic data are usually noisy, so we used short period filtering (1.5 to 2.5 Hz, used throughout unless otherwise stated) to increase the SNR. After filtering, two strong signals P'dP' are evident in the data (Fig. 2A). The strongest signals, arriving at ~2280 s, are mainly P'SurfP', with some possible energy from P3KP (24). The other clear signals arrive 150 to 160 s before the theoretical arrival time of P'SurfP'. This 150- to 160-s lead time is consistent with the time interval between P'•660•P' and P'SurfP', implying that the earlier distinct signals could be P'•660•P'. If we take the smoothed envelope of this short-period seismogram (Fig. 2A, bottom trace), both the earlier signals and P'SurfP' have spindle shapes, which suggests scattering as their origin (12, 30), although the SNR of P'•660•P' is not high.

Fig. 2 P'•660•P' observations from individual stations and arrays.

(A) P'•660•P' observed at station LPAZ for the 20 June 2003 Mw 7.0 South America earthquake. The top trace is broadband velocity waveform data. The middle and bottom traces are high-frequency (1.5 to 2.5 Hz) filtered data and its smoothed envelope, respectively. (B) Smoothed envelope observed at the IL array from the 24 November 2008 Mw 7.3 Okhotsk Sea earthquake. Black lines represent the smoothed-envelope seismograms of the 19 constituent stations, and the red line is their stacked trace. (C) Stacked smoothed envelopes at AS, IL, and YK arrays from 11 events. Each trace represents stacked data of an array from one event, equivalent to the red line in (B). The three dashed lines are the predicted arrivals of P'•660•P' (red), P3KP (orange), and P'SurfP' (blue), respectively. The dotted green line shows the arrival of P'•410•P'.

In order to confirm the nature of the probable back-scattered seismic signals, we analyzed seismograms at the IL array from a different earthquake, beneath the Okhotsk Sea (Fig. 2B). Both P'SurfP' and P'•660•P' are clear in the stacked result, and P3KP arrives between them (Fig. 2B). The arrival time (~2130 s) and ensemble shape of P'•660•P' envelopes are similar to these at LPAZ (La Paz). Additionally, all the 19 smoothed envelope seismograms are very similar to each other, arguing for coherent and robust seismic phases rather than local noise. The stacked result yields a reliable observation, with a SNR greater than 2.

We collected more data of the three arrays (AS, IL, and YK) from 10 more earthquakes (table S1) and repeated the stacking process (Fig. 2C, stacked seismograms, and fig. S2, detailed data). We examined the National Earthquake Information Center (NEIC) earthquake catalog (31) and ruled out the possibility of contaminating signals from any aftershock in the time window of concern. With angular epicentral distances from podal distance (0°) to ~40°, the arrivals of P'•660•P' are almost constant ~160 s before P'SurfP' arrivals (Fig. 2C). A constant travel time, independent of angular epicentral distance, is the specific character of asymmetrical back-scattering (out of great-circle plane scattering) (24), in contrast to other commonly observed seismic phases traveling on a great-circle path (for example, Fig. 2C, P3KP). Thus, we confirmed that the signal is P'•660•P'. However, we did not see the P'•410•P' signal from the 410-km boundary.

An advantage of array data is that directional informational can be resolved for seismic signals in the form of wave slowness (28), which is valuable for identifying asymmetrical back-scattering (24, 25). We used the larger-aperture (~20 km) YK and AS arrays to investigate the slowness of P'•660•P'. Smaller-aperture arrays (such as the 10-km IL array) might have insufficient resolution. Signals arrive with azimuths of about ±120° off the great circle paths and slowness >1.9 s/degree (fig. S3). These results are consistent with the slowness of P'•660•P' at angular epicentral distance around 30° (24, 25) because of the asymmetrical back-scattering of PKP waves.

This type of asymmetric back-scattering can be generated by either volumetric heterogeneities near the 660-km boundary or the topography of the 660-km discontinuity (fig. S4). The arrival times of observed P'•660•P' restrict the heterogeneities to a thin layer near the 660-km discontinuity. Either deeper or shallower scatters would cause substantial advance or delay of its travel time, contrary to our observations (fig. S5). Previous studies (32) did not observe the strong precursors or coda waves of P'P', which would imply volumetric heterogeneities, so rough topography seems to be the more likely candidate.

The topographic variations of the 660-km discontinuity are constrained from the amplitudes of P'•660•P'. We used the power law P = C2Dk−3, where C2D indicates how rough the interface is for the two-dimensional (2D) power spectral density (PSD) of the 660-km interface and ray theory to model asymmetrical scattering synthetics (33). In order to reduce the uncertainties, we used P'SurfP' as a reference phase to investigate P'•660•P' because their ray paths are very close to each other.

As the second largest deep-focus event ever recorded (34), the 9 June 1994 moment magnitude (Mw) 8.2 Bolivia deep earthquake provides an opportunity to estimate the topography of the 660-km discontinuity by modeling the ratios of amplitudes of P'•660•P' to P'SurfP' (Fig. 3). Sixteen broadband stations in Bolivia and Peru recorded strong P'•660•P' and much stronger P'SurfP', with angular epicentral distances from 2° to 10° (Fig. 3C, black lines). To obtain the theoretical ratio of P'•660•P' to P'SurfP', the synthesized envelope functions of P'•660•P' (fig. S4A), P'•410•P', and P'SurfP' were formed by using the Preliminary Reference Earth Model (PREM) (27). The smoothed and squared envelope of the first 50 s direct P wave from the teleseismic station SJG (San Juan) was taken as an empirical source time function (fig. S6) and convolved with the squared synthetic envelope functions. The energy levels of background noise were estimated from the data by using the average of squared envelope in the time window 2000 to 2100 s (Fig. 3C, black dashed line) and then accounted for in the synthetic envelopes.

Fig. 3 Maps, synthetic and observed envelopes of P'•660•P' at near podal distances from the Mw 8.2 Bolivia earthquake.

(A) Map of the earthquake (red star) and seismic stations (green triangles) used in (C). (B) The sampled region of P'•660•P' (gray area) for this earthquake (red star). The gray area starts from the PKP caustic distance of 141.4° on the 660-km discontinuity and gradually vanishes because of the decaying amplitude of PKP with increasing distance. (C) Observed and synthetic smoothed envelopes of high frequency (1.5 to 2.5 Hz) P'•660•P', P'•410•P', and P'SurfP'. Smoothed envelopes of velocity seismograms of observations (corresponding to black lines) and synthetics (colored lines) are plotted with normalized amplitude. The black dashed lines indicate the zero baselines of the seismograms. C2D for the 660-km discontinuity topography models is 100 m for all the synthetics. P'SurfP' is generated by volumetric heterogeneities with exponential autocorrelation function (autocorrelation length L = 13 km and perturbations root mean square = 5%) in the top 10 km crust (24) and free surface topography (a power law PSD with C2D = 0.3 m). Taking P'SurfP' as reference, synthesized P'•660•P' with C2D = 100 m seem to match the observations well. P'•410•P' is invisible or very weak in observations.

P'•660•P' synthetics (Fig. 3C, red lines) with C2D = 100 m matched the data well in terms of both the onset and shapes, and we used this value as a reasonable estimation of the 660-km discontinuity topography for this sampled region. We attempted only to estimate the order of magnitude of C2D rather than a precise measurement (fig. S7) because substantial uncertainties are involved owing to poorly constrained aspects of Earth’s structure (33, 35). This order of magnitude of 100 m is much larger than the global average C2D of the free surface (C2D = 0.3 m) (Fig. 1A). We also computed P'•660•P' from alternative models of volumetric heterogeneities around the 660-km boundary, and this heterogeneous layer, if present, must be thin (<250 km) (fig. S8). The durations of P'SurfP' in the observations were longer than in the synthetics. This might be due to the presence of P3KP arrivals in observations, which we did not account for in the synthetics.

The 410-km global discontinuity is associated with the olivine-to-wadsleyite phase transition. Similar to P'•660•P', we would expect to see the asymmetrical scattering wave P'•410•P' if the 410-km is sharp and rough at small lateral scales. However, we do not discern P'•410•P' on individual seismograms. In order to enhance its observability, we stacked all the smoothed envelopes and were still not able to identify P'•410•P'. This indicates weaker 410-km discontinuity small-scale topography than that of the 660-km discontinuity. One interpretation of this observation is very weak small-scale topography on a sharp 410-km discontinuity (Fig. 3C). This interpretation implies that the 410-km discontinuity is a pure phase-transition boundary, which only has topographic variations at large and intermediate scales mostly caused by smoothly changing thermal structures. An alternative interpretation is that a broad transition width of the 410-km discontinuity (22, 36, 37) could substantially reduce the reflection coefficient of short-period seismic waves and decrease the P'•410•P' amplitude as the transition width approaches the wavelength of seismic waves. However, because other reports support a sharp rather than broad 410-km discontinuity (8, 38), the 410-km boundary would seem to be predominantly due to the phase transition.

P'SurfP' travels twice as far in the upper mantle and crust as P'•660•P', so uncertainty in the attenuation model of the upper mantle might bias our topography estimation. Furthermore, our synthesized P'SurfP' also had large uncertainty because of poorly quantified small-scale heterogeneities in the crust (24). Upper-mantle and crustal structures affect our estimation of small-scale 660-km discontinuity topography. Thus, we used another phase—PKiKP, which reflects from the inner core boundary (ICB)—as a reference phase to estimate the topography of the 660-km discontinuity. PKiKP is often observed at high frequency, and its ray path is similar to that of P'•660•P'. Although P'•660•P' travels four times through the lower mantle whereas PKiKP does so only twice, the lower mantle has much lower attenuation than the upper mantle (27, 39). By using PKiKP as a reference and PREM to describe the attenuation of both P'•660•P' and PKiKP, we estimated 660-km discontinuity topography by comparing envelope synthetics to observations. After applying a high-frequency (3 to 4 Hz) filter, PKiKP is clear, and the SNR of P'•660•P' is larger than 4 on the IL array for the Sea of Okhotsk earthquake (Fig. 4A and table S1). Unlike the Bolivia earthquake (Fig. 3C), P'SurfP' is even weaker than P'•660•P', probably because of attenuation in the upper mantle. We compared the measured and synthetic ratios of P'•660•P' to PKiKP (Fig. 4B) and estimated a C2D larger than 5000 m, which is several orders of magnitude larger than the globally averaged topography of the free surface. Factors such as uncertainty in the ICB reflection coefficient or topography may cause amplitude fluctuations of PKiKP (40), thus biasing the estimation of 660-km topographic variation. However, the P'•660•P'/PKiKP ratios measured at the IL array from other three events (Fig. 4B) show similar results of C2D > 1000 m. This C2D is much higher than the C2D = 100 m estimated from the Bolivia earthquake (Fig. 3C). This substantial difference could be due to large uncertainty in the estimations of C2D and/or geographical difference in the small-scale 660-km topography (fig. S9). Nonetheless, we observed strong topography of the 660-km discontinuity in both cases.

Fig. 4 P'•660•P', PKiKP, and their amplitude ratios.

(A) Smoothed envelopes of velocity seismograms of high-frequency (3 to 4 Hz) PKiKP (top) and P'•660•P' (bottom) recorded on the IL array from the 14 August 2012 Sea of Okhotsk earthquake. The red traces are the stacked results. (B) P'•660•P'/PKiKP amplitude ratios for four events. The envelope functions of P'•660•P' (fig. S4) and PKiKP were calculated by using ray theory. The red lines represent predicted P'•660•P'/PKiKP amplitude ratios for four topography models with different coefficients C2D (10, 100, 1000, and 5000 m). The measured P'•660•P'/PKiKP amplitude ratios (circles) are higher than the red line with C2D = 1000 m and well distributed around the line of C2D = 5000 m.

The P'•660•P' in this study samples not only currently occurring deep subduction zones (such as Fig. 3B, Japan) but also regions without any known major subduction (figs. S10 and S11). In the former setting, a slab stagnant at the 660-km boundary (6) could accumulate a large amount of chemical heterogeneities, which cause 660-km small-scale topography and volumetric heterogeneities. In the latter case, accumulated oceanic crusts from ancient slabs could be buoyant above the 660-km boundary and form a garnetite layer filled with chemical heterogeneities (41).

Small-scale topography of the 660-km boundary, or a less likely thin layer of volumetric heterogeneities, would best be explained with a chemical origin. A gravitationally stable garnetite layer above the 660-km interface, due to oceanic crust accumulation from currently subducting and/or ancient slabs, is possible. This scenario was argued by Ringwood (42) and further discussed by Irifune and Tsuchiya (41). Some regions lack small-scale topography of the 660-km interface, implying a globally discontinuous chemical layer. Our observations support simulations that describe subducting slabs as transient features of the transition zone, which eventually penetrate into the lower mantle. They also support a picture of partially blocked upper- to lower-mantle circulation (43), which effectively alternates between layered and whole-mantle convection. This type of model is more complicated than either the layered or whole-mantle end-member case, but adding a time-dependent degree of material exchange over Earth’s history may help unify geochemical, geodynamical, seismological, and petrological observations of the mantle.

Supplementary Materials

Materials and Methods

Figs. S1 to S11

Table S1

References (4464)

References and Notes

  1. The National Earthquake Information Center (NEIC) earthquake catalog is provided by the USGS;
  2. Materials and methods are available as supplementary materials.
Acknowledgments: The authors thank three anonymous reviewers for their constructive suggestions to improve the manuscript. W.W. is grateful to F. Simons for useful discussions. Funding: S.N. was supported by Chinese Academy of Sciences grant XDB18000000. W.W. and S.N. were supported by National Basic Research Program of China (973 Program) through grant 2014CB845901 and National Natural Science Foundation of China grant 41590854. W.W. and J.C.E.I. acknowledge support from NSF (EAR1644399 and 1736046) and the use of computing facilities provided by the Princeton Institute for Computational Science and Engineering. Author contributions: S.N. conceived the study. S.N. and J.C.E.I. supervised the study. All authors contributed to the analysis of data. W.W. processed the data, computed the synthetic seismograms, and produced figures. The manuscript was drafted by W.W. and edited by J.C.E.I. and S.N. Competing interests: All authors declare no conflicts of interest. Data and materials availability: Seismic data are archived at the Incorporated Research Institutions for Seismology (IRIS). All seismic data have been collected through the Data Management Center of IRIS, using BREQ FAST. The bathymetry/topography data are available from General Bathymetric Chart of the Oceans (

Correction (19 February 2019): The Funding section of the Acknowledgments has been corrected.

Stay Connected to Science

Navigate This Article