Summary
The origin of the flattening of the seafloor bathymetry and geoid at old ages is still a matter of debate. Is this departure from the square-root-of-age law induced by dynamic uplift of the lithosphere linked to the flow in the mantle or by cooling of the lithosphere different from that predicted for an infinite half-space? The amplitude of the dynamic uplift should strongly decrease at short wavelengths and should not be correlated with the seafloor age on a local scale. Therefore, in order to eliminate the ‘dynamic signal’, the cooling of the lithosphere is studied here by means of depth–age and geoid–age relationships computed from high-pass filtered bathymetric, geoid and age fields, or from local dtopography/d
and dgeoid/dage slopes. The geoid–age and depth–age relationships obtained here are in good agreement with those previously deduced from global non-filtered data. At all observation scales, a significant departure from the half-space solution is found. However, the data are not sufficiently accurate to discriminate between the Plate and CHABLIS (Constant Heat flow At the Base of the Lithospheric ISotherm) models. Our results support the existence of small-scale convection supplying heat at the base of the old oceanic lithosphere.
Introduction
The seafloor bathymetry is known to depart at old ages from the subsidence expected for a purely conductive half-space cooling of the oceanic lithosphere. Two main processes may explain this depth ‘anomaly’: (1) small-scale advection of heat or shear heating below the lithosphere may slow down the cooling of the lithosphere, or (2) large-scale convection currents and mass anomalies in the mantle may sustain a positive dynamic topography on old seafloor. It is essential for our understanding of mantle convection to be able to constrain these phenomena, that is, the intensity of the small-scale convection or shear heating at the base of the lithosphere and the amplitude of dynamic topography.
Cooling of the oceanic lithosphere
The flattening of topography with age observed on averaged subsidence profiles may be explained by many models in which the cooling of the lithosphere is hampered by shear heating (Schubert 1976) or advection of heat at the base of the lithosphere. Parametric models such as the Plate model (Parsons & Sclater 1977) or the constant basal heat flow model (CHABLIS model) (Doin & Fleitout 1996) assume a given onset time (at 70 or 0 Ma, respectively) and a given evolution with age of the heat flow advected or the shear heating at the base of the lithosphere. They may be considered as mathematical analogues of a wide range of physical models (e.g. Schubert 1976, Davaille & Jaupart 1994, Fleitout & Yuen 1984a, Buck 1985). The onset time and vigour of small-scale convection in the latter models strongly depend on the value of the viscosity beneath the lithosphere and on the activation energy in the rheological law (Yuen & Fleitout 1984; Davaille & Jaupart 1994; Doin 1997; Dumoulin 1999; Marquart 1999), both of which are expected to vary regionally with small changes in asthenospheric temperature and volatile content.
Some geophysicists consider that the lithosphere cools as though it were on top of an infinite conductive half-space (e.g. Heestand & Crough 1981; Davies 1988; Morgan & Smith 1992; Carlson & Johnson 1994; Davies & Pribac 1993; Hohertz & Carlson 1998). In this view, the Earth's mantle loses heat only by convection on the scale of plate tectonics, except in limited areas where localized plume uprisings occur. The departure of the seafloor topography and heat flow from the half-space predictions is considered to be either insignificant (that is, this departure is only observed on small areas of old seafloor heavily covered by hotspots and hotspot swells) or due to broad topography uplifts on old oceanic areas linked with dynamic overpressure in the mantle.
Seismologists have observed that the seismic lithosphere continues to cool after reaching an age of 80 Myr until up to 150 Myr (Woodward & Master 1991; Nataf & Ricard 1996; Trampert & Woodhouse 1996; Zhang & Lay 1999). This has often been used as an argument against the Plate model or in favour of the half-space cooling model. However, lithospheric cooling beyond 80 Myr is also in agreement with the CHABLIS model.
Retrieving the cooling signal
As argued in the next section, the analysis of the depth–age or geoid–age curve from high-pass filtered data gives an unbiased evaluation of the departure from the half-space model predictions, that is, it is almost free of the influence of dynamic topography, and of the artefact introduced by the numerous hotspots covering old oceanic seafloor. Previous attempts to eliminate the dynamic signal through a high-pass filtering process comprise studies of the geoid jump across transform faults (Sandwell & Schubert 1982; Cazenave 1983; Driscoll & Parsons 1988; Freedman & Parsons 1990) and the evaluation of the geoid–age curve on globally filtered data (Doin 1996). In the present paper, we retrieve the cooling signal from the topography and geoid data, on which we apply two types of high-pass filtering processes (Fleitout & Doin 1996). After a general presentation of the inversion methods and of the available data, the topography as a function of age is evaluated using a method similar to that used in Doin (1996). Another type of inversion based on local 
and 
slopes is then presented. Finally, the implications concerning the cooling of the oceanic lithosphere are discussed.
THE SOURCES OF DYNAMIC TOPOGRAPHY AND GEOID, AND THE SIGNAL LINKED TO SMALL-SCALE CONVECTION
The sources of dynamic topography and geoid
A ‘dynamic’ topography and geoid distort the determination of the signal linked to the cooling of the oceanic lithosphere. What we call the ‘dynamic’ topography or geoid is the topography or geoid not directly linked to lithospheric mass anomalies, that is, induced by deep mass anomalies and by plate motions. Plate motions and deep mass anomalies are known to be the main sources of the large-wavelength geoid anomalies as the lithospheric contribution to the large-wavelength geoid is only of the order of 10 m (Turcotte & McAdoo 1979; Doin 1996). The problem is that the computed dynamic topography and the observed large-wavelength geoid show a high correlation with the age of the ocean floor at large wavelengths. In particular, Colin & Fleitout (1990) observed that the degree 2 geoid (mostly dynamic) and the functions such as age or age squared were highly correlated. Therefore, any estimate of the dynamic topography is affected by the age correction applied to take into account the cooling of the lithosphere. Vice versa, the estimated topography or geoid induced by the lithospheric cooling with age is affected by the dynamic signal (Cazenave & Lago 1991). Colin & Fleitout (1990) attempted to separate dynamic and lithospheric effects using an inverse method made possible by the low correlation between age and geoid for harmonics larger than 2. This inversion yielded almost no dynamic topography and a large flattening of the topography at old ages. These results need to be confirmed using other types of analyses.
In order to retrieve the lithospheric cooling signal the ideal would be to subtract the modelled ‘dynamic’ topography and geoid from the observed signal. By adjusting a sufficient number of model parameters, a good fit to the observed large-wavelength geoid with both one-layer and two-layer dynamic models can be found (e.g. Hager 1985; Corrieu 1994). Sedimentation and transgression in India (Fleitout & Singh 1988) and in North America (Gurnis 1990) place bounds on temporal variations in dynamic topography over these continents. They appear to be small (less than 500 m). Most one-layer models that fit the geoid well predict a very strong dynamic topography (amplitude greater than 1000 m) (e.g. Hager 1985), unless ‘hidden’ masses linked to the phase transitions at the top of the lower mantle compensate lower mantle mass anomalies (Thoraval 1995; Le Stunff & Ricard 1997; Cadek & Fleitout 2000). The various models proposed imply widely different dynamic topographies and it is therefore impossible to predict the ‘dynamic’ geoid and topography with sufficient accuracy to retrieve the lithospheric signal by a simple subtraction of the ‘dynamic’ signal.
We argue here that, at short wavelength, the dynamic topography and geoid signals have a moderate amplitude and are probably uncorrelated with the age of the lithosphere. The amplitude of the topography and geoid due to deep masses strongly decreases with wavelength. For a mass anomaly at 800 km depth and a realistic viscosity distribution, the predicted signal is divided by at least 3 when the wavelength decreases from spherical harmonic degrees l=2 to l=10 and it becomes negligible for l=20 (Ricard 1984). Moreover, the deep mass anomalies might be related to the oceanic age at large wavelengths (former subducted slabs are more numerous below old oceans), but should not be linked to the short-wavelength structure of the moving plates. Plate motions also induce a large-wavelength signal correlated to the age of the oceanic lithosphere. In a model without lateral viscosity variations, the dynamic topography and geoid linked to plate motions simply relate to the horizontal divergence of surface plate velocities (∇HvS), that is, they show a high where surface plate velocities converge (subduction zones) and a low where they diverge (ridges) (Schubert 1978; Ricard 1984). Quantitatively for a homogeneous earth and a given degree l, the spectra of the topography and the geoid vary with ∇HvS/l and ∇HvS/l(2l+1), respectively, and both peak at degrees 4 and 5 (Ricard 1986). An important layering of the flow due to the increase in viscosity from the asthenosphere to the lower mantle steepens the amplitude decrease of both fields with l. The spectral decay of the expected dynamic topography thus depends somewhat upon the input parameters in mantle dynamics models. However, as an example, Cadek & Fleitout (2000) modelled the partially layered flow induced by mantle mass anomalies (derived from seismic tomography) and plate motions. Models that accurately fit the observed geoid yield a decrease in dynamic topography from 174–149 m at degree 2 to 39–26 m at degree 5 and to 7–11 m at degree 9, which correspond to 25, 8 and 3 per cent, respectively, of the amplitude of the oceanic topography at the same degree. Therefore, the topography and geoid linked to plate motions should have a limited amplitude and should not be correlated with age at short wavelengths.
Signal associated with small-scale convection
This paper uses the topography and geoid as functions of age in the wavelength range between 500 and 10 000 km to assess the presence of small-scale convection bringing heat to the base of the lithosphere. It is therefore important to discuss the nature of this ‘small-scale’ convection and its associated signal on topography and geoid. Furthermore, we will discuss how plumes rising from the upper boundary layer may perturb our analysis.
The structure of convection in the Earth's mantle is inferred from experimental or numerical simulations. Among the studies using a temperature-dependent viscosity fluid and aimed at understanding the thermal evolution of the lithosphere and its equilibrium thickness, there are numerical and laboratory experiments on a fluid cooled from above (Buck 1985; Davaille & Jaupart 1993), on a fluid heated from within (Grasset & Parmentier 1998) and on a bottom-heated convective box (Doin 1997; Dumoulin 1999). All have heat transfer at the base of the lithosphere controlled by ‘small-scale’ destabilization of the base of the lithosphere and proportional to the viscosity in the asthenosphere to the power 1/3. The role of plumes in the experiments with a bottom boundary layer is to increase locally the temperature in the asthenosphere, and hence to decrease the viscosity and to increase the heat transfer by ‘small-scale’ destabilization of the base of the lithosphere.
The wavelength of these instabilities is inferred to be of the order of 200 km (Dumoulin 1999). For a lateral temperature variation of 250 °C over a thickness of 50 km (which is an upper bound), the expected topographic signal amounts to some tens of metres and the geoid to a few tens of centimetres (Fleitout & Moriceau 1991). If the instabilities form random blobs, the signal would be undetectable. However, if they form downgoing sheets parallel to the shear direction below the plate, the geoid signal may be visible as parallel lines such as the lines observed in the Pacific ocean (e.g. Haxby & Weissel 1986) with a wavelength of about 200 km (Marquart 1999). However, Fleitout & Moriceau (1991) argued that the amplitude of the observed associated topography is too large to reflect small-scale instabilities at the base of the lithosphere. As ‘direct’ evidence for present-day small-scale convection is ambiguous, we will concentrate here on its time-integrated effect on the cooling of the lithosphere. This effect can be approximated by a simple heat transfer at the base of the lithosphere at a scale larger than the typical distance between instabilities. This approximation may not be valid very close to transform faults, but is reasonable at scales between 500 and 10 000 km, which we consider here.
Hotspots probably represent the main and perhaps the only localized rising currents in the Earth's mantle. They will induce geoid and topographic signals in three ways: (1) the deep part of the rising current will induce topography and geoid anomalies; (2) the cushion of hot material spreading below the lithosphere will contribute to the geoid and topographic swells; and (3) the small-scale convection enhanced by a low viscosity will increase the amount of heat brought up by small-scale convection and thin the lithosphere. The signal linked to the first two sources of mass anomalies is not locally correlated with age and will therefore not be detected in our study. The third process (increased heat transfer at the base of the lithosphere) is part of the process that we wish to study. Indeed, heat transfer is certainly not uniform over the whole oceanic lithosphere and the estimates that we obtain are averages of this heat transfer. In particular, our analysis is not in contradiction to the observation that the subsidence rate is largest in areas far away from hotspots (Carlson & Johnson 1994). At a short scale, heat transfer by small-scale convection may be related to the age of the lithosphere, in the case of a hotspot interacting with a transform fault (Dalloubeix & Fleitout 1989). There, the hot material carried by the plume may stay predominantly on the young side of the transform fault, enhancing the heat transfer and increasing the topography and geoid difference between the two sides of the fault. If this phenomenon is significant on a global scale, our analysis somewhat underestimates the amount of heat carried by small-scale convection.
Local versus average subsidence rates
Large local and regional variations in the seafloor subsidence along swaths perpendicular to the ridge are observed (e.g. Marty & Cazenave 1989; Hayes & Kane 1994). When they are interpreted as variations in the lithospheric cooling parameters, they correspond to large variations of either the heat input in the case of the CHABLIS model, or the mantle temperature in the case of the Plate and half-space models. However, ‘noise’ (for example, variations in oceanic crust thickness, or highs and lows of dynamic topography) can also produce strong variations in the subsidence rate (Dumoulin 2000). Furthermore, the meaning of individual bathymetric profiles along swaths in terms of a single lithospheric cooling history is unclear: areas along a swath may be more than 2000 km apart and thus may have experienced very different cooling histories. In our analysis we will instead interpret statistical averages for a given age range of strongly fluctuating subsidence rates.
GENERAL PRESENTATION OF THE INVERSIONS AND THE DATA TREATMENT
To retrieve the seafloor subsidence corresponding to the cooling of the oceanic lithosphere, we need to invert the relation
where FT(age) is the thermal subsidence predicted by one of the cooling models. The dynamic topography due to lower mantle masses, surface plate velocities and swells is part of the noise here. As argued above it should not correlate with age at medium and short wavelengths and then should not slant the determination of FT(age). The other important and annoying source of noise is the presence of seamounts and volcanic structures. This will be discussed later.
We also want to invert the relation
where FN describes the geoid–age function as predicted by one of the cooling models. Here the noise is mainly due to lower mantle masses and to mantle return flow linked with plate motion. At the largest wavelengths the noise is about 10 times larger than the amplitude expected for FN(age), and the adjusted coefficients are not reliable. However, the noise amplitude should strongly decrease with wavelength (in 1/l2). High-pass filtering of (2) is therefore necessary.
Simple parametrizations of the three lithospheric cooling models are presented in Appendix A. The functions FT and FN are simply expressed as the sum of the expected response due to the cooling of an infinite half-space and a flattening term, representing the departure from the half-space model due to heat carried to the base of the lithosphere.
Surface observables
The oceanic topography data is extracted either from the ETOPO-5 database on a 5′×5′ grid or from the data set built from satellite altimetry and ship depth soundings by Smith & Sandwell (1997) (referred in this paper as SMITH97). It is then sampled again on a coarser 1°×1° grid using the normal mode estimator from GMT (Wessel & Smith 1991) that computes the maximum likelihood value in each 1°×1° cell. This procedure reduces the influence of small-dimension topographic highs or lows, and is therefore more relevant than pure sampling on 1°×1° nodes or computing the average over each 1°×1° cell. This topographic field is then corrected for sedimentation as in Colin & Fleitout (1990). The sediment thickness data are available between 71°S and 71°N. Regions shallower than −1500 m are later excluded from the data set as probable continental slopes or anomalous volcanic structures. The global geoid data sets are OSU91A (Rapp 1991) and EGM96 (Lemoine 1998). The marine geoid data set has been computed between 52°S and 52°N on a 7′×7′ grid using the data from the SEASAT, repeat GEOSAT, ERS-1 and ERS-2 altimetric missions (Doin 1995). The mapping technique is described in Doin (1995) and enables one to compute directly mean sea-level heights. The marine geoid data set is also resampled on a coarser 1°×1° grid with the normal mode estimator from GMT. The age data set is taken from Mueller (1993).
INVERSIONS OF THE DEPTH–AGE AND GEOID–AGE CURVES ON GLOBAL DATA SETS
Principle
We follow here the inversion procedure described by Doin (1996). Each field appearing in eq. (1) or eq. (2) (or more exactly in eqs A1–A3 and A7–A9 in Appendix A) is first high-pass filtered by removing its first 0 to n−1 spherical harmonic degrees. The coefficients in eqs (A1)–(A3) and (A7)–(A9) are then estimated by multilinear least-squares regression. The misfit in this regression is evaluated for every 1°×1° cell where the topography (or geoid) and age field are defined and is weighted by the surface covered by each cell. Because the number of independent ‘measurements’ is unknown, errors on the coefficients cannot be evaluated.
Topography
The coefficients p0, p1, p2, b and c are computed using
or
or
where the subscript n−30 refers to the spherical harmonic degrees kept in the analysis.
Two treatments must in this case be applied to the topography data (ETOPO-5 or SMITH97) before the high-pass filtering process. (1) Seamounts and volcanic structures, more numerous over old seafloor than over young seafloor, seriously perturb the ‘averaged’ depth–age relationship. They must be removed prior to the high-pass filtering. We therefore follow the procedure of Lewitt & Sandwell (1996) and remove the influence of small-scale structures by applying to the topography the GMT modal filter with a diameter of 1000 km (Fig. 1) or 1500 km. Spherical harmonic degrees larger than 38 or 25, respectively, are strongly damped by this procedure. Numerous volcanic features, for example the Hawaiian chain and the Ninetyeast ridge, disappear from the filtered topography maps. In Fig. 2 the mean seafloor depth and standard deviation are plotted as a function of age. ‘Seamount-free’ topographic maps show deeper seafloor at old ages and smaller standard deviations than raw data. However, mid-oceanic ridges are also less pronounced. For homogeneity, the same modal filter is applied to the age, square root of age, and age squared data sets. (2) The high-pass filtering must be performed on data sets (i.e. topography, age, square root of age and age squared) known everywhere on the Earth's surface. Therefore, the topography and age fields are extrapolated to cover the whole surface of the Earth prior to the high-pass filtering. We assign to nodes with unknown field values the average of neighbouring defined values. This procedure avoids creating jumps at the borders of the undefined areas. These jumps would result in large undulations in high-pass filtered maps, contaminating the real data around undefined areas.
Topography of the seafloor (ETOPO-5) with small-scale volcanic structures such as seamounts removed by the low-pass modal filter of GMT. The diameter of the area for which the mode is evaluated is 1000 km. Contours are drawn every 500 m from −6000 to −2000 m.
Topography of the seafloor (ETOPO-5) with small-scale volcanic structures such as seamounts removed by the low-pass modal filter of GMT. The diameter of the area for which the mode is evaluated is 1000 km. Contours are drawn every 500 m from −6000 to −2000 m.
Mean and standard deviation of the topography filtered by the low-pass modal filter of GMT as a function of age. The legend indicates the diameter of the area for which the mode is evaluated.
Mean and standard deviation of the topography filtered by the low-pass modal filter of GMT as a function of age. The legend indicates the diameter of the area for which the mode is evaluated.
Geoid
Δ N0, ΔN1, ΔN2, β and γ are evaluated using OSU91A or EGM96 and for various threshold spherical harmonic degrees n between 6 and 13 by multilinear least-squares regression of
or
or
The ‘age’, ‘age3/2’, and ‘age2’ fields have been extrapolated prior to the high-pass filtering process as mentioned above. Note that the flattening coefficients b (eq. 4), c (eq. 5), β (eq. 7) and γ (eq. 8) are proportional to the heat flow carried to the base of the lithosphere (Appendix A).
RESULTS
Topography
The coefficients p0, b and p1, and c and p2 computed using eqs (3)–(5) are shown for various threshold harmonic degrees (n=0 to n=10) in Figs 3(a) and (b). The half-space lithospheric cooling model implies that the subsidence rate is constant with age. However, the computed value for this subsidence rate, p0 (about −220 m 
Results from the global inversion of FT(age) and FN(age) on high-pass filtered data sets. In panels (a) and (b) the topography flattening (coefficients b or c) is plotted as a function of p1 or p2 (eqs 4 and 5) for threshold spherical harmonic degrees 0≤n<11 (full symbols, annotated by n). p0 values (eq. 3) are also shown with b or c set to 0 (open symbols). A modal low-pass filter with a diameter of 1000 km (triangles) or 1500 km (squares) is applied to ETOPO-5 (black symbols) or SMITH97 (grey symbols) data. The dashed lines show the relationships defined by eqs (9) and (10). The point obtained by Colin & Fleitout (1990) on unfiltered data is also shown for comparison. The thick arrows represent the trend expected as n increases if the flattening of the depth–age curve is due to dynamic topography. In panels (c) and (d), the geoid flattening (coefficients β or γ) is shown as a function of ΔN1 or ΔN2 (eqs 7 and 8) for threshold spherical harmonic degrees 5<n<14 (full symbols, annotated by n). The values of ΔN0 (eq. 6) are also shown with the flattening coefficient set to 0 (open symbols). Black symbols: OSU91A; grey symbols: EGM96.
Results from the global inversion of FT(age) and FN(age) on high-pass filtered data sets. In panels (a) and (b) the topography flattening (coefficients b or c) is plotted as a function of p1 or p2 (eqs 4 and 5) for threshold spherical harmonic degrees 0≤n<11 (full symbols, annotated by n). p0 values (eq. 3) are also shown with b or c set to 0 (open symbols). A modal low-pass filter with a diameter of 1000 km (triangles) or 1500 km (squares) is applied to ETOPO-5 (black symbols) or SMITH97 (grey symbols) data. The dashed lines show the relationships defined by eqs (9) and (10). The point obtained by Colin & Fleitout (1990) on unfiltered data is also shown for comparison. The thick arrows represent the trend expected as n increases if the flattening of the depth–age curve is due to dynamic topography. In panels (c) and (d), the geoid flattening (coefficients β or γ) is shown as a function of ΔN1 or ΔN2 (eqs 7 and 8) for threshold spherical harmonic degrees 5<n<14 (full symbols, annotated by n). The values of ΔN0 (eq. 6) are also shown with the flattening coefficient set to 0 (open symbols). Black symbols: OSU91A; grey symbols: EGM96.
−1), is much slower than the subsidence rate (about −340 m 
−1) observed on young ocean floor. Furthermore, when a flattening term is included in FT(age) (b or c), it shows strong scatter, but is always non-zero (between 8 and 20 m Myr−1 for b and between 0.035 and 0.085 m Myr−2 for c). On average the estimates of b and p1 and of c and p2 are close to those obtained using unfiltered bathymetry (Colin & Fleitout 1990). (In the latter case, seamounts are iteratively removed from the data set by applying a cut-off on the residual topography.)
We find faster subsidence rates in the square root of age p1 or p2 for larger estimates of the flattening coefficient b or c (Figs 3a and b). This is not surprising since the errors in b and p1 (in c and p2, respectively) have a correlation coefficient of 0.96 (0.82, respectively) due to the strong correlation between 
n−30 and 
n −30 (
2n −30, respectively). The relation
or
respectively, holds when a curve in p1
+ bt or p2
+ ct2, respectively, is fitted by a curve in p0
for t between 0 and tmax. The same relations hold for the unevenly distributed age of the ocean floor if tmax is set to about 140 Myr. The straight lines b(p1) (eq. 9) and c(p2) (eq. 10) are shown in Figs 3(a) and (b) for a given p0=−215 m 
−1. All points (p1, b) and (p2, c) plotting along these lines define similar subsidence curves FT(age); noise fluctuations (introduced by volcanic plateaus, seamounts, the data treatment processes, etc.) superimposed on the same cooling signal will induce displacements mainly along these lines. Thus, the results obtained with both ETOPO-5 and SMITH97 data and both modal filters (Figs 3a and b) reflect similar cooling relationships.
On the other hand, if the oceanic lithosphere was cooling as a half-space, as n increases (that is, as dynamic topography is filtered out), the subsidence rate should speed up from −215 m Myr−1 to the subsidence rate observed at young ages (about −340 m 
−1) and the flattening coefficient should tend to zero: the points in Figs 3(a) and (b) should depart from lines (9) and (10), scatter along the trend shown by the thick arrow, and for n>8 plot around b=0 or c=0. Obviously, the computed values (b, p1) and (c, p2) do not follow this trend but fall next to the lines eqs (9) and (10). No trend can be detected as n increases. Therefore, the flattening of the topography versus age relationship is clear in unfiltered topography, in topography low-pass filtered by modal evaluation and in high-pass filtered topography, in disagreement with the (half-space+dynamic topography) model.
Geoid
The coefficients ΔN0, ΔN1 and β, and ΔN2 and γ defined in eqs (6)–(8) are shown in Figs 3(c) and (d) for 5<n<14. When included in the analysis, the flattening coefficients β and γ are always non-zero, between 1.0 and 1.8 cm Myr−3/2 for β and between 0.05 and 0.08 cm Myr−2 for γ. Therefore, even at wavelengths where the amplitude of the geoid due to lower mantle mass anomalies and plate motions is low, the flattening of the geoid–age curve is obvious.
INVERSION USING LOCAL TOPOGRAPHY OR GEOID SLOPES
Principle
The idea here is the same as in the studies of the geoid jump across transform faults (Sandwell & Schubert 1982; Cazenave 1983; Driscoll & Parsons 1988; Freedman & Parsons 1990), with a rougher data analysis but a global coverage. The local slopes of the topography with square root of age and of the geoid with age are studied as a function of age. Both are expected to be constant if the half-space cooling model is valid. A departure from the half-space constant value shows that the cooling of the lithosphere has been slowed down by heat advection at its base. We choose to compute the local slopes on 5°×5° to 17°×17° wide areas—slopes computed on smaller areas will be more affected by high-frequency noise, data accuracy limits and small-scale convection effects, whereas slopes computed on larger areas are less numerous and smooth data over areas with crustal age differences that may be greater than 35 Myr.
Topography analysis
Taking the derivatives of eq. (1) (or more precisely of eqs (A1)–(A3) in Appendix A) leads to
and
Geoid analysis
The altimetric geoid must be high-pass filtered prior to the computation of geoid slopes with age to reduce further the noise in the inversion of eq. (2). The same filter f is therefore applied to the age, age3/2 and age2 fields. Eqs (A7)–(A9) then become
and
The methods for computing the topography and geoid slopes and their associated standard deviation, and for evaluating the coefficients in eqs (11)–(16) are described in Appendix B.
Topography slopes
There are two main sources of noise in topography slopes. First, areas covered with seamounts lead to ill-defined topographic slopes in the square root of age with very large standard deviations. However, the noise introduced by these seamounts could be negative as well as positive and therefore should not systematically slant the computed FT(age) relationships that represent a data average. Second, the ETOPO-5 data set has some drawbacks (Smith 1993) that yield a low accuracy in the computation of local topography slopes. In particular, it was gridded by digitizing contours at most 200 m apart. Therefore, the error on each data point can be of the order of 100 m. However, topographic slopes computed using the SMITH97 data set do not on average differ much from those computed with ETOPO-5. Third, as southern oceans are not well covered by ship depth soundings (Smith 1993), only areas north of 30°S have been retained in the inversion.
Along a 500 km profile with a spreading velocity of 5 cm Myr−1 and an age of 45 Myr, the expected topography decrease with age, according to the half-space cooling model, would be 245 m. On older seafloor or for oceans with a fast spreading rate this topography decrease is less important. Therefore, the best-defined slopes should be close to ridges, in oceans with a slow spreading rate, and across fracture zones. In areas where the seafloor topography does not correlate well with age, we obtain a large error on the topographic slopes, and the weight of these slopes in the inversion is small (see Appendix B). It therefore reduces the influence of topography anomalies that are not related to the progressive thermal cooling of the lithosphere. As a result, only a small fraction of the seafloor has a significant weight in the inversion. In Fig. 4 the areas where the topographic slopes are well defined are displayed. They are mainly along the mid-Atlantic, Pacific and Indian ridges plus some Pacific fracture zones. Therefore, the computed average subsidence curves will mostly be representative of oceans of age <90 Myr.
Symbols represent 5°×5° areas where the topographic slope in square root of age p is best defined (standard error on the topographic slope lower than 250 m Myr−1/2). The greyscale is proportional to the value of the topographic slope in m Myr−1/2. Squares: p<−440; hexagons: −440<p<−360; circles: −360<p<−280; diamonds: −280<p<−200; upward triangles: −200<p<−120; downward triangles: −120<p<−40; stars: −40<p<40.
Symbols represent 5°×5° areas where the topographic slope in square root of age p is best defined (standard error on the topographic slope lower than 250 m Myr−1/2). The greyscale is proportional to the value of the topographic slope in m Myr−1/2. Squares: p<−440; hexagons: −440<p<−360; circles: −360<p<−280; diamonds: −280<p<−200; upward triangles: −200<p<−120; downward triangles: −120<p<−40; stars: −40<p<40.
The histograms of the computed topographic slopes in square root of age are shown in Fig. 5(a) for 10 Myr age bins and for 5°×5° wide areas. The mean and standard deviation of these histograms are shown in Fig. 5(b). Variability is large around the mean and may represent a different thermal history for each area, but it may also be linked to volcanic structures and variations in crustal thickness. The 95 per cent confidence interval of the mean (Fig. 5b) gives the expected accuracy of the computed mean representative of the averaged cooling behaviour of the lithosphere. The average subsidence rate is significantly faster at 5 Myr (−285±24 m 
(a) Histograms of the computed local topographic slopes for each 10 Myr age bin. The amplitude of the curve at a given age is proportional to the frequency of measurement of this value, weighted by the inverse of its squared standard error. Grey curves: SMITH97; black curves: ETOPO-5. (b) Average, standard deviation (dotted errors bars) and 95 per cent interval of confidence for the average (solid error bars) compared with the predictions of the CHABLIS (dashed curve) and Plate (dot-dashed curve) models. Grey symbols: SMITH97; black symbols: ETOPO-5.
(a) Histograms of the computed local topographic slopes for each 10 Myr age bin. The amplitude of the curve at a given age is proportional to the frequency of measurement of this value, weighted by the inverse of its squared standard error. Grey curves: SMITH97; black curves: ETOPO-5. (b) Average, standard deviation (dotted errors bars) and 95 per cent interval of confidence for the average (solid error bars) compared with the predictions of the CHABLIS (dashed curve) and Plate (dot-dashed curve) models. Grey symbols: SMITH97; black symbols: ETOPO-5.
−1) than at 105 Myr (36±101 m 
−1). Note that according to the lithospheric half-space cooling model the topographic slope should remain constant. This result remains consistent when larger areas (up to 17°×17°) are used for computing the slopes.
Results of inversions using local topographic slopes
The flattening coefficients of the subsidence curve, b and c, are plotted against the coefficients in square root of age, p1 and p2 (see eqs 11–13 and Figs 6a and b). The points with their error ellipses correspond to different widths of the regression areas, from 17° to 5°, and to two different ways of computing local topographic slopes. Their accuracy decreases (the error ellipse widens) with the number of computed topographic slopes, that is, with the width of the regression area. They show a large scatter. However, the flattening coefficient b or c is always significantly larger than zero (between 14 and 22 m Myr−1 for b and between 0.055 and 0.075 m Myr−2 for c). If the lithosphere were cooling as an infinite half-space, the flattening coefficient should vary around zero, because the effects of dynamic topography on local topographic slopes are negligible. The values of p1 and b and of p2 and c computed using local topographic slopes are on average close to those found using unfiltered topography data (Colin & Fleitout 1990) or those found by inversion of high-pass-filtered bathymetric maps. Therefore, dynamic topography cannot be responsible for the flattening of the seafloor subsidence at old ages.
Results of the inversions of depth–age curves using local topography slopes (from ETOPO-5) computed over 5°×5° to 17°×17° areas. The flattening coefficients b (a) and c (b) are represented as functions p1 (a) and p2 (b) with their 1σ error ellipses (see eqs 12 and 13). p0 (a and b) is also shown with open symbols and flattening coefficients set to zero. The inversions are performed using local topography slopes with age (squares) or with square root of age (circles). Filled diamonds: coefficients inverted using unfiltered topography by Colin & Fleitout (1990).
Results of the inversions of depth–age curves using local topography slopes (from ETOPO-5) computed over 5°×5° to 17°×17° areas. The flattening coefficients b (a) and c (b) are represented as functions p1 (a) and p2 (b) with their 1σ error ellipses (see eqs 12 and 13). p0 (a and b) is also shown with open symbols and flattening coefficients set to zero. The inversions are performed using local topography slopes with age (squares) or with square root of age (circles). Filled diamonds: coefficients inverted using unfiltered topography by Colin & Fleitout (1990).
The Plate model assumes that the lithosphere first cools as a half-space in young oceanic areas (this cooling is quantified by p2
) and then at 70 Myr stops cooling. Accordingly, the values of p2 computed here (between −240 and −420 m 
−1, Figs 3b and 6b) bracket the subsidence rate observed on young ocean floor (i.e. −320 m 
−1) (e.g. Parsons & Sclater 1977; Schroeder 1984). In the CHABLIS model, the lithosphere near the ridge cools slower than expected from a conductive half-space. This results in values more negative for p1 (between −320 and −550 m 
−1, Figs 3a and 6a) than for p2. These lower values for p1 are reasonable as the thermal expansion coefficient might be higher (i.e. 3.8×10−5 K−1) than commonly estimated (i.e. 3.2×10−5 K−1) (Doin & Fleitout 1996).
Computed geoid slopes
The geoid slopes with age df(geoid)/df(age) are computed using either the EGM96 or the DOIN97 data set and different filters f defined as follows: spherical harmonic degrees larger than m+4 are kept, degrees m+3 to m are progressively smoothed out, and degrees m−1 to 2 are set to zero, with (1) m=4, (2) m=6, (3) m=8, (4) m=10 and (5) m=12. Results are consistent between the various filters, except that the average computed geoid slope for age less than 15 Myr regularly increases from −14 to −8 cm Myr−1 as m increases. Thus, the smoothing effect of the filter artificially reduces the amplitude of the geoid slopes for age less than 15 Myr. In Fig. 7, the areas where the geoid slopes are well defined are displayed. They are mainly along the mid-Atlantic, South Pacific and North Indian ridges plus some Pacific fracture zones and some small oceanic basins. We show in Fig. 8(a) for 10 Myr age bins the histograms of the geoid slopes that are computed for the five filters and then added together. This tends to smooth the effects of chancy correlation between the large-wavelength geoid and the age field, and the effects of some artefacts introduced by the filter. The average value of the computed geoid slopes in age significantly increases with age from −11±2 cm Myr−1 at 5 Myr to 1±4 cm Myr−1 at 125 Myr. The mean, standard deviations and 95 per cent confidence interval of these histograms are shown in Fig. 8(b). Again, variability is large around the mean and may represent various thermal histories and also the effects of ‘noise’ (seamounts, low-degree dynamic geoid, etc.) and errors in the age field. The geoid slopes computed here can be compared to the geoid jump across fracture zones (Fig. 8c) calculated along single altimetric profiles and fitted with thermo-elastic models of the lithosphere.
Symbols represent 5°×5° areas where the geoid slope in age ΔN is best defined (standard error on the geoid slope lower than 7 cm Myr−1). The greyscale is proportional to the weighted average of the geoid slopes obtained for the five filters (in cm Myr−1). Squares: ΔN<−22; hexagons: −22<ΔN<−18; circles: −18<ΔN<−14; diamonds: −14<ΔN<−10; upward triangles: −10<ΔN<−6; downward triangles: −6<ΔN<−2; stars: −2<p<6.
Symbols represent 5°×5° areas where the geoid slope in age ΔN is best defined (standard error on the geoid slope lower than 7 cm Myr−1). The greyscale is proportional to the weighted average of the geoid slopes obtained for the five filters (in cm Myr−1). Squares: ΔN<−22; hexagons: −22<ΔN<−18; circles: −18<ΔN<−14; diamonds: −14<ΔN<−10; upward triangles: −10<ΔN<−6; downward triangles: −6<ΔN<−2; stars: −2<p<6.
(a) Histograms of the computed geoid slopes (for the five filters) with age for each 10 Myr age bin. The amplitude of the curve at a given age is proportional to the frequency of measurement of this value, weighted by the inverse of its squared standard error. Grey curves: EGM96; black curves: DOIN95. (b) Average (squares), standard deviation (dotted errors bars) and 95 per cent interval of confidence for the average (solid error bars) compared to the predictions of the CHABLIS (dashed curves) and Plate (dash-dotted curves) models. Grey symbols: EGM96; black symbols: DOIN95. (c) Geoid jump across transform faults divided by the age offset between the two sides of the faults versus age, compared with the predictions of the Plate model GDH1 (Stein & Stein 1992) and the CHABLIS model (Doin & Fleitout 1996).
(a) Histograms of the computed geoid slopes (for the five filters) with age for each 10 Myr age bin. The amplitude of the curve at a given age is proportional to the frequency of measurement of this value, weighted by the inverse of its squared standard error. Grey curves: EGM96; black curves: DOIN95. (b) Average (squares), standard deviation (dotted errors bars) and 95 per cent interval of confidence for the average (solid error bars) compared to the predictions of the CHABLIS (dashed curves) and Plate (dash-dotted curves) models. Grey symbols: EGM96; black symbols: DOIN95. (c) Geoid jump across transform faults divided by the age offset between the two sides of the faults versus age, compared with the predictions of the Plate model GDH1 (Stein & Stein 1992) and the CHABLIS model (Doin & Fleitout 1996).
Results of inversions using local geoid slopes
The coefficients of FT(age) adjusted using local geoid slopes are plotted in Figs 9(a) and (b) for two filters, (3) and (5). Again, the flattening coefficients β and γ (defined by eqs 11–13) are non-zero, between 0.85 and 1.5 cm Myr−3/2 for β, and between 0.035 and 0.065 cm Myr−2 for γ. The computed values of ΔN1 and β, and ΔN2 and γ are also in good agreement with the coefficients found using global filtered data sets and threshold degrees n greater than 5 this study; Doin 1996). As the flattening of the geoid–age curve is similar and significant at both observation scales, it cannot originate from spurious correlations between age and geoid anomalies due to lower mantle masses and plate motions. It can thus be attributed to the cooling of the lithosphere itself. The adjusted geoid–age curves follow the general trend of the predicted geoid–age curves for a plate or a constant basal heat flow model.
Results of the inversions of geoid–age curves using local geoid slopes computed over 5°×5° to 17°×17° areas. The flattening coefficients β (a) and γ (b) are represented as functions of ΔN1 (a) and ΔN2 (b) with their 1σ error ellipses (see eqs 15 and 16). ΔN0 are also shown with open symbols and flattening coefficients set to zero. Squares and circles correspond to the results of inversions using local geoid slopes with age using the high-pass filter 3 or 5. For comparison, the coefficients inverted on global data sets are also shown by filled diamonds, (a) for n>5, this study, and (b) after Doin (1996).
Results of the inversions of geoid–age curves using local geoid slopes computed over 5°×5° to 17°×17° areas. The flattening coefficients β (a) and γ (b) are represented as functions of ΔN1 (a) and ΔN2 (b) with their 1σ error ellipses (see eqs 15 and 16). ΔN0 are also shown with open symbols and flattening coefficients set to zero. Squares and circles correspond to the results of inversions using local geoid slopes with age using the high-pass filter 3 or 5. For comparison, the coefficients inverted on global data sets are also shown by filled diamonds, (a) for n>5, this study, and (b) after Doin (1996).
CONCLUSIONS
We carried out two types of inversion using as input data either high-pass-filtered topography and geoid or local slopes of the topography with square root of age and local slopes of the geoid with age. These data should be little affected by the dynamic effect of deep mantle masses and surface velocities. They present complementary information because of their various frequency contents. They also have different drawbacks and advantages. The inversion using local slopes is at the limit of resolution of the various data sets used here, particularly the age data set from Mueller (1993) and the ETOPO5 bathymetry data set, whereas these data sets are clearly accurate enough for the wavelength considered in the global inversion (larger than 1000 km). Moreover, the global inversions average data evenly over the Earth's oceans, whereas the inversions using local slopes average data only where these slopes are well defined, that is, in small portions of the Earth's oceans, mostly in areas with ages <80 Myr. The inversion using local slopes provides realistic estimates of errors in the flattening coefficients. Furthermore, the departure from the half-space model is clearly visible early (before 100 Myr) in the dtopo/d
versus age or dgeoid/dage versus age curves. Therefore, the small surface coverage of old seafloor heavily covered by hotspots swell (>100 Myr) does not lead to doubts about the credibility of the observed flattening.
Note that the inversions here consisted of finding an average cooling behaviour (or subsidence depth–age curve) for the oceanic lithosphere, although local variations in subsidence rates are great. These variations are expected not only because of lateral variations in the cooling history of the lithosphere (variations in asthenospheric heat flow or in temperature) but also because of heterogeneities in crustal thickness and volcanic structures. These effects should be cancelled out by the averaging process.
The flattening of the depth–age and geoid–age relationships in oceanic basins is observed at various scales, from local observations on 500×500 km areas to global observations for wavelengths greater than 4000 km. The inverted subsidence curves at short and medium wavelengths are consistent with those previously noted by many authors (e.g. Stein & Stein 1992) using unfiltered topography data. When subsidence curves are parametrized in p1
(CHABLIS model) or in p2
2 (Plate model), all types of analyses yield results consistent with p1∼−440 m Myr−1/2 and b∼16 m Myr−1, and p2∼−310 m Myr−1/2 and c∼0.050 m Myr−2. The inverted geoid curves at short wavelength are also consistent with geoid jumps across fracture zones e.g. Cazenave 1983) or observed geoid slopes (e.g. Sandwell & Schubert 1980). When geoid–age curves are parametrized in ΔN1t+βt3/2 or in ΔN1t+γt2, all types of analyses yield ΔN1∼−22 cm Myr−1 and β∼1.2 cm Myr−3/2, ΔN2∼−14 cm Myr−1 and γ∼0.055 cm Myr−1/2. One can therefore define an average cooling behaviour of the lithosphere independent of the observation scale. This study is complementary to previous studies, in particular to studies concerning geoid jumps across transform faults. Its advantage is in analysing global geoid and topography sets simultaneously, sweeping a large wavelength spectrum. However, note that as geoid and topography data are noisy, all analyses yield only a rough estimate of the average subsidence–age and geoid–age curves due to cooling.
We conclude that hotspot swells and geoid anomalies or dynamic topography due to the return flow in the mantle or to the effect of lower mantle masses cannot be responsible for the deviation of the observed bathymetry from the prediction of the half-space cooling model. Humler (1999) suggested that the topography flattening might be partly due to a thicker Cretaceous than present-day seafloor crust (by about 200–400 m). However, this crustal thickening should not affect the geoid anomalies. The flattening of both topography and geoid strongly support models where heat is carried to the base of the lithosphere, slowing down the cooling of the lithosphere. However, the geoid and bathymetric data are very noisy, too noisy to be used alone to discriminate between the Plate and CHABLIS models. In order to explain the flattening of the geoid–age and topography–age relationships, depending on which model is used, the heat carried by small-scale convection or by shear flow beneath the oceanic lithosphere must range from 35 to 50 mW m−2 see Appendix A and Doin & Fleitout 1996). As the flattening does not seem to depend upon the absolute velocity of the plates, we favour the existence of small-scale convection. This constrains the mode of heat transport and the viscosity profile in the Earth's mantle.
APPENDIX A:PARAMETRIZATION OF LITHOSPHERIC COOLING MODELS
To simplify the inversion of eqs (1) and (2), the lithospheric cooling models FT(age) and FN(age) are linearized as follows. If the oceanic lithosphere cools by pure conduction as a half-space medium, the seafloor subsidence is proportional to the square root of age, and the geoid decreases linearly with age (Davis & Lister 1974). Doin & Fleitout (1996) showed that the subsidence departure from the square root of age law (or ‘flattening’f) is expected to be proportional to the age if a constant heat flow, qm, is advected by convection at the base of the lithosphere (f≈b age). For the Plate model this departure can be expressed in a first approximation by a term in age squared (f≈c age2). FT(age) can then be written as follows. For the half-space model:
for the CHABLIS model:
and for the Plate model:
with
and
where α is the thermal expansion coefficient, ΔT is the temperature drop across the lithosphere, κ is the thermal diffusivity, Cp is the heat capacity, k is the thermal conductivity, and ρm and ρw are the mantle and water volumic masses. Note that the flattening coefficients b and c vary with the heat carried to the base of the lithosphere, qm.
Similarly, the flattening of the geoid–age curve is expressed roughly by adding to the ‘reference’ linear law valid for the half-space model,
a term either β age3/2 representing an early development of small-scale convection (for a constant basal heat flow model),
or a term γ age2 more appropriate if small-scale instabilities develop later (for a Plate model),
(see Fig. A1). In these expressions,
Comparison of the predictions of the CHABLIS and Plate models (thick dashed and dot-dashed lines, respectively) with the analytic functions FN–CH and FN–PL. Time and geoid height are non-dimensionalized using the conductive timescale through the equilibrium thickness of the lithosphere and 2πGρmΔTακ/g0, respectively. For the same set of parameters, the geoid decrease, that is, the cooling, is much slower for the CHABLIS model than for the Plate model. When the parameters are adjusted to fit oceanic data, the maximum non-dimensional time (≈0.33 for the CHABLIS model and ≈0.4 for the Plate model) corresponds to about 140 Myr.
Comparison of the predictions of the CHABLIS and Plate models (thick dashed and dot-dashed lines, respectively) with the analytic functions FN–CH and FN–PL. Time and geoid height are non-dimensionalized using the conductive timescale through the equilibrium thickness of the lithosphere and 2πGρmΔTακ/g0, respectively. For the same set of parameters, the geoid decrease, that is, the cooling, is much slower for the CHABLIS model than for the Plate model. When the parameters are adjusted to fit oceanic data, the maximum non-dimensional time (≈0.33 for the CHABLIS model and ≈0.4 for the Plate model) corresponds to about 140 Myr.
and
Again, the flattening coefficients β and γ vary with the heat carried to the base of the lithosphere, qm.
APPENDIX B:LEAST-SQUARES REGRESSION
The equations (11)–(16) are in the form Yi=a+bXi. They are inverted by a non-linear least-squares procedure (Tarantola & Valette 1982) that takes into account errors on both Yi and Xi. The parameters are a, b and Xi, with a priori values of a0, b0, and Xi0, and standard errors σa, σb and σXi. The data are Yi, with their standard errors σYi. This procedure yields an estimate of errors that is considered to be reliable if the slopes Yi form a set of independent data points.
The local slopes Yi=dtopo/d
or Yi=df(geoid)/df(age) are computed on non-overlapping areas covering surfaces from 5°×5° to 17°×17°. In each given area i, Yi is the coefficient of linear least-squares regression between the observed topography (y) and square root of age (x) fields or between the observed geoid (y) and age (x) fields; σYi is the standard deviation of the error on Yi:
where var(x)i, var(y)i and cov(x, y)i are the variances of x and y and the covariance of x and y, respectively, computed on a given area i with a node every 1°×1°.
The average values Xi0=2 
i and Xi0=4 
i3/2 that appear in eqs (12) and (13) and their standard deviations σXi are computed over each given area i. There is no straightforward derivation here of Xi0=df(age3/2)/df(age) and Xi0=df(age2)/df(age), which appear in eqs (15) and (16). They are therefore calculated by a linear regression of the corresponding fields in each given area i as is done for Yi. σXi is also computed in the same way as σYi.
Acknowledgments
M-PD thanks D. Sandwell, who made her stay at the Scripps Institution of Oceanography possible. This work was supported by the French Ministry of Foreign Affairs and by a NATO grant. We thank the reviewers for their helpful comments on the manuscript.
