Glacial cycles drive variations in the production of oceanic crust

See allHide authors and affiliations

Science  13 Mar 2015:
Vol. 347, Issue 6227, pp. 1237-1240
DOI: 10.1126/science.1261508

Connecting orbit to the ocean floor

The amount of magma erupted at mid-ocean ridges can be modified by periodic ice ages that alter sea level. Crowley et al. analyzed high-resolution ocean depth data across the Australian-Antarctic ocean ridge (see the Perspective by Conrad). The results revealed 23-, 41-, and 100-thousand-year periodicity. These periods are similar to the well-known Milankovitch cycles associated with ice ages that are triggered by changes in Earth's orbit. Decreasing sea levels decrease the overlying pressure, thereby increasing the amount of erupted magma. The cyclic nature of glaciations and sea level creates a series of spaced topographic highs along the sea floor. Thus, Earth's atmosphere and mantle are coupled on a glacial time scale.

Science, this issue p. 1237; see also p. 1204


Glacial cycles redistribute water between oceans and continents, causing pressure changes in the upper mantle, with consequences for the melting of Earth’s interior. Using Plio-Pleistocene sea-level variations as a forcing function, theoretical models of mid-ocean ridge dynamics that include melt transport predict temporal variations in crustal thickness of hundreds of meters. New bathymetry from the Australian-Antarctic ridge shows statistically significant spectral energy near the Milankovitch periods of 23, 41, and 100 thousand years, which is consistent with model predictions. These results suggest that abyssal hills, one of the most common bathymetric features on Earth, record the magmatic response to changes in sea level. The models and data support a link between glacial cycles at the surface and mantle melting at depth, recorded in the bathymetric fabric of the sea floor.

The bathymetry of the sea floor has strikingly regular variations around intermediate and fast-spreading ocean ridges. Parallel to the ridge are long, linear features with quasi-regular spacing called abyssal hills (1). High-resolution mapping of the sea floor over the past few decades (24) has shown that these hills are among the most common topographic features of the planet, populating the sea floor over ~50,000 km of ridge length. Hypothesized models for these features include extensional faulting parallel to the ridge (3), variations in the magmatic budget of ridge volcanoes (5), and variation in mantle melting under ridges owing to sea-level change associated with glacial cycles (6). This latter model stems from the fact that glacial-interglacial variations transfer ~5 × 1019 kg of water between the oceans and the continents. This mass redistribution translates to sea-level variations of ~100 m and modifies the lithostatic pressure beneath the entire ocean. Because mantle melting beneath ridges is driven by depressurization, ocean ridge volcanism should respond to sea-level changes, potentially leading to changes in the thickness and elevation of ocean crust.

Plate spreading at mid-ocean ridges draws mantle flow upward beneath the ridge; rising parcels of mantle experience decreasing pressure and hence decreasing melting point, causing partial melting. Mantle upwelling rates are ~3 cm/year on average, whereas sea-level change during the last deglaciation was at a mean rate of 1 cm/year over 10 thousand years (ky). Because water has one third the density of rock, sea-level changes would modify the depressurization rate associated with upwelling by ±10%, with corresponding effects on the rate of melt production. Mantle upwelling rate scales with the mid-ocean ridge spreading rate, but the rate of sea-level change over the global mid-ocean ridge system is roughly uniform. On this basis, previous workers inferred that the relative effect of sea-level change should scale inversely with spreading rate, reaching a maximum at the slowest rates (6). An elaboration of this model with parameterized melt transport gave a similar scaling (7).

To test these qualitative inferences, we investigated the crustal response to sea-level change using a model that computes mantle flow, thermal structure, melting, and pathways of melt transport. The model is based on canonical statements of conservation of mass, momentum, and energy for partially molten mantle (8, 9) and has previously been used to simulate mid-ocean ridge dynamics with homogeneous (10) and heterogeneous (11) mantle composition. It predicts time scales of melt transport that are consistent with those estimated from 230Th disequilibium in young lavas (12). In the present work, the model is used to predict crustal thickness time series arising from changes in sea level (Fig. 1) (13).

Fig. 1 Simulated bathymetric relief driven by Plio-Pleistocene sea-level variation.

(A) Imposed sea-level variation [black, from (14)] and predicted bathymetric relief (color) for the past 1.25 million years from simulations at three half-spreading rates U0 and three permeability levels K0. Isostatic compensation is assumed to scale the amplitude of crustal thickness variation by 6/23 to give bathymetric relief. Permeability in the simulations is computed by applying K(x,z) = K0(ϕ/ϕ0)3 m2 to the porosity field ϕ(x,z), where ϕ0 = 0.01 is a reference porosity. Light blue, dark blue, and red lines correspond to log10K0 = –(13, 12.5, 12), respectively. (B) Power spectral density estimates for each time series, made by using the multitaper method with seven tapers. Axes are logarithmic.

A suite of nine model runs for three permeability scales and three spreading rates was driven over a 5-million-year period by using a Plio-Pleistocene sea-level reconstruction (14). Crustal curves from simulations with larger permeability and faster spreading rate contain relatively more high-frequency content than those of lower permeability and slower-spreading-rate runs (Fig. 1). Our model results contradict the previous scaling arguments (6, 7) in not showing a simple decrease in the sea-level effect on ridge magmatism with increasing spreading rate.

To better understand these numerical results, we carried out an analysis of leading-order processes using a reduced complexity model. This model provides a solution for crustal thickness response to changes in sea level, approximating the results of the full numerical model, but with greater transparency. Assuming that all melt produced by sea-level change is focused to the ridge axis, we obtain a magmatic flux in units of kilograms per year per meter along the ridge ofEmbedded Image (1)where ρwm is the density ratio of sea water to mantle rock, Π is the adiabatic productivity of upwelling mantle (in kilograms of melt per cubic meter of mantle per meter of upwelling), xl(z) is the half-width of the partially molten region beneath the mid-ocean ridge at a depth z, zm is the maximum depth of silicate melting beneath the ridge, and Embedded Image is the rate of sea-level change τ years before time t (13).

This formulation reveals why our numerical model results (Fig. 1) contradict earlier work (6, 7). Whereas earlier work noted that variations in crustal thickness are inversely proportional to spreading rate, CSL = MSL/(U0ρc), our model shows that mass flux is proportional to the width of the partially molten region beneath the ridge. This width can be expressed as xl(z) = U0R(z)/(4κ), where U0 is the half-spreading rate, κ is the thermal diffusivity, and R(z) accounts for depth-dependent influences on melting that are independent of spreading rate (fig. S1). The competing influences associated with the volume of mantle from which melt is extracted and the rate at which new crust is formed means that sensitivity to sea-level variation does not simply decrease with increasing spreading rate.

Instead, the magnitude of the crustal response depends on the time scale of sea-level forcing relative to the time required to deliver melt from depth to the surface. Melt delivery times τ are computed in the reduced model by using a one-dimensional melt column formulation and decrease with higher pemeability and faster spreading rate (13, 15). The same response occurs in the numerical model; in both cases, τ decreases with increasing spreading rate because the background melting rate, dynamic melt fraction, and permeability of the melting region all increase.

To quantify crustal response as a function of time scale, we use the amplitude ratio of crustal to sea-level variation, called admittance, computed at discrete frequencies by applying sinusoidal forcing. Admittance curves for both the numerical (Fig. 2A) and reduced (Fig. 2B) models show a distinct maximum that shifts toward higher frequencies and larger magnitudes with shorter τ. When the period of sea-level forcing is short relative to the characteristic transport time τm = τ(zm), additional melt produced at depth (falling sea-level phase) does not have time to reach the surface before a negative perturbation to melt production occurs (rising sea-level phase); positive and negative perturbations cancel, and crustal variation is small. When forcing periods are long relative to τm, melt perturbations reach the surface but are again small because melt production scales with the rate-of-change of sea level. Forcing periods near τm give maximum admittance because of a combination of large perturbation of melting rates and sufficient time to reach the surface (Fig. 2C). These results suggest that ridges are tuned according to melt-transport rates to respond most strongly to certain frequencies of sea-level variability.

Fig. 2 Crustal thickness admittance, computed for a sinusoidal variation in sea level with period TSL.

(A and B) Admittance curves derived from (A) numerical simulations and (B) the reduced model (13). (C) A plot of depth z versus the integrand from the reduced model of magma production owing to sea-level variation, Embedded Image (Eq. 1 and text following). The model is evaluated for U0 = 4 cm/year, K0 = 10−13 m2, and three values of sea-level oscillation period TSL.

The correspondence of the results from the numerical and reduced models provides a sound basis for investigating the potential effects of sea-level change on sea-floor bathymetry. Variations in melt production lead to variations in crustal thickness, and through isostatic compensation, such thickness variations should produce changes in bathymetry identifiable in high-resolution surveys. The prominent spectral peaks of late Pleistocene sea-level variation at the approximately 1/100 ky−1 ice age, 1/41 ky−1 obliquity, and 1/23 ky−1 precession frequencies (16) therefore translate into a prediction for a bathymetric response that depends on permeability and spreading rate.

Our model results suggest that the best chance to detect a sea-level response between 1/100 ky−1 to 1/20 ky−1 frequencies is at intermediate spreading ridges. Slow spreading ridges show little precession response, an obliquity response that is sensitive to uncertainties in permeability, and the effects of intense normal faulting. Such faulting causes rift valleys with larger relief than expected from sea level–induced melting variations. The sea-level signal should be less polluted by tectonic effects at fast spreading ridges but may have peak admittances at frequencies higher than 1/20 ky−1 that would obscure the responses at predicted frequencies. For example, the numerical simulation with the fastest spreading and highest pemeability has peak spectral energy at frequencies above precession (Fig. 1B).

At intermediate half-spreading rates of 3 cm/year, 40-ky periods lead to predicted bathymetric variations with a wavelength of 1200 m on each side of the ridge. Such fine-scale variations can be obscured in global topographic databases that grid data from multiple cruises and may have offsets in navigation or depth. To investigate the model predictions, a modern data set with uniform navigation and data reduction from a single survey is preferred. Such data are available for two areas of the Australian-Antarctica ridge that were surveyed by the icebreaker Araon of the Korean Polar Research Institute in 2011 and 2013 (Fig. 3).

Fig. 3 Bathymetry at a section of the Australian-Antarctic Ridge.

A region of consistent bathymetry is indicated between the black lines (top right) and is shown in profile (bottom left, blue) after converting off-axis distance to an estimate of time. Time is zero at the approximate ridge center. Also shown is bathymetry after filtering frequencies outside of 1/150 ky and 1/10 ky (green), and simulated bathymetry (black, for U0 = 3.3 cm/year and K0 = 10−13 m2). Spectral estimates (bottom right) are shown for the unfiltered bathymetry (blue) and model results (black), where the latter are offset upward by an order of magnitude for visual clarity. Data availability is uneven across the ridge, and spectral estimates are for the longer, southern flank. Unlike in Fig. 1B, spectral estimates are prewhitened in order to improve the detectability of spectral peaks (supplementary materials). Vertical dashed lines indicate frequencies associated with 100-ky late-Pleistocene ice ages, obliquity, and precession. Axes are logarithmic. Statistical significance is indicated by the black bar at the top right; spectral peaks rising further than the distance between the mean background continuum (corresponding to the black dot) and 95th percentile (top of black bar) are significant.

Analysis was undertaken by identifying a region whose abyssal hill variability is relatively undisturbed by localized anomalies, averaging off-axis variability into a single bathymetric line and converting off-axis distance into an estimate of elapsed time by using a plate motion solution (17). Spectral analysis of the associated bathymetry time series is performed by using the multitaper procedure (18) and shows spectral peaks that are significant at an ~95% confidence level near the predicted ice age, obliquity, and precession frequencies (Fig. 3). Although absolute ages are uncertain because we lack sea-floor magnetic reversal data, spectral analysis only requires constraining the relative passage of time. The 2σ uncertainties associated with relative Australian-Antarctic plate motion are ±4% (17), implying, for example, that the 1/41 ky−1 obliquity signal resides in a band from 1/39 ky−1 to 1/43 ky−1, a width that is smaller than our spectral resolution.

Another check on model-data consistency is to compare magnitudes of variability. Surface bathymetry will be roughly 6/23rds of the total variation in crustal thickness because of the relative density differences of crust-water and crust-mantle, assuming conditions of crustal isostasy. The closest match between simulation results and observations, in terms of the distribution of spectral energy, is achieved by specifying a permeability at 1% porosity of K0 = 10−13 m2 (Fig. 3). The standard deviation of the simulated bathymetry is 36 m, after multiplying crustal thickness by 6/23 and filtering (13). To minimize the contribution from non-sea-level–induced variations in the observed bathymetry, it is useful to filter frequencies outside of those between 1/150 ky−1 and 1/10 ky−1. The standard deviation of the filtered observations is 44 m, where the slightly larger value is consistent with changes in sea level being an important but not exclusive driver of changes in crustal thickness.

Analysis of bathymetry in another area of the Australian-Antarctic ridge 400 km to the southeast (fig. S2) shows a significant spectral peak at the obliquity frequency and indication of a peak near 1/100 ky−1, but no peak near the precession frequencies. Predicted and observed bathymetry is also similar, with standard deviations of 33 and 34 m, respectively, after accounting for fractional surface expression and filtering. Absence of a precession peak may result from spectral estimates being more sensitive to elapsed time errors at higher frequencies (19), where such errors may be introduced through extensional faulting or asymmetric spreading. Detection could also be obscured by the previously noted influence of faulting (3, 20) as well as off-axis volcanism or sediment infilling of abyssal troughs. Detection of significant spectral peaks at predicted frequencies at two locations of the Australian-Antarctic ridge nonetheless constitutes strong evidence for modulation of crustal production by variations in sea level.

Our numerical and analytical results show a complex relationship between spreading rate and amplitudes of crustal thickness variations associated with changes in sea level. Perturbations to the background melt production and delivery depend on the frequency content of the sea-level signal, as a result of the dynamics of magma transport. Reference mantle permeability and ridge-spreading rate are key controls on this frequency dependence. This result could be useful: The spreading rate can be accurately determined for a ridge, but parameters associated with magma dynamics are far less certain, such as the amplitude and scaling of permeability. Uncertainty associated with spectral estimates of bathymetry and sea-level estimates need to be better characterized, but together, these may provide a constraint on the admittance and hence dynamical parameters of a ridge.

Although results from the high-resolution bathymetry are promising, much remains to be done to further test the hypothesis advanced here. Crustal thickness is not an instantaneous response to melt delivery from the mantle but also reflects crustal processes that may introduce temporal and spatial averaging. Where long-lived magma chambers are present, for example, there may also be a crustal time-averaging that depends on spreading rate. In addition, faulting at all spreading rates is an observed and important phenomenon, and sea-floor bathymetry reflects the combined effects of magma output and crustal faulting (3, 20). Deconvolving the relative roles of such processes will be important. High-resolution surveys in targeted regions will provide the opportunity for a more complete and rigorous analysis than is presently possible.

Supplementary Materials

Figs. S1 to S3

Table S1

References (2129)

Data Files KR1 and KR2

References and Notes

  1. Materials and methods are available as supplementary materials on Science Online.
  2. Acknowledgments: The research leading to these results has received funding from the European Research Council (ERC) under the European Union’s Seventh Framework Programme (FP7/2007–2013)/ERC grant agreement 279925 and from the U.S. National Science Foundation under grant 1338832. It is a part of project PP13040 and PE14050 funded by the Korea Polar Research Institute. J.W.C. thanks J. Mitrovica, and R.F.K. thanks the Leverhulme Trust for additional support. Numerical models were run at Oxford’s Advanced Research Computing facility. Bathymetry data are included in the supplementary materials.
View Abstract

Stay Connected to Science

Navigate This Article