Investigating Mantle's Density Resolution Using the
Neighborhood Algorithm

Sébastien Rousset and Barbara Romanowicz


Unlike travel times or waveform data, normal mode data are directly sensitive to density. However, the sensitivity kernels for density are much smaller than those for velocities, so the controversy about the possibility of resolution of the mantle's density is still vivid, especially since the publication of model SPRD6 (Ishii and Tromp, 1999). Several authors (Resovsky and Ritzwoller, 1999, Romanowicz, 2001, Kuo and Romanowicz, 2002) objected that density cannot yet be constrained and the controversy is still going.

However the inversion processes used by previous studies all rely on least-square inversions and require the use of a starting model, the choice of which is critical for the reliability of the results. Unlike this simple inversion scheme whose result is one "best" model, stochastic methods sample the parameter space, and their result is a set of models whose statistical properties reflect the likelihood function. In this study, in order to investigate the resolution, we use the Neighborhood Algorithm ( Sambridge 1999a,b). The first step of the algorithm generates a set of models that samples the parameter space preferentially where the fit is better. The second step approximates the posterior probability density (ppd) in the parameter space using the previously generated set and quantitative information is extracted from this approximate ppd by a Bayesian approach.

Data Set and Parameterization

Our data set consists in a set of splitting coefficients (Giardini and al., 1988) inverted from normal mode spectra. These coefficients are linearly related to the aspherical structure of the Earth considered as a perturbation $\delta x$ of the elastic parameter $x$ of our reference model PREM, integrated over depth:

C_{st} = \int_{0}^{a}(\frac{\delta x}{x})_{s,t}(z) K^{x}_{s}(z).dz
\end{displaymath} (37.1)

where s, t are degree and order in spherical harmonic expansion and a the Earth's radius. The sensitivity kernel $K^x_s$ is calculated for the reference model and depends only on the degree s of the expansion.

A first data set used 63 well constrained modes (He and Tromp, 1996 and Resovsky and Ritzwoller, 1998), corrected for the contribution of the crust; tests proved that the choice of the crustal model has little impact on the correction. Another data set adds a large number of upper mantle modes (Widmer, 2002); these modes are less well constrained but may improve the resolution in the upper mantle.

The Model Parameter Space

We search the degree 2 of the spherical harmonic expansion of the perturbation, with five coefficients $C_{2,0}$, $C_{2,1}$, $S_{2,1}$, $C_{2,2}$ and $S_{2,2}$. After trying various parameterizations for the radial variations of the perturbations, we selected a 7 (8 when more modes are included) cubic splines parameterization that gives naturally smooth variations.

To further reduce the number of parameters, only the shear velocity structure is fully discretized by splines. The bulk velocity and density structures are scaled to shear velocity structure with 3 scaling coefficients: lower mantle (bottom 2 splines), the middle mantle (3 splines) and one for the upper mantle (2 or 3 splines). Note that we are allowing lateral variations in the scaling coefficient as the scaling coefficients can be different for different spherical harmonic coefficients. Finally, the contribution of the topography of the CMB to the splitting of the modes is also included. Parameter space dimension is then 14 (resp. 15).


Sampling the Parameter Space

The neighborhood algorithm sampler (Sambridge 1999a) produces a set of models that sample the parameter space preferentially where the fit is good. Figure 37.1 shows 300 of the $C_{2,0}$ models for an exploration with 7 splines, $Vp$ and $\rho$ being scaled as described previously. While most models appear similar to the best model, a secondary minimum of misfit also appears that would be missed by a least square inversion. Note how the density in the uppermost mantle is poorly constrained in this case with a broad variety of values.

Figure 37.1: 300 good (low misfit), randomly selected C2,0 models. The black line in front is the SPRD6 model, color/grey is given by the misfit between synthetic coefficients and data.

Quantitative Results from the Bayesian Integrals

The second program, the appraiser (Sambridge 1999b) uses a set of models to approximate the posterior probability density in the parameter space, and computes Bayesian integrals over this approximate ppd, allowing a quantitative use of the set of models. It is then possible to get the marginal probability for each parameter - or quantities based on theses parameters. Such quantities have been estimated and are shown in Figure 37.2. The spherical harmonic coefficient shown is $S_{2,2}$, for a parameterization with 8 splines, scaling and topography of CMB. Note how the scaling values in the middle mantle (depth range 670 to 1800 km) are more poorly constrained than values in the lowermost and upper mantle. While the average shear velocity anomaly at these depths noted $ < V_s > mm$ is well constrained, the scaling coefficient for bulk velocity anomaly $V_p / V_s$ (for $d(ln V_p)/d(ln V_s)$ and density $Rho / V_s$ can take a broad range of values.


Further work will include stability tests for the method. Preliminary tests show that the confidence intervals may not be completely reliable (!), especially for some parameters for which convergence of the Bayesian integrals is slow and/or currently insufficient. The addition of the modes of Widmer (Widmer, 2002) seems to help to constrain the perturbation in the upper mantle but the effect on other values remains to be investigated.


Special thanks to Malcolm Sambridge for making the NA software package available. This package, its very helpful online help and a short description of the algorithm can be found at Figure 37.2 was made using a plotting utility of the package.


Giardini, D., X.D. Li and J. Woodhouse, Splitting functions of long-period normal modes of the Earth, Jour. Geophys. Res., 93 B11, 1988.

He, X. and J. Tromp, Normal-mode constraints on the structure of the Earth, J. Geophys. Res., 102, 17981-17994, 1996.

Ishii, M. and J. Tromp, Normal-Mode and Free-Air Gravity Constraints on Lateral Variations in Velocity and Density of Earth's Mantle, Science, 285, 1231-1236, 1999.

Kuo, C. and B. Romanowicz, On the resolution of density anomalies in the Earth's mantle using spectral fitting of normal mode data, Geophys. J. Inter., in press

Romanowicz, B, Can we resolve 3D density heterogeneity in the lower mantle ?, Geophys. Res. Lett., 28, 15, 1107-1110, 2001.

Resovsky, J. and M. Ritzwoller, Regularization uncertainty in density model estimated from normal mode data, Geophys. Res. Lett., 26, 15, 2319-2322, 1999.

Resovsky, J. and M. Ritzwoller, New and refined constraints on 3D Earth structure from normal modes below 3mHz., J. Geophys. Res., 103, 783-810, 1998.

Sambridge, M., Geophysical inversion with a neighborhood algorithm I - Searching a parameter space, Geophys. J. Int., 138, 479-494, 1999.

Sambridge, M., Geophysical inversion with a neighborhood algorithm II - Appraising the ensemble, Geophys. J. Int., 138, 727-746, 1999.

Widmer-Schnidrig, R., Application of regionalized multiplet stripping to retrieval of aspherical structure constraints, Geophys. J. Int., 148, 1-18, 2002.

Figure 37.2: 2D marginal probabilities for the $S_{2,2}$ anomalies. Average of the shear velocity anomaly is noted <Vs>, scaling coefficient from Vs to density (resp. Vp) rho/vs (resp. Vp/Vs). lm stands for a value in the lower mantle, mm middle and um upper mantle.

Berkeley Seismological Laboratory
215 McCone Hall, UC Berkeley, Berkeley, CA 94720-4760
Questions or comments? Send e-mail:
© 2004, The Regents of the University of California