## 1 Forty years of gravity observations

Since the first artificial satellite orbited around the Earth (1957), the tracking of satellite orbits has been used to constrain the lateral variations of the gravity field of the Earth [12]. In this last 40 years, our knowledge of the Earth's gravity field has improved significantly in both accuracy and resolution [20]. The accuracy of the lowest degrees of the field has been increased by the addition of new and better data (a larger number of satellite orbits, some dedicated, with variable inclinations). At the same time, the resolution benefited from adding the constraints of satellite altimetry over ocean [2] and from expending the geographical coverage and the quality of direct gravity measurements at the Earth's surface. Combined models have been published up to degree and order 360 and in this paper we will use the Egm96 model [13].

However, the gravity community of scientists had been longing for years to have a direct satellite dedicated to the gravity field until these last years when, at last, various satellites have been launched or planned (Champ2000 [21], Grace2002 [22] and soon Goce2006 [29]). Of course, a direct measurement of gravity is impossible in orbits, but the derivatives of the gravity can be measured by monitoring the distances between two or more probes orbiting around the Earth. These probes can be two active satellites (Grace), an active satellite and a passive orbiting target, or two probes inside a single satellite (accelerometric concept, Champ).

The knowledge of satellite-only gravity models is necessary for various reasons. First, the direct surface gravity data that controls most of the present-day gravity models at, say, degree larger than 50, is very geographically uneven and of variable quality. The missing data are sometimes estimated from the topography itself, which bias further study of gravity–topography relationships. Second, the altimetric data assumes that the ocean topography is an equipotential surface, although surface currents are driven by the dynamic topography of the sea surface (i.e. the topography in excess of the equipotential). Third, the quality of the gravity models is now such that the time variations of the gravity of astronomical (tides), meteorological (annual perturbations), hydrological (variations of water table levels), and secular sources (post glacial rebound) must be taken into account.

The combined data of Champ and Grace have recently allowed geodesists to describe the gravity field using satellite only data up to degree ∼160, although it is probably safer to limit the expansion around degrees 100–120. This is already a striking progress, as the previous satellite-only solutions were probably only accurate up to degree 20–30. This decade should see a huge further progress in gravity observations.

As usual in the geodynamic community, we do not refer the geoid to a best-fitting ellipsoid, but to the shape that the Earth should have if gravity and rotation were in equilibrium [18]. This non-hydrostatic geoid and the corresponding free-air gravity only differ at even degrees and order 0 (practically, only at degrees 2 and 4) from those used by geodesists. The degree-2 order-0 density induces a gravity anomaly seen in the non-hydrostatic geoid of the geodynamicists, but would apparently have no effect when the best-fitting ellipsoid of the geodesists is used.

Here, we want to discuss a few simple ideas about the source of the gravity field and the implications for the solid-Earth geophysics. We only consider the global gravity field and all studies with local gravity measurements or with oceanic altimetry are beyond the scope of this paper. The gravity field is depicted in Fig. 1 as the geoid surface (top) and the free-air gravity (bottom) anomalies. The relationships between these quantities will be given in Section 2. The geoid undulations emphasize the long-wavelength components of the free-air anomalies. An unexpected feature of the geoid is its lack of correlation with the usual pattern of plate tectonics or surface topography.

## 2 A few basic equations

The gravitational potential outside the volume including all mass sources is simply harmonic and can be written as:

$$V(r,\theta ,\varphi )=\sum _{l=0}^{\infty}\sum _{m=0}^{l}{\left(\frac{a}{r}\right)}^{l+1}{V}_{lm}(a){Y}_{lm}(\theta ,\varphi )$$ | (1) |

Two possible mappings of the gravity data can be done, either on the form of a geoid anomaly (see Fig. 1, top):

$$N(\theta ,\varphi )=\frac{V(a,\theta ,\varphi )}{{g}_{0}}$$ | (2) |

$$\delta {g}_{lm}=(l-1)\frac{{V}_{lm}(a)}{a}$$ | (3) |

The natural use of the spherical harmonic functions in the gravity representation suggests similarly to expand the topographies ${h}_{i}(\theta ,\varphi )$ (e.g., surface, Moho, core–mantle undulations...) and the internal density $\rho (r,\theta ,\varphi )$ on the same basis:

$$\rho (r,\theta ,\varphi )=\sum _{l,m}{\rho}_{lm}(r){Y}_{lm}(\theta ,\varphi )$$ | (4) |

$${h}_{i}(\theta ,\varphi )=\sum _{l,m}{h}_{lm}^{i}(r){Y}_{lm}(\theta ,\varphi )$$ | (5) |

The coefficients ${V}_{lm}(a)$ are related to the density variations by

$${V}_{lm}(a)=\frac{G}{(2l+1)a}\underset{E}{\int}\rho (r,\theta ,\varphi ){\left(\frac{r}{a}\right)}^{l}{Y}_{lm}(\theta ,\varphi )\phantom{\rule{0.2em}{0ex}}\mathrm{d}V$$ | (6) |

By introducing the spherical harmonic expansion (4) into (6) and taking into account the density jumps $\mathrm{\Delta}{\rho}_{i}$ (here assumed uniform, and corresponding to the difference between the densities above and below the discontinuity) associated to each interface i located around the mean radius ${a}_{i}$, one gets:

$${V}_{l}^{m}(a)\simeq \frac{4\pi G}{(2l+1){a}^{l+1}}(\phantom{\rule{0.2em}{0ex}}\underset{0}{\overset{a}{\int}}{r}^{l+2}{\rho}_{lm}(r)\phantom{\rule{0.2em}{0ex}}\mathrm{d}r-\sum _{i}\mathrm{\Delta}{\rho}_{i}{h}_{lm}^{i}{a}_{i}^{l+2})$$ | (7) |

An obvious but important implication of Eq. (7) is that the geoid undulations generated by the Earth topography alone are:

$${N}_{l}^{m}=\frac{3}{2l+1}\frac{{\rho}_{c}}{\overline{\rho}}{h}_{lm}$$ | (8) |

Compensation means that as soon as density anomalies are present, they induce stresses that deflect all the density interfaces. It is only the total mass anomalies (i.e. density anomalies and interface deflections) that can be safely used in gravity modeling. Deflections of interfaces are generally computed assuming perfect compensation under each point (Airy model) or using elastic models of the lithosphere for shallow mass anomalies (e.g., [31]) and using viscous models for deep heterogeneities related to mantle convection. We do not want to develop here a complete model of viscous dynamic compensation that can be found in [26] and [23] (see also [7] in Cartesian coordinates), and we only summarize the results of these papers.

The general idea is to consider a model of the Earth where the rheology and the reference density are known, then from this mechanical model to compute the interface deformations (surface, CMB...) imposed by the presence of internal density loads. Only mass distributions (internal masses plus interface deformations) that are solutions of the mechanical problem are physically acceptable and it is only within these solutions that the 3D density structure of the Earth can be found. In other words, the models of Moho topography, CMB topography and mantle 3D tomography taken separately cannot be trusted, unless when taken together they satisfy the requirement of mechanical equilibrium.

From the interface deflections and from the internal loads the resulting gravity can be easily computed following (7). The results of this modeling are as follows:

- • the presence of an internal mass of degree l induces an interface topography that reaches exponentially an equilibrium shape after a time constant of order $2\eta l/({\rho}_{m}{g}_{0}a)$. This is the constant of post-glacial rebound, a few thousands of years for degrees ∼10. This is very short compared to most geological timescales. Except over the zones previously covered by the last glaciation, the long-wavelength topography is in steady equilibrium with internal loads;
- • density anomalies close to the surface (between the radii $a-d$ and d) induce a topography that verifies the usual isostatic rule:

with $\epsilon (0)=0$, i.e. the departure from isostasy goes to zero with d. Under a stiff viscous lithosphere, $\epsilon (d/a)$ varies as ${(d/a)}^{2}$. Isostasy is even more closely verified when the elasticity of the lithosphere is also taken into account. In that case, the right-hand side of Eq. (9) varies as ${(d/a)}^{4}$ [30];$${\rho}_{c}{h}_{lm}+\underset{a-d}{\overset{a}{\int}}{\rho}_{lm}(r)\phantom{\rule{0.2em}{0ex}}\mathrm{d}r=\epsilon \left(\frac{d}{a}\right)$$ (9) - • in agreement with the previous point, masses located directly on interfaces (surface or Cmb) are exactly compensated and do not generate a gravity field at first order. This implies that we cannot understand the geoid by only considering compensated mass anomalies very close to the surface or very close to the Cmb;
- • the resulting gravity has generally a sign opposite to that of the deep-density anomaly; e.g., a topographic high is associated with a crustal root (internal deficit of density), but a geoid high. A positive correlation between internal densities and geoid can however be obtained when the deep masses are close to a significant viscosity increase. This point will be discussed in Section 4.

$${N}_{l}^{m}=\frac{3}{2}\frac{l+2}{2l+1}\frac{{\rho}_{c}}{\overline{\rho}}\frac{d}{a}{h}_{lm}=-\frac{3}{2}\frac{l+2}{2l+1}\frac{\delta {\rho}_{lm}}{\overline{\rho}}\frac{{d}^{2}}{a}$$ | (10) |

The fact that the gravity perturbations induced by the surface topography and the internal loads tend to cancel each other strengthens the need to use the exact gravity equation in order to model the short-wavelength components. For the short-wavelength anomalies, we therefore represent the mass anomalies by parallelepipedic cells of uniform density (a volume limited by two longitudes, two latitudes and two radii). To compute the crustal field, we use for example a grid of 1024 longitudes and 512 latitudes. The computation of Eq. (6) for parallelepipedic units can be very precisely performed up to high order by recurrence relations [11] (here we use ${l}_{\mathrm{max}}=250$). Notice however that although we compute exactly the gravity field from the crustal model, this model remains based on the isostatic rule that in principle is only valid to first order. For the lowest harmonic degrees that we will ascribe to the mantle and, in part, to the lithosphere, we can safely use the first-order expression (7).

## 3 Crustal and lithospheric sources of the gravity field

The most direct candidates for the sources of the gravity field are of course related to the density heterogeneities associated with the surface topography and the associated Moho undulation. However, we know that at least one major structure of the Earth's topography, namely the oceanic ridges, is not related to the crust, but to the cooling and contraction of the oceanic lithosphere (e.g., [30]).

Assuming that the topography is locally compensated (Airy compensation), we can estimate the Moho depth from the isostatic rule. We define the densities of ice, water, oceanic or continental crust by ${\rho}_{\mathrm{i}}$, ${\rho}_{\mathrm{w}}$, ${\rho}_{\mathrm{c}}$ and their thicknesses by ${h}_{\mathrm{i}}$, ${h}_{\mathrm{w}}$ and ${h}_{\mathrm{c}}$, respectively. These layers overlay a lithosphere of density ${\rho}_{L}$ with thickness L and a mantle of density ${\rho}_{\mathrm{m}}$. We define H as the external topography of the Earth, being zero over oceans and following the surface of the crust or the ice sheets on continents. Isostasy implies:

$$({\rho}_{\mathrm{m}}-{\rho}_{\mathrm{i}}){h}_{\mathrm{i}}+({\rho}_{\mathrm{m}}-{\rho}_{\mathrm{w}}){h}_{\mathrm{w}}+({\rho}_{\mathrm{m}}-{\rho}_{\mathrm{c}}){h}_{\mathrm{c}}+({\rho}_{\mathrm{m}}-{\rho}_{L})L-{\rho}_{\mathrm{m}}H=Ct{e}_{1}$$ | (11) |

To build an isostatic model of crustal thickness, ${h}_{\mathrm{c}}(\theta ,\varphi )$, we assume that the lithospheric thickness is only a function of the age of the sea floor, $L(\mathit{age})$. Therefore, using a compilation of sea-floor ages [19], we computed ${h}_{\mathit{age}}$, the excess topography of the sea floor with respect to its depth far from the ridge. Our empirical fit, ${h}_{\mathit{age}}$, reaches 2360 m over ridges, varies as expected with $\sqrt{\mathit{age}}$ at young lithospheric ages and flattens to zero near 120 Myr and under continents.

We assume that all this age-dependent excess sea-floor topography is related to the lithospheric thickening that we compute from:

$$-({\rho}_{\mathrm{m}}-{\rho}_{\mathrm{w}}){h}_{\mathit{age}}+({\rho}_{\mathrm{m}}-{\rho}_{L})L(\mathit{age})=Ct{e}_{2}$$ | (12) |

The resulting geoid and free-air gravity computed to degree and order 250 are depicted in Fig. 3 (using ${\rho}_{\mathrm{i}}=900\text{}\text{kg}\phantom{\rule{0.2em}{0ex}}{\text{m}}^{\mathrm{-3}}$, ${\rho}_{\mathrm{w}}=1000\text{}\text{kg}\phantom{\rule{0.2em}{0ex}}{\text{m}}^{\mathrm{-3}}$, ${\rho}_{\mathrm{c}}=2800\text{}\text{kg}\phantom{\rule{0.2em}{0ex}}{\text{m}}^{\mathrm{-3}}$, ${\rho}_{L}=3270\text{}\text{kg}\phantom{\rule{0.2em}{0ex}}{\text{m}}^{\mathrm{-3}}$ and ${\rho}_{\mathrm{m}}=3200\text{}\text{kg}\phantom{\rule{0.2em}{0ex}}{\text{m}}^{\mathrm{-3}}$). This degree of resolution should be attained within the next 10 years by satellite-only solutions. The predicted geoid has a rather low amplitude, about one third of that observed (compare with Fig. 1, top), and little resemblance to the observed pattern. The free-air gravity map, on the other hand, is highly correlated with the observations (compare with Fig. 1, bottom). Over trenches, our model that forces isostasy cannot explain the very local signal of non-isostatic origin.

We quantify the misfit between model and observation in Fig. 4 where we plot the degree correlations, the degree amplitude ratio between model and observation and the variance reduction of the observed gravity. The thick lines correspond to our total-crust and lithospheric model. The thin dotted lines only include the crustal model. As expected from Section 2, for the lowest degrees, the isostatically compensated crustal signal is much lower in amplitude (middle panel) than the observations and not even correlated (top panel). However, for degrees larger than ∼15, a very signification correlation (top panel) is observed that leads to a $\sim 45\mathrm{\%}$ variance reduction of the signal. Including the oceanic-lithosphere thickening improves moderately but uniformly the correlations and the variance reductions, particularly between degrees 8–40 (compare the thin and thick lines).

The elastic support of mass anomalies that exists at short wavelength [31] has not been considered in this modelling of gravity of crustal and lithospheric density anomalies. This does not seem to affect the fit of our synthetic gravity with observations, except maybe by a slight reduction of the variance for degrees 150–250 (Fig. 4, bottom panel). The presence of elasticity affects the compensating topography and the gravity by a term varying like ${(k{E}_{\mathrm{e}})}^{4}$ (where ${E}_{\mathrm{e}}$ is the elastic thickness of the lithosphere and k the wavenumber of the topography [30,31]). Replacing k by $(l+1/2)/a$, the good fit depicted on Fig. 4 implies that the average elastic thickness of the Earth is not larger than ∼30 km. Any attempt to explain the Earth's gravity for degrees larger than 200–250 will need to account for the lithospheric elasticity.

## 4 The long-wavelength mantle sources

The previous exercise clearly shows that the lowest wavelengths of the gravity field cannot be related to near-surface features: shallow anomalies would be compensated and would not generate much gravity anomaly and, at any rate, no obvious near-surface feature (geological structure, ocean–continent distribution...) correlates with the geoid. We must therefore look toward the 3D structure of the mantle to identify the large-scale density heterogeneities. Undulations of the CMB also affect the geoid, but are consequences of the general mechanical equilibrium. The liquid core cannot sustain large mass anomalies and the small inner core density heterogeneities, if they exist, are too far away to reasonably have a large influence on the surface gravity.

Starting in the mid-1980s (e.g., [6]), mantle tomography, i.e. the determination of the 3D structure of the mantle through seismic-wave propagation modeling has made enormous progress. Like the geoid observations, the quality of the models has improved both in precision and in resolution. To search for deep-mantle mass anomalies, we make use of the synthetic S model (‘Smean model’ of [1]), which is a weighted average of previously published models [9,17,28]. This model may not be a very good model for any of the dataset that have been used to construct it, but it should somewhat emphasize the common structures found by various tomographic approaches.

Although the relationship that relates seismic velocities and densities may not be so simple, we will consider these velocities as a proxy for densities, where faster velocities correspond to higher densities. The correlation between the observed gravity field and the Smean tomographic model is depicted in Fig. 5. This correlation is plotted as a function of degree (horizontal line) and depth (vertical axis). This correlation indicates very clearly that at long wavelength, the gravity is positively correlated with fast mantle zones from the transition zone to about ∼1000 km deep (red colors) and negatively correlated with fast mantle zones in most of the lower mantle and in the first ∼200 km (blue color).

The anti-correlation (blue color) is in fact in agreement with the isostatic rule: the dense lithosphere (fast velocities) should depress the surface topography and the resulting gravity anomalies should have the sign of the surface depression and be negative. It is in fact the depth range of positive correlations, across the upper–lower mantle interface that needs to be explained.

Since 1984 [23,26], an explanation for the change of sign of the correlation between mantle density and geoid has been proposed. In the case when the mantle viscosity increases rapidly with depth by a factor in the order of ∼30, a dense anomaly is supported by the stiffer underlying mantle, which results in a smaller surface deflection, larger CMB deflection and a positive correlation with the geoid. Away from a sharp viscosity increase, the induced surface deflection is larger, and imposes its sign to the gravity anomaly. The general theory also includes the dynamic topography induced by the flow at the CMB.

The geoid models based on this approach are well correlated with the observed geoid (e.g., [8,10]). Unfortunately, this approach does not constrain tightly the exact amount of viscosity increase that, according to various authors, ranges between one or two orders of magnitude and is localized across the transition zone, at 670-km depth, across the top of the lower mantle, with a discontinuous jump or a continuous increase. Rather than using a given tomographic model to derive the gravity contributions of the mantle heterogeneities, we can also use a somewhat different approach, maybe more speculative, but that brings the present structure of the mantle within the general paradigm of plate tectonics.

The general agreement on convection of fluids heated largely by internal radioactivity is that the regions of lithospheric downwellings (the subduction zones) should be underlain by cold descending plumes (the slabs). These structures should dominate the mantle heterogeneity structure. This view results in a testable consequence that the present-day mantle structure should, at first order, be the result of past subductions [27].

Starting from a compilation of plate reconstruction models spanning the last 200 Myr [16], we compute the position that subducted slabs could have in the mantle. This very simple model assumes that each piece of slab sinks vertically and that the slab excess density is conserved through the whole mantle [15,25]. The only free parameter of the model is the viscosity increase at 670 km depth that slows down the sinking velocity of slabs in the lower mantle.

The quality of this geodynamical model based on paleo-plate reconstructions lies in its robustness. Indeed there is only one parameter to adjust, the viscosity jump at 670 km depth. The age relationship for the slab density of the lithosphere at subduction is well known (e.g., [30]).

Using a very simple viscosity profile (a 100-km-thick lithosphere with a viscosity 13 times that of the upper mantle, and a lower mantle 30 times more viscous, the resulting density model is well correlated with the available tomographic models, and provides a remarkable fit to the lowest degrees of the Earth's gravity field [25]. We perform the computation up to degree 10 only as our slab model is too simplistic to be trusted at very high degrees. The computed geoid of Fig. 6 is strikingly similar to the observed one in Fig. 1. The correlations depicted in Fig. 7 (dotted lines) confirm that the fit of the degrees 2 and 3 is close to perfect, but the first 10 degrees are indeed highly correlated with observations. This confirms the general geodynamic findings that slabs are the major structures of the convective mantle and that the structure of the whole mantle records the last 200 Myr of plate tectonics [15]. Further complexities in the modelling can be added, related to the details of the rheology in the transition zone, to the lateral viscosity variations and to the nature of the lithospheric rheology, but the gain in fitting the observed long-wavelength gravity remains small.

## 5 Conclusion: present-day understanding of the static gravity field

When the contributions of our crust, lithosphere and mantle sources are added together, the resulting gravity field is correlated at more than 99% with the observed gravity (see Fig. 7, solid line). The residual geoid (Fig. 8), has an amplitude reduced by ∼3 with respect of the original geoid. The total variance reduction is however not so good in the range $l=6\text{\u2013}15$, where the two end-member models (deep mantle sources versus shallow crust and lithospheric sources) are both inaccurate.

Certainly the lithosphere under continents is not uniform and density heterogeneities are observed by surface wave tomography up to 200–250-km depth (e.g., [5]). However, the density anomaly due to the subcontinental lithosphere should be very low, as no major anomaly of the residual geoid seems to correlate to the ocean–continent distribution or even to the position of cratons [14].

The fact that the mantle contribution of the geoid is highly correlated with observations but does not lead to a major variance reduction except for the very first degrees could be improved. We select the viscosity profile in order to achieve the best total variance reduction, which is strongly dominated by the signal at degrees 2 and 3. Fitting the total gravity rather than the total geoid would certainly select a viscosity profile yielding a lower total variance reduction for the geoid, but a more uniform variance reduction over the first 10 degrees.

Although an isostatic density model provides a satisfactory modeling of the gravity field, it is clear that the processes of mechanical compensation only require isostasy approximatively. The existence of deviatoric stresses (viscous or elastic), even at long wavelength, forbids the existence of a strict Archimedean equilibrium. The fact that the new generation of geoid models will be obtained with accuracy and without the potential bias of using surface data will impose to re-discuss the small departure of the shallow masses from isostasy. In particular, it will be necessary to take into account the elastic support of the lithosphere to extend the gravity modeling toward higher degrees.