Research Article

Normal-Mode and Free-Air Gravity Constraints on Lateral Variations in Velocity and Density of Earth's Mantle

See allHide authors and affiliations

Science  20 Aug 1999:
Vol. 285, Issue 5431, pp. 1231-1236
DOI: 10.1126/science.285.5431.1231


With the use of a large collection of free-oscillation data and additional constraints imposed by the free-air gravity anomaly, lateral variations in shear velocity, compressional velocity, and density within the mantle; dynamic topography on the free surface; and topography on the 660-km discontinuity and the core-mantle boundary were determined. The velocity models are consistent with existing models based on travel-time and waveform inversions. In the lowermost mantle, near the core-mantle boundary, denser than average material is found beneath regions of upwellings centered on the Pacific Ocean and Africa that are characterized by slow shear velocities. These anomalies suggest the existence of compositional heterogeneity near the core-mantle boundary.

Seismologists have produced numerous shear (S) and compressional (P) velocity models of the mantle based mainly on body-wave data (1–4). Unfortunately, they have been unable to develop a whole-mantle density model because short-period body waves are not directly sensitive to density and long-period normal-mode observations have been too scarce and of insufficient quality. There are density models based on current and past plate motions in which the sources of heterogeneity are cold subducted slabs (5), and an upper mantle density model, containing strong high-degree components, has been determined with surface wave data (6).

A density model of the mantle is important to geophysics because mantle-flow calculations cannot be performed without it (7, 8). However, a whole-mantle density model is usually obtained by scaling an S velocity model, which is justifiable if the heterogeneity has a single cause that affects both, for example, variations in temperature.

In recent years, mainly as a result of the 9 June 1994 Bolivia earthquake, the quantity and quality of normal-mode, or free-oscillation, data have improved substantially (9), allowing us to constrain the long-wavelength structure of Earth. We used this data set, combined with further constraints imposed by Earth's gravity field, to determine lateral variations inS velocity, P velocity, and density within the mantle, in addition to dynamic topography on the free surface and topography on the 660-km discontinuity (660) and the core-mantle boundary (CMB).


The sensitivity of a normal mode to lateral heterogeneity is usually visualized in the form of a splitting function (10), which represents a local radial average of Earth's three-dimensional (3D) heterogeneity. Each mode has its own unique splitting function, which at high degrees corresponds to a surface-wave phase velocity map. An isolated mode of degree l “sees” even-degree structure up to degree 2l, and some clusters of modes have sensitivity to Earth's odd-degree heterogeneity (11). A normal-mode splitting function σ depends on colatitude θ and longitude φ and may be expanded in spherical harmonicsY s t of degree s and ordert as (12)Embedded Image(1)The splitting-function coefficients c s t are linearly related to 3D relative variations in S velocity (δβ/β)s t, Pvelocity (δα/α)s t, density (δρ/ρ)s t, and topography on discontinuities (δh/h)s t by (12)Embedded Image Embedded ImageThe integration over radius r is from the CMB with radius b (3480 km) to Earth's surface with radius a (6371 km), and the summation is over all discontinuities d. The sensitivity kernelsK β, K α,K ρ, and Kd are given in terms of the eigenfunctions of the modes that define the splitting function (12). The number of unknown parameters can be reduced by assuming that lateral variations in P velocity and density are proportional to lateral variations in Svelocity. In the past, this has been the common approach (12, 13), with typical scaling values of δ ln α/δ ln β = 0.55 and δ ln ρ/δ ln β = 0.2 to 0.4 (8, 14).

Lateral variations in density, combined with flow-induced topography on discontinuities, determine Earth's gravity anomaly. Undulations on boundaries are determined by the density contrast and radial stresses imposed by mantle flow, which depend on lateral variations in density and the relative viscosity of the layers above and below the discontinuity. Frequently, the objective of mantle-flow calculations is to find a viscosity profile and velocity-to-density scaling that reconcile a given velocity model with the observed gravity anomaly. Here, we determined lateral variations in density and boundary topography simultaneously; this circumvents the need for a viscosity model when fitting the gravity anomaly. The spherical harmonic coefficients of the free-air gravity anomaly,f s t (15), are linearly related to those describing lateral variations in density, (δρ/ρ)s t, and boundary topography, (δh/h)s t, by (16)Embedded Image(2)The kernels K ρ andK d describe the sensitivity of the gravity anomaly to 3D density and boundary topography, respectively.

Velocity and Density Models

The splitting functions were corrected for the effects of Earth's crust with the recent crustal model Crust5.1 (17). We assumed an isostatically compensated crust in the calculation of the gravity field. Our “dynamic topography,” therefore, is the nonisostatic topography of the free surface that is produced by density anomalies within the entire mantle, including the lithosphere.

In the inversion, we allowed for lateral variations inS velocity, P velocity, density, dynamic free-surface topography, and topography on the 660 and the CMB (18). Initially, we inverted for topography on the 410-km discontinuity (410) as well but found that undulations on this boundary are poorly constrained by the data (19). Therefore, topographic variations on 410 are those determined in a recent travel-time study (20). The starting model for the inversion included S velocity model SKS12WM13 (1) andP velocity model P16B30 (3). The density component of the starting model was obtained by scaling velocity model SKS12WM13 by a factor of 0.2. There was no starting dynamic topography on the free surface or the CMB, except for its excess ellipticity, which was determined by very long baseline interferometry (21). Starting topography on 660 was taken from Gu et al. (20). This starting model explains 74% of the variance in the mode data but has a relatively high χ2/N of 6.6, where N = 2850 is the number of normal-mode data (22). Model SPRD6, obtained by an inversion of normal-mode and gravity data, explains 92% of the variance in the mode data and 96% of the variance in the gravity data, which correspond to χ2/N = 2.0 and χ2/N f = 0.07, whereN f = 27 is the number of free-air gravity coefficients. The fit to the free-air gravity anomaly is excellent; however, it is well known that the gravity anomaly can be relatively easily fit because of its highly nonunique dependence on density and boundary topography. Therefore, we did not allow free-air gravity to change models of density and boundary topography substantially from models determined with only normal-mode data.

We compared the S model obtained from the normal-mode inversion with a recent model based on body-wave data, SKS12WM13 (1) (Fig. 1). In the shallow mantle, continents are characterized by fast velocity anomalies, and mid-ocean ridges correspond with slow velocity anomalies. There are strong heterogeneities in the lowermost mantle, with a distinct ring of high-velocity anomalies around the Pacific Ocean. The strong negative anomalies in the lowermost mantle underneath the Pacific Ocean and Africa are interpreted as mantle upwellings. The correlation between the two models is high in the upper mantle and near the CMB, but a lower correlation is found in the mid-mantle, which is relatively poorly constrained by travel-time and waveform data.

Figure 1

Comparison of theS part of SPRD6 and SKS12WM13. (A) Relative perturbations in S velocity at (from the top) 100 (± 4.5%) 600 (± 1.0%), 1300 (± 1.0%), 1800 (± 0.8%), 2300 (± 1.0%), and 2850 (± 3.0%) km, where the values in parentheses indicate the scale for the maps at each depth. TheS part of model SPRD6 is shown on the left, and model SKS12WM13 (1) is shown on the right. Blue colors indicate regions of higher than average velocities, and red colors indicate regions of slower than average velocities. Coastlines have been superimposed for reference. Because the maximum harmonic degree in SPRD6 is 6, we truncated degree 12 model SKS12WM13 at the same degree. (B) Root-mean-square (RMS) amplitudes of the S part of SPRD6 and SKS12WM13 as a function of depth. The red line is SPRD6, and the black line is SKS12WM13. (C) Correlation between the S model of SPRD6 and SKS12WM13 as a function of depth. For this number of model parameters, two maps are correlated at the 95% confidence level if the correlation is greater than 0.24.

We compared the P part of model SPRD6 withP velocity model P16B30 (3) (Fig. 2). The elongated velocity anomalies at mid-mantle depths underneath the Americas and Eurasia are characteristic of several recent travel-time models (23). The two P velocity models are well correlated at the top of the mantle and around 2500-km depth but show lower correlation in the mid-mantle and near the CMB. The P velocity part of SPRD6 agrees well with a recent model obtained by Kárason and van der Hilst (24).

Figure 2

Comparison of the P part of SPRD6 and P16B30. (A) Relative perturbations inP velocity at six discrete depths throughout the mantle, as in Fig. 1. The scale for the maps at each depth is (from the top) ±2.5, ±0.8, ±0.5, ±0.5, ±0.5, and ±0.7%. TheP part of model SPRD6 is shown on the left, and model P16B30 (3), truncated at degree 6, is shown on the right. The color scheme is the same as in Fig. 1. (B) RMS amplitudes of theP part of SPRD6 and P16B30 as a function of depth. The red line is SPRD6, and the black line is P16B30. (C) Correlation between the P model of SPRD6 and P16B30 as a function of depth. For this number of model parameters, two maps are correlated at the 95% confidence level if the correlation is greater than 0.24.

The laterally heterogeneous even-degree whole-mantle density model (Fig. 3 and Table 1) is generally poorly correlated with either theS or P velocity models. Near the CMB, high-density anomalies appear in areas of mantle uplift beneath the central Pacific and Africa. To further illustrate the nature of the anomalies near the CMB, we inverted for lateral variations in S velocity, bulk sound velocity, and density. Our results confirm earlier observations of an anticorrelation between shear and bulk sound velocity (25) (Fig. 4).

Figure 3

Relative perturbations in even-degree density at six discrete depths as in Figs. 1 and 2. The scale for each map (from the top) is ±1.0, ±0.5, ±0.5, ±0.5, ±0.8, and ±1.0%, respectively. Blue regions denote higher than average density, and red regions denote lower than average density. The maximum and minimum values of the density model are summarized in Table 1.

Figure 4

Shear velocity, bulk sound velocity, and density models at 2800-km depth. (A) Map views of shear velocity (top), bulk sound velocity (middle), and density (bottom) from an inversion with these model parameters (model SBRD6). The color scheme is the same as in the previous figures. The scale for the maps is ±2.5, ±2.0, and ±1.5% for S velocity, bulk sound velocity, and density, respectively. Note the strong anticorrelation between shear and bulk sound velocity (correlation coefficient of −0.63). (B) Same as in (A), but showing only the degree 2 component. The scales for the maps are ±0.7, ±0.4, and ±0.5%, respectively.

Table 1

The maximum and minimum values of the density model with the associated error at the given depth. The average value of the model is always zero because of the absence of a degree 0 coefficient. The errors are uniformly distributed over the sphere, which reflects the fact that modes are sensitive to structure everywhere on the globe. The error values are much smaller than the extrema of the model.

View this table:

The root-mean-square (RMS) amplitudes, or the power of lateral heterogeneity, of the S, P, and density models exhibit an increase near the surface and the CMB, with the exception of the P model, which is relatively flat in the lower mantle (Fig. 5A). The RMS amplitude of density is comparable to that of P in the lower mantle and exceeds it below 2000-km depth. As observed in recent travel-time models (3, 4), the S and Pparts of model SPRD6 are well correlated in the upper mantle but not as well in the mid-mantle (Fig. 5B). In the lower mantle, the correlation peaks around 2500-km depth and drops toward the CMB. Compared with the correlation between P and S, density does not correlate as well with either of the two velocity models. It is substantially decorrelated from the velocity models in the transition zone (400- to 700-km depth).

Figure 5

RMS amplitude and correlation. (A)S, P, and density RMS amplitudes of model SPRD6 as a function of depth. The red curve is the S model, the blue curve is the P model, and the green curve is the RMS amplitude of the density model. (B) Correlations between theS, P, and density parts of SPRD6 as a function of depth. Correlation between S and P is shown by the blue line, that between S and density is shown by the red line, and that between P and density is shown by the green line. The black line is the zero line. For this number of model parameters, two maps are correlated at the 95% confidence level if the correlation is greater than 0.24.

The last part of model SPRD6 consists of topography on the free surface, the 660, and the CMB (Fig. 6). The 660 part of model SPRD6 is compared with two recent travel-time models (20, 26), and our CMB model is compared with a seismological (27) and a geodynamical model (28).

Figure 6

Topography of the free surface, the 660, and the CMB. (A) Map of dynamic free-surface topography. Blue colors indicate areas of depression, and red colors indicate areas of elevation. The scale for the map is ±1.0 km. (B) Maps of 660 topography. The color scheme is the same as in (A), but the scale is ±15.0 km. The top map is part of model SPRD6, the middle map is from Gu et al. (20), and the bottom map is from Flanagan and Shearer (26). The latter two maps are truncated at harmonic degree 6 to ease visual comparison. (C) Maps of CMB topography plotted on a scale of ±5.0 km. The top map is from the SPRD6 inversion, the middle map is from Morelli and Dziewonski (27), and the bottom map is from Forte et al. (28). The bottom map is based on S velocity model S.F1.K/WM13 (8) and a depth-dependent density scaling proposed by Karato (39). The Morelli and Dziewonski map is expanded in spherical harmonics up to degree 4, and the Forte et al. map is truncated at degree 6.

Robustness of the Models

First, we examined the dependence of our models, especially the density model, on the starting model used in the inversion. Density models obtained from different starting models (with different scaling values) are consistent with one another (supplementary Figs. 1 and 2) (29), and the strong positive density anomalies near the CMB persist regardless of the starting model. The density model can also be influenced by the parameterization of the inversion due to the effects of damping. We performed inversions with three parameterizations: (i)S velocity, P velocity, and density; (ii)S velocity, bulk sound velocity, and density; and (iii) shear modulus, bulk modulus, and density. The density models obtained on the basis of these three inversions are compatible with one another, indicating that the effect of damping on parameterization is limited (supplementary Fig. 3). We also investigated results of inversions for density with and without boundary topography (supplementary Fig. 4). Introduction of topography reduced the amplitude of density heterogeneity near a boundary, but it did not substantially influence its pattern.

Next, we examined the resolution of the inversion. Resolution is dependent on the damping scheme (30) and the relative weighting between the free-air gravity anomaly and the mode data (31). We experimented with a variety of damping schemes, for example, size, gradient, and second derivative damping, as well as a range of weighting of the gravity anomaly, and found that the pattern of density heterogeneity did not change substantially. Furthermore, the number of resolved parameters for the density model, determined by the trace of the appropriate portions of the resolution matrix, is consistently greater than that of the S or Pvelocity models, regardless of the weighting of the gravity anomaly (this included an inversion with only normal-mode data).

Resolution in the lower mantle was investigated with input models peaked at 2700-km depth (Fig. 7A). In this test, density is the best resolved of the three models. TheS velocity is also well resolved, but the recovery of theP model is poor. We also investigated the components of Earth structure that contribute to produce an output model with a spike at 2800-km depth, known as a Backus-Gilbert resolution kernel (32) (Fig. 7B). Resolution tests (supplementary Fig. 5), Backus-Gilbert kernels (supplementary Fig. 6), and checkerboard resolution tests (supplementary Fig. 7) consistently show that the density structure is well resolved by the data.

Figure 7

Resolution tests. (A) The input is a sharp peak in S, P, and density heterogeneity at 2700-km depth at spherical harmonic degree 2 and order 1. (Left) Map views of models at 2700-km depth. The input models are shown in the first column, and recovered models are shown in the second column. The models are (from the top)S velocity, P velocity, density, and topography on the free surface, the 660, and the CMB. There is no initial topography on any of the boundaries, and the output shows that there is virtually no leakage of volumetric structure to topography on boundaries. (Right) Cross sections of the same resolution test for S velocity (top), P velocity (middle), and density (bottom). The models in the first column are input models; the output models in the second column show smearing of structure in the radial direction. (B) Backus-Gilbert resolution kernel (32) for structure at 2800-km depth for S velocity (S),P velocity (P), and density (R) heterogeneity at spherical harmonic degree 2 and order 1. The red curve represents the S velocity, the blue curve represents theP velocity, and the green curve represents densityR. The black line is the zero line. A Backus-Gilbert test indicates where structure in the inverted model has its origin in the real Earth. For example, ideally, the Backus-Gilbert resolution kernels in this figure should be spiked at 2800-km depth. Instead, we see that degree 2, order 1 structure at 2800 km in our models is due to Earth's structure between roughly 2400 and 2880 km and that there is a small contribution from heterogeneity near the surface.

The pattern of density heterogeneity that we obtained from the mode data alone is robust but has a range of amplitudes that satisfy the data equally well. Addition of the free-air gravity anomaly to the data set narrowed the amplitude range and slightly modified the pattern in the mid-mantle. It also introduced a trade-off between density heterogeneity and boundary topography. Of the three boundary topographic models, which include the free surface, the 660, and the CMB, the free surface is the least constrained. There is a wide range of acceptable topographic amplitudes with small-amplitude changes to the 3D density model and little change in the gravity fit. To obtain reasonable topography on the free surface, we damped the amplitude of our dynamic topography to a size (±1 km) between the observed dynamic topography (33) and that predicted from mantle-flow calculations (34). This resulted in dynamic topography resolution tests that failed to recover any input models. Topography on 660 is constrained better than the free surface, but the data are insufficient to give a reliable model by themselves; therefore, the 660 model was damped toward the model of Gu et al. (20) (supplementary Fig. 8A). The data do have sufficient sensitivity to undulations of the CMB, although there is a range of permissible amplitudes (supplementary Fig. 8B). We required the CMB to have an overall amplitude of ±5 km to be compatible with the amplitude determined by Morelli and Dziewonski (27).

The fit to the mode data was monitored as the number of degrees of freedom was increased from an S-velocity–only inversion to model SPRD6, with independent lateral variations in S,P, density, and boundary topography (supplementary Table 1and supplementary Fig. 9). The effects of uncertainties in the data on model structure were investigated by adding various levels of random noise to the splitting functions (supplementary Fig. 10). Even when random noise with an RMS of 30% was added to the data, the resulting models are highly consistent with model SPRD6.

Implications of the Models for Mantle Dynamics

We have demonstrated that normal-mode data and the free-air gravity anomaly can constrain independent large-scale lateral variations in S velocity, P velocity, and density. Our velocity models are consistent with existing travel-time models (1–4). The main differences are in the mid-mantle where travel-time and waveform models are least constrained.

Our whole-mantle density model has an RMS amplitude that is almost constant throughout the mantle, with a slight increase near the surface and the CMB. In the transition zone, the density model is substantially decorrelated from the velocity models, providing direct evidence that lateral variations in this region may not be entirely due to variations in temperature (35). The most prominent features of our density model appear near the CMB, where high-density anomalies are found in regions of mantle uplift beneath the Pacific and Africa. This observation is consistent with convection simulations where dense material piles up beneath upwellings because it is too heavy to be entrained in the uplift (36). Introduction of heavy core material into the mantle beneath upwellings can also cause such density anomalies (37). These anomalies, combined with the regional anticorrelation between shear velocity and density, suggest that compositional as well as thermal heterogeneity is required to explain lateral variations in the lowermost mantle, especially near the CMB.

  • * To whom correspondence should be addressed. E-mail: ishii{at}


View Abstract

Stay Connected to Science

Navigate This Article