The oceanic age-depth relationship is a fundamental consequence of Earth’s convective system. The new analysis presented here integrates four major data sets to assess the global age-depth distribution of the oceanic lithosphere. First, the distribution of depths along the modern mid-oceanic-ridge axis was obtained using an updated digitized divergent plate-boundary data set for the main ocean basins. Mid-ocean ridges vary in depth by ∼8 km, from −6453 m to +1998 m, with a median depth of −3000 m, modal depth of −2670 m (with a second mode at −2770 m), and area-weighted mean depth of −3044 m. Ridge depth has some correlation with spreading rate up to about 50 mm/yr if the hotspot-impacted North Atlantic and Sheba Ridges are excluded. Second, zero-age crustal thickness, based on a compilation of ∼1000 estimates, was used to test the relationship of oceanic crustal thickness with axial depth. There is some correlation, but it is not a compact correlation, implying that other processes are contributing to both axial depth and crustal thickness variations. Third, a revised sediment-load correction is presented based on assessment of 9,260,185 records from 1263 deep-sea drilling sites with wet bulk density measurements. This can be extrapolated to thicknesses in excess of 8000 m, with an exponential porosity-depth relationship with an initial porosity of ϕ0 = 0.61. Dispersion of integrated wet bulk density measurement-derived sediment-load corrections for the 1263 sites implies ∼±100 m uncertainty in the sediment-load correction at 1300 m sediment thickness. Fourth, the age-depth relationship of oceanic lithosphere was derived based on sediment-load corrected depths of ∼90,000 explicitly dated, globally distributed magnetic reversal picks restricted to the main ocean basins and filtered to exclude large igneous provinces and trench-impacted locations. A best-fit global age-median sediment-corrected depth relationship was derived using a plate model with zero-age ridge depth of −2309 m, temperature difference of 1304 °C, and asymptotic plate thickness of 122 km. The dispersion of sediment-load corrected depths as a function of age is comparable to that observed along the mid-oceanic-ridge system, suggesting that zero-age depths contribute to dispersion at all ages. This is borne out for oceanic lithosphere at relatively young ages (ca. <70 Ma), but the dispersion at older ages appears to be dominated by superposition of hotspots and other large-scale perturbations. In addition, a comparison was made with various vintages of age-grid–derived age-depth distributions, as well as previous age-depth models. The implied age-depth relationship using older versions of the age grid deviates significantly from the age-depth relationship obtained using the magnetic reversal picks alone. There is considerable disparity between pick ages and corresponding grid cell ages for these data sets. In contrast, the most recent age grid (2016) more closely matches the age-depth relationship derived here and shows significantly less disparity with magnetic reversal pick ages. A qualitative suggestion is made that the deviation of the age-depth relationship from a half-space cooling to exponential model reflects the advective contribution of asthenospheric flow associated with asthenospheric tractions contributing to plate-driving forces.
The bathymetry of the oceanic lithosphere is the most obvious surficial evidence on Earth of mantle convection directly coupled to plate tectonics. The variation of depth with age of the oceanic lithosphere has played an important role in our understanding of lithospheric creation and cooling. Absent this mantle convective signature, and hence the mantle dynamic contribution to oceanic topography, seafloor bathymetry would be nearly flat, reflecting the relative homogeneity of oceanic crustal thickness. Instead, there is ∼3 km of dynamic topographic relief from young to old oceanic lithosphere upon which are superimposed variations due to hotspot magmatism and shorter-wavelength, non-plate-cooling, dynamic topography contributions. Thus, establishing the age-depth relationship of the oceanic lithosphere potentially provides an important framework for identifying anomalous topography of the oceanic basement (Hoggard et al., 2017). The interested reader is directed to Hoggard et al. (2017) for a review of previous analyses and a discussion of various approaches to identifying anomalous ocean bathymetry relative to various age-depth relationships.
Simple theoretical grounds lead to an expectation that conductive boundary layer cooling of the oceanic lithosphere will give rise to depth varying as square root of age (Turcotte and Schubert, 1982). Parsons and Sclater (1977) showed that square root of age also arises from a plate cooling model (McKenzie, 1967). Early analyses of the age-depth distribution of the oceans focused on data from particular areas (Nagihara et al., 1996; Parsons and Sclater, 1977; Sclater et al., 1980; Stein and Stein, 1992). Recently, several analyses have used global data sets to try to determine the global age-depth relationship (Crosby and McKenzie, 2009; Crosby et al., 2006; Hoggard et al., 2017; Marty and Cazenave, 1989; van der Meer et al., 2017). Most existing analyses have separated the age-depth history into an interval of younger ages that follow a square root age relationship, and an older interval, beyond which depth increased less rapidly, fit by an exponential function conforming with a plate model. More recently, Hoggard et al. (2017), using ∼2300 widely dispersed locations constrained in terms of sediment load and about half by crustal thickness corrections, and with ages extracted mostly from Müller et al. (2008), derived best-fit half-space and plate-thickening models for their data for all ages. In their analysis, the plate model yielded a better fit to their data and thus they emphasized their best-fit plate model to describe the global age-depth distribution. The general approach of Hoggard et al. (2017) will be followed here, except that the crustal thickness correction will not be included, because the distribution of crustal thickness measurements is so limited, and according to Hoggard et al. (2017, p. 2337): “(t)he median size of the crustal correction is 200 m,” which is small compared with both the sediment-load correction and the dispersion of depth at all ages.
This equation has two independent parameters, zr and which can be varied to find a best-fit solution.
Equations 1 and 2 model the water-loaded depth as a function of age starting from a zero-age depth (zr) and subsidence as a function of age. Thus, in order to assess the depth-age relationship, it is necessary to (1) establish the correspondence between age and depth for each point used in the analysis; and (2) correct for effects of sediment loading. In addition, it is worth being aware of the spatial distribution of other sources of anomalous bathymetry, such as trenches, hotspot tracks and associated flexural moats, oceanic plateaus, and short-wavelength dynamic topography, among other sources, that may perturb the age-depth relationship. To better interpret the results, it is also useful to have data on the zero-age depth distribution that can be directly compared with model-derived estimates.
AGE OF OCEANIC LITHOSPHERE
Establishing the age of oceanic lithosphere is obviously a critical component of each age-depth analysis. It is important to recognize that oceanic lithosphere age is derived from four data types: (1) zero-age lithosphere along axes of present-day mid-oceanic ridges; (2) direct sampling of upper oceanic lithosphere by dredging or submersibles; (3) direct sampling by deep-sea drilling; and (4) locations of magnetic reversal picks derived from ship or airborne magnetometer surveys correlated to a magnetic reversal time scale (Fig. 1). Oceanic age maps, termed age grids, use plate-kinematic models derived by fitting data from magnetometer surveys together with oceanic transform/fracture zone systems to interpolate and extrapolate ages between and beyond known ages to fill in ages of virtually all oceanic lithosphere (Seton et al., 2012). Most previous analyses of the age-depth relationship have used various maps of ocean age distribution to estimate the ages corresponding to sites for which depth and sediment loading have been determined. It is generally not possible to know the basis upon which a given map location or grid cell age assignment has been made, and thus it is difficult to know the reliability of the derived age corresponding to a given location. A characteristic of the more recent age-grid data sets is that finite rotations for each plate pair have been estimated for specific magnetic reversal times that are then used to compute interval stage rotations to interpolate and extrapolate intervening ages under an assumption of constant spreading rate. The interpolation and extrapolation intervals can range up to 10 or more million years. For example, during the Cretaceous normal polarity superchron, which has an estimated duration of 37.6 m.y. between 120.6 Ma and 83.0 Ma (Gee and Kent, 2007), changes in fracture zone geometries imply changes in rotation parameters, but without any constraint on timing of changes in rotations, resulting in potentially large, and effectively unconstrained uncertainties of all ages within the Cretaceous normal polarity superchron. Müller et al. (2008a) published estimates of error based primarily on propagated interpolation and extrapolation rotation uncertainties, which are obviously helpful, but it is still difficult to know the full reliability of the existing age grid, and hence the reliability of age-depth models that depend on these sources for their constraint on oceanic lithosphere ages. This topic will be revisited below.
Of the three direct sources of oceanic lithosphere ages older than 0 Ma, by far the most data are provided by the distribution of magnetic anomaly picks. The distribution of magnetic reversal picks has been assembled in the Global Seafloor Fabric and Magnetic Lineation (GSFML) database (Seton et al., 2014; http://www.soest.hawaii.edu/PT/GSFML/; downloaded 10 May 2017). The GSFML database contains 101,418 total records (Data Repository Table S11). Of these, 2328 are duplicate records, 3120 are from back-arc basins, and 2347 are from regions with no sediment thickness data because they are located in regions outside the sediment thickness grid of Whittaker et al. (2013), and 1045 lie within 200 km of trenches or have been identified on already-subducted oceanic lithosphere. In addition, 680 picks correspond with grid cells defined as nonoceanic crust in the Müller et al. (2016) age grid, with values of NaN in the grid, whereas 876 correspond with nonoceanic cells in the Müller et al. (2013) age grid.
It is worth emphasizing that these are the reversal picks that have been used to determine the finite rotations that underlie the interpolation and extrapolation used to fill the various generations of the age grid (Müller et al., 1993, 2008a, 2013, 2016). These are thus the data upon which the age grid is predicated, and therefore they are essentially the only explicitly dated locations for oceanic lithosphere ages older than 0 Ma. This statement ignores the fewer than ∼1300 deep-sea drilling sites that either date basement or place a younger age bound on basement age. In a later section of this paper, a comparison is made between age-depth values derived from only the picks and those in the age grids of Müller et al. (2008a, 2013, 2016). For now, one reason for suspecting that the results might be different is provided by comparing the locations with known ages using data from reversal pick locations against the corresponding grid cell ages based on the most recent Müller et al. (2016) age grid (Fig. 2). There are 361 discrete reversal pick ages in the GSFML data set with age-grid ages distributed along the line of correlation but dispersed about the 1:1 line. This suggests that using the age-grid ages may contribute to dispersion of depths that may impact the resulting age-depth relationship. For this reason, the analysis presented here will focus solely on the age-depth relationship as determined from the locations of magnetic reversal picks, and hence with explicit ages.
Each magnetic reversal pick location reflects the position on Earth’s surface where a particular change in the magnetic reversal record has been identified (Fig. 1; Seton et al., 2014). The dating of that point in the magnetic reversal record then derives from correlation of the magnetic reversal sequence to age in millions of years. For the purposes of the present paper, ages of the magnetic reversals are based on the time scale of Gee and Kent (2007). This analysis is therefore similar to that of Johnson and Carlson (1992), who used ∼50 directly dated oceanic basement locations penetrated by Deep Sea Drilling Project (DSDP) and Ocean Drilling Program (ODP) boreholes where depth, age, and sediment thickness were independently constrained, but it expands that data set by about a factor of 2000.
METHODOLOGY AND APPROACH TO AGE-DEPTH ANALYSIS
In order to assess the age-depth distribution of the oceans, this paper presents: (1) an analysis of the modern-day mid-ocean-ridge axial depth distribution to examine the spectrum of zero-age depths together with an assessment of the relationship of axial depth to spreading rate and zero-age crustal thickness; (2) a new assessment of the sediment-load correction of oceanic basement; and (3) sediment-load–corrected depth as a function of age, where ages are constrained explicitly by the locations with magnetic anomaly picks. Several additional aspects are discussed, including: (a) Is there a simple relationship between the distribution of depths at zero-age and the spectrum of depths as a function of age? If so, what is the underlying reason for that relationship? (b) How does the magnetic reversal pick age-depth relationship compare with the age-depth relationship derived over the same geographic extent using age-grid ages?
MID-OCEAN-RIDGE AXIAL DEPTHS, SPREADING RATE, AND ZERO-AGE CRUSTAL THICKNESS
New oceanic crust and lithosphere are produced at mid-ocean ridges. The age-depth relationship described by Equations 1 and 2 defines a zero-age depth from which subsidence begins. It is useful to independently constrain the zero-age depth distribution so as to better interpret the processes controlling the age-depth distribution. In order to assess the zero-age depth distribution, the locations of the mid-ocean-ridge axis in the main oceanic basins were significantly updated from that of Bird (2003), shown here in Figure 1, by combining three data sources. (1) The youngest magnetic anomaly picks in the GSFML database were segregated and plotted. This database includes 136 picks at 0 Ma, 3749 picks of the 780 ka Brunhes-Matuyama reversal, and 3596 picks of the 2.581 Ma C2An.1ny boundary, which are particularly useful along slow spreading ridges. These reversals straddle and closely parallel the mid-ocean-ridge axis. (2) Shallow earthquake epicenters derived from the Global Centroid-Moment-Tensor (CMT) Catalog (http://www.globalcmt.org/CMTfiles.html/jan76_dec13.ndk) were used to help delineate the boundary where no reversals are present. (3) Finally, bathymetry and digitizing were done using Google Earth starting from the Bird plate boundaries (Bird, 2003). With the exception of Iceland, shallow continuations of, for example, the Sheba, Gakkel, and Gulf Rise ridges (Fig. 3) were not digitized. In addition, the Red Sea ridge was not included, because there are insufficient magnetic anomaly picks along this ridge axis in the GSFML data set and an absence of clear morphology with which to confidently trace the ridge axis along the entire length of this basin. The Red Sea is therefore not a part of the depth-age analysis.
The digitized ridge axis is a continuous polyline, including transforms, for each of the 15 main ocean basin divergent plate-boundary segments between triple junctions or to the points where they are buried by sediments in the Laptev Sea, head of the Gulf of Aden, and Gulf of California (Fig. 3). Hence, each polyline segment corresponds with ridge and transform segments separating individual plate pairs along their mutual boundary. The relatively coarsely digitized ridge axis was then sampled at a spatial interval of 0.5 arcmin, resulting in 54,921 points that define the ridge axis (Fig. 3; Data Repository Table S2 [see footnote 1]). The depth and source identifier (SID) were extracted from SRTM30_Plus (Becker et al., 2009) at each of these points. Transforms were excluded from the depth extraction algorithmically using MORVEL instantaneous rotations of the adjacent plate pairs (DeMets et al., 2010), which were used to determine whether sequential points along the polyline defined great circles with a declination that was within 50° of the great circle about the corresponding instantaneous rotation axis. Only adjacent points following approximately great circle trajectories about the relative plate-motion rotation were extracted. The total effective ridge length is 50,795 km long, along which ∼2.69 km2/yr of lithosphere is being produced, based on the MORVEL rotations (DeMets et al., 2010). The main ocean basin length weighted mean spreading rate for these ridges is 52.9 mm/yr.
Figure 4 shows the distribution of zero-age axial depths using a 10 m vertical bin in which the total area within each height bin was tabulated. The area of each 0.5 arcmin SRTM30_Plus (Becker et al., 2009) grid cell was used rather than frequency because SRTM30_Plus is a rectilinear grid that inherently more densely samples high latitudes relative to low latitudes, thereby potentially skewing a simple frequency tabulation. The SRTM30_Plus data set distinguishes between bathymetric cells that are constrained by direct ship soundings (presence of a nonzero SID) versus those interpolated using satellite gravity interpolations (SID = 0). This allows an assessment of the impact of only sampling ridge hypsometry from soundings or gravity interpolation or all points along the ridge axis (Fig. 4). Using this discriminator, 33,337 depths were constrained by ship tracks, and 20,584 were interpolated using gravity. As is clear from this comparison, there is no significant difference between the ship and gravity-interpolated hypsometries. For the purposes of the present paper, all points along the ridge axis were used.
The median is −3000 m, the mode at 1 m vertical resolution is −2841 m, and with 10 m binning, the mode is −2660 m, with a secondary mode at −2770 m that corresponds to the mode of the shiptrack-only data. The area-weighted mean (−3044 m) depths of the global main ocean basin mid-ocean-ridge system are reasonably tightly clustered, with a ±1σ range about the median of +539/–738 m, and a ±2σ range of +1726/–1568 m. The deepest point is at −6453 m, as constrained by ship-derived bathymetry at the eastern ridge-transform intersection of the Romanche transform (Fig. 3), and the highest point is at +1998 m in Iceland. As shown by the ±1σ values and deviations of the median and weighted mean depths from the mode, the distribution is skewed on the deep side, but it has a long tail on the shallow end (Fig. 4). The skewness largely derives from contributions from those ridge segments characterized by deeper median depths (see Data Repository Figs. S1A–S1T [see footnote 1]), including the Gakkel, SW Indian, and South and Central Mid-Atlantic Ridges, as well as the deep tail on the Australia-Antarctic Ridge associated with the Australia-Antarctic Discordance (Fig. 3). The Pacific-Antarctic Ridge has a narrow depth range with a well-defined mode at −2460 m that counterbalances some of the skewness from the deeper ridges. Only the North Atlantic and Sheba Ridges are characterized by broad shallow distributions, but they do not balance out the contributions by the deeper ridges.
The relationship between ridge length-weighted mean spreading rate and median, mode, and area-weighted mean depths of the global main ocean basin ridges and the 15 individual plate-boundary segments is shown in Figure 5A. This plot is similar to one in Gale et al. (2014), except that Figure 5A plots length-weighted mean rate for the entire boundary segment rather than the rate at the midpoint of some smaller subdivision of each plate-boundary segment. This is a more representative rate because it accounts for the nonlinearity of spreading rate as a function of colatitude. In addition, this plot includes the modal and median depths, not just the mean, as the median and mode tend to be more robust representations of populations of data characterized by skewed distributions. The fastest spreading segments, the Nazca-Pacific (148.2 mm/yr) and Cocos-Pacific (104.1 mm/yr), comprising the East Pacific Rise, are tightly clustered with modal depths of ∼–2800 m. The slowest spreading segments, along the Southwest Indian Ridge between Nubia and Antarctica (13.9 mm/yr), between Somalia and Antarctica (14.2 mm/yr), and along the Antarctic-America Ridge (17.3 mm/yr) southwest of the Bouvet triple junction, are associated with the deepest median (−3619 m, −3669 m, and −3867 m, respectively) depths. Modal and weighted-mean depths show a similar pattern. However, the next two slowest spreading ridges, the North Mid-Atlantic Ridge between North America (NA) and Eurasia (EU) and Sheba Ridge between Arabia (AR) and Somalia (SO), have median and weighted-mean depths shallower than the global values. Both of these ridges display the largest variations in ridge depths, presumably reflecting their proximity to large hotspots. Both also extend above sea level. If one ignores the Iceland- and Afar-impacted North Mid-Atlantic Ridge and Sheba Ridge, then there is some support for a shallowing in ridge depth with increasing spreading rate. However, the Chile Rise and Juan de Fuca Ridge have comparable spreading rates and yet have mean and median depths that are different by >600 m, although their modal depths are quite similar. Finally, the Gulf Rise is comparable with the Juan de Fuca and Pacific-Antarctic Ridges, indicating little systematic variability for rates above 50 mm/yr.
An alternative assessment uses the local spreading rate at each point along the ridge axis and plots the distribution of depths as a function of spreading rate independent of plate-boundary segment, as shown in Figure 5B. If the hotspot-impacted North Mid-Atlantic Ridge and Sheba Ridge (20 mm/yr) are excluded, then there is a reasonably clear shallowing in mid-oceanic-ridge axial depth with increasing spreading rate up to ∼70 mm/yr. At faster rates, there is no clear trend.
The dispersion of depths at zero-age has been suggested to be linked to variations in crustal thickness that in turn link back to variations in mantle melting associated with mantle potential temperatures (Behn and Grove, 2015; Klein and Langmuir, 1987; among many others workers). This link has been tested locally using ridge-parallel seismic reflection and refraction surveys to map crustal thickness in relation to bathymetry (Fig. 6A). It is thus reasonable to determine whether the dispersion of zero-age depths seen globally (Fig. 4) reflects a comparable dispersion in zero-age crustal thickness. In order to test this possibility, crustal thickness estimates derived from seismic surveys at zero or very young (younger than ca. 2 Ma) ages were compiled and compared with corresponding (here this means spatially closest) ridge-axial depths. In total, 1108 measured zero-age crustal thicknesses were compiled from the literature (see citations above) where a link to a specific location along the ridge axis could be established or reasonably estimated with an average angular deviation from the ridge axis of 0.05° and a maximum of 0.4° (Fig. 6A; Data Repository Table S3 [see footnote 1]). The average crustal thickness as sampled by this very limited data set is 5.5 ± 1.0 km (1σ). These thickness estimates were then linked to various ridge properties, including axial depth and spreading rate. Figure 6B shows axial depth versus crustal thickness. There is a general tendency for thicker crust to be associated with shallower ridges and thinner crust to be associated with deeper ridge axes, but the relationship is not tight. This suggests that there are other contributions beyond simple thermal and crustal isostasy. It is important to emphasize that this data set represents a sampling of ∼0.5% of the total ridge axial length. It is likely that zero-age depth variations reflect similar contributions from variations in crustal thickness and mantle temperature along with spatially varying contributions from mantle dynamic support associated with radial flow under ridges, and radial and horizontal flow associated with hotspots (Detrick et al., 2002; among many other potential contributors). In order to help visualize this, the residuals about the best-fit line relating crustal thickness to axial depth shown in Figure 6B were used to define ±1σ and ±2σ deviations, corresponding to +293 m/–362 m and +794 m/–874 m, respectively, of the difference between axial depths predicted by this best-fit relationship and the actual axial depths. These depth anomalies are shown on Figure 6A. From this, it is clear that there is a tendency for the too-deep anomalies to cluster in the region of the Australia-Antarctic Discordance along the Southeast Indian Ridge, and along the Southwest Indian Ridge in regions with relative lows in dynamic topography (Forte et al., 2015). In contrast, many of the too-shallow anomalies are located in the vicinities of hotpots, including Iceland, the Azores, and the Galapagos (Fig. 6A), with relative highs in dynamic topography (Forte et al., 2015). Figure 6C shows the relationship between crustal thickness and spreading rate. For the most part, there are not enough data to warrant drawing any firm conclusion, but there does not appear to be a simple correlation between spreading rate and crustal thickness.
The depth data used for this analysis are based on SRTM30_Plus (Becker et al., 2009), which has a spatial resolution of 0.5 arcmin. The correlation of age and water depth requires effects associated with sediment loading to be removed to recover the sediment-load–corrected basement depth. The marine sediment thickness data set of Whittaker et al. (2013) synthesizes the current state of knowledge. This data set has a 5 arcmin spatial resolution but is limited to regions south of 80°N latitude and north of 70°S latitude, with some additional areas with no sediment thickness data, particularly in the western Pacific. For the present purposes, these sediment thickness data are accepted as provided by Whittaker et al. (2013).
Uncertainty in the sediment-load correction can be explicitly assessed for sediment thicknesses of 1300 m or less by comparing the actual IODP drill-site mean sediment column density and corresponding sediment-load correction relative to that predicted by Equation 8, as shown in Figure 10. The largest misfits are associated with actual mean column densities, which significantly exceed those predicted by Equation 5. These reflect the presence of dense sediments typically associated with relatively rare, sulfide-rich, black smoker deposits. Overall, the uncertainty derived from this comparison is +52 m/–109 m at ±2σ. This level of uncertainty is small compared with the dispersion of sediment-load–corrected depths and therefore is ignored in the age-depth analysis.
Age versus Sediment-Load–Corrected Depth Analysis of the Main Ocean Basins
The sediment-load–corrected depth versus age distribution of the main ocean basins based on 88,907 picks from the GSFML data is shown in Figure 11. These data have been filtered to also exclude (1) picks from oceanic plateaus, as delineated by Coffin and Eldholm (1994); (2) digitized and available data from large igneous province data sets (http://www-udc.ig.utexas.edu/external/plates/data/LIPS/Data/LIPS.2011.dat, downloaded 30 August 2018; Fig. 6A); (3) data within 200 km of trenches based on Bird (2003) (Fig. 6A); and (4) data from regions with sediment thickness of 0 m (Whittaker et al., 2013) for ages older than 20.2 Ma or thickness greater than 8000 m. The gap in the magnetic reversal data between 83 Ma and 120.6 Ma reflects the Cretaceous normal polarity superchron, during which there are no reversals and therefore no picks. In order to more easily visualize these data, the present analysis used estimates of statistical measures of the depth distribution of the picks as a function of age over the main ocean basins with or without hotspot-impacted ridges. As with the mid-ocean-ridge depth analysis, this assessment focused on the median, ±1σ, and ±2σ ranges and also tabulated the mean and modal depth at each of the 361 unique magnetic reversal pick ages represented in the GSFML database (Data Repository Table S1 [see footnote 1]). Because the depth distribution at each age is distinctly skewed with long tails toward shallower depths, the median and modal depths are more robust representations of the populations than the mean.
Best-fitting solutions using Equations 1 through 3 were derived from subsets of these 361 estimates filtered by the number of picks contributing to the depth distribution at each age. The number (N) of picks per pick age ranged from 1 to 3706, resulting in variably resolved estimates of the characteristic statistics. The best-fitting solutions will focus on picks with N ≥ 40 because this value provides a robust number of reversal ages (≥190) and reasonably tight distributions of depth with age.
The GSFML data set includes 133 picks along the present main ocean basin ridge axis (i.e., with age = 0 Ma), which allow a comparison to be made between measures based on sampling the full length, as described above, and a small subset based on these picks that is provided in Table 1. In general, the differences in the median, mode, and mean are small (∼200 m), particularly considering that 80 picks are from the East Pacific Rise, 42 are from the deeper-than-average Southwest Indian Ridge, and the remainder are from the Gulf Rise in the Gulf of California. This comparison provides some confidence that even this quite limited sampling set may yield reasonable estimates of the global distribution of depths as a function of age.
Figure 12 shows the median depths as a function of age for all 361 pick ages and emphasizes those with N ≥ 40 (filled black dots). Best-fit half-space (Eq. 1) and plate (Eq. 2) models constrained by fitting pick ages with N ≥ 40 are shown as well (Table 2). Although a half-space subsidence model can fit the median depth data with zr = –2741 m and B (product of all parameters in front of the square root of t) = 272, with a corresponding ΔT = 941 °C, this fit is worse as measured by the misfit, M, relative to the best-fit plate model. It is worth noting that if the maximum age of the half-space solution is limited to <84 Ma, then zr = –2388 m, and B = 354, with a corresponding ΔT = 1224 °C, assuming κ = 1 × 10−6 m2/s. These latter results are close to those of Parsons and Sclater (1977). The best-fit plate model has zr = –2309 m, zp = 122 km, and ΔT = 1304 °C. Based on the relative fit and closer match of ΔT to best estimates of mean mantle melting temperatures under mid-ocean ridges (Behn and Grove, 2015), the plate model provides a better overall fit to the GSFML pick data. These values can be compared with the best estimates derived by Hoggard et al. (2017), who reported zr = –2380 m, zp = 146 km, and ΔT = 1253 °C, also shown on Figure 12. It is clear that the median depths for ages older than the Cretaceous normal polarity superchron based on the GSFML data are shallower than the data used to constrain the Hoggard et al. (2017) fit, resulting in a warmer ΔT and significantly thinner zp.
The age versus depth data compiled by Hoggard et al. (2017) scatter around their best-fit solution by more than ±1 km. Given the factor of ∼40 increase in data employed in the current analysis, as well as the assessment of zero-age depth, it is possible to characterize the depth distribution at each pick age and thereby analyze the subsidence behavior of these different cohorts. This provides a basis for assessing the origins of this dispersion and for suggesting possible explanations. Figure 13A plots the modal and mean depths as functions of age for N ≥ 40 at each pick age. In general, the modal depths are reasonably tightly clustered around the best-fit plate solution (see Table 2). The area-weighted mean depths are more scattered at younger ages but have a much more restricted depth distribution at ages older than the Cretaceous normal polarity superchron. As shown by Figure 13A and Table 2, the best-fit solutions for the median, mean, and modes are quite similar and largely differ as a result of different zero-age ridge depths (Table 2). This is not a surprising result, and it suggests that the population of depths evolves reasonably coherently with age. Figure 13B shows the ±1σ deviations and best-fit plate solutions for these depth cohorts as a function of age. The −1σ deviations are tightly clustered about the best-fit solution at all ages, and, in fact, this is the fit with the smallest M of all main ocean basin plate-model solutions (Table 2). In contrast, the +1σ deviations show much greater dispersion. Again, this is not surprising, because the depth distributions are skewed, as with the mid-oceanic-ridge axial depths, with long tails toward the shallow end, and hence small changes in frequency can result in large changes in the depth corresponding with the +1σ deviation. Nonetheless, the best-fit plate-model solutions for +1σ, median, and −1σ deviations are also quite similar. As shown by the values in Table 2, the differences between the best-fit solutions are deeper, cooler, and thicker for the −1σ solution, and shallower, hotter, and thinner for the +1σ best-fit solution relative to the median. If picks associated with the Iceland-impacted North Mid-Atlantic Ridge and Afar-impacted Sheba Ridge are excluded, the same relations hold but with a tighter distribution (Table 2), as would be expected. Based on these latter best-fit plate-model solutions, this would suggest that there is a ±1σ range of melting temperatures along mid-oceanic ridges of ±∼80 °C, giving rise to approximately ±500 m of difference in zero-age depth and ∼±8 km variation in oceanic lithosphere thickness. This range of zr is comparable with the ±1σ range derived from the axial depth analysis presented above. Thus, one possible explanation for the dispersion of depths as a function of age is that each cohort starts with a different initial condition, i.e., zr, ΔT, and zp, and that these differences give rise to different age-depth histories. Figure 14 maps the spatial distribution of depth anomalies of 93,041 picks. Here, the categories were defined by differencing the sediment-load–corrected depth of each pick from the best-fit plate model of the median depths for those pick ages with N ≥ 40 and determining depth differences corresponding to ±1σ and ±2σ values. The categories correspond to +876 m/–386 m and +1938 m/–883 m for ±1σ and ±2σ, respectively. The results shown in Figure 14 can be interpreted as providing some support for the above stated hypothesis. The densely picked North Atlantic and region adjacent to the Azores are both characterized by broad swaths of anomalously shallow depths. Regions including the Australia-Antarctic Discordance, south end of the Central Mid-Atlantic Ridge, and area adjacent to the Southwest Indian Ridge are characterized by broad swaths of anomalously deep pick depths, extending out to oceanic crust younger than the Cretaceous normal polarity superchron. These would account for the coherence of the age-depth evolution linking effectively zero-age depths with older picks, as shown in Figure 13B. However, there are also many regions characterized by anomalous ridges bounded by regions of “normal” or even picks with anomalous depths of the opposite sign, as, for example, along the Pacific-Antarctic Ridge and East Pacific Rise. Plate kinematics indicate that the Southwest Indian Ridge has been characterized by slow spreading rates since ca. 49 Ma (Rowan and Rowley, 2016), and yet there are many picks that are anomalously shallow.
Figure 14 shows that there is a general correlation of anomalously shallow picks with regions impacted by both modern and ancient hotspot activity. Obviously, in the North Atlantic around Iceland; in the Central Atlantic around the Azores, Cape Verde, and Canaries; in the South Atlantic around the Walvis and Rio Grande; in the Indian Ocean around the Crozet, Mascarene, Seychelles, and Kerguelen Plateau; and in the Pacific, surrounding the Ontong-Java Plateau, Darwin Rise, and Southwest Pacific Swell, values are anomalously shallow. The shallowing of many of these regions likely postdated the production of the oceanic lithosphere and thus would reflect superposition of other processes, including thermal rejuvenation (Crough, 1978), within-plate magmatic thickening (Abrams et al., 1993), and dynamic topography, such as, for example, in the too-shallow Southwest Pacific Swell (Rowley et al., 2016) or too-deep Argentine Basin (Forte et al., 2015). From Figure 14, it appears that both ridge-dominated and off-axis components contributed about equally. Thus, while there do appear to be cohorts that subside coherently, as apparent at relatively young ages in Figures 13A and 13B, this is not the sole explanation of these patterns.
The hotspot-influenced Gulf of Aden and particularly the North Atlantic, where the emphasis on the latter reflects the large number of magnetic reversal picks in the North Atlantic basin, significantly increase the scatter of almost all measures of the age-depth relationship. Table 2 provides results of fitting in which data from these regions are excluded from the analysis. Figure 15 shows results for the median and ±1σ depths for this subset of the data. The median depths are tightly clustered about the best-fit plate model, which is also true of the −1σ depths. However, overall removing the contributions from these hotspot-impacted regions does not make a significant difference to the results, and thus the remaining analysis will use all of the main ocean basin results for purposes of comparison.
Magnetic Reversal Picks Compared to Age-Grid Age-Depth Relations
Several recent analyses of the oceanic age-depth relationship have used some version of the Müller et al. (1997, 2008) or Seton et al. (2012) age grids to derive the correlation between age and depth (Crosby and McKenzie, 2009; Crosby et al., 2006; van der Meer et al., 2017). Hoggard et al. (2017) assigned ages to sites for which they determined sediment-loading and crustal corrections using an augmented Müller et al. (2008) age-grid data set. This general approach is obviously appealing because the age grids are global and thus allow virtually any point on the oceanic lithosphere to be assigned an age. However, Figure 2 suggests that there is considerable spread in grid cell ages for the locations where the age is independently constrained from discrete magnetic reversal picks. It is thus worth comparing the age-depth relationships as derived from the explicitly dated magnetic reversal picks relative to those derived from several of these age grids. In order to make these as comparable as possible, the age grids were sampled over the same broad basins where both magnetic reversals and sediment thickness estimates exist (gray shaded regions of Fig. 3). Oceanic lithosphere located in regions without sediment thickness data were excluded. Data were also filtered to exclude regions corresponding with large igneous provinces and trenches (Fig. 6A), as for the pick data sets. The equatorial Atlantic and Red Sea, where there are no reversals, and ages are entirely extrapolated, were also excluded. All back-arc and marginal basins were excluded. The ocean basin outlines are crude in that they overlap regions with continental crust, but, when integrated with information in the age grids, particularly the distribution of nonoceanic lithosphere, it is possible to overlay these to segregate data to be included or not from the analysis.
The age-grid–derived age-depth data were binned at 0.1 m.y. intervals for the Müller et al. (2013, 2016) versions of the age grid and a 1.0 m.y. binning for the Müller et al. (2008b) version, as shown in Figure 16. For each version of the age grid, the sampling reflected the resolution of the age grid. The Müller et al. (2008, 2013) age grids have 6 arcmin resolution, whereas the Müller et al. (2016) age grid has 2 arcmin resolution. The corresponding depths were derived from SRTM_Plus 0.5 arcmin data (Becker et al., 2009), regridded at the same resolution and same grid registration as the corresponding age-grid data source. The 5-arcmin-resolution sediment thickness data (Whittaker et al., 2013) were not regridded, and instead the spatial correlation between the various grids was used to extract the sediment thickness data corresponding to each age-grid cell. In addition, the data were filtered to remove both oceanic plateaus and trenches. Figure 16 shows that the age grids of Müller et al. (2008) and Müller et al. (2013) yield effectively identical results, but that there has been a significant change in the implied age-depth relationship between Müller et al. (2008a, 2013) and Müller et al. (2016), such that there is now much closer correspondence with estimates derived from explicitly dated oceanic lithosphere using the GSFML data set. Figure 17 shows the spatial distribution of age differences between the ages derived from the picks and those from two of the age grids, Müller et al. (2013), shown in Figure 17A, and Müller et al. (2016), shown in Figure 17B.
Figure 16 highlights the fact that there are significant differences when using the gridded age data sets alone. For example, at very young ages (younger than 20 Ma), the GSFML-derived estimates correspond closely the with −1σ deviation of the Müller et al. (2016) age-depth distribution. For ages between ca. 0 to ca. 50 Ma, the GSFML-derived estimates and the median depths derived with all of the age grids diverge by >500 m. For ages from ca. 55 Ma to ca. 75 Ma, the age-depth relationships derived from the Müller et al. (2008, 2013, 2016) age grids converge with that derived using the GSFML data. At ages older than ca. 75 Ma, the Müller et al. (2008, 2013) age-depth relationship diverges sharply from both the best-fit GSFML curve and associated data, and also from the Müller et al. (2016) age-grid–derived age-depth relationship. The Müller et al. (2016) medians agree quite well with the GSFML relationship from ca. 50 Ma to ca. 135 Ma. At ages older than ca. 135 Ma, these values sharply diverge from the GSFML medians. It is perhaps worth noting here that the older age-grid data sets have been most heavily used to derive ages for the most recent age-depth analyses (Crosby et al., 2006; Hoggard et al., 2017; van der Meer et al., 2017), and yet they appear to yield the most divergent results when compared with the age-depth relationship reported here based on the GSFML data. It thus appears that these recent age-depth analyses may have been significantly impacted by the age-depth correlation derived from these older grids.
The differences between these various versions of the age-grid age-depth relationship relative to the picks can be understood by comparing the spatial patterns of differences in ages between the pick age and age-grid ages. The differences in age at the pick locations between the pick ages and the corresponding grid ages for all 99,090 nonduplicated picks, including back-arc basins based on the 2013 (Müller et al., 2013) and 2016 (Müller et al., 2016) age grids, are shown in Figures 17A and 17B, respectively. The comparison of the 2013 age-grid ages (Fig. 17A) relative to the pick ages shows that much of the Cretaceous record in the Pacific, where ∼50% of all older oceanic lithosphere exists, is characterized by blues, indicating that the age-grid ages are older than the pick ages. This will lead to a shallowing of the corresponding age-depth relationship (Fig. 16) if the interpolated and extrapolated ages in the intervening grid are similarly biased older. This contrasts with the 2016 age grid (Fig. 17B) over the same region, where there is no clear bias, and hence a closer match between the age-depth relations.
Focusing on the Müller et al. (2013) age-grid age-depth relationship, as this grid is comparable to the widely used 2008 age grid, the depth differences with respect to the best-fit plate model for the median depths of the GSFML data are shown in Figure 18. The same filtering criteria were applied to this map as to the other analyses. Thus, only grid cells with ages and sediment thickness data in the main ocean basins were included. As shown by the histogram, depth differences have a skewed frequency distribution, with the mode at −164 m, and with a long tail to the too-shallow end. Since a Mollweide projection is equal area, the counts are proportional to area. The regions that are shallow relative to the best-fit model are readily identified as mostly associated with hotspot-related magmatism and thus likely impacted by magmatic thickening of the crust (Abrams et al., 1993). If this is true, then it is likely that the intrusion of sills and eruption of off-axis flows would obscure the primary magnetic anomalies, leaving only those areas unaffected to retain their seafloor spreading–related magnetic reversals and primary age-depth relationships and hence contribute to the GSFML assessment. Those regions where reversals have been obscured would then require extrapolation or interpolation of ages, thereby contributing shallower depths to the age-grid–derived age-depth relationship, biasing the overall results. This then would account for the significant contrast between the GSFML age-depth distribution and the 2008 and 2013 age-grid–derived age-depth distributions. This argues in favor of the GSFML-derived relations, which are more likely capturing the real age-depth relationship. Obviously, filtering out all of the anomalous regions will yield an age-depth relationship comparable to that derived directly from the picks.
Comparison with Other Age-Depth Models
Figure 19 compares the present model for the age-depth relationship of the main ocean basins with models of Stein and Stein (1992), Johnson and Carlson (1992), Crosby and McKenzie (2009), van der Meer et al. (2017), and Hoggard et al. (2017). The van der Meer (2017) age-depth model is included because it has been used to assess ocean basin volume over time, which, given its disparity with all other models and data, would suggest caution in using their results for those or other purposes. The best-fit model derived here is in the middle of all the other models, perhaps suggesting something of a consensus approximation to the actual average age-depth relationship. The directly constrained age-depth data provide no support for the shallowing of depths extending through the Cretaceous normal polarity superchron to about ca. 160 Ma as derived by Crosby and McKenzie (2009). Of course, the majority of this interval corresponds with the Cretaceous normal polarity superchron, and therefore there are no magnetic anomaly picks with which to directly address that part of the history. However, depth-age data derived from deep-sea drilling results in the Atlantic, Indian, and Pacific Oceans from Johnson and Carlson (1992) are not compatible with the deflection modeled by Crosby and McKenzie (2009). In addition, the match between the Müller et al. (2016) age-grid–derived age-depth relationship and the best-fit GSFML-derived plate model limit any potential support for this interpretation.
One of the interesting observations based on a comparison of the best-fit model estimate of the zero-age depth with that derived by sampling along the mid-ocean-ridge system is that these are not at all similar. The median depth of the main ocean basin mid-ocean ridge is −3000 m, the modal depth is −2670 m, with a secondary mode at −2770 m, and the mean depth is −3044 m. The corresponding best-fit models depths are −2309 m, −2487 m, and −2220 m. Only the modal depths are close; the other two are different by ∼700 m or more. The disparity in depths for ages older than ca. 2–1 Ma, as shown in Figures 12 and 13A, is generally much smaller than ∼200 m, suggesting that this large misfit is only true for very young ages. This would be compatible with axial valleys being depressed relative to thermal isostasy along ridges and then recovered by ascent up the rift valley walls. The obvious counter is that along ridges with axial highs, such as the East Pacific Rise or Pacific-Antarctic Ridge, the median, modal, and mean depths are also deeper by ∼400 m compared to those predicted by plate or half-space fitting. This suggests that it is not simply the morphology, rift valley versus axial high, but a more general retardation near the ridge axis. This suggests that perhaps high-temperature hydrothermal heat loss, which is generally focused along the ridge axis, is locally reducing the thermal buoyancy and/or increasing mantle viscosity and flexural rigidity sufficiently to suppress ridge axial regions relative to that predicted by fits to data that are only a few million years older.
The axial depth distribution of mid-oceanic ridges is characterized globally by a reasonably tight clustering with a ±1σ range of ∼±500 m about the median depth. There is similar dispersion of depths at all ages that might suggest that these are related. However, an analysis of the spatial distribution of the depth dispersion suggests that inheritance from the ridge axial depth may only account for perhaps 50% of this dispersion. The other half appears to be related to superposition of other sources of anomalous topography over different regions and ages. Thus, there does not appear to be a simple correlation to axial depth. It is not clear why the roughness along the ridge should decrease with age, such that other sources begin to dominate at older ages. Presumably, some of this reflects the fact that sediment thickness, in general, increases with age, smoothing the seafloor roughness. The lack of detailed knowledge of small-scale sediment thicknesses results in a smoothly varying mapping of sediment thickness variations (Whittaker et al., 2013), which in turn results in considerable underestimation of subsediment seafloor roughness. If this is the case, then the larger, longer-wavelength, secondary sources, such as hotspot magmatism and thermal rejuvenation, together with other sources of dynamic topography (Fig. 18), dominate over the much shorter-wavelength variability characteristic of the mid-ocean-ridge axis (Fig. 3).
Klein and Langmuir (1987) provided a clear rationale as to why ridge axial depth and crustal thickness should covary. Their hypothesis was explicitly tested by Tolstoy et al. (1993) and by all subsequent experiments focused on determining zero-age crustal thickness. Gale et al (2014) summarized existing petrological and geochemical data derived from modern mid-oceanic-ridge basalts. A comparison of the relationships of Fe8.0/(Si8.0–41) (see Gale et al., 2014, their fig. 22) against axial depth and spreading rate with Figure 6 shows closely similar relationships, with some correlation to axial depth and none to spreading rate. They argued that variations in mantle potential temperature dominate, with secondary effects due to mantle heterogeneity and addition of plume influences that together largely explain their results. They derived a full range of ∼200 °C in mantle melting temperature to explain their results, which is broadly compatible with the ∼140 °C range based on estimates of the best-fit range at ±1σ of ΔT derived here. These in turn are reasonably matched by the melting model of Behn and Grove (2015).
The most recent age grid (Müller et al., 2016) represents a substantial improvement in the correlation between magnetic reversal pick and age-grid ages. The greater disparity between these in older versions of the age grid (Müller et al., 2008a, 2013; Seton et al., 2012), which have been used in many previous age-depth analyses, implies that mismapped ages may have significantly impacted those analyses. This suggests some caution is needed in addressing questions that are dependent upon details of the age grid.
Perhaps future generations of the age grid can employ one of approaches that have been adopted by Sandwell and colleagues in their generation of various global topography data sets (Becker et al., 2009). They have either used an odd versus even value for depths to indicate whether a given grid cell is constrained by ship or gravity interpolation, respectively, or generated a corresponding grid of SID values in which SID = 0 indicates gravity interpolation. For example, a simple solution would be to constrain subsequent versions of the age grid explicitly using the GSFML anomaly data and interpolate and extrapolate away from those constraints using plate-motion models. A useful component would be to tag each grid cell with an indicator of the source and quality of the age control, perhaps by adding a last digit to the ages to indicate whether it is directly constrained by magnetic anomaly picks, interpolated from nearby anomalies, or extrapolated from plate motions. In this way, future users could know explicitly the quality of the constraints on the age-grid reconstructions.
There has been much discussion of flattening of ocean bathymetry with increasing age relative to the half-space cooling model. A half-space cooling model fits the age-depth data well to ca. 83 Ma, beyond which the rate of subsidence clearly slows relative to the expectation of that model. Similarly, analyses of lithospheric thickness as a function of age (Steinberger and Becker, 2016) are compatible with a half-space model to ca. 110 Ma; older than that, the rate of thickness increase slows as predicted by the plate model. From an Occam’s razor perspective, the half-space model would represent the simplest hypothesis, as cooling and thickening are simply governed by physical properties of lithospheric materials and mantle potential temperature at the ridge axis (Jaupart et al., 2007). The plate model requires specification of a bottom boundary condition, the choice of which results in different heat-flow and age-depth histories at younger ages and thus more complexity. The failure of the half-space model at older ages and the compatibility of these various data with a plate model beg the questions as to the mechanism that is responsible for providing the extra heat at older ages. This is a classic problem that has been extensively discussed in the literature, including, but not limited to, onset of secondary convection (Richter and Parsons, 1975) and lithospheric reheating (Nagihara et al., 1996). Rowley et al. (2016) recently presented kinematic and mantle dynamic modeling results interpreted in terms of the East Pacific Rise axis overlying an actively upwelling convective roll in the mantle that has locked the longitudinal position of the East Pacific Rise in place on a time scale of at least 50 m.y. and likely longer. This is documented by long ridge residence times along the East Pacific Rise that are only exceeded in the southern Mid-Atlantic Ridge.
This recent work strongly supports mantle tractions at the base of the lithosphere as a major contributor to driving plate motions (Bird et al., 2008; Forte, 2011; Ghosh and Holt, 2012; Glišović and Forte, 2014; Rowley et al., 2013, 2016). This coupling requires advection of hotter asthenosphere away from ridges, thereby adding a plate-like heat source that will be most apparent at older ages. where the temperature contrast will be largest, and where the time scale for conductive heating from below will be longest. The resulting combination of decreasing the rate of cooling relative to that expected by half-space cooling and thermal erosion of the bases of old plates may be the most parsimonious explanation for the transition from half-space to plate-like cooling and thickening of plates. Such a mechanism that links active asthenospheric flow with basal tractions obviates the need for convective destabilization of old oceanic lithosphere (Parsons and McKenzie, 1978) and instead links the age-depth pattern directly to large-scale mantle convection processes driving plate motions. In detail, this leads to testable predictions that should allow this suggestion to be validated or falsified. Rather than a set the age of onset of the transition from half-space to plate cooling, as expected for growth of instabilities, this model predicts that this transition age should vary spatially depending in detail on the differential velocity histories of the asthenosphere and overlying lithosphere. Models capable of integrating this differential evolution in terms of the basal heat-flow boundary condition are not generally available, but the recent work of Glišović and Forte (2014, 2017) holds promise for such testing in the near future.
A revised age-depth relationship was derived here using the depth distribution of explicitly dated oceanic lithosphere using locations of ∼90,000 magnetic reversals compiled in the GSFML database. These data were fit with a plate model characterized by a zero-age ridge depth of −2309 m, asymptotic plate thickness of 122 km, and ΔT of 1304 °C. The ±1σ dispersion of data about the median, and the modal and area-weighted mean depths were also fit with plate-model solutions characterized by cooler, thicker, and deeper to warmer, thinner, and shallower conditions, as described by ΔT, zp, and zr among these different cohorts. This suggests that the starting zero-age depth is a contributor to the dispersion of depths at all ages. However, this is not the only factor, particularly at older ages, where superposition of other processes appears to contribute to the dispersion of depths with age. These factors make it difficult to identify the various potential perturbations in depths due to subsequent superposition of dynamic topography, lithospheric reheating, or thickening associated with non-plate-margin magmatism with amplitudes less than about ±500 m or more. Recent support for asthenospheric flow independent of plates giving rise to lithospheric tractions as a contributor to plate-driving forces implies that the advection of asthenosphere away from ridges may be a contributor to the transition from half-space to plate-like cooling.
I thank the three anonymous reviewers for the constructive critiquing of previous versions of this manuscript. I also thank Mark Hoggard for discussions and for helpful suggestions for improving various aspects of the presentation. This research was supported by the contributors to the Paleogeographic Atlas Project at The University of Chicago.