Research Article

Long-Term Sea-Level Fluctuations Driven by Ocean Basin Dynamics

See allHide authors and affiliations

Science  07 Mar 2008:
Vol. 319, Issue 5868, pp. 1357-1362
DOI: 10.1126/science.1151540


Earth's long-term sea-level history is characterized by widespread continental flooding in the Cretaceous period (∼145 to 65 million years ago), followed by gradual regression of inland seas. However, published estimates of the Late Cretaceous sea-level high differ by half an order of magnitude, from ∼40 to ∼250 meters above the present level. The low estimate is based on the stratigraphy of the New Jersey margin. By assimilating marine geophysical data into reconstructions of ancient ocean basins, we model a Late Cretaceous sea level that is 170 (85 to 270) meters higher than it is today. We use a mantle convection model to suggest that New Jersey subsided by 105 to 180 meters in the past 70 million years because of North America's westward passage over the subducted Farallon plate. This mechanism reconciles New Jersey margin–based sea-level estimates with ocean basin reconstructions.

Long-term (107 to 108 years) global (eustatic) sea-level fluctuations have been a driving force of biogeography, climate change, and organic evolution (1), yet they are poorly understood compared with short- to medium-term fluctuations (103 to 106 years) (2). This is due to the diversity of potential driving factors involved, including changes in mid-ocean ridge length, spreading rates, oceanic area, sedimentation, mantle convection, superplumes, large igneous province emplacement, and ice volume (2). It is well documented that many continents, in particular North America, Europe, and Africa, were inundated by shallow seas peaking in the Late Cretaceous 80 million years ago (Ma), followed by a gradual sea-level fall (3), but the associated magnitude of global sea-level change is very controversial.

Subsidence of the sea floor away from mid-ocean ridges is caused by the thickening and cooling of a thermal boundary layer at Earth's surface (4). Most published investigations of long-term sea-level change rely on the analysis of the age-area and associated depth-area distribution of presently preserved ocean floor to derive an estimate of ocean depth and volume change through time, as pioneered by Pitman (5). Rowley (6) reviewed such reconstructions and concluded that substantial ocean basin volume change and associated sea-level change have not occurred since the breakup of Pangaea. In contrast, Kominz (7) partially reconstructed past mid-ocean ridges (isochrons) for the last 80 million years and estimated the sea level to have been 230 +135/–185 m at 80 Ma, and Xu et al. (8), following a similar methodology, concluded that crustal production rates have decreased by 20 to 30% since 65 Ma, accompanied by a lowering of sea level between 125 and 250 m. Haq et al. (9) estimated the magnitude of long-term sea-level change from the present elevation of ancient marine deposits after reconstructing their subsidence history and adopting a sea-level tiepoint of 242 m above the present level for 80 Ma from Harrison (10), primarily on the basis of Kominz's (7) mid-ocean ridge volume analysis. In contrast, long-term sea-level estimates from stratigraphic sections of the New Jersey margin of North America yield a substantially lower 40-m sea-level high at 80 Ma (2), building on the work of Watts and co-workers (11, 12), who suggested that global sea-level signals can be isolated from the North American east coast sedimentary record.

We present a comprehensive reconstruction of the global age-area and depth-area distribution of ocean floor, including remnants of subducted crust, since the Early Cretaceous (140 Ma), to compute the effects of changes in crustal production, sediment thickness, and ocean-basin depth and area on sea-level fluctuations through time. In addition, we use a mantle convection model to test Miller et al.'s (2) sea-level curve and its premise that the New Jersey shelf has not been affected by tectonic processes other than thermal subsidence since the Cretaceous. This assumption allowed Miller et al. (2) to argue that the regional sea-level curve from New Jersey margin stratigraphy represents global (eustatic) sea-level variations.

Reconstructing vanished oceans. We establish the locations and geometry of mid-ocean ridges through time on the basis of marine magnetic anomaly identifications; geological information such as paleomagnetic data from terranes and microcontinents, especially in the Tethys Ocean (13); mid-oceanic ridge subduction events; and the rules of plate tectonics (14). On the basis of a global set of tectonic plate rotations (15, 16), we construct a set of refined sea-floor isochrons (Fig. 1) following the interpolation technique outlined by Müller et al. (15, 16) but including a multitude of additional data (see Supporting Online Material).

Fig. 1.

Age-area distribution of the ocean floor at (A) 140 Ma, (B) 100 Ma, (C) 50 Ma, and (D) the present day. See fig. S1 for paleo-age-area distribution in 10-My intervals. PAC, Pacific Plate; FAR, Farallon Plate; PHX, Phoenix Plate; IZA, Izanagi Plate; KU, Kula Plate; AUS, Australian Plate; PC, Proto-Caribbean; EB, Enderby Basin; AB, Amerasian Basin; C, Cocos Plate. Active mid-ocean ridges are represented as white lines, and subduction locations as black lines and triangle symbols.

In areas where one flank of a mid-ocean ridge system has been subducted, we assume spreading symmetry to reconstruct the subducted flank. This is a reasonable assumption given that globally, the maximum cumulative spreading asymmetry has been found to be less than 10% (17). The now entirely subducted ocean floor north of Africa, India, and Australia (the Meso-Tethys, Fig. 1) was reconstructed with the model of Heine et al. (18). Sea-floor isochrons for the Proto-Caribbean were constructed on the basis of the Early Cretaceous divergence history of the Americas, followed by the eastward insertion of the Farallon Plate into the Caribbean, leading to subduction of Proto-Caribbean ocean crust in the Late Cretaceous (19). Our Cretaceous reconstruction of the Panthalassa Ocean (Fig. 1) is based on the premise that the Pacific Plate formed at 170 Ma in the Late Jurassic as a triangle, originating from a triple junction between the Farallon, Phoenix, and Izanagi plates (20), which are reconstructed on the basis of preserved magnetic lineations and fracture zones (21, 22) and the assumption of spreading symmetry, with the exception of areas where ridge jumps are documented. We reconstruct the now entirely subducted Izanagi Plate and its conjugate Pacific Ocean floor (Fig. 1) following Whittaker et al. (23). Locations of plate boundaries within the Pacific/Panthalassa Ocean relative to surrounding continents are reconstructed with a plate circuit through Antarctica after 83.5 Ma and with Pacific and African plate motion relative to hotspots independently determined for earlier times (15, 16).

Plate boundary reconstructions between the Meso-Tethys and Panthalassa oceans are based on modeled motion vectors between the Izanagi, Phoenix, Australian, and Meso-Tethys plates. We use the combined model of finite rotations, sea-floor isochrons, mid-ocean ridge locations, and outlines of boundaries between continental and oceanic crust to construct a set of oceanic paleo-age grids from 140 Ma to the present (Fig. 1 and fig. S1), with accompanying age uncertainty grids (fig. S2). The mean error of our modeled age-area distribution of today's ocean floor is about 2.8 million years. We scale our error analysis for ocean basin reconstructions such that the mean reconstruction error of now subducted ocean floor is about 10 million years, namely, a fourfold increase in the combined uncertainties reflecting errors in mid-ocean ridge location and geometry.

Oceanic depth through time. An analysis of crustal production through time in 5-million-year intervals (fig. S3) illustrates that the most severe fluctuations in spreading-ridge crustal production are recorded in the Pacific Ocean. The largest drop in global crustal production, of about 50%, is recorded from 65 to 60 Ma, corresponding to the subduction of the Izanagi-Pacific spreading ridge (23) (Fig. 1 and fig. S1). This event, together with subsequent ridge subduction episodes along the Americas (fig. S1), leads to a mean Pacific sea-floor age difference of 35 million years (My) between 120 Ma and today (excluding back-arc basins), whereas all oceanic crust outside the Pacific Ocean taken together has experienced a total mean age difference of about 10 My (Fig. 2A). This demonstrates that dramatic changes in the age-area distribution of ocean floor have been driven primarily by cyclic changes in mid-ocean ridge creation, evolution, and destruction in the Pacific Ocean. It also shows that the use of the present distribution of the ocean-floor ages to obtain paleo-age distributions without reconstructing subducted sea floor (6, 24, 25) is not valid.

Fig. 2.

(A) Mean oceanic crustal age of the Pacific versus Indian, Atlantic, and Arctic ocean basins including all back-arc basins. The Pacific Ocean mean age increased from 30 to 65 My between 120 Ma and the present, whereas the mean age of all other ocean basins changed relatively little in comparison (a maximum of 10 My), challenging the notion that the Pacific Ocean mean age has remained constant through time (25). (B) Mean oceanic crustal age in 1-My intervals from 140 Ma to the present with error envelope computed from the age error grids (see fig. S2 for age error grids in 10-My intervals). (C) Mean oceanic basement depth since 140 Ma with uncertainties derived from age errors based on the GHD-1 age-depth model (26).

The maximum global difference in mean oceanic crustal age between 120 Ma and the present is about 26 My, ranging from 38.4 ± 4.1 My at 120 Ma to 64.1 ±1.4 My at present (Fig. 2B). This is in contrast to the 6-My difference predicted by Cogne et al. (25), who erroneously assumed that the Pacific age-area distribution has remained constant since the breakup of Pangaea. Using a published age-depth relation (GDH-1) (26), we compute the depth-area distribution of the ocean basins (Fig. 2C). We choose GDH-1 for converting age to depth because this relation is based on sediment-corrected depths without excluding data from hotspot swells and sea-mounts. Therefore, GDH-1 provides a good average fit to sediment-unloaded oceanic basement depths (26) and is preferable for predicting the average oceanic basement depth through time, including thermally rejuvenated lithosphere, as compared to models [e.g., (27, 28)] that reflect ocean-depth changes related to plate aging through time only.

Three additional factors play an important role in controlling global ocean basin depth through time, namely, the generation of oceanic large igneous provinces (LIPs), oceanic sedimentation (2), and changes in oceanic crustal area. We use Schubert and Sandwell's (29) method to determine the average elevation of individual major oceanic plateaus (table S1) relative to the surrounding sea floor, based on the difference between the modal depths within two polygons, one outlining the perimeter of a given plateau from a LIP database (30) and another polygon including surrounding ocean floor. Plateaus are added to our basement depth grids at their emplacement time based on a revised compilation of eruption ages (31). We estimate the uncertainty in LIPs-related mean oceanic depth changes to be on the order of up to ±5 m. An additional effect not considered here is the potential uplift of the sea floor some time before LIP eruption when the body of hot material (plume head) arrives in the upper mantle, followed by a collapse during eruption. From estimates for the size (about 500 km in diameter) and density anomaly (about 1%) of the plume head, we infer about a 2-m sea-level rise before eruption of a LIP, with 5 m probably an upper limit (corresponding to ∼700 km in diameter).

As the ocean floor ages, its sediment cover thickens, but abyssal sediment thickness is also latitude dependent, as illustrated in a polynomial surface fit of global sediment thickness (32) as a function of oceanic crustal age and latitude (fig. S4). We use this model to estimate the distribution of total oceanic sediment thickness and its isostatic compensation (33) through time, associated with an estimated mean sediment thickness uncertainty of about ±15 m for the Tertiary and ±25 m for the Cretaceous; however, this uncertainty is not well constrained. We create a set of paleo-bathymetry maps by adding major oceanic plateaus (table S1) and sediment thickness to our reconstructed basement depth maps (Fig. 3). These maps allow us to compute oceanic crustal area and mean depth through time.

Fig. 3.

Reconstructed bathymetry through time by combining basement depths derived from paleo-age grids (Fig. 1 and fig. S1) and the GDH-1 age-depth model (26) with estimates of sediment thickness through time (figs. S4 and S5) and inclusion of all major oceanic plateaus (table S1).

The emplacement of LIPs alone has resulted in a total sea-level rise of nearly 100 m since 140 Ma, with major pulses at about 120, 110, 90, 60, and 40 Ma (Fig. 4 and fig. S6). The sediment thickness effect on sea level through time is estimated to be an 18-m drop in sea level from 140 to 82 Ma, followed by a sea-level rise of 61 m (fig. S6). Consequently, the modeled increase in mean basement depth of 350 m from 120 Ma to the present (Fig. 2C) is alleviated by oceanic plateau emplacement and sediment thickness fluctuations to a net depth increase of about 260 (+100, –85) m since 118 Ma (fig. S6).

Fig. 4.

(A) Predicted relative sea-level change (red) and error envelope (reflecting basement age-depth errors) (light pink) since 140 Ma based on changes in the average depth and area of the ocean basins using the GDH-1 plate model (26) age-depth conversion (Database S1) compared with the sea-level curves by Kominz (7) (magenta), Haq et al. (9) (dark blue), adjusted to the Berggren et al. time scale (2), Haq and Al-Qahtani (43) (dark green), Miller et al. (2)(black), Watts and Steckler (11) (pink), Watts and Thorne (12) (light green), Pitman (5) (light blue), and the two curves by Xu et al. (8) (light and dark orange). The Haq et al. (9), Haq and Al-Qahtani (43), and Miller et al. (2) curves have been smoothed by a 10-My cosine-arch low-pass filter (in bold), representing long-term sea-level change. The original curves from Haq et al. (9) and from Miller et al. (2) are shown as thin lines. (B) Sea-level effect of large igneous provinces (LIPs) (100 m cumulatively). (C) Total crustal production through time since 140 Ma with major tectonic events highlighted in blue. (D) Miller et al. (2) sea-level curve, uncorrected (black) and corrected for post–70 Ma modeled tectonic subsidence of the New Jersey margin based on the s20rts seismic tomography model (39) for the pure backward advection (p.b.a) (purple) and modified backward advection (m.b.a) (light blue) compared with long-term sea-level variation from our computation of the volume of the ocean basins through time (red), as well as with the models by Haq and co-workers (9, 43) and Watts and co-workers (11, 12).

Long-term sea-level change. Global sea-level change during the past 33 My was dominated by large ice-sheet growth and decay (glacio-eustacy) (2). As our sea-level curve does not include the effects of glacial ice, we shift it by the 54 m that would be added to sea level if all present ice sheets melted (2) (Fig. 4), resulting in a Cretaceous high of 304 m at 119 Ma. Changes in sea-level Δsl as they would appear to an observer on land are related to changes in mean ocean basement depth Δh via isostatic compensation (34, 35): Δsl = (ρm ρw)/ρm × Δhm is mantle density, 3300 kg m–3; ρw is water density), yielding an apparent sea-level drop of about 235 m since a maximum at 118 Ma (fig. S6). However, this estimate, based purely on mean oceanic depth change, does not account for changes in oceanic area. Our reconstructions suggest that since 140 Ma, the global oceanic area has decreased by 3.3%, corresponding to a loss of about 1 million km2. This is due to continental crustal stretching during rifting, preceding continental breakup and ocean basin formation. Assuming no change in global oceanic water volume, the modeled decrease in oceanic area through time partly counteracts the aging and associated deepening of the global oceans after 80 Ma and attenuates the modeled Late Cretaceous sea-level high to about 170 m (fig. S6). Continental margin thermal subsidence following continental stretching is compensated mostly by sedimentation that keeps continental shelves relatively close to sea level (36), resulting in small changes in ocean basin depth from passive margin aging, which is not considered here.

A sea-level curve (Fig. 4) that takes into account both changes in depth and area of the global ocean basins reduces the modeled sea-level high at 118 Ma to about 150 m (Fig. 4). We obtain a total sea-level rise of about 70 m between 140 and 110 Ma; this is a better match with most Cretaceous continental flooding estimates [summarized by (2)] than our alternative sea-level curve not considering changes in oceanic area, which implies a sea-level rise of less than 30 m for the same time period (fig. S6). We model a long-term sea-level fall measuring about 100 m from 80 to 30 Ma, driven by the aging and deepening of the global ocean floor during this time (Fig. 2 and fig. S6). Our Late Cretaceous sea-level maximum of 170 m (90 to 265 m) occurs at 82 Ma, close to the benchmark time of 80 Ma (2, 7), supporting estimates from continental inundation as summarized by Miller (2) and close to Bond's (37) estimate of 200 m as well as that by Watts and co-workers (11, 12), but differing considerably from Miller et al.'s (2) best sea-level estimate of ∼40 m at 80 Ma (Fig. 4), based on New Jersey margin stratigraphy, raising the question of what causes this substantial disagreement.

New Jersey margin subsidence. We propose that the large 130-m discrepancy at 80 Ma between our sea-level curve and that from Miller et al. (2) (Fig. 4) reflects the progressive, mantle-driven tectonic subsidence of the New Jersey margin over the negatively buoyant Farallon slab imaged in the lower mantle (Fig. 5), while global sea level was falling after the Late Cretaceous highstand. Our reconstructions of subduction trench positions through time (Fig. 5A) suggest that the east coast of North America has moved over the Farallon slab (Fig. 5B) during the past 60 My. We test this hypothesis using regional results from a well-established global mantle convection modeling approach (38) to compute time-dependent dynamic surface topography by advecting present mantle density anomalies, derived from seismic mantle tomography, backward in the mantle flow field, which in turn is computed from mantle density anomalies and given surface plate velocities. We model the flow field of the mantle through time using a spectral method based on spherical harmonic expansion of surface plate velocities and internal density heterogeneities at each depth level (38). Dynamic topography is computed beneath air with a free upper boundary and the same viscosity structure as in (38) with a high-viscosity lithosphere (2.4 × 1022 Pa·s). We compare regional time-dependent dynamic topography since 70 Ma for the New Jersey coast of North America on the basis of mantle density anomalies derived for the present day from three commonly used seismic tomography models, s20rts, smean, and ngrand (3941), using a conversion factor of 0.25 from relative seismic velocity to relative density variations below 220 km and disregarding velocity variations above 220 km. In addition to the pure backward advection previously used (38), we also employ a modified procedure to reconstruct past mantle density anomalies. It differs in two ways: (i) Backward-advected upwellings are always continued up to 220 km at each time step, and thus the restoration of past mantle upwellings is improved; (ii) backward-advected downwellings are removed in the uppermost 220 km, and thus the effect of reconstructed ocean floor, which is already explicitly included, is not counted twice.

Fig. 5.

(A) Reconstructed subduction zone locations in present-day coordinates with coastlines in black and present-day plate boundaries in gray with ages of trench positions color coded. Legend numbers indicate reconstructed ages (in million years) in an absolute mantle reference frame (see text). Horizontal bold blue line and gray circles outline location of mantle seismic tomography profile below. (B) West-east mantle seismic tomography cross-section (39) centered on the New Jersey margin [see (A) for profile location], imaging the subducted Farallon slab in the lower mantle underneath the margin. Our kinematic model predicts North America's east coast to be underlain by subducted slabs about 80 to 140 million years old, coinciding with the location of the imaged Farallon slab (3941). (C) Predicted dynamic topography of the New Jersey margin from 70 Ma to the present based on the s20rts, smean, and ngrand seismic tomography models for the pure (p.b.a) and modified (m.b.a) backward advection cases (3941) (table S2). (D) Modeled dynamic topography of eastern North America at 70, 50, 25 Ma and at present based on the s20rts (p.b.a) seismic tomography model, in the reference frame of the North American plate. The black dot outlines the location of the New Jersey margin.

All models result in New Jersey tectonic subsidence during the past 70 million years, with amplitudes of 105 to 180 m (s20rts), 330 to 360 m (smean), and 255 to 385 m (ngrand), illustrating differences in the seismic tomography models and resulting model uncertainties (Fig. 5C). Mantle flow models tend to overpredict dynamic topography compared to estimates based on observed topography (42). However, all models agree with the qualitative consideration of the New Jersey margin overriding the Farallon slab, with s20rts providing the most reasonable estimates. Without taking dynamic topography into consideration, Miller et al.'s (2) model leads to an underestimate of the post–Late Cretaceous global sea-level drop. Therefore, we add the s20rts-based modeled New Jersey margin mantle-driven subsidence to Miller et al.'s (2) sea-level curve (Fig. 4D) and assume that before 70 Ma, prior to overriding the Farallon slab, the New Jersey margin was not affected by changes in dynamic topography. Miller et al.'s (2) corrected sea-level curves based on the s20rts seismic tomography models are both largely within the error envelope of our Cretaceous sea-level curve based on ocean basin volume analysis (Fig. 4D), implying a Late Cretaceous sea-level high between 142 and 217 m (based on s20rts) at 80 Ma. In particular, the modified backward advection s20rts–based tectonic subsidence correction brings Miller et al.'s (2) curve in good alignment with the Cretaceous portion of our ocean basin volume–derived curve (Fig. 4D), whereas the pure backward advection (38) s20rts-based corrected Miller et al. (2) curve is similar in amplitude to Haq et al.'s curve (9, 43).

Our results support the idea that mantle convection–driven dynamic topography has played an important role in the subsidence history of the New Jersey margin, predicting a total post–70 Ma tectonic subsidence between 105 and 180 m at average rates of 1.5 to 2.5 m/My. For times after 110 Ma, our ocean basin volume–derived sea-level curve yields lower estimates than those by Haq and co-workers (9, 43). The discrepancies may reflect global effects of dynamic topography, not modeled here, because the mantle backward advection that we use cannot realistically restore all global mantle density anomalies in the past. However, a combination of a global ocean basin volume analysis with modern geodynamic models provides a powerful tool for discriminating eustasy from regional, time-dependent sea-level variations caused by mantle convection. Our global sea-level curve calibration provides an improved framework for sequence stratigraphy, resource exploration, and models for long-term climate change.

Supporting Online Material

Materials and Methods

Figs. S1 to S6

Tables S1 and S2

Database S1

References and Notes

View Abstract

Navigate This Article