## Abstract

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 in*S* 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).

## Theory

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 2*l*, 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 harmonics*Y*
_{s}
^{t} of degree *s* and order*t* as (12)(1)The splitting-function coefficients *c*
_{s}
^{t} are linearly related to 3D relative variations in *S* velocity (δβ/β)_{s}
^{t}, *P*velocity (δα/α)_{s}
^{t}, density (δρ/ρ)_{s}
^{t}, and topography on discontinuities (δ*h*/*h*)_{s}
^{t} by (12)
The 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 kernels*K*
_{β}, *K*
_{α},*K*
_{ρ}, and *K _{d}
* 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

*S*velocity. 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)(2)The kernels *K*
_{ρ}
^{′} and*K*
_{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 in*S* 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) and*P* 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, where*N*
_{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.

We compared the *P* part of model SPRD6 with*P* 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).

The laterally heterogeneous even-degree whole-mantle density model (Fig. 3 and Table 1) is generally poorly correlated with either the*S* 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).

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 *P*parts 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).

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).

## 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 *P*velocity 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. The*S* velocity is also well resolved, but the recovery of the*P* 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.

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}seismology.harvard.edu