## Abstract

A digital bathymetric map of the oceans with a horizontal resolution of 1 to 12 kilometers was derived by combining available depth soundings with high-resolution marine gravity information from the Geosat and ERS-1 spacecraft. Previous global bathymetric maps lacked features such as the 1600-kilometer-long Foundation Seamounts chain in the South Pacific. This map shows relations among the distributions of depth, sea floor area, and sea floor age that do not fit the predictions of deterministic models of subsidence due to lithosphere cooling but may be explained by a stochastic model in which randomly distributed reheating events warm the lithosphere and raise the ocean floor.

Knowledge of ocean floor topography data is essential for understanding physical oceanography, marine biology, chemistry, and geology. Currents, tides, mixing, and upwelling of nutrient-rich water are all influenced by topography. Seamounts may be particularly important in mixing and tidal dissipation (1), and deep water fisheries on seamount flanks have become economically significant (2). Seamounts, oceanic plateaus, and other geologic structures associated with intraplate volcanism, plate boundary processes, and the cooling and subsidence of the oceanic lithosphere should all be manifest in accurate bathymetric maps.

Conventional sea floor mapping is a tedious process. Ships have measured depth with single-beam echo sounders since the 1950s, but these data are sparsely distributed (hundreds of kilometers between surveys) and may have large errors in navigation and digitization (3). More accurate multibeam swath-mapping systems came into use on some ships in the 1980s, but in the deep ocean, these were deployed primarily along mid-ocean ridges (4). Some surveys are classified as secret in military archives (5) or remain proprietary for economic or political reasons. Global bathymetric mapping requires some means of combining these heterogeneous soundings and estimating depths in the regions where survey data are sparse. Traditionally, bathymetric contours have been drawn by hand so that intuition (or prejudice) fills the gaps in coverage. The contours may then be digitized and interpolated to produce gridded estimates. The last global syntheses were made in the late 1970s and early 1980s, yielding the fifth edition of the General Bathymetric Charts of the Oceans (GEBCO) (6) and the Earth Topography 5-arc-min grid (ETOPO-5) (7).

Recent developments allow another approach to this problem. International cooperation has yielded access to a greater variety of sounding data (4, 6, 8), and automated quality control (3, 9) and archiving (10) methods have been devised. In addition, the ERS-1 and Geosat spacecrafts have surveyed the gravity field over nearly all of the world's ocean areas (11, 12). Over the 15- to 200-km wavelength band, marine gravity anomalies are caused primarily by topographic variations on the ocean floor; thus, in principle, satellite gravity data can be used to infer some aspects of the ocean's depths. However, the topography/gravity ratio varies from one region to another because of changes in sediment thickness and other factors, so that the estimation of topography from gravity is not straightforward and requires accurate depth soundings for calibration. Here, we report our efforts to combine quality-controlled ship depth soundings with interpolation guided by satellite-derived gravity data to yield a high-resolution grid of sea floor topography.

We assembled digital depth soundings from the U.S. National Geophysical Data Center (NGDC) (8), the Scripps Institution of Oceanography (SIO) (13), and two databases derived from data originally archived at the Lamont Doherty Earth Observatory (LDEO) (3, 9, 14). Although the four sources have many cruises in common, each has unique strengths (15). We also obtained recent survey data directly from various investigators (16). Global gridded gravity anomaly data were derived from satellite altimeter measurements of sea surface slope (11,17). The accuracy of the gravity grid is 3 to 7 mgal (1 milli-Galileo = 10^{–}
^{5} m/s^{2}) with a resolution limit of 20 to 25 km, depending on factors such as local sea state and proximity to areas of high mesoscale ocean variability (11, 18). Land elevation and shoreline data (19) furnished additional constraints in nearshore areas.

**Method and limitations**. Dixon *et al.*(20) have summarized the basic theory for estimating sea floor topography from gravity anomalies. Models of the isostatic compensation of sea floor topography furnish a spectral transfer function that predicts the gravity anomaly expected from sea floor topography (21). This transfer function is isotropic and depends on mean depth, crustal density and thickness, and elastic lithosphere thickness, and although its inverse provides a theoretical basis for estimating sea floor topography from observed gravity anomalies, there are a number of complications that require careful treatment (20, 22-25). (i) The estimation must be restricted to a limited wavelength band, because the gravity-to-topography transfer function becomes singular at wavelengths much longer than the flexural wavelength of the lithosphere and also at wavelengths much shorter than 2π × the mean depth, because of isostatic compensation and upward continuation, respectively. (ii) At short wavelengths, the transfer function depends on well-constrained parameters (mean depth and crustal density), but at longer wavelengths, it also depends on more variable and less certain parameters (elastic lithosphere thickness or crustal thickness). (iii) Sedimentation preferentially fills bathymetric lows and can eventually bury the original topography, adding a spatially dependent and nonlinear aspect to the transfer function. (iv) The transfer function is two dimensional and so requires complete coverage of gravity data. (v) The transfer function is nonlinear (22) in areas of high-amplitude topography, especially where summits approach the sea surface.

Our method uses a Gaussian filter with a half-amplitude at 160-km wavelength to smoothly separate the topography into two bands: a long-wavelength regional topography, where the transfer function requires an assumed elastic thickness, and a short-wavelength local topography, where the transfer function is independent of the elastic thickness. We restricted gravity-derived estimation to the shorter wavelength band (26). We regionally calibrated the topography/gravity ratio using the local topography at points constrained by soundings, to adjust the transfer function for the effects of regionally varying sediment thickness. Our earlier method (23) was improved by adding a constraint-propagation step: grid cells constrained by data were set to the median of data values in the cell, and then a finite-difference interpolation routine (27) was used to perturb neighboring estimated values toward the constrained values. Thus, in well-surveyed areas, the accuracy and resolution depended only on the grid spacing (17) and the quality of the constraints and not on the gravity data. Constrained and estimated values were separately encoded so that they could be distinguished. During our process, inaccurate ship soundings became evident that were not detectable in initial quality control (3); these soundings were deleted from the databases, and a new estimate was made. This approach differed from another solution (24) in several respects (28).

**Results, verification, and assessment**. Our topography (29) reveals all of the intermediate- and large-scale structures of the ocean basins (Fig. 1). Incised canyons are seen in the continental margins. Spreading ridges stand out as broad highs with an axial valley along the Mid-Atlantic Ridge and an axial high along the East Pacific Rise and the Pacific-Antarctic Rise. Fracture zones reflect the direction of opening of the Atlantic basin, whereas in the Pacific they record a more complex history of major plate reorganizations. Numerous seamounts, some in linear chains, display a variety of patterns of volcanism.

We tested the precision and resolution of our estimates in a worst case scenario: a high-relief area lying more than 160 km from most soundings. After making a topography estimate, we obtained depths measured by the research vessel (R/V) *Atalante*(30) in a remote area of the South Pacific near the Foundation Seamounts (31-33). *Atalante* found seamounts with summits less than 1 km from the surface and a 6500-m-deep trough in this area, where the ocean floor typically lies at 4-km depth. The root-mean- square (rms) amplitude of these topographic variations was 791 m, whereas the rms difference between our estimates and the observed values was 250 m; thus, the estimates recovered nearly 70% of the signal. Cross-spectral coherency (Fig. 2) shows high correlation between our estimates and the observed depths at all wavelengths greater than 25 km, with a slight decrease at wavelengths greater than 100 km, where the estimation relied on interpolated soundings more than altimetry. The low coherency at wavelengths less than 25 km means that in the estimated topography, two narrow objects may blur into one object if they are closer than 12.5 km apart and objects much narrower than 12.5 km will be poorly resolved.

The Foundation Seamounts were unknown until they were revealed by satellite altimetry (33). We compared the*Atalante* depths with ETOPO-5 depths to assess how much information was missing in older maps; the (rms) difference between the two was 580 m, more than 70% of the signal. Uncharted seamounts were a significant source of topographic variation, and information from satellite gravity can reduce the error in estimated topographic variation by more than half.

**Hypsometry and thermal subsidence models**. Various estimates of the distribution of ocean floor area with depth, called hypsometric curves, are shown in Fig. 3A. Previous curves (34) were calculated in 1-km intervals. The ETOPO-5 data cannot yield a more detailed curve because of biases toward multiples of 100, 200, and 500 m, the contours that were digitized to produce ETOPO-5. Our solution yielded a smooth curve at 50-m intervals. Viewed in 1-km intervals, our solution had more area in the 3- to 4-km range and less in the 5- to 6-km range than was seen previously, reflecting the increased number of seamounts mapped by satellite altimetry (35).

We interpolated ages (36) and sediment thicknesses (37) to our depth grids, excluded anomalous or problematic areas (38) but not thermal swells, and isostatically corrected the depths in the remaining areas for the sediment load using an average sediment density (39). Grouping the areas of the remaining data into intervals of 50 m of depth and 1 million years (My) of age shows the variation of depth with age (Fig.4A). These exclusions and corrections affected the hypsometry (Fig. 3B); the sediment correction slightly increased the area deeper than 5.2 km, whereas the exclusions reduced the area elsewhere. The area reduction was roughly independent of depth at depths less than 3.5 km. There is an overall decrease of area with age, as expected if spreading ridges generate area at a constant rate, whereas subduction zones consume area independently of age (Fig.5) (40). Superimposed on this are steplike changes that primarily reflect variations in spreading rate but are also caused by our exclusions of some areas (41). It appears that a large area of sea floor was generated during the 20 to 33 million years ago (Ma) period, implying more rapid sea floor generation with a consequent rise in sea level (42) and a change in global carbon chemistry (43).

From the distribution of area with depth within each age interval, we found the modal depth, the median depth, the first and third quartile depths, and the interquartile range (IQR) (44) (Fig. 4A). The mode and quantiles deepened with age from 0 to 55 Ma and then abruptly flatten and become nearly constant at ages older than 70 Ma. Comparing these statistics with the area distribution (Fig. 4A) shows that in every age interval the distribution of depths is skewed, with more area above the mode than below it. Seamounts would cause this. The IQR was about 0.5 km for the first 27 My, increased to 0.75 km at 80 Ma, and stayed nearly constant at older ages. These statistics are not reliably determined at ages older than 130 Ma because the available sample of area is small.

Thermal models of the cooling and subsidence of oceanic lithosphere (45) predict that depth *z* should increase monotonically with the age of the sea floor *t* as the lithospheric plate moves away from the hot mid-ocean ridge where it formed. Because these models treat *z* deterministically, they will not explain why the variation we observe in *z* increases with *t*, and we therefore compare them with the median*z* at each *t*. The simplest model is the “boundary layer” (BL) model (46), which predicts that*z* should increase with *t* according to(1)where *r* = 2.5 km is the depth at zero age and *a* = 0.35 km/(My)^{1/2} is the subsidence rate. This model gives a good fit to the median depths (Fig.4B) in the 0 to 55 Ma age range, but not beyond that, where the observed depths flatten. A model that imitates BL for a time and then flattens is the “plate model,” which predicts that(2)in which *h* is the overall subsidence from zero to infinite age and τ is a characteristic time. For*t* < τ, Eq. 2 gives approximately the same behavior as Eq. 1 with *a* = 4*h*/π^{3}τ, whereas for larger*t*, the *z*(*t*) curve will flatten. A plate model with *r* and *a* as in BL and τ = 62.8 My, we call model “PS” after its authors (45); a model with *r* = 2.6 km, *a* = 0.365 km/(My)^{1/2}, and τ = 36 My was called “Global Depth and Heat Flow 1” or GDH1 by its authors (47). Model PS fits the median depths in the same range as model BL (Fig. 4B) because it has the same initial behavior; it also begins to flatten at about the right time. However, the flattening is not abrupt enough to fit the data. Model GDH1 achieves greater flattening by using a smaller τ than model PS, but to have approximately the right overall subsidence*h*, GDH1 increases *a* to compensate for the decreased τ; this causes it to predict too much initial subsidence and a too early onset of flattening while still failing to follow the abrupt flattening observed. Hypsometric curves predicted by models PS and GDH1 (48) (Fig. 3C) are not a good match to the observations; they imitate the variations in the area-age distribution too closely, they predict less area than is observed at depths shallower than 5000 m, and they predict that large areas of the ocean floor should lie at the asymptotic depth of the model.

We propose a stochastic reheating model to explain the hypsometry, the abrupt flattening of median depth with age, and the increased variability of depth at older ages. It is based on the following observations: (i) different regions in the oceans subside like BL until different flattening ages, some as old as 170 Ma (49); (ii) depth and heat flow in some old basins behave like BL but with an effective thermal age *t*
_{e,}which is younger than the (paleomagnetically determined) actual age*t*
_{m} (50); (iii) effective elastic lithosphere thicknesses under some intraplate volcanoes suggest that hot spots produce a “thermal rejuvenation,” reducing*t*
_{e} (51); (iv) the depths and subsidence of hot spot swells suggest that many swells reduce*t*
_{e} to the same value (52); and (v) on old lithosphere, it is hard to find areas that have not been near hot spots (53). If *N* sites of reheating, each affecting an area of diameter *D*, are uniformly and randomly distributed in an area *A*
_{h}, and the plates making up *A*
_{h} move over these sites with velocity*V*, then in a steady-state situation with*A*
_{h}, *D*, *N*, and *V*all independent of time, the probability density function*p*(*t*) for *t*, the time since any point's last reheating, is(3)where(4)is the mean time between reheating events (54). If reheating events always reset*t*
_{e} to a fixed value *t*
_{0}, whereas lithosphere younger than *t*
_{0} is not affected by reheating, and if the sea floor we observe is a random sample from such a system, then *X*, the chance that sea floor of (actual) age *t*
_{m} has not been reheated since*t*
_{0}, is(5)and the probability density function for*t*
_{e} given an observed *t*
_{m}is(6)This model predicts that (i) when *t*
_{m} ≤*t*
_{0}, *t*
_{e} is deterministic and *t*
_{e} = *t*
_{m}; (ii) as *t*
_{m} exceeds *t*
_{0}, a fraction (1 − *X*) of the area has*t*
_{e} < *t*
_{m}; (iii) when*t*
_{m} > *t*
_{0} + α log 2, the median *t*
_{e} is less than*t*
_{m}; and (iv) at large *t*
_{m},*X* → 0 and a steady state prevails in which*t*
_{e} is independent of *t*
_{m}. In the steady state, the mean effective age is μ_{t} = *t*
_{0} + α and the standard deviation of the effective ages is σ_{t}= α.

If *z* follows the BL cooling rule with*t*
_{e} in Eq. 1, then the probability density for*z* given *t*
_{m} is(7)where(8)and *z*
_{0} and*z*
_{m} are the depths corresponding to*t*
_{0} and *t*
_{m} in Eq. 1. The change of variables from *t*
_{e} to *z* is nonlinear so the mean and standard deviation of *z* cannot be found from μ_{t} and σ_{t}; however, the quantiles (44) can be transformed:*z*
_{Q} = *r* +*a*
always, and*t*
_{Q} = *t*
_{0} − α log(1 − *Q*) in the steady state.

The predictions formed by our model with *t*
_{0} = 45 My and α = 25 My combined with a Gaussian noise process having a zero mean and an IQR of 0.5 km at all ages are shown in Figs. 3C and 4C. We added this noise to simulate the variability of the data at ages younger than *t*
_{0}, where the model is deterministic. Taking *A*
_{h} = 115 × 10^{12} m^{2}, the area of the sea floor older than 45 Ma, if *D* and *V* are roughly 1000 km and 80 km/My, respectively, then α = 25 My gives *N* = 57 reheating sites in this area, and α*V*, the characteristic distance between reheating sites, is about 2000 km. We were able to obtain reasonable fits to the data in Figs. 3C and 4C with other values as well (for example, *t*
_{0} = 40 and α = 30), but we caution against overfitting the data; the Gaussian noise, being symmetric, does not realistically simulate the skewed variability produced by seamounts.

If one merely desired to explain the hypsometry, any reheating process that converted *t*
_{m} to *t*
_{e}such that the distribution of area with *t*
_{e} was as shown by the curve in Fig. 5 would suffice. This curve is obtained from the hypsometric curve by converting *z* to*t*
_{e} with the inverse of Eq. 1. It shows that*t*
_{e} is usually less than*t*
_{m} and that areas with *t*
_{e}older than 130 My are rare. Reheating that is concentrated primarily in one area (55) or at one age (56) might explain some features of the area-*t*
_{e} distribution in Fig. 5, but the age independence of the depth distribution at old age (Fig. 4A) suggests a model with reheating events randomly distributed in space and time.

Boundary layer cooling is the physically expected behavior unless some other thermal process intervenes. The plate model assumes that a fixed temperature is maintained at a fixed depth everywhere and so predicts that subsidence curves should flatten at the same age everywhere; they do not (49). Our model indicates that from time to time some areas will probably be reheated (57) and uplift and rejuvenated subsidence may occur anywhere at any time after*t*
_{0}. The depth at any individual site older than*t*
_{0} may be between*z*
_{0} and *z*
_{m}, but the aggregate properties of a large sample of sites should have order statistics (44) predicted by the model; these statistics can flatten more abruptly than the plate model subsidence curve if α is small compared with τ. Old areas are rare in comparison with young areas (Fig. 5) and constitute a small sample. Thus, the fact that some sites on old sea floor have depths lying on the PS plate model curve (45, 50) does not disprove our model. The heat flow at these sites argues for reheating and against the plate model (50).

↵* To whom correspondence should be addressed. E-mail: walter{at}amos.grdl.noaa.gov