Research ArticleGEOCHEMISTRY

Oceanic crustal carbon cycle drives 26-million-year atmospheric carbon dioxide periodicities

See allHide authors and affiliations

Science Advances  14 Feb 2018:
Vol. 4, no. 2, eaaq0500
DOI: 10.1126/sciadv.aaq0500


Atmospheric carbon dioxide (CO2) data for the last 420 million years (My) show long-term fluctuations related to supercontinent cycles as well as shorter cycles at 26 to 32 My whose origin is unknown. Periodicities of 26 to 30 My occur in diverse geological phenomena including mass extinctions, flood basalt volcanism, ocean anoxic events, deposition of massive evaporites, sequence boundaries, and orogenic events and have previously been linked to an extraterrestrial mechanism. The vast oceanic crustal carbon reservoir is an alternative potential driving force of climate fluctuations at these time scales, with hydrothermal crustal carbon uptake occurring mostly in young crust with a strong dependence on ocean bottom water temperature. We combine a global plate model and oceanic paleo-age grids with estimates of paleo-ocean bottom water temperatures to track the evolution of the oceanic crustal carbon reservoir over the past 230 My. We show that seafloor spreading rates as well as the storage, subduction, and emission of oceanic crustal and mantle CO2 fluctuate with a period of 26 My. A connection with seafloor spreading rates and equivalent cycles in subduction zone rollback suggests that these periodicities are driven by the dynamics of subduction zone migration. The oceanic crust-mantle carbon cycle is thus a previously overlooked mechanism that connects plate tectonic pulsing with fluctuations in atmospheric carbon and surface environments.


The long-term carbon cycle, driven by successions of supercontinent amalgamation, stability, and dispersal, controls the concentration of CO2 in the atmosphere on time scales of hundreds of millions of years, regulating the Earth’s surface temperature and driving its climate between icehouse and greenhouse extremes (1, 2). Continental rifts were recently identified as one of the main drivers of elevating atmospheric CO2 during supercontinent fragmentation over long time scales [half-periods of 50 million years (My) and longer] (3). A recent, statistically robust synthesis of CO2 data for the last 420 My (4) reveals additional 26-to 32-My CO2 cycles based on our spectral analysis of these data, but they have not been previously analyzed and their origin is unknown. A time-series analysis of a diverse range of Mesozoic-Cenozoic geological events, particularly extinctions, was independently found to display a similar, dominant ~26- to 30-My periodicity, which was thought to reflect an extraterrestrial forcing, possibly cosmic showers (5), related to a ~30- to 42-My half-period (6) of the solar system’s oscillation about the plane of the Milky Way galaxy (7). In contrast, it has been argued that there is no evidence or theoretical justification for well-defined periodic bursts of asteroids or comets when the solar system crosses the galactic plane, and that the periodicity of the oscillation of the solar system about the galactic plane is too variable to produce narrow peaks (8). Instead, terrestrial causes are generally favored to explain periodicities of extinctions and related changes in surface environments (9, 10); however, a tectonic mechanism that may drive 26- to 30-My periodicities in surface environments and, in particular, atmospheric CO2 fluctuations has yet to be identified. Subduction and seafloor spreading together constitute the main plate tectonic driving forces via slab pull, slab suction, and ridge push (11), and plate boundary processes are known to play a large role in driving the global carbon cycle (12). These processes include degassing of CO2 along mid-ocean ridges (13), storage of CO2 in the ocean crust due to hydrothermal alteration of basalt (14, 15), storage of CO2 in the oceanic mantle due to serpentinization of mantle peridotites along the trench-parallel outer rise where the subducting plate bends and is faulted, and venting of arc volcanoes where a portion of subducted CO2 reenters the atmosphere (Fig. 1) (12).

Fig. 1 Schematic diagram showing oceanic crustal carbon cycle.

Values indicate maximum and minimum carbon fluxes for the oceanic crustal carbon reservoir since 230 Ma. Downward-pointing arrows indicate that the carbon is sequestered into the crust or mantle. Upward-pointing arrows indicate carbon flux into the atmosphere. The ocean crust acts as either a sink or source of atmospheric carbon depending on its changing capacity to hold CO2 through time, subject to fluctuations in the age-area distribution of ocean crust and changes in bottom water temperature. The percentage of subducted carbon degassing into the atmosphere is not well known, and we show end-member fluxes based on 35 and 65% of subducted carbon escaping into the atmosphere.

Oceanic sediments, including shelf and deep-sea carbonates, turbidites, and pelagic clay, contain variable amounts of carbon and potentially represent up to one-third of the total oceanic crustal and mantle carbon (12). However, the relative volumetric contributions of these sediment types to the carbon cycle, either today or through time, have not yet been quantified, although it is known that marine carbonates and organic-rich shales were more prevalent during the Cretaceous than today (16, 17). It is estimated that 45 to 65% of subducting carbon is released back into the atmosphere via arc volcanism and diffuse outgassing (12), whereas estimates from Dasgupta and Hirschmann (18) are somewhat smaller (~30 to 35%).

To date, the contributions of these processes to the global carbon cycle have only been analyzed for the present (12), implying that plate tectonics may exert a strong control on the carbon cycle. To test this hypothesis, the time dependence of these carbon fluxes must be reconstructed. This provides the context and motivation for our effort to quantify the role of the interplay of seafloor spreading and subduction in modulating atmospheric CO2 (4) since the breakup of Pangaea (Fig. 2A). Here, we focus on the crustal and mantle components that represent the bulk of the oceanic carbon storage.

Fig. 2 Atmospheric CO2 and tectonic model components.

(A) Atmospheric CO2 through time with 68% confidence intervals (4). (B) Mid-ocean ridge length (orange) and subduction zone length (light blue) through time based on the study of Müller et al. (19). (C) Seafloor spreading rates (red) and convergence rates (dark blue) from the study of Müller et al. (19). (D) Power spectra with SEs of the data shown in (A) to (C), following the same color scheme.

We use a recently published global tectonic plate model with continuously closing plate boundaries and associated paleo-oceanic age grids (19) as a tectonic framework for reconstructing the main oceanic crustal components of the global carbon cycle. Degassing of mid-ocean ridge basalt emits CO2 into the atmosphere (13), whereas the ocean crust stores carbon via hydrothermal alteration of basalt on young mid-ocean ridge flanks (14, 15). Oceanic crust is recognized as a significant sink for carbon, which accumulates at an estimated rate of 3.4 × 1012 mol C/year (14) as a consequence of long-term alteration of the crust by low-temperature basement fluids. This alteration causes a gradual loss of macro-porosity and a CO2 enrichment of the upper 300 m of the crust by precipitation of calcite cements in veins and voids from seawater-derived solutions (14, 20, 21). Using U-Pb dating of vein calcite, Coogan et al. (22) showed that >80% of the carbonate forms within 20 My of crustal accretion, with the remainder continuing to precipitate in crust aged 20 to 50 My. Precipitation of calcite is further controlled by ocean bottom seawater temperature, which determines reaction kinetics (20); this explains why the CO2 content of Cretaceous crust is significantly elevated relative to Late Cenozoic crust (20, 23).


Tectonic model

Our tectonic model (19) implies a doubling of mid-ocean ridge lengths from ~40,000 km since the Jurassic breakup of Pangaea around 200 million years ago (Ma) (Fig. 2B) to ~80,000 km during the Early Cretaceous period (145 to 100 Ma), subsequently dropping again to between ~60,000 and 70,000 km, whereas the global system of subduction zones is relatively stable through time, largely oscillating between 60,000 and 70,000 km. The increase in mid-ocean ridge length during the stepwise disintegration of the supercontinent Pangaea is extremely well constrained because of the detailed mapping of rifts and passive margins through time (24), with the geometry and length of mid-ocean ridges through time constrained by about 100,000 compiled marine magnetic anomaly identifications from all ocean basins (25). Destroyed back-arc basins and subducted mid-ocean ridges have been restored based on combined evidence from preserved magnetic lineations, fracture zones, geological data from overriding plates, including accreted terranes, and the rules of plate tectonics that demand continuity of all plate boundaries (19).

Fluctuating seafloor spreading rates display a slow long-term decline from ~8 to 5 cm/year (Fig. 2C), whereas plate convergence rates rise from ~5 to ~10 cm/year from the Jurassic to the Early Cretaceous, reflecting the immense increase in mid-ocean ridge length during this period. The tectonic cycles embedded in this plate model are highlighted by a power spectral analysis (Fig. 2D) of the main tectonic elements (Fig. 2, B and C), which reveals that changes in mid-ocean ridge and subduction zone lengths are dominated by relatively long 50- to 100-My cyclicities. In contrast, seafloor spreading rates, a key parameter driving fluctuations in the vast oceanic crustal carbon store, are characterized by cycles with a dominant period of 26 My (Fig. 2C), matching those seen in observed atmospheric CO2 cycles (Fig. 2A). Wright et al. (26) carried out a detailed analysis of how seafloor spreading rates based on our plate model in the Pacific Ocean depend on the geological time scale used, a parameter regarded as the main uncertainty involved in computing spreading rates. They found that the amplitudes of spreading rate cyclicities may depend on the time scale, but only at time scales of up to several million years, not tens of millions of years. Their analysis demonstrates that the time-dependent patterns of spreading rates at wavelengths analyzed here (that is, ~25- to 30-My cycles) are robust and time scale–independent.

Oceanic crustal carbon model

Our spectral analysis of tectonic model elements related to carbon cycle components through time suggests a profound connection between seafloor spreading cycles, atmospheric CO2, and surface environments at periods of ~26 My. Coogan and Dosso (27) recently speculated on a correlation between the time-dependent composition of altered oceanic crust and global surface environmental changes on similar time scales, because this corresponds to the time scale of the majority of chemical exchange on young mid-ocean ridge flanks.

Using all known inventories of crustal CO2 from ocean drilling sites (Fig. 3A) (20), we calculate the volumetric distribution of crustal carbon sequestered in present-day and now subducted seafloor back to 230 Ma (fig. S1) following the proposition that oceanic crustal CO2 content is driven by ocean bottom water temperature (Fig. 3B) (22, 23, 27) and age. Our robust, weighted log-linear relationship between crustal CO2 content, crustal age, and bottom water temperature (Fig. 3C) captures the rapid increase of crustal CO2 in 0- to 20-My-old ocean crust to about 2 weight % (wt %) CO2 when bottom water is relatively cold, but reaching values around 3 wt % when the bottom water temperature is as high as 15° to 20°C during hothouse climates (Fig. 3C). The maximum age of calcite precipitation in ocean crust was estimated to be 50 My by Coogan and Gillis (23). Our CO2–age–bottom water temperature relationship complies with this expectation for cool bottom water temperature scenarios but shows a slow continuing increase of crustal CO2 in aging ocean crust under hot bottom water conditions (15° to 20°C) (Fig. 3C) as a consequence of the sparse CO2 data (Fig. 2A). The Cretaceous Atlantic sites that display unusually high CO2 values (Fig. 3A) represent outliers, suggesting that frequently intercalated sediment in these cores may have led to precipitation of excess vein and cementing carbonate with a sedimentary origin within the oceanic crust (20).

Fig. 3 Oceanic crustal CO2 dependence on crustal age and bottom water temperature.

(A) Oceanic crustal CO2 dependence on crustal age and oceanic bottom water temperature (20). (B) Estimated oceanic bottom water temperature through time and geological periods (see Materials and Methods), categorized into five temperature regimes. (C) Relationship between oceanic crustal CO2, crustal age, and bottom water temperature, based on best-fit log-linear relationships (see Materials and Methods).

To apply this relationship to reconstructions of the ocean basins (19), we use paleo-ocean bottom water temperature estimates covering a time period from ~400 Ma to the present day because at 230 Ma, our starting time, ocean floor as old as ~180 My exists (fig. S2). We apply our CO2-age-temperature relationship to produce oceanic crustal CO2 grids at 1-My intervals (Fig. 4 and fig. S1) based on a combined set of oceanic paleo-age and paleo-ocean bottom water temperature grids (Fig. 4 and figs. S2 and S3) at 0.1° resolution. We compute the carbon flux due to mid-ocean ridge degassing by combining present-day flux estimates (13) with our estimates of mid-ocean ridge length (Fig. 2B), and we estimate the oceanic mantle carbon content by scaling published present-day estimates using our modeled changes in subduction zone length (Fig. 2B), following a uniformitarian approach.

Fig. 4 Global model grids.

Modeled oceanic crustal CO2 content through time (left), bottom water temperature at the time of crustal accretion for any given parcel of ocean crust through time (middle), and paleo-age of the ocean crust (right) (19). Note that the increased crustal CO2 content at 50 and 100 Ma, relative to today, is driven mainly by increased ocean bottom water temperatures during crustal formation. Crustal CO2 content decreases from 100 to 200 Ma because of the presence of large areas of crust formed during the Late Carboniferous–Early Permian ice age associated with cool ocean bottom water temperatures (Fig. 3B) in the Panthalassic Ocean.

Oceanic crustal carbon storage and fluxes

Mean upper oceanic crustal CO2 content, representative of the upper 300 m of ocean crust, is 1.85 wt % at present day and reaches a maximum of 2.25 wt % at 180 Ma in the Early Jurassic, reflecting the high bottom water temperatures during the latest Permian and Triassic hothouse climates (Fig. 3B)—crust formed during these periods dominated the seafloor at that time. Together with fluctuations in ocean floor area, and considering that the upper 300 m only represents two-thirds of the total oceanic crustal CO2 (20), this results in a total oceanic crustal store of 2000 to 2400 million tons of carbon, following the method of Kelemen and Manning (12) (Fig. 5A). The ocean crust became relatively deprived of its ability to hold CO2 during the Early Cretaceous because of the gradual subduction of old, CO2-rich Triassic ocean crust (fig. S3). The first time derivative (Fig. 5B) of the crustal CO2 storage curve (Fig. 5A), filtered using a 5-My cosine arch filter to suppress high-frequency noise, portrays how the capacity of the ocean crust to hold CO2 changes through time. These fluctuations are driven by changes in ridge length, seafloor spreading rates, convergence rates, and the resulting age-area distribution of CO2 storage. As a consequence, the ocean crust acts alternatively as a net source or a net sink of CO2, reflecting its time dependence on bottom water temperature and tectonic forcing. Subduction-induced flux of carbon into the atmosphere is dependent on the length of subduction zones, convergence rates, the CO2 content of the downgoing crust, and mantle and mass-balance calculations, suggesting that between 45 and 65% of the subducting CO2 is emitted into the atmosphere (12), although mass-balance calculations by others (18) have suggested somewhat lower values [see summary by Kelemen and Manning (12)]. We not only initially follow the simple assumption that half of the subducting carbon escapes into the atmosphere via volcanism and diffuse degassing (Fig. 5B) but also explore the effect on the total carbon budget analyzed here if this percentage is as low as 35% or as high as 65%. The degassing fraction of subducting carbon affects the coherence of our total modeled atmospheric carbon flux through time with atmospheric CO2 fluctuations derived from geological observations (Fig. 5D) (4). The carbon cycle components included in our model capture the main components of the oceanic crustal and mantle carbon cycle (Fig. 1), but the only two elements that fluctuate substantially through time are the crustal carbon subduction flux and crustal carbon storage (Fig. 5B).

Fig. 5 Model outputs.

(A) Oceanic crustal carbon content through time (blue) with SD (gray). (B) Carbon flux through time showing change in oceanic crustal carbon storage (blue), mid-ocean ridge degassing (purple), subduction flux into the atmosphere (red), and downgoing subduction carbon destined either for the deep mantle or the lithosphere (cyan) (12) assuming a simple 50-50 split. This ratio is not well known and may be anywhere between 1:3 and 3:1. Oceanic mantle carbon subduction flux is based on scaling today’s flux (12) with subduction zone length (Fig. 2B) through time (orange). (C) High-pass–filtered (cosine arch filter with 60-My width) atmospheric CO2 (black curve with gray error envelope) and modeled total carbon flux (green) including the relevant components shown in (B) (that is, all but the subducted carbon flux partitioned between the lithosphere and the convecting mantle); light red bars indicate bandwidths of the main periods of correlation between ~26-My period peaks in atmospheric CO2 and modeled carbon flux. (D) Spectral coherence of unfiltered atmospheric CO2 and modeled carbon flux time series (black, with 1 SD error bars; see text for discussion) peaks at 26 and 16 My. The three curves illustrate the coherence based on assuming that 35% (red), 50% (black), or 65% (blue) of subducted carbon degasses into the atmosphere. Note that the coherence at the 26-My period increases with decreasing subducted atmospheric carbon flux, whereas the remainder of the coherence plot is largely unaffected. (E) Power spectra of globally averaged trench migration speeds faster than 30 mm/year of eight global plate models with alternative reference frames [no net rotation model in light blue; models based on paleomagnetic data in dark blue, green, and dark yellow; and all other models based on hotspot tracks—see figure 4 in the study of Williams et al. (28) for details], revealing a dominant 26-My periodicity in trench migration.

Spectral analysis

The high-pass–filtered (see Materials and Methods) time series of our modeled carbon flux into the atmosphere displays a good overall correlation between peaks and troughs at periods between 25 and 30 My (Fig. 5C). There is spectral coherence (coherence > 0.5; see Materials and Methods) between these two signals (Fig. 5C) at 26 and 16 My. Exploring the effect of assuming that 65%, 50%, or 35% of the subducting carbon escapes into the atmosphere, we find that the coherence between modeled oceanic and atmospheric carbon emissions at a wavelength of 26 My varies between 0.62, 0.67, and 0.73. That is, the coherence at 26 My increases with a decreasing portion of subducting carbon escaping into the atmosphere (Fig. 5D).

Because coherence is a measure of the phase consistency between signals, two signals can be coherent while signal power is relatively low, and this is the case for the coherence peak at 16 My, as the observed CO2 data do not exhibit a high-amplitude rhythm at this wavelength (Fig. 5D). In contrast, the observed 26-My correlation coincides with power spectral peaks in the CO2 data and the carbon subduction flux, in turn correlating with fluctuations in seafloor spreading rates at the same periodicity (Fig. 2, C and D) and suggesting an origin related to changes in plate driving forces.


Considering that slab pull and slab suction have been identified as the dominant plate driving forces (11), we examine a set of plate kinematic models for periodicities in subduction processes. An analysis of subduction trench migration behavior since the Cretaceous period, including eight alternative absolute plate motion models (28), suggests that the driving force in the 26-My periodicity originates from the episodicity in slab rollback and advance. This episodicity is particularly well expressed in the global length of subduction zones rolling back at relatively fast speeds of more than 30 mm/year. Power spectra of this quantity for eight plate models (28) reveal a dominant 26-My cycle (Fig. 5E), illustrating a rhythm in subduction dynamics that reflects processes governing the episodicity of back-arc basin formation (29) and the cyclicity of slab penetration into the lower mantle (30). A recent review of cyclicities in large igneous province emplacement found a shift from a 64.5-My cyclicity to a 28- to 35-My cyclicity during the Jurassic period (31), similar to the cycles found in our oceanic crustal and mantle carbon model and associated seafloor spreading and trench migration cycles. Together with a ~30-My period of magnetic field reversals (32), this possibly reflects episodic release of plume material at the core-mantle boundary driven by a “top-down” episodicity of seafloor spreading rates and subduction dynamics. Deep Earth pulsing at a period of ~30 My is further supported by 30-My cycles in Mesozoic-Cenozoic volcanism found in the circum-Arctic region (33).

The connection between seafloor spreading rates, cyclicity in carbon subduction, and volcanic degassing along arcs described here suggests a compelling connection between tectonic processes, atmospheric CO2 cycles, and surface environmental conditions, as reflected in the 26-My cyclicity of ocean anoxic events, deposition of massive evaporites, and sequence boundaries (5). The coherence of our modeled total carbon flux with reconstructed atmospheric CO2 fluctuations at a wavelength of 26 My is maximized when assuming that only 35% of the subducting carbon is released back into the atmosphere (Fig. 5D), as opposed to 50 or 65%, leading to the suggestion that lower estimates for this fraction are more reasonable than higher-end estimates. The coupled plate-mantle processes driving these cycles on a global scale remain to be discovered, but their unraveling promises to have profound implications for modeling of the entire Earth system. Future generations of fully dynamic plate-mantle models may be capable of modeling the physical processes driving cyclicities in the interplay between plate driving forces, creation, and destruction of ocean crust as well as the equivalent periodicities in the global carbon cycle. Future extensions of the oceanic carbon cycle described here should consider the role of deep-sea sediments to store and release carbon through time, as well as a combination of oceanic and terrestrial carbon cycle components in a time-dependent global carbon model.


Plate tectonic parameters

Mid-ocean ridge length, subduction zone length, seafloor spreading rates, and convergence rates orthogonal to subduction zone azimuths through time were computed in 1-My intervals using the pyGPlates python library (, based on the plate model by Müller et al. (19).

Oceanic crustal CO2 and carbon through time

We converted crustal CO2 weight percent to volume of carbon by applying the methodology established by Kelemen and Manning (12) for the present-day oceanic crustal carbon reservoir to the geological past, based on our crustal paleo-CO2 grids (Figs. 2 and 5, A and B, and fig. S1).

Paleo-ocean bottom water temperatures

Using the oceanic crustal CO2 and temperature data of Gillis and Coogan (20), we established a robust, weighted log-linear relationship between crustal CO2 content, bottom water temperature at the time of its formation, and crustal age (19). We used the Generic Mapping Tools (GMT) trend2d function (34), which iteratively reweights input data to reduce the influence of outliers. To model paleo-oceanic crustal CO2 (fig. S1) as a function of crustal age and ocean bottom water temperature, we combined paleo-ocean bottom water temperature grids (fig. S3) with paleo-oceanic crustal age grids (fig. S2) (19). We assigned a bottom water paleo-temperature estimate to any given parcel of ocean crust for the time of its formation, using the Mg-temperature record of benthic foraminifera (35) for 65 to 0 Ma, the δ18O temperature record of benthic foraminifera for 110 to 70 Ma (36) and 120 to 115 Ma (37), and the δ18O temperature record of benthic foraminifera (38) for 160 to 130 Ma. For the Jurassic, we used a combination of δ18O temperature record of oysters, belemnites and trigoniids (39), and belemnites and demersal fishes (40). Because of the lack of data on bottom temperatures for Middle Devonian to Late Triassic oceans, our temperatures were estimates based on a comparison of published paleo-climate scenarios (41, 42) and calculations of seawater temperatures from δ18O of apatite fossils (4346) for these periods, following the well-established bottom water temperatures for similar climate modes in the Cenozoic to construct a complete estimated bottom water temperature record from 400 Ma to the present day (Fig. 3B).

Spectral data analysis

We applied a high-pass filter to our modeled carbon emission time series, as well as the observed CO2 data, using a cosine arch filter with a total width of 60 My to remove low-frequency parts of the signals, focusing on the correlation between relatively high-frequency fluctuations (tens of millions of years) between the two signals (Fig. 5C), because the processes that likely drive longer-term fluctuations (4) were not considered here. The two time series display a good overall visual correlation between peaks and troughs at periods between 25 and 30 My (Fig. 5C). To quantify the frequency-dependent correlation between modeled oceanic carbon emissions and atmospheric CO2 fluctuations, we computed the spectral coherence between the unfiltered signals with GMT spectrum1d (34), using the Welch’s method of ensemble averaging of multiple overlapping windows 128 My long and SE estimates (47). The spectral coherence between these two signals (Fig. 5D) ranging from 0 to 1 measured the normalized correlation between two power spectra, with values exceeding 0.5 considered coherent (47).


Supplementary material for this article is available at

fig. S1. Modeled oceanic crustal CO2 content through time.

fig. S2. Paleo-age of the ocean crust.

fig. S3. Modeled paleo-ocean bottom water temperature.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank J. Cannon for his development of pyGPlates workflows. We are grateful to F. A. Capitanio, an anonymous reviewer, and E. Kevin Furlong for constructive reviews. All figures were created using GMT software. Funding: This research was supported by the Australian Research Council grant DP130101946 and the Alfred P. Sloan Foundation through the Deep Carbon Observatory. Author contributions: R.D.M. and A.D. contributed equally to the conceptual framework, analysis, and manuscript preparation. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors. All data sets and digital grids are available for download at
View Abstract

Stay Connected to Science Advances

Navigate This Article