Research Article

Global Sea Floor Topography from Satellite Altimetry and Ship Depth Soundings

See allHide authors and affiliations

Science  26 Sep 1997:
Vol. 277, Issue 5334, pp. 1956-1962
DOI: 10.1126/science.277.5334.1956


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/s2) 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.

Figure 1

Color shaded-relief image of the modeled topography, illuminated from the northwest.

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.

Figure 2

Coherency between our depth estimates and theAtalante measurements as a function of wavelength, with one sigma error bars. Coherency values of 1, 0, and 0.5 correspond to perfect correlation, complete absence of correlation, and equal magnitudes of correlated and uncorrelated components, respectively.

The Foundation Seamounts were unknown until they were revealed by satellite altimetry (33). We compared theAtalante 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).

Figure 3

Hypsometric diagrams showing sea floor area per 50-m interval of depth, as a function of depth. (A) Using all data. Our solution (solid line) smoothly resolves 50-m intervals, whereas ETOPO-5 (7) (gray bars) is biased toward contour values. Previous estimates (34) (dashed line) gave only 1-km intervals. Our solution for 1-km intervals (dotted line) has less area in the 5- to 6-km bin and more area in the 3 to 4-km bin than previously found, because the satellite gravity reveals more seamounts. (B) Solution from all data [solid line, same as in (A)] and from a restricted data set corrected for sediment thickness (gray area). (C) Observed hypsometry of the restricted data set [gray area, same as in (B)] and predictions of three model curves: plate models PS (45) (dashed line) and GDH1 (47) (dotted line) and our model (solid line).

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).

Figure 4

Sediment-corrected depth as a function of sea floor age. (A) Area of the sea floor in 50-m depth and 1-My age intervals (colors, scale at bottom). Upper and lower quartiles of depth (thin lines), median depth (heavy line), and modal depth (dots) are shown in each age interval from 0 to 130 Ma, where they are reliably determined. (B) Observed median depth [heavy line, same as in (A)] and depth predicted by deterministic models PS (45) (thin solid line), GDH1 (47) (dashed line), and BL (46) (dotted line). (C) Synthetic data predicted by our model give an area distribution (colors, scale at bottom) and a median and quartiles of depth at each age (thin lines); the observed median [heavy line, same as in (A) and (B)] is shown for comparison.

Figure 5

Sea floor area of the sediment-corrected data set as a function of age. The gray pattern shows the relation obtained from actual sea floor age t m. The solid line shows synthetic relation obtained from hypsometric curve by converting actual depth to thermal age t e with the use of the BL model.

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 medianz at each t. The simplest model is the “boundary layer” (BL) model (46), which predicts thatz should increase with t according toEmbedded Image(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 thatEmbedded Image(2)in which h is the overall subsidence from zero to infinite age and τ is a characteristic time. Fort < τ, Eq. 2 gives approximately the same behavior as Eq. 1 with a = 4h3τ, whereas for largert, 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 subsidenceh, 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 aget m (50); (iii) effective elastic lithosphere thicknesses under some intraplate volcanoes suggest that hot spots produce a “thermal rejuvenation,” reducingt e (51); (iv) the depths and subsidence of hot spot swells suggest that many swells reducet 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 velocityV, then in a steady-state situation withA h, D, N, and Vall independent of time, the probability density functionp(t) for t, the time since any point's last reheating, isEmbedded Image(3)whereEmbedded Image(4)is the mean time between reheating events (54). If reheating events always resett 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 sincet 0, isEmbedded Image(5)and the probability density function fort e given an observed t misEmbedded Image(6)This model predicts that (i) when t mt 0, t e is deterministic and t e = t m; (ii) as t m exceeds t 0, a fraction (1 − X) of the area hast e < t m; (iii) whent m > t 0 + α log 2, the median t e is less thant m; and (iv) at large t m,X → 0 and a steady state prevails in whicht 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 witht e in Eq. 1, then the probability density forz given t m isEmbedded Image(7)whereEmbedded Image(8)and z 0 andz m are the depths corresponding tot 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 Embedded Image always, andt 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 × 1012 m2, 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 esuch 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 tot e with the inverse of Eq. 1. It shows thatt e is usually less thant m and that areas with t eolder 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 aftert 0. The depth at any individual site older thant 0 may be betweenz 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}


Stay Connected to Science

Navigate This Article