next up previous contents
Home: Berkeley Seismological Laboratory
Next: Towards Forward Modeling of Up: Ongoing Research Projects Previous: Progress in modeling deep   Contents

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 likehood function. In this study, in order to investigate the resolution, we use the Neighbourhood Algorithm ( Sambridge 1999a,b) which uses forward modeling to sample the parameter space preferentially where the fit is better. A second part of the algorithm can use this sampling set to retrieve quantitative information on the distribution of models.

Dataset and Model Parameterization

Theoretical background and dataset

We use a 2 steps inversion, as was done by Ishii and Tromp ; this method is simpler, but the coefficients measured may not correspond to a unique Earth model. Some studies (Kuo and Romanowicz, 2002) prefer a 1-step inversion, but this method is not practicable with the neighbourhood algorithm since the forward modeling of synthetic seismograms for each model would require too much computation.

During the first step, splitting coefficients (Giardini and al., 1988) of normal mode were obtained by non-linear inversion. Our dataset consist in splitting coefficients of 63 modes (He and Tromp, 1996 and Resovsky and Ritwoller, 1998), corrected for the crust contribution and the effects of rotation. We have excluded a few modes for which 2 measurements were available but incompatible.

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})_{st}(z)\,K^{x}_{s}(z).dz
\end{displaymath} (35.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.

In this study, we limit ourselves to degree 2 perturbations and invert the 5 degree 2 splitting coefficients $C_{20}$, $C_{21}$, $S_{21}$, $C_{22}$ and $S_{22}$ of the selected modes for the corresponding degree 2 structure of $V_p$, $V_s$ and $\rho$ in the mantle.

Inversion settings

The neighbourhood algorithm (NA) is a stochastic method comparable to Monte-Carlo method, genetic algorithm or simulated annealing (Moosegard and Sambridge, 2002). We are using the first part (Sambridge 1999a) of the algorithm to sample the parameter space preferentially in good fitting areas, using the geometrical construction of Voronoi cells (closest neighbourhood) to fill the space. At each iteration, $n_s$ new models are randomly generated in the $n_r$ best fitting cells, the cells being defined by all previously generated models. After a sufficient number of iterations, we obtain a large (about 150 000 for our problem) set of models with more models in good fitting areas. The parameters $n_s$ and $n_r$ must be tuned for each problem ; an important value is the parameter space dimension which largely governs the choice of $n_s$ and $n_r$ and seems to be the main limitation of the use of the algorithm : the concentration of models in good fitting areas failed for high dimension, at least in reasonable computation time.

Model parameterization

The model is naturally laterally parameterized by spherical harmonics, but we have the choice of the radial parameterization : layered models (equivalent to gate functions), splines functions, polynomials... Several trials show that 6 to 8 functions are a sufficient number to retrieve feature while keeping the total number of parameters low enough for use with the NA. Layered models oscillate and forcing them to be smooth causes instability and introduces arbitrary constraints that we would like to avoid. The splines-based parameterization that we selected allow smoother models.


The most acceptable models are obtained when adding a spline-based model to the contribution of tomographic input models (Figure 35.1). The $V_s$ perturbation is well retrieved while the $V_p$ perturbation doesn't agree well with the input. The $\rho$ perturbation is still oscillating, although smoother acceptable models exist.

Figure 35.1: Best $V_s$ model (black) obtained when using the tomographic (grey) model as input and 7 splines, for spherical harmonic coefficient $C_{20}$
\epsfig{file=rousset02_model.epsi, width=4cm}\end{center}\end{figure}

One noteworhty feature of all parameterizations is the good fit between the $V_s$ model obtained in our inversion of splitting coefficients and the $V_s$ tomographic model, even when this model is not used as input, the difference residing mostly in the uppermost mantle where simple parameterization like ours is more likely to fail.

Figure 35.2 show the distribution of models for the $C_{20}$ contribution of the share velocity perturbation. This distribution is stable and satisfactory, giving good confidence in the uniqueness of the $V_s$ solution, but is obtained using tomographic models as input and the stability if the result with respect to the choice of the input model remains to be tested.

Further work will make use of the second part of the algorithm (Sambridge, 1999b) which use Bayesian integrals to analyze quantitatively the set of model obtained during the first part (or any other method generating a large set of models).


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

We thank the National Science Foundation for support of this research.


Giardini, D., X.D. Li and J. Woodhouse, Splitting functions of long-period normal modes of the Earth, J. 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

Li, X.D., D. Giardini and J.H. Woodhouse, Large-Scale Three-Dimensional Even-Degree Structure of the Earth from Splitting of Long Period Normal Modes, J. Geophys. Res., 91 B1, 551-577, 1991.

Moosegard, K. and M. Sambridge, Monte Carlo Analysis of inverse problems, Inverse Problem, 18 R29-R54, 2002.

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, Regularisation 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 neighbourhood algorithm I - Searching a parameter space, Geophys. J. Int., 138, 479-494, 1999.

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

Figure 35.2: Projections of the set of $V_s$ models obtained as output of the NA program ; Vs[x] is the value of the coefficient for the $x^{th}$ spline, 1 being the deepest and 7 the most shallow spline.

next up previous contents
Berkeley Seismological Laboratory
215 McCone Hall, UC Berkeley, Berkeley, CA 94 720-4760
Questions or comments? Send e-mail:
© 2002, The Regents of the University of California.