Oceanic upper mantle structure from experimental scaling of V_S and density at different depths

Summary

Isotropic shear-wave velocities from the Ekström/Dziewonski S20 global tomographic model and residual ‘crust-free’ gravity are used to determine density–velocity relations for various mantle layers by inversion. The residual gravity field was obtained by subtracting the crustal density inhomogeneities from the observed surface gravity. The residual field is dominated by upper mantle signatures. As the effect of dynamic topography is also reduced from the gravity field before-hand, the gravity–seismic inversion is almost independent of the viscosity distribution within the Earth. There exists a reliable correlation between residual gravity and shear-wave velocities in the upper mantle under oceans but not under continents. This result supports the idea that lateral density variations in the oceanic upper mantle are mostly driven by the temperature distribution, while under the continents temperature- and composition-induced density anomalies are superimposed. The inversion is therefore restricted to the oceanic area and was done by a least-squares adjustment in the spatial domain.

The density/velocity ratios obtained for the subcrustal layer (down to 75km) and for the deep part of the oceanic upper mantle (below 225km) are very close to the ones predicted from mineral physics analyses. However, for the oceanic layer from 75 to 225km depth, a significantly lower density/velocity ratio is obtained. This may be explained by a compositional layering: dry, partially depleted peridotite residuum underlain by ‘damp’ fertile peridotite and separated by the G discontinuity at a depth of about 70km. For the mantle transition zone and for the lower mantle the results obtained are less reliable, but relative minima in the density/velocity ratio are found for the layers around the phase transition discontinuities and for the depth interval from about 900 to 1100km. Applying the adjusted velocity-to-density relations to the seismic tomography model, a mantle density distribution under the oceans is found which explains about 85 per cent of the longest-wave ‘crust-free’ gravity field (degrees 1–5) and 70 per cent of the overall field up to degree 20, which was the limit of this investigation.

density, Earth's interior, gravity, seismic tomography, upper mantle

1 Introduction

The determination of the Earth's mantle density structure is important for a number of geophysical problems: density differentiation in the Earth's mantle is the driving source for mantle convection and for plate tectonic processes at the Earth's surface; and from the point of view of mineral physics the knowledge of the densities within any Earth layer gives additional constraints to determine properties of the material.

Unfortunately, the determination of the density structure of the Earth's mantle has an ambiguous solution when using only surface gravity data, which depend directly on the density distribution but as a result of an integration over the whole Earth's body. The usual way to overcome this difficulty is to combine gravity data with other geophysical data sets to produce a joint solution that fits all data sets and thus contains fewer degrees of freedom. Seismic tomography models are usually used for these purposes. In global studies, this method has been applied since the pioneer works of Hager & O'Connel (1981), Ricard et al.(1984), Richards & Hager (1984), Forte & Peltier (1987), King & Masters (1992), Corrieu et al.(1994) and many others. Three principal difficulties have to be solved when jointly analysing gravity and seismic mantle tomography data. First, the observed geoid and gravity variations reflect not only internal mantle density inhomogeneities but also the crustal structure and dynamic topography that is produced by flows deforming the Earth's surface. Therefore, in order to separate the sources it is necessary to model the Earth's crust and simultaneously Earth dynamics, which depend on the viscosity distribution inside the Earth and on the type of mantle flow (whole mantle or layered convection), a topic still under discussion. Second, the long-wavelength observed gravity field and geoid have, at least for the harmonics up to 20, almost no correlation with lithospheric inhomogeneities and discontinuities because, due to mass compensation in the crust and upper mantle, the total gravitational effect is close to zero. As a result, most previous global inversion studies did not provide sufficient resolution for the upper mantle, either in the horizontal or in the vertical direction. Third, density variations, on the other hand, are only partly reflected in seismic velocity variations, which are primarily related to the temperature distribution rather than chemo-physical material properties, making the velocity-to-density scaling an inherent problem. This is especially important for the continental upper mantle.

A joint interpretation of gravity and seismic tomography data therefore implies a solution of the velocity-to-density scaling problem. It is assumed that the isotropic part of the velocity variations determined by seismic methods is mainly associated with temperature changes (e.g. Karato 1993), although the degree of this relation depends on the actual temperature–pressure state: the velocities vary much faster in the vicinity of the melting point, which changes with the composition of the material (Murase & Fukuyama 1980; Sato et al. 1989). The relation between temperature and density as described by thermal expansion is almost linear over the entire temperature range, but density is more sensitive than velocity to changes in the composition, for example those arising from the degree of upper mantle depletion (Jordan 1978). When deriving the relation between velocity and density perturbations from theoretical mineral physics equations, one has to take into account all these factors, which leads obviously to a complex mantle density determination based on seismic tomography data (e.g. Karato & Jung 1998). When determining the velocity-to-density scaling using seismic tomography and gravity data, one can easily draw conclusions about the composition and temperature at various depths afterwards.

From an analysis of recent studies about the derivation of velocity-to-density scaling factors from seismic mantle tomography and gravity data as carried out by, for example, King & Masters (1992), Ricard et al.(1993), Forte et al.(1993, 1995), Corrieu et al.(1994), Mitrovica & Forte (1997), Thoraval & Richards (1997), Wen & Anderson (1997a) and Pari & Peltier (1998), the following general remarks can be made.

• 1

  The observed geoid, which includes a significant dynamic component, is usually taken to constrain the investigations. It is thus necessary to determine simultaneously the viscosity and density structure or to adopt an a priori model for one of them as well as to select one of the mantle flow models.

• 2

  Most of the observed geoid signal (about 70 per cent) is concentrated at the lowest-degree harmonics, where the effect of the lower mantle is dominant. The signal from the upper mantle is very small because upper mantle density inhomogeneities are to a large extent compensated by opposite crustal density inhomogeneities, for example under oceans by seafloor topography. Consequently, when using the observed geoid as a constraint, the resulting scaling coefficients for the upper mantle are determined with a significant error and tend to be smaller than those determined in experiments and predicted by mineral physics equations (Sobolev et al. 1996; Karato & Jung 1998). This is a result of the missing correlation between seismic tomography models of the upper mantle and the long-wavelength geoid.

• 3

  The inversion is usually performed in terms of spherical harmonic coefficients, which causes difficulties in discriminating between regimes that are different in the space domain. Under continents, a correlation between density and velocity variations should not necessarily be expected because of the possible superposition of lateral thermal and material composition changes. Forte et al.(1995) divided the gravity and seismic velocity fields into two components, namely correlated and not correlated with the continent–ocean distribution, and made separate inversions for these two components. The first solution therefore characterizes the main differences between subcontinental and suboceanic mantle, while the second one is attributed by the authors to be primarily thermal in origin. However, within continents there may be significant variations of the mantle composition (for example between young and old continental lithosphere), and therefore the fields not correlated with the continent–ocean distribution may contain signals of a mixed origin.

The main problem in the determination of the upper mantle's velocity-to-density scaling coefficients and density structure is to free the observed gravity field from the crustal signal, which due to isostatic compensation largely masks the gravitational signal from upper mantle density inhomogeneities. Large-scale crustal structures are isostatically compensated, but the mechanisms and depths of compensation are still under debate. By determining the upper mantle density structure our knowledge of isostatic compensation can be significantly improved.

In this study, based upon the most recent global compilation of the thickness of crustal layers and density data, crust5.1 (Mooney et al. 1998), and topography/bathymetry data, residual gravity has been calculated by subtracting the ‘known’ crustal lithosphere effects from the observed surface field (Kaban et al. 1999). By doing this, the effect of dynamic topography as part of the total perturbations of the near-surface boundaries is no longer present in the residual ‘crust-free’ gravity field and geoid, which now reveal only the impact of subcrustal density inhomogeneities. The same is true for the effect of postglacial rebounds. This makes it possible to model the mantle density structure almost independently of its viscosity distribution.

The variations in the residual gravity field are, due to unreduced compensating masses in the upper mantle, much larger than in the initial (observed) field, and they are now obviously correlated with the seismic velocity distribution in the upper mantle (Fig.1).

Figure 1

Two profiles of upper mantle V_S seismic velocity perturbations (per cent), observed gravity (dashed line), and residual ‘crust-free’ gravity variations (solid lines) along latitudes 22.5°S and 42.5°S. Seismic velocities are from the S20 tomographic model (Ekström & Dziewonski 1998), observed gravity is from the GRIM4-S4 model (Schwintzer et al. 1997), and residual gravity from the field shown in Fig.2. The gravity fields are truncated after degree 20 to be of the same resolution as the tomographic model.

The joint analysis of the residual gravity field and seismic tomography data then allows us to investigate in detail the density structure of the upper mantle by determining the most appropriate velocity-to-density scaling coefficients. These results are more reliable than those based upon the observed and unreduced gravity field or geoid.

The present study concentrates on oceanic areas. Only under oceans can one expect to explain most of the lateral density inhomogeneities within a specific layer in the upper mantle by temperature variations, and consequently to recover the density structure using seismic tomography and gravity data. Furthermore, the oceanic crust is more homogeneous than the continental crust, and therefore the residual ‘crust-free’ gravity field over oceanic areas is less affected by shortcomings in the crustal model.

When investigating velocity-to-density relations, the horizontal layering of the oceanic upper mantle is an essential issue. It has been proposed by Gaherty et al.(1996, 1999) and by Hirth & Kohlstedt (1996) that a boundary, named the G discontinuity, is located at a depth of about 70km under oceans, and that it is a compositional boundary placed at the depth of partial melting during the production of the overlying oceanic crust. The resulting compositional layering—dry, partially depleted peridotite residuum underlain by ‘damp’ fertile peridotite—persists as the plate translates away from the ridge. The existence of this boundary leads to principal geodynamic consequences because the damp peridotite will be two orders of magnitude less viscous than the overlying ‘dry’ harzburgite owing to the water-enhanced mobility in olivine (Hirth & Kohlstedt 1996). Evidence for the G discontinuity comes from seismic data, which show a significant decrease of V_S under oceans at depths between 70 and 200km (e.g. Gaherty et al. 1996). The determination of the density structure under oceans could provide valuable information regarding this discontinuity. If the assumed layering in the oceanic mantle is real, the velocity-to-density scaling will be remarkably different for different depths, although the limited resolution of the tomography data does not allow the exact depth of this discontinuity nor the type of transition across the layers to be fixed.

There is also debate concerning the ocean spreading model. According to the well-known cooling oceanic lithosphere model, major temperature differences occur at depths above 130km (e.g. Hager 1983; Cazenave et al. 1986). Existing seismic tomography models, however, indicate that the velocity anomalies correlated with the mid-ocean-ridge structures extend at least to a depth of 200km (Gaherty et al. 1996; Ekström & Dziewonski 1998). Moreover, the most pronounced anomaly occurs at depths between 120 and 200km (cf. Fig.1). This would be an argument to increase drastically the thickness of the cooling plate up to 210km (Wen & Anderson 1997b) or to modify the classic cooling oceanic lithosphere model towards an active spreading model resulting from a super-plume. It is therefore important to investigate the suboceanic density structure by an experimental determination of the velocity-to-density scaling at various depths in order to support one or the other model assumption. If the velocity anomaly at larger depths were to correspond to a comparable density anomaly, it would mean that the real dynamics differ from the classical cooling oceanic lithosphere model. Alternatively, the strong and deep velocity anomaly may not correspond to an equivalent density anomaly, because of the above-mentioned perturbations due to the effect of the layer below the G discontinuity, and thus the classical model would not be contradicted.

The main purpose of this paper is to construct a mantle density model that is consistent with the residual ‘crust-free’ gravity field and the most recent V_S tomography model presented by Ekström & Dziewonski (1998). Experimental velocity-to-density scaling factors will be determined for the principal mantle layers. The experimental velocity-to-density relation over a given depth range should characterize the general composition (variation of scaling factors with depth) and thermal state (lateral density variations) of the specified layer. The study concentrates on suboceanic upper mantle structure.

2 Gravity, crust and seismic tomography: preparatory considerations

2.1 The residual ‘crust-free’ Earth gravity field

The gravitational effect of topography, bathymetry, sediments and consolidated crust down to the Moho relative to a horizontally homogeneous reference density model was calculated as a 3-D direct gravity problem solution, taking into account density anomalies in the lateral and vertical directions. This crustal field was then subtracted from the observed one to give the residual ‘crust-free’ gravity field (Fig.2), as described in Kaban et al.(1999). The reference density model corresponds to a ‘normal’ old (200Ma) oceanic lithosphere. Its choice affects only the mean value of the residual field and is of no importance for this study. The observed gravity field is computed in terms of free-air gravity disturbances instead of gravity anomalies in order to account for the indirect effect of the undulation of the geoid (Heiskanen & Moritz 1967; Chapman & Bodine 1979).

Figure 2

Residual (crust-free) gravity field (mgal) truncated after degree 20 (Kaban et al. 1999) and centred by subtracting the overall mean value.

When the crustal field is subtracted, the amplitudes in the residual field are increased compared with the observed gravity field, as can be deduced from Fig.1. The global root mean square (rms) of the residual gravity disturbances exceeds the rms of the initial field by a factor of 4, that is, 85 compared with 22mgal. This is true for both oceanic and continental areas. All fields are restricted to a spectral resolution up to degree/order 20. The signal magnification is due to the fact that long-wavelength crustal anomalies are to a large extent compensated in the upper mantle below the Moho (Kaban et al. 1999), resulting in the observation of only a small net gravitational signal at the Earth's surface. The crustal reduction can be compared with a Bouguer correction extending from the surface throughout the whole crust. The increase in the rms by a factor of 4 means that the topography/bathymetry and crustal model reflect the upper mantle density structure and form the basis for interpretation, rather than the observed gravity. It is also possible to use just the topography and crustal data, to compute a residual topography as in Kaban et al.(1999). Contrary to the case of residual gravity, however, we have to apply some additional initial suppositions, for example about isostasy, and to take into account dynamic effects such as postglacial rebounds when interpreting the residual topography.

The following components contribute to the overall error budget of the residual ‘crust-free’ gravity disturbances: errors in

• 1

  observed gravity and topography;

• 2

  crust model with

• (2.1) depth to Moho,

• (2.2) thickness and density of sediments, and

• (2.3) density of consolidated crust; and

• 3

  adopted mean density of the upper mantle.

In order to quantify statistically the propagation of these error sources to the residual gravity disturbances, the relation

 

1

(Strang van Hees 1993) is used, connecting an anomalous mass of density δρ and thickness δh with surface gravity δg in terms of spherical harmonics of degree l and order m. G denotes the gravitational constant. Eq (1) gives the first-order gravitational effect of near-surface anomalous mass distributions.

It has to be kept in mind that this study is restricted to long-wavelength features up to degree/order 20 in the spectral domain, corresponding to a minimum wavelength of 2000km in the spatial domain or a 10°×10° grid resolution; that is, insufficiently recovered small-scale anomalous features (e.g. Ninetyeast Ridge) in the crustal and gravity models are filtered out and grid-value errors are reduced due to averaging.

The accuracy of topography/bathymetry data is hard to estimate, but, based on Topographic Science Working Group (1988) and Smith (1993), an error standard deviation of 50m for block-averaged grid values may be realistic. Introducing constant density contrasts of 2.7 and 1.7gcm⁻³ with respect to atmosphere and water, respectively, application of eq. (1) yields error standard deviations in gravity of 3.4 and 5.4mgal for oceanic and continental areas, respectively. Therefore, the error of 50m is distributed over the harmonics up to degree and order l,m=20 (δh_l,m=2.5m), and the resulting δg_l,m are then accumulated over l,m to form a rms value representing the standard deviation in computed gravity. In the same way, the following error sources are propagated to gravity.

The Moho depths have an error standard deviation for normal oceanic areas and well-studied continental areas of about 500m. With a mean density contrast of 0.4gcm⁻³ between crust and upper mantle, this results in an error contribution of 8.0mgal. For oceans, the Moho depth error estimate of 500m is based on White et al.(1992) and Mooney et al.(1998). For well-studied continental areas, about 100 independent seismic Moho depth determinations over one 10°×10° grid cell are available for averaging: thus a 500m error estimate is realistic for the US, South Canada, Australia and most parts of Eurasia, even if a single determination has an error of several kilometres.

A change in the adopted mean mantle density of 0.05gcm⁻³ yields an additional rms in gravity of 5.0 and 14mgal for oceanic and continental areas, respectively, applying an overall rms of 2.5km for oceanic areas and 7km for continental areas in Moho depths. This error source also changes the mean gravity over oceans with respect to continents due to the different level in mean Moho depths, but this is of only minor concern for this study, which concentrates on oceanic fields only.

As the density of sediments approaches consolidated crust densities at depths of greater than 5km (e.g. Beyer et al. 1985; Kennet et al. 1994), errors in the sedimentary thickness are not the driving error source. The density distribution within sedimentary layers is assumed to be known to an accuracy of about 15 per cent (Kaban & Mooney 2001). An error standard deviation of 0.1gcm⁻³ in the density contrasts of distribution of the modelled sediments leads to an error contribution in the computed gravitational effect of 1.6mgal for an average sedimentary thickness of 400m in the normal ocean, and 8mgal for an average sedimentary thickness of 2km on the continents.

Error standard deviations in the mean density of consolidated crust of 0.02gcm⁻³ (Carlson & Herrick 1990) for oceanic and 0.03gcm⁻³ for continental crust over an average crustal thickness of 7km (oceans) and 35km (continents) result in gravitational error standard deviations of 5.6 and 42mgal, respectively. Uncertainties in the crustal density distribution are therefore the dominating error source for continental areas when reducing the crustal gravitational field from the observed field. The error standard deviation for continental crust density is derived from the uncertainty given in Christensen & Mooney (1995) for each particular crustal layer. This value is reduced when averaging the densities over the three layers of consolidated crust given in the crust5.1 model.

Combining all error sources and considering that the observed long-wavelength gravity field has an accuracy of better than 1mgal (Schwintzer et al. 1997), the estimated standard deviations in the residual gravity disturbances amount to 12mgal for most parts of the oceans and 46mgal for well-studied continental regions. This means that the signal-to-noise ratio is about 7 over oceans and less than 2 over continents, which is an additional argument for basing this study only on oceanic fields.

These error propagation results give a rough idea of the overall accuracy in residual long-wavelength gravity disturbances after having reduced the gravitational effect of the Earth's crust based on the most recent crustal model. The errors increase for regions with weak data coverage and localized anomalous crustal features revealing significant deviations from the general assumptions made above. To gain insight into the geographical distribution of the errors in the residual gravity disturbances, a global distribution of location-dependent estimates for the error sources is introduced and propagated to errors in the computed crustal gravitational field.

The following deviations from the average values assumed above are applied to give a regionally differentiated error distribution. According to White et al.(1992), Moho depth uncertainties are, for example, smaller for normal crust in the Pacific and larger for very young oceanic crust and the Indic. For anomalous oceanic crust (for example thickened and melt-affected crust, plateaux, marginal basins) an error of 10 per cent of its thickness and an increased density error (0.03gcm⁻³) are adopted. As in crust5.1, the unresolved structure of the plateaux close to New Zealand are not considered in this error analysis.

On continents, errors in Moho depths are increased to 1km under tectonically stable areas and to 2km under large-scale orogens if the areas are only sparsely covered with seismic data (part of Africa, the southern part of Asia, part of South America, part of Canada). If no data are available in these continental areas, the crustal model crust5.1 gives average values derived from well-studied structures of a similar kind. Consequently, in these cases the rms of the Moho depths (from 3.3km for Archean shields to 7.3km for Alpine-type orogens) and densities (∼0.055gcm⁻³) for a particular structure are taken as the error standard deviation. For sediments, an error of 15 per cent of its gravitational effect is always assumed.

All these sources of crustal mass and density uncertainties, together with the error values given above for topography/bathymetry and upper mantle mean density, are then propagated to yield the error pattern in the residual ‘crust-free’ gravity field in terms of error standard deviations, as shown in Fig.3.

Figure 3

Error standard deviations (mgal) of residual gravity disturbances (Fig.2) due to crust reduction.

On continents, the values range from 30 to 120mgal, almost exceeding the signal in Fig.2. For oceans, the error distribution looks much more homogeneous and has a level that is low enough for the purpose of this study. The errors vary from 8mgal over normal oceanic crust to 40mgal over anomalous oceanic crust. The indirect effect of the large uncertainties in the continental crust model when computing residual gravity disturbances over oceans has proved to be insignificant.

2.2 Seismic tomography model

Ekström & Dziewonski (1998) recently developed a high-resolution tomography model (S20) of shear-wave velocity (V_S) variations which is split into isotropic and anisotropic parts. In previous studies in which tomography models were used together with gravity to determine velocity-to-density scaling factors, it was assumed that the effect of anisotropy was small and that velocity perturbations only reflected the spatially inhomogeneous temperature and composition field. Ekström & Dziewonski (1998), however, showed that the effect of anisotropy in the upper mantle under oceans produces velocity anomalies that are at least as large as those due to temperature effects. The effect of anisotropy also explains the disagreement in velocities for the upper mantle among different tomography models, for example S20 versus S12 of Su et al.(1994), which are based on different groups of data with varying sensitivity to vertically and horizontally polarized shear-wave velocities. Previous attempts to determine the anisotropy of the mantle on a global scale were limited to a few low-degree harmonic coefficients (e.g. Montagner & Tanimoto 1991). The model S20 provides a unique possibility to study the upper mantle structure on a global scale with a considerably higher resolution. It is also important for this study that the S20 tomography model includes a correction for the impact of the Earth's crust using basically the same data (Mooney et al. 1998) as were used for the calculation of the residual gravity field. A consistent reduction is especially critical for the upper mantle because of the significant interference between crustal and mantle effects in both seismic waves and gravity.

For the purpose of this study, mantle density distribution by velocity–gravity inversion, only the isotropic part of the S20 seismic model is used. The Ekström & Dziewonski (1998) model has a spherical harmonic resolution up to degree and order 20 and a vertical resolution of about 50km in the upper mantle. Velocity variations with depth along constant latitudes and horizontal variations at constant depths resulting from the model are given in Figs1, 4 and 5, respectively.

Figure 4

V_S perturbations (per cent) for the spherical surface at 50km depth, as given by the S20 tomographic model (Ekström & Dziewonski 1998).

Figure 5

V_S perturbations (per cent) vertically averaged between 75 and 225km depth, as resulting from the S20 tomographic model (Ekström & Dziewonski 1998).

2.3 Correlation between residual gravity and seismic velocity signatures

The development and interpretation of tomographic models of the Earth's mantle have usually proceeded under the assumption that the seismic velocity variations represent a spatially heterogeneous temperature field associated with mantle convection. Comparing the velocity pattern in Figs4 and 5 with the pattern of the residual gravity field in Fig.2 reveals similarities for oceanic areas, whereas under continents the similarities are not as obvious. This statement is quantitatively supported by the values of the correlation coefficients between the residual gravity field (integrated over all depths) and the isotropic seismic anomalies extracted from the tomographic model at specific depths, as given in Fig.6. High correlation coefficients exceeding 0.6 occur for depths between the Moho and 225km under oceans. Under continents, the correlation is significantly smaller and almost zero in the subcrustal layer. We can conclude that under the oceans in the uppermost mantle the lateral variations of the velocity and density are principally controlled by the same factor, namely the temperature field. Under continents, density variations arise from temperature and chemical composition, and in many cases compensate each other. Kaban et al.(1999) showed, based on an analysis of residual topography and gravity, that there are almost no places under the continents where the average density of the upper mantle is greater than that under old ‘normal’ ocean. All seismic tomography models show a remarkable difference between the continental upper mantle (especially cold and therefore high-velocity cratons) and oceanic upper mantle. Owing to the missing correlation with the gravity field, this supports the hypothesis of Jordan (1978) that density anomalies under continents are largely controlled by the presence of a low-density peridotite depleted in its basaltic constituents.

Figure 6

Correlation coefficients between the residual ‘crust-free’ gravity field (Fig.2) and V_S perturbations over depth. Seismic velocities are taken from the S20 tomographic model (Ekström & Dziewonski 1998).

Fig.6 shows the unusual behaviour of the correlation coefficient above the transition zones: a negative correlation peak is found for the layer on top of the 670km boundary, while the correlation coefficient is close to zero above the 410km discontinuity. This could be the result of temperature anomalies, mantle flows and phase transition boundary deformations, and their composite influence on overall gravity and seismic velocities.

3 Experimental velocity-to-density scaling results

3.1 Inversion technique

An initial density model is produced by a simple linear conversion of the seismic tomography velocity disturbances into density variations: dρ=c_j⁰dV_S/V_S. The gravitational impact of each layer j is then propagated to the Earth's surface for use in the inversion. The choice of the initial scaling factors (c_j⁰) is arbitrary and does not change the results because these factors are rescaled in a least-squares adjustment. The inversion to solve for the unknown scaling factors is performed in the space domain, with 10°×10° equal-angular grid values computed from the spherical harmonic representation of the fields. This type of inversion rather than an inversion in the spectral domain with spherical harmonic coefficients is preferred because it facilitates the restriction of the inversion to selected areas. Computational tests have proved, that, if the area is large enough, an inversion in the spectral domain, namely in terms of spherical harmonic coefficients which are appropriately filtered by a convolution with an area function, gives essentially the same results as working directly in the spatial domain. The scaling coefficients a_j (to rescale the initially adopted values c_j⁰) are estimated in a least-squares adjustment starting with the linear observation equations

 

2

supplemented by the pseudo-observation equations with respect to the unknowns a_j:

 

3

where dgⁱ_res is the residual gravity disturbance at point i, dg_jⁱ the gravity disturbance due to the initial density (c_j⁰ applied) distribution in layer j propagated to the surface point i, a_j the unknown scaling coefficient for layer j, and cosφⁱ a formal weight factor to account for the convergence of the meridians in the equal-angular grid. The coefficients β_jk control the selection and strength of regularization by limiting the difference of the scaling coefficients between layers j and k (damping). Both the residual gravity and layer-wise propagated gravity are centred to zero mean over the area under consideration; p denotes the weighting and ε the residuum. The least-squares principle, that is minimization of the function

 

4

leads to the normal equation system

 

5

where the elements of the normal matrix B are b_jk=Σ_icosφⁱdg_jⁱdg_kⁱ + s_jk, where s_jk = −β_jk (j ≠ k) and s_jj = Σ_lβ_jl (j = k), and the components of the right-hand-side vector c are Σ_icosφⁱdgⁱ_resdg_kⁱ. System (5) is solved for the vector â by inversion of the normal matrix B. The solved-for parameters â_j are a posteriori transformed to the convenient dimensionless ratio

 

6

by computing an average density ρ¯_j for each layer from PREM (Dziewonski & Anderson 1981). The density–velocity scaling coefficient standard deviations s_ĉj resulting from the fit in the least-squares adjustment are then

 

7

where q_jj is a diagonal element of B⁻¹, and f denotes the degree of freedom [number of equations (2) plus (3) minus number of unknowns].

The pseudo-observation equations (3) constrain the variations of the scaling factors from one layer to another and control the parametrization: introducing large values for the weight β_jk joins originally separate layers into one, based on the results of first tentative adjustments and on the knowledge about the vertical stratification of the mantle, whereas adopting β_jk=0 decouples the associated unknowns a_j and a_k in the adjustment.

Using the seismic tomography spherical harmonic model, the horizontal velocity distribution was computed for depths down to 2850km with a step size of 50km, each field representing the state within a 50-km thick layer bounded by spherical shells. The top of the uppermost layer is undulating and defined by the Moho discontinuity, the same surface as used for the determination of the residual gravity field.

As mentioned above, the dynamic effect of the inner inhomogeneities reflected in dynamic topography at the Earth's surface need not be taken into account, but the gravitational effect of dynamic deformations at the core–mantle boundary (CMB), or at any internal boundary in the case of a layered convection, is still present in the residual gravity. Any internal boundary perturbation should, unlike CMB deformations, already be reflected in the seismic tomography model, and therefore should enter into the results when adjusting for the empirical velocity-to-density scaling factors. The magnitude of the CMB dynamic topography impact on surface gravity was estimated using the kernel technique with a computer program that was kindly provided by Yanick Ricard (Ricard et al. 1993). It was modified to implement the latest innovations, for example the effects of a compressible mantle and of gravity changes with depth (Corrieu et al. 1995; Panasyuk et al. 1996). For density variations located in the upper mantle and in the transition zone, the induced dynamic topography at the CMB leads to variations in gravity at the surface of not more than 3mgal for any mantle viscosity model and can be neglected. For density variations located in the lower mantle, the induced dynamic CMB effects may reach amplitudes of 18mgal, but only at the lowest harmonics of degree 1 and 2 for which viscosity differences throughout the Earth are less important. Four recent viscosity models, discussed in Corrieu et al.(1994), Forte et al.(1995), Milne et al.(1998) and Pari & Peltier (1998), were applied and the scaling coefficients estimated taking into account dynamic CMB gravity effects. It is found that the scattering of the coefficients obtained for the different viscosity profiles is significantly smaller than the estimated coefficients' standard deviations even for the lower mantle layers. Therefore there is no need to deal with viscosity profiles for the purpose of this study.

3.2 Velocity-to-density scaling for the oceanic upper mantle

For a first tentative inversion, the upper mantle down to 375km was divided into seven layers of 50km thickness each, the highest possible resolution with the S20 model of Ekström & Dziewonski (1998). The mantle transition zone was represented by one layer down to 775km, and the lower mantle by two layers separated at an artificial boundary at a depth of 1625km. The adjusted scaling coefficients and their standard deviations for each layer are shown in Fig.7. A stable solution could only be obtained after introducing non-zero values for β_jk in eq. (3), which are equal for all layer combinations jk. The largest velocity-to-density scaling factor was found for the uppermost layer; that is, the subcrustal layer extending down to 75km. For the depth interval 75–225km, the three velocity-to-density scaling coefficients are relatively small, while for the depth range 225–375km the amplitudes become almost equal to the amplitude in the uppermost mantle. The total contribution of the transition zone and the lower mantle to the residual ‘crust-free’ gravity obtained after the inversion is not more than 15 per cent, independent of the model parametrization, and therefore a higher resolution of the middle and lower mantle was suppressed in this first step.

Figure 7

The scaling coefficients d ln ρ/d ln V_S and their standard deviations for the suboceanic mantle obtained from a tentative inversion of the residual gravity field and seismic model of Ekström & Dziewonski (1998) with maximum vertical resolution in the upper mantle (50km in accordance with the resolution of the seismic model). Small error bars show the error contribution from crust reduction (Fig.3). The transition zone and the lower mantle are represented by three layers: 375–775km, 775–1625km and 1625–2875km depth. The inversion is performed for the fields derived from the spherical harmonic coefficients up to degree and order 20.

It is important to understand what is the main source of the standard deviations of the estimated scaling coefficients. Since uncertainties of the mantle velocity data are not specified, a simple test with the ‘crust-free’ gravity was performed. Random noise values were generated at each data point of the 10°×10° grid with an amplitude according to the map in Fig.3. These noise values were added to the ‘crust-free’ gravity, and the scaling coefficients were computed in the same way as with the initial field. After several repetitions of these calculations (about 100), the obtained rms of the scattering in the scaling factors became stable, revealing the error due to the uncertainties in the crustal correction. In Fig.7 these rms values are shown together with the error bars estimated from eq. (7). It is evident that the error contributions arising from uncertainties in the crustal model are much smaller than the overall standard deviations after adjustment. This means that the seismic tomography data or simplifications in the linear model (eq.2) relating gravity and seismic velocities remain up to now the main error source in this kind of study.

To increase the accuracy and reliability of the inversion for the upper mantle and to avoid the necessity of stabilization, which might bias the solution in view of the significant variations of the scaling factors with depth, several of the initially 50km thick layers were joined to form thicker ones, taking into account the main features of the mantle structure and the result of the first tentative inversion. The mantle division and the velocity statistics underlying the following adjustments are summarized in Table1. The subcrustal layer down to 75km produces the most significant effect and is therefore kept isolated. The second upper mantle layer reaches down to 225km and is characterized by abnormally low velocities under oceans, in particular under the Pacific (Figs1 and 5), and also by small scaling coefficients (Fig.7). The third main layer includes the lower part of the upper mantle down to 375km. The scaling factors of the first solution differ only insignificantly within the latter two layers, so the concatenation is justified. The division within the transition zone and lower mantle is kept the same at this stage of the inversion.

Table 1

Division of the mantle into layers, and S model seismic velocity statistics.

The scaling coefficients and their errors as obtained for the most general (six layers) vertical mantle structure are shown in Fig.8. The adjustment is performed for the fields derived from the spherical harmonics up to degree and order 20, including and excluding the degree 1 and 2 terms, respectively. The gravity signal contribution from the deep layers is concentrated at the lowest harmonics. Therefore, in order to reduce the impacts of the CMB, the lower mantle and the transition zone in the inversion for the upper mantle, the degree 1 and 2 terms were filtered out a priori. However, filtering out the very long-wavelength spectral content leads, if unconstrained, to unacceptably large standard deviations for the scaling coefficients of deeper layers, and each, even minor regularization, results in very low coefficients for these layers compared with an inversion including the degree 1 and 2 terms.

Figure 8

The scaling coefficients d ln ρ/d ln V_S and their standard deviations for the suboceanic mantle obtained from the inversion of the residual gravity field and the seismic model of Ekström & Dziewonski (1998) for the most general (six layers, cf.Table1) vertical mantle structure. Small error bars show the error contribution from crust reduction (Fig.3). The inversion is performed for the fields derived from the spherical harmonic coefficients up to degree and order 20, including and excluding the degree 1 and 2 terms, respectively. For purposes of illustration, the horizontal axis in the graph is truncated at 1800km depth, but the inversion was performed for the whole mantle.

Therefore, the inversion was split into two steps: first, an inversion to determine the upper mantle scaling factors without taking into account the degree 1 and 2 terms, and introducing a stabilization only for the transition zone and lower mantle, leaving the coefficients for the upper mantle independent; and second, an inversion of the complete field for the remaining deeper layers while fixing the scaling coefficients for the upper mantle to the values obtained from the first stage. By this separation of the inversion process into two steps, the determination of the scaling factors for the upper and lower mantle is largely decoupled.

It turns out that the strength of stabilization applied to the transition zone and lower mantle in the first step of the solution (Fig.8) affects only marginally the solution for the upper mantle, where all β_jk in eq. (3) are equal to zero for j=1 to 3. The relative differences in the scaling coefficients for the upper mantle layers as resulting from the step 1 adjustment can be considered to be significant if compared with their standard deviations. Once again, the standard deviations far exceed the error contributions from the uncertainties in the crustal correction. It can be seen from Fig.8 that the scaling coefficients for the subcrustal layer and for the third layer (225–375km) are significantly larger than the coefficient for the second layer (75–225km). The values obtained for c_l and c₃ (0.31 and 0.23, respectively) in the adjustment with the degree 1 and 2 terms excluded are close to the ones derived from standard mineral physics equations and experimental studies (e.g. Sobolev et al. 1996; Karato & Jung 1998), but the velocity-to-density scaling factor obtained for the second layer has a value of only 0.05.

3.3 ‘Normal’ versus ‘anomalous’ ocean

‘Normal’ ocean may be defined as that part of the world's ocean where, in general, the lithosphere fits to the cooling oceanic lithosphere model; in other words, the bathymetry can be described as a function of age, and the thickness of the crust is close to a constant value. The rest of the ocean is affected by processes others than normal ocean spreading, mostly by recent or former mantle plumes. The origin of the plumes is still under debate. In order to focus on this open problem, the velocity-to-density relationships for ‘normal’ and ‘anomalous’ oceanic regions are estimated separately and then compared.

In addition to the classification provided by Mooney et al.(1998), the residual topography obtained in Kaban et al.(1999) is used as a criterion to discriminate between the ‘normal’ and ‘anomalous’ parts of the oceanic lithosphere. The bathymetry resulting from the cooling oceanic lithosphere model has been computed using the formulation of Hager (1983) and was removed from the observed field. Thus, choosing those areas of the world's ocean where the residual topography is small in amplitude means that areas affected by dynamics other than ocean spreading are largely excluded. Of course, the real oceanic lithosphere could differ from the schematic boundary layer model, but it is important that the parameters of this model fit statistically the distribution of ocean bathymetry over the mid-ocean ridges. The area of ‘normal’ ocean where the amplitude of residual topography does not exceed the adopted value of 400m is depicted in Fig.9: almost the whole of the western Pacific and the North Atlantic are considered to be ‘anomalous’.

Figure 9

Geographical distribution of ‘normal’ ocean (shaded area); that is, the area where observed bathymetry agrees with the bathymetry predicted by the cooling oceanic lithosphere model within 400m. Shading indicates the age of the ocean floor in Ma.

The velocity-to-density scaling coefficients obtained in separate adjustments for the ‘normal’ and ‘anomalous’ parts of the world's ocean are shown in Fig.10. For the subcrustal layer down to 75km, the scaling factor under the ‘anomalous’ ocean is somewhat smaller than that under the ‘normal’ ocean, but this statement is not definitive because the error intervals of the two coefficients overlap. For the layer between 75 and 225km the scaling factors are nearly the same for ‘normal’ and ‘anomalous’ ocean, and agree with the solution obtained when estimating only one parameter for the whole ocean. The low scaling coefficient and the large standard deviation obtained under ‘normal’ ocean for the layer below 225km is due to an almost missing signal in both gravity and velocity.

Figure 10

The scaling coefficients d ln ρ/d ln V_S adjusted for ‘normal’ and ‘anomalous’ oceanic areas, respectively, compared with the result for the whole oceanic area (same as in Fig.8), excluding degree 1 and 2 terms in the inversion. Only the scaling coefficients and standard deviations for the upper mantle are shown, but the inversion was performed for the whole mantle.

3.4 Non-linear inversion for the 75ȁ225Ȁkm layer

Because of the extremely low velocity-to-density scaling factor obtained for the oceanic upper mantle in the 75–225km deep layer, we investigate whether the velocity-to-density relationship there is essentially non-linear, that is varying with the amplitude of the velocity perturbation. This variation may be due to changes in the derivative of velocity with respect to temperature when the temperature is close to the melting point. In this case, the velocity-to-density scaling factor decreases with increasing amplitude of (negative) velocity perturbation. On the other hand, a low value of the scaling factor may be due to a superposition of thermal and lateral compositional changes that affect the velocity and density in different ways.

We do not intend to perform a comprehensive analysis of this problem because of the many uncertain parameters controlling a real velocity-to-density relation. Solving for too many parameters would lead to an unstable solution when attempting an inversion with the presently available data sets. Instead, the general type of a possibly non-linear dependence will be determined in an experimental approach: a combination of a linear and a non-linear relation described by a third-order polynomial is introduced (Fig.11). The sum of the linear and non-linear parts is then used for the velocity-to-density conversion within the 75–225km deep layer. Two coefficients, one to scale the linear and one to scale the non-linear part, will be determined as unknowns in the gravity–velocity inversion, such that one can discriminate between the two extreme cases, namely a large density/velocity ratio for low velocities, and small ratio for high velocities and vice versa, which can both be approximated by the composite function with appropriate scaling. Compared with the previous inversion for the whole mantle, the parametrization of the other layers is unchanged. The inversion, as described in Section 3.2 and represented in Fig.8 (thick line), is repeated, but this time also adjusting for the additional factor to scale the non-linear function. The third-order polynomial in Fig.11 is centred at +0.6 per cent on the velocity axis, which is the average of the maximum and minimum perturbations within this layer. The scaling coefficients obtained in the inversion for the remaining layers of the oceanic mantle fit well within their standard deviations to the ones obtained in the previous inversion with linear density–velocity relations throughout the mantle.

Figure 11

Non-linear velocity-to-density relation composed of a linear part (dot-dashed curve) and a third-order polynomial (dashed curve), both with a scaling factor to be adjusted. The solid line shows the experimental d ln ρ/d ln V_S-curve from the combination of the a posteriori scaled linear and non-linear functions for the 75–225km deep layer determined for the whole oceanic area. The thin lines indicate the error band of the relation computed from the standard deviations of the scaling coefficients. The other mantle layers are treated as in the inversion represented by Fig.8.

Thed ln ρ/d ln V_S curve and its error band as resulting from the adjustment are shown in Fig.11. Most of the error is due to the non-linear scaling coefficient, which has an error of 50 per cent. All attempts to determine the non-linear part separately for the ‘normal’ and ‘anomalous’ parts of the ocean failed because of an increase in the relative error of the non-linear coefficient due to the decrease of the size of the area under inversion. From the experimentally determined shape of the density–velocity relation (Fig.11) one can deduce that low-amplitude velocity anomalies (−1 to +1.5 per cent) transform into close-to-zero density changes, while for negative velocity anomalies of larger amplitude the density–velocity relation is almost linear with a ratio of d ln ρ/d ln V_S equal to 0.10–0.13.

3.5 Transition zone and lower mantle

After fixing the upper mantle scaling coefficients obtained in the first-step inversion (Fig.8), a second-step inversion is performed to obtain the solution for the transition zone and lower mantle. Since the input of the transition zone and lower mantle into the ‘crust-free’ gravity is relatively low, a stable solution can only be obtained by applying some damping (cf.eq.3). The strength of the damping depends on the number of unknowns (layers): when increasing the number of layers, one has to increase the β factors in eq. (3). Two cases are investigated: only a few layers with weak damping, yielding scaling coefficients averaged over a wide range; and many layers with strong damping, giving an increased resolution with depth but also an increased correlation between the solved-for parameters. The transition zone is, after previously having treated it as one layer (Fig.8), divided into two sublayers according to Table1, the upper one containing the 410km and the lower one the 670km discontinuity. The lower mantle was investigated in a two- and four-layer stratification, respectively, (see Table1). In addition, we perform an inversion for a homogeneous division of the mantle into 100km thick layers down from 400km depth. In the last case the damping coefficients β are 3.5 times larger than in the previous case. The inversions are performed for the whole ocean taking into account, contrary to the first-step inversion, degree 1 and 2 terms of the spherical harmonic representation of the gravitational and seismic fields.

The results of the inversion for the transition zone and lower mantle are shown in Fig.12. All the solutions demonstrate principally the same features. For the layers affected by the 410 and 670km discontinuities, the experimental scaling factors obtained are very small, presumably because they reflect a composite of distinct sources. For the lower mantle the scaling coefficient values vary remarkably: from a velocity-to-density scaling factor close to zero at depths between about 900 and 1000km, to a maximum value at depths between 1800 and 2300km.

Figure 12

The scaling coefficients d ln ρ/d ln V_S obtained for the lower mantle and the transition zone under oceans. The inversion is performed including all terms from degree 1 to 20. The scaling coefficients for the three upper mantle layers were determined in the previous step (Fig.8) and fixed. Both solutions are damped, but for the 25-layer solution the damping is approximately 3.5 times stronger than that for the six layers. Small error bars show the error contribution from crust reduction (Fig.3). For the 25-layer solution all error bars are approximately 1.5 times larger than for the six-layer solution.

The results for the transition zone and the lower mantle are considered as preliminary because the signal amplitude of about 30–40mgal of these layers in the residual gravity field hardly exceeds the noise even at long wavelengths. It is only possible to compare the relative trends of the obtained scaling factors; the absolute values depend strongly on the applied damping.

4 Discussion and conclusions

4.1 Experimental scaling and mantle structure

There exists a reliable correlation between shear-wave velocities (isotropic part) in the uppermost mantle and the residual ‘crust-free’ gravity field under oceans. This correlation is remarkably smaller or even missed in the subcrustal layer under continents. Supposing that lateral velocity variations in a specific layer are mainly due to temperature heterogeneities, while upper mantle density (and gravity) is sensitive to changes in temperature as well as in material composition, one can conclude that the lateral density changes in the oceanic uppermost mantle are mainly due to temperature heterogeneities while under continents temperature- and composition-induced density anomalies are superimposed and to a large extent compensate each other. The same conclusion was reached by Jordan (1978) based on mineral physics analysis.

It was found that the density/velocity ratio varies significantly with depth for the upper mantle under oceans. The value d ln ρ/d ln V_S=0.31 in the subcrustal layer down to about 75km is close to the one predicted for dry peridotites at appropriate conditions and assuming that the velocity variations are due to temperature changes (e.g. Sobolev et al. 1996). The same is true for the lower part of the upper mantle between about 225 and 375km depth: the value d ln ρ/d ln V_S=0.23 agrees well with the one predicted by Karato (1993). However, the density/velocity ratio for the layer from 75 to 225km depth is significantly lower than that in the overlying subcrustal layer, and than could be expected from ‘normal’ mineral physics equations. This result can be explained if one takes into account that, in the vicinity of the melting point, the V_S velocities vary much faster than density, which only depends on the thermal expansion coefficient. Moreover, the decrease in seismic velocity by about −6.5 per cent in the oceanic asthenosphere does not necessarily require the existence of partial melt, but can be explained if the temperature exceeds 90–95 per cent of the homologous one (Sato et al. 1989; Murase & Fukuyama 1980). Approaching the melting point at larger depths is only realistic if the homologous temperature at 75–225km depth is lower than in the subcrustal layer due to a different chemical composition at these depths.

Karato & Jung (1998) and Gaherty et al.(1999) give a perfect explanation of this phenomenon by introducing the G discontinuity at about 60–80km depth under oceans (Fig.13). This discontinuity can be attributed to a change in water content due to partial melting: the presence of water significantly reduces seismic wave velocities through anelastic relaxation (Karato & Jung 1998). The same idea was proposed by Hirth & Kohlstedt (1996), who suggested that the base of an oceanic plate is defined by a compositional rather than thermal boundary layer, or, at least, that the base is strongly influenced by a compositional boundary. Hirth & Kohlstedt (1996) estimated that the damp peridotite will be two orders of magnitude less viscous than the overlying dry harzburgite owing to the water-enhanced defect mobility in olivine. A discontinuity in mantle viscosity can develop at a depth of about 60–80km as a result of the extraction of water from olivine during the mid-ocean-ridge basalt melting process. The results of this study support these ideas, although the determination of the exact position of the G discontinuity under oceans is not possible because the vertical resolution of the used data sets is not high enough.

Figure 13

Schematic cross-section of oceanic upper mantle showing the compositional layering and the G discontinuity and the experimental scaling for the upper mantle. Gaherty et al.(1999) argued that decompression, melting and extraction of basaltic magma occurs in a narrow melt separation zone beneath the ridge crest. Any volatiles will enter the melt phase, resulting in a dry layer of depleted peridotite residuum overlying water-undersaturated (damp) mantle. This compositional boundary is preserved as the plate ages. The position of the 1200°C isotherm is estimated according to Hager (1983)

The suboceanic mantle density structure inferred from the seismic observations and from the experimental velocity-to-density scaling also agrees with the conventional idea (e.g. Hager 1983; Cazenave et al. 1986) that most of the temperature changes occur above 100km, while the large velocity low at depths down to 225km is generated by an only moderate temperature increase. At greater depths below about 225km, the pressure (and consequently homologous temperature) is high enough even for ‘wet’ peridotite, so that the velocity-to-density scaling factor turns out to be almost the same as for the subcrustal layer (Figs8 and 10).

The results, derived without any assumptions on isostasy, confirm a high level of isostatic compensation of the ocean-floor topography and Moho variations by subcrustal density anomalies associated mainly with temperature variations under oceans. The classical ‘plate’ cooling lithosphere model, which is a completely isostatic one, agrees well with the obtained results, at least within the accuracy provided by the seismic tomography model. Consequently, due to compensation, there exists only a very reduced sensitivity of long-wavelength observed surface gravity (i.e. free-air gravity disturbances prior to crustal reduction) to the obtained scaling coefficients for the upper mantle. For the deeper layers, starting with the transition zone, however, the input from observed free-air gravity disturbances to constrain the determination of velocity-to-density scaling becomes important.

The non-linear velocity-to-density scaling functions obtained for the 75–225km deep layer (Fig.11) reveal that a ratio d ln ρ/d ln V_S equal to 0.10–0.13 applies for dV_S/V_S<−1.5 per cent, and a ratio close to zero for smaller variations of V_S. The velocity interval −1.5 to + 1.5 per cent at depths from 75 to 225km under oceans in the Ekström & Dziewonski (1998) model largely coincides with the ‘anomalous’ parts of the oceans, which deviate significantly from the ideal cooling lithosphere model. Thus, the near-zero scaling coefficients could be evidence for changes in mantle composition in these parts of the ocean, possibly due to changes in mineral water content, affecting seismic velocities but not density. Another possible explanation could be that stronger chemical composition changes like those under continents occur, but the effect on density (and by this gravity) is compensated by thermal variations reflected in the velocity signature. For the other parts of the suboceanic mantle in the 75–225 deep layer characterized by high-amplitude negative velocity anomalies, a stable near-linear density–velocity relation is obtained but the gradient is 2–3 times smaller than the one obtained for the subcrustal layer above 75km. This amount of decrease in scaling is more reasonable when discussing the impact of a vertical compositional stratification of the oceanic mantle according to the G discontinuity hypothesis, than the extremely low velocity-to-density scaling factor of 0.05 obtained for this layer when adjusting only one linear coefficient.

From the results obtained for the scaling coefficients under ‘anomalous’ ocean one can conclude that the low in the residual gravity over the southwestern Pacific originates from the lower part of the upper mantle below 225km depth. This conclusion agrees with the southwestern Pacific superswell model suggested by McNutt (1998), which predicts for the subcrustal layer only very small temperature differences and thus unresolvable anomalies in both heat flow and effective elastic plate thickness.

The results obtained for the transition zone and for the lower mantle are less reliable because their signals in residual gravity are about as large as the crustal reduction errors. Without considering the absolute values of the scaling factors obtained for these layers, relative minima are nevertheless found for the layers affected by the transition zone discontinuities and for the depth interval 900–1100km. The minima near the phase transition zones are probably due to a complex relation between temperature, composition and phase state of the material at these depths. For the lower mantle, Forte & Perry (2000) also found that the velocity-to-density scaling factor decreases in its upper part. One possible explanation could be a discontinuity at a depth of about 900–1000km (e.g. Wen & Anderson 1997a; Ritzwoller & Lavely 1995), which may cause a similar effect to in the transition zone, although presently available data do not allow a definitive conclusion.

4.2 Removal of the velocity-to-density scaled seismic model from the residual ‘crust-free’ gravity field

The estimated experimental velocity-to-density scaling factors will be used to convert velocity into density anomalies, which are then propagated by forward computation to yield the disturbing mantle gravity field for the whole Earth. The scaling factors were derived only for oceanic areas where the correlation between residual gravity disturbances and mantle seismic velocities was high, and it was assumed that this correlation is due to the fact that, under these areas, upper mantle temperature anomalies dominate the compositional ones. Therefore, when subtracting the disturbing mantle gravity field from the residual ‘crust-free’ field, one can in principle, apart from modelling and data errors, identify in the remaining gravity field the signals coming from significant lateral chemical composition changes, especially under the continents.

The main problem with this procedure arises from the abnormally low experimental scaling factor obtained for the oceanic mantle between 75 and 225km depth, which is considered to be specific for oceanic areas. If the reasoning is correct, that the low scaling factor results from a high mineral water content, then the relation between velocity perturbations and temperature-induced density anomalies within this layer could be completely different under continents: the layer of partially depleted mantle is probably much thicker under the continents than under the oceans, and the temperature under the major part of the continental area (especially under cratons) is much lower, not approaching the homologous temperature as under oceans. Hence the density/velocity gradient within this layer is probably much larger under continents than under oceans.

For the upper mantle under continents (above 250km), it was therefore decided to use the velocity-to-density conversion factor based on a mineral physics approach, d ln ρ/d ln V_S= 0.21–0.32 (Karato & Jung 1998; Sobolev et al. 1996), to reproduce the effect of temperature changes, which are responsible for most of the variations in V_S. As this approach is somewhat arbitrary and the non-linear scaling obtained for the suboceanic layer cannot simply be propagated to positive velocity variations under continents, the results obtained with this approach for the subcontinental mantle are considered to be preliminary, and only the largest remaining gravity perturbations with an amplitude of several hundred mgal are worth discussing.

The velocity variations of the S20 tomographic model (Ekström & Dziewonski 1998) are then transformed into density variations, and propagated by forward computation to gravity perturbations at the Earth's surface. This field is removed from the ‘crust-free’ residual gravity disturbances (Fig.2) to yield the crust- and (thermal) mantle-free gravity field as shown in Fig.14. This field varies from −350 to +150mgal. The mean value over the oceans was constrained to be equal to zero, and over oceans the remaining features turn out to be of short wavelength with a root mean square (rms) of 40mgal. The most pronounced remaining signals over oceans in Fig.14 are the gravity lows around New Zealand and Iceland and the gravity highs at the northwestern border of the Pacific.

Figure 14

Residual gravity (mgal) after the removal of mantle gravitational effects applying the experimentally determined velocity-to-density scaling coefficients d ln ρ/d ln V_S to the S20 mantle tomographic model of Ekström & Dziewonski (1998). The scaling coefficients result from a gravity–seismic inversion restricted to oceanic areas. Different velocity-to-density scalings for the middle part of the upper mantle (75–225km depth) were applied for the suboceanic and subcontinental mantle.

With regard to continental areas, huge gravity lows of long wavelengths appear for cratons. The amplitudes reach −350mgal, which is considered to exceed the modelling errors in the crust (Fig.3) and (thermal) mantle gravity reduction. This result fits Jordan's (1978) idea of a chemical boundary layer under cratons whose mantle component consists of a low-density peridotite depleted in its basaltic constituents. This compositional low in mantle density is to a large extent compensated by a high in density because of low temperatures. The present analysis reveals the density constituent that is due to the chemical composition.

To summarize the results of the experimental velocity-to-density scaling and subsequent gravity field reduction, the signal variances per degree of the spherical harmonic expansion of the gravity field residuals over oceanic areas for the various stages of the reduction process are given in Figs15 and 16. Before expanding the grid values into spherical harmonics, the field was multiplied by a smoothed (truncated after degree 20) ocean–continent function to obtain the spherical harmonic coefficients related to the area of inversion only. The reduction in power is obvious, probably because the upper mantle lateral density structure under oceans is to a large extent determined by the thermal state.

Figure 15

Signal degree variances of the spherical harmonic expansion of the gravity field residuals for the various stages of the reduction process. The gravity field residuals are filtered to represent oceanic areas only. The solid line denotes the power of the initial ‘crust-free’ gravity field, the thick dashed line represents the field after having reduced the gravitational effects of upper mantle density inhomogeneities, and the thin dashed line represents the final stage after removal of whole mantle effects.

Figure 16

Fraction of residual power versus initial power over degree of fields' truncation, expressed by the ratios of the accumulated signal degree variances (cf.Fig.15) of the gravity field residuals and the initial ‘crust-free’ gravity field (ocean areas only) from degree 2 onwards. The solid line is obtained after removal of the effects of the upper mantle, and the dashed line after removal of the effects of the whole mantle.

The rms of the gravity field residuals over oceanic areas is reduced from the initial residual ‘crust-free’ field's 76mgal to 49mgal, taking into account only the upper mantle layers, and further down to 40mgal when additionally reducing the gravitational impact of the upper/lower mantle transition zone and the lower mantle down to the CMB. The remaining power after crust and mantle reduction still exceeds the power of the observed surface gravity field by a factor of 2, even in the area of normal ocean where the crustal model errors are small. This means that the seismic signal contains a significant component that is not correlated with gravitational signatures. The investigation of the origin of this component is left for future studies.

The upper mantle reduction is significant for all terms throughout the spectrum almost up to degree 20, except for the degree 1 terms. The increase in power for the degree 1 terms when reducing the upper mantle may be an artefact introduced by the convolution with the ocean–continent function, or may result from the fact that the upper mantle scaling factors are derived from an inversion in which the lowest harmonics were filtered out beforehand. The decrease in residual power when reducing also the deeper layers is best visible in the terms up to degree 7 (Fig.15). The spectrum of the remaining signal turns out to be flat and may represent the noise level of this investigation combining all error sources. Fig.16 shows the ratios of the accumulated degree variances of the residuals versus the initial ‘crust-free’ gravity disturbances over oceans from degree 2 onwards. One can deduce from Fig.16 that about 85 per cent of the power in the longest-wave ‘crust-free’ gravity field (degrees 1–5) and 70 per cent of the overall power up to degree 20, which was the limit of this investigation, are explained by the inverted seismic velocity model.

Acknowledgments

This work was supported by GeoForschungsZentrum Potsdam (GFZ), in particular by Prof. Ch. Reigber. The authors are grateful to Dr S. Sobolev of GFZ for fruitful discussions and to anonymous reviewers for useful comments. The coloured plots were produced using the freely available and maintained plot software Generic Mapping Tools (GMT 3.3), see Wessel & Smith (1995).