Research ArticleCLIMATOLOGY

Reconciling atmospheric CO2, weathering, and calcite compensation depth across the Cenozoic

See allHide authors and affiliations

Science Advances  22 Jan 2021:
Vol. 7, no. 4, eabd4876
DOI: 10.1126/sciadv.abd4876


The Cenozoic era (66 to 0 million years) is marked by long-term aberrations in carbon cycling and large climatic shifts, some of which challenge the current understanding of carbon cycle dynamics. Here, we investigate possible mechanisms responsible for the observed long-term trends by using a novel approach that features a full-fledged ocean carbonate chemistry model. Using a compilation of pCO2, pH, and calcite compensation depth (CCD) observational evidence and a suite of simulations, we reconcile long-term Cenozoic climate and CCD trends. We show that the CCD response was decoupled from changes in silicate and carbonate weathering rates, challenging the continental uplift hypothesis. The two dominant mechanisms for decoupling are shelf-basin carbonate burial fractionation combined with proliferation of pelagic calcifiers. The temperature effect on remineralization rates of marine organic matter also plays a critical role in controlling the carbon cycle dynamics, especially during the warmer periods of the Cenozoic.


Global climate and the global carbon cycle have undergone substantial changes since the Early Eocene [∼50 million years (Ma)], the period of peak Cenozoic deep ocean temperatures (13). Warming dominated during the Paleocene (~58 Ma), and the average temperature was steadily rising (high-latitude surface and deep ocean warmed by ~4°C) until it reached its maximum around 6 Ma later (Fig. 1A) (1, 2), as inferred by a negative excursion in stable oxygen isotope ratios (δ18O). The temperature increase during the Late Paleocene and Early Eocene (LPEE; 58 to 52 Ma) was accompanied by a long-term decrease in stable carbon isotope ratios (δ13C) in both surface and deep ocean (Fig. 1B), indicating a drop in the net organic carbon output from the ocean and atmosphere and reordering in carbon cycling (4, 5). Concomitant with the temperature rise was a global deepening of the calcite compensation depth (CCD) of at least 500 m (4, 5).

Fig. 1 Compilation of Cenozoic temperature, δ13C, pCO2, and calcite compensation depth (CCD) history for different ocean basins.

(A) The reconstructed deep ocean temperature across the Cenozoic is based on the δ18O data from (46). The prescribed temperature change in the model is based on the δ18O-temperature relationship described in (47): Tt = 16.9 − 4.0 × (δ18Ot − δ18Osw), where δ18Ot is the observed data and δ18Osw is the δ18O of the seawater at a given time in the past (B) δ13C compilation at different depths [benthic, (46); surface, (48); bulk, (49)] and (C) atmospheric CO2 reconstruction based on paleo-proxies from multiple sources: dataset 1 (50), dataset 2 (51), dataset 3 (52), dataset 4 (52), dataset 5 (53), dataset 6 (54), and expressed in parts per million by volume. (D) Pacific, Atlantic, and Indian Ocean CCD expressed in meters below sea level. Equatorial Pacific A compiled from (55); Equatorial Pacific B and Pacific A from (13); Pacific B, North Atlantic, South Athlantic, and Indian A from (56); and Indian B from (14). Hyperthermals are excluded from the original Indian Ocean CCD data (green line) (14).

Since the Cenozoic temperature maximum, the Earth system has experienced gradual cooling, which eventually resulted in major glaciation and ice cap formation at high latitudes, as reflected in δ18O (1, 3). Evidence suggests that the primary cause of the long-term cooling was a decreasing concentration of atmospheric CO2 (Fig. 1C) (6, 7). However, the reason behind the CO2 decline remains enigmatic as some studies suggest a generally constant CO2 degassing rate over the same time period while others suggest a decrease (810). Lower CO2 concentrations in the atmosphere result in colder temperatures and a weakened hydrological cycle, which led to decelerated weathering rates lowering ocean alkalinity (and carbonate saturation). This ultimately induces shoaling of the carbonate compensation depth. The CCD is often operationally defined as the depth at which sediments bear less than 5 or 10 weight % CaCO3 (11), and as such, its evolution (temporal and spatial) can be traced over the geologic past by inspecting the CaCO3 content of the sediment cores (1214).

Contrary to expectations, however, the deep-sea carbonate records indicate that the Pacific CCD deepened over 1 km, from ∼3.0 to 3.5 km in the earliest Eocene to about 4.6 km at present (12, 13), creating a carbon cycle conundrum. Similar trends are also observed in the Atlantic and Indian basins (Fig. 1D). As the carbonate preservation in the sediments is mostly a function of the carbonate saturation state of the ocean, the position of the CCD reflects changes in carbonate saturation (11), which is, in turn, controlled by the balance between sinks and sources of carbon, which govern atmospheric CO2. Thus, the variable position of the paleo-CCD over time carries a signal of the combined carbon cycle dynamics of the past. Tracing the CCD evolution across the Cenozoic and identifying mechanisms responsible for its fluctuations are therefore important in deconvolving past changes in atmospheric CO2, weathering, and deep-sea carbonate burial (5).

Here, we present the first modeling study that investigates carbon cycling across the Cenozoic with a full-fledged ocean carbonate chemistry model and predictive organic carbon burial that also incorporates a marine sediment component, crucial for the CCD reconstruction. We use an expanded version of the LOSCAR model [Long-term Ocean-Atmosphere Sediment CArbon cycle Reservoir; (15)], called LOSCAR-P, which combines the long-term phosphorus (P) and carbon cycle (16) coupled to a modified version of the GEOCARB III model (17) (see the Supplementary Materials).


To simulate and reconcile the observed pCO2 and δ13C trends during the Cenozoic (Fig. 2, A and B), we include a temperature-dependent mechanism capable of decoupling high production (due to high nutrient concentration) in the surface ocean from high organic carbon export. Temperature exerts two competing effects on organic carbon export and remineralization. On the one hand, higher temperatures elicit increased primary production and stronger carbon export (18). On the other hand, high temperatures instigate a higher degradation rate of organic matter sinking through the water column, albeit resulting in a lower amount of organic matter being exported to the ocean floor (19). Matsumoto (20) quantified the opposing effects that temperature has on the organic carbon pump and, consequently, on pCO2. The study showed that rising temperatures promote a boost in remineralization at a rate higher than the rise of the export production (the opposite is true for declining temperatures), essentially decoupling export production rates from the organic carbon pump. Hence, it is possible to increase the nutrient supply to surface waters via nutrient redistribution due to higher regeneration rates and higher deep-ocean nutrient inventories while simultaneously decreasing organic carbon burial in the deep ocean. We modeled this behavior by assuming a temperature-dependent Martin curve (see Materials and Methods and the Supplementary Materials). This temperature dependency is critical for reproducing the observed pCO2 and δ13C trends in the model (as illustrated in fig. S18, A and B, when T dependency is omitted).

Fig. 2 The reference model simulation (black line with closed black circles) and error envelope (orange area) plotted against the paleo-proxy data displayed and described in Fig. 1.

(A) Atmospheric pCO2 preditcted by the model in parts per million by volume. (B) Model predicted sediment bulk carbonate δ13C versus the observed sediment data [open blue circles, (49)]. ‰, per mil. (C) Model predicted evolution of the Pacific CCD against the Pacific CCD data range displayed in Fig. 1. The orange area represents model results for a 95% confidence interval of the Martin curve temperature dependence and the uncertainty associated with the GEOCARB weathering parameters (57). The Martin curve temperature dependence is as follows: b(t) = 0.062 × T(t) + 0.303( ± 2σ), where b(t) is the Martin curve attenuation coefficient at time t, and T is the median water temperature of the upper 500 m of the water column at time t. σ (0.16, calculated in this study) is the SD of residuals (disagreement between the linear regression model and the dataset), and 2σ is the 95% confidence interval.

On multimillion year time scales, the continental CaCO3 fluxes must be balanced by the burial of CaCO3 in marine sediments (11). Hence, when combined with the δ13C record and pCO2 proxies, the CCD provides an additional constraint on global weathering rates and ocean carbonate saturation. There have been a number of studies centered around investigation of the carbon cycle dynamics and climate across the Cenozoic, e.g., (8, 17, 2129), yet each of the studies is missing one or more critical components of the long-term ocean-atmosphere carbon cycle. For example, (8, 17, 2228) did not model the CCD; or they modeled the CCD without considering balancing the long-term ocean-atmosphere carbon cycle (29). The only previous study modeling the CCD across the entire Cenozoic, which reconciles the weathering fluxes and the CCD trends (albeit lacking a true sediment model), concluded that the changes in sea level and low weatherability of the Paleogene played a major role in controlling the CCD (21). In theory, such a scenario allows for decoupling of weathering rates from the CCD response and permits the multimillion year long periods in which terrestrial weathering steadily declines while the CCD deepens (e.g., during the past ∼50 Ma). According to our results, this mechanism may also have played a role in controlling the CCD across this time interval (see below).

To account for the effects of the eustatic sea-level fluctuations in our model, the relative shelf to deep CaCO3 burial was made dependent on the reconstructed sea-level curve (see the Supplementary Materials). Our results indicate that the change in shelf-basin carbonate deposition due to sea-level change is responsible for the overall CCD deepening across the Cenozoic but that the sea-level rise itself is likely not the main mechanism. We propose that latitudinal expansion of carbonate platforms and proliferation of calcifying organisms played a substantial role in controlling the CCD over long time scales.

The peak Cenozoic carbonate mass accumulation rate over the ocean shelf areas (including both shelves and slopes) is observed during the Paleocene and the Eocene (30, 31). The mean rate of carbonate accumulation over the shelves during this time period was about two to three times higher compared to the mean rate during the rest of the Cenozoic (30). A larger total Paleogene carbonate accumulation on the shelves is not only attributable to the longitudinal carbonate platform expansion, but it was also a result of the latitudinal expansion of carbonate platforms. While the longitudinal expansion is a result of a higher sea-level stand and a larger shelf area, the latitudinal expansion is a result of the Paleogene extreme warmth (30). In the modern ocean, carbonate accumulation is confined between ∼30° south and north of the equator, whereas during the early Cenozoic, the carbonate platforms expanded and occupied a wider band around the equator (∼45° north and south). A more uniform climate during warm periods, such as during the Late Paleocene and the Early Eocene, likely led to a lower CaCO3 saturation gradient between the poles and the equator, which promoted the latitudinal expansion of carbonate platforms (30, 31).

Coeval with high carbonate accumulation on the shelves was a proliferation (increase in relative species diversity/richness) of the calcifying taxa (coccolithophorids and globigerinina) during the Paleogene (32, 33). The Cenozoic diversity maximum coincides with the warm Paleocene-Eocene epoch, during which the species richness was two to three times higher than during the colder periods (33). As the planet cooled during the Late Eocene and into the Oligocene, the species richness exhibited a decline (33).

The combination of the Paleocene-Eocene sea-level rise and temperature rise (which peaked around 52 Ma) allowed calcareous organisms to migrate to shallow habitats, expanding the carbonate platforms both longitudinally and latitudinally (30, 31). An expanded shelf area during high sea stands (such as during the LPEE) combined with a larger number of CaCO3-producing organisms over the shelf areas (33) [which we postulate was due to increased nutrient supply caused by accelerated weathering and higher phosphate remineralization rates (34, 35)] can accommodate excess alkalinity inputs (caused by higher pCO2 and temperature). The above-described mechanism largely decouples the CCD from global weathering, resulting in CCD shoaling during the Early Eocene. As the planet cooled and seas regressed following the Eocene climatic optimum, proliferation of calcareous species and migration to the open ocean reversed the above trend and sent the CCD to greater depths (Fig. 2C). We modeled the proliferation and CaCO3 expansion effect by amplifying the shelf-deep CaCO3 partitioning factor (fsh; shelf deposition increase relative to modern) so that the carbonate burial over shelves was ∼3 times greater before the Eocene-Oligocene boundary, consistent with the evidence provided above (see the Supplementary Materials).

Our modeling framework demonstrates that the global CCD deepening over the Cenozoic does not necessarily require an increase in global chemical erosion rates, a mechanism often invoked to reconcile the observed trends across this time period (13, 36). Our results are also consistent with the current understanding of the long-term carbon cycling and feedbacks between weathering and atmospheric pCO2 (5). Declining atmospheric CO2 concentration (6), combined with decreasing (or roughly constant) CO2 degassing rates over the same time period (810), requires a slowdown in (or roughly constant) weathering rates because of the pCO2-silicate weathering feedback that stabilizes atmospheric CO2 on geological time scales (8). Other mechanisms, such as changes in weatherability, could also have played a role (see Materials and Methods).

The model predicted changes in silicate and carbonate weathering fluxes are also consistent with phosphorus accumulation rates, which appear to be declining between 50 and ∼25 Ma (consistent with the total P burial trend in our model; see fig. S17F) (37) with notable peaks at around 5 and 15 Ma. However, the long-term Cenozoic decline in phosphorus accumulation rate was interpreted as a contradiction (37), because the trend appeared to be in conflict with the Himalayan uplift theory (36), which states that the Himalayan orogenesis must have caused an increase in chemical erosion across the Cenozoic. If the weathering rates did not increase during this time interval (see fig. S17), for which there is evidence (38), then there is no inconsistency between P accumulation rates and the Cenozoic weathering trends, which further validates our model results.


Our results show that including complete ocean carbonate chemistry, prognostic organic carbon/phosphorus burial rates, and δ13C at various depths, combined with the ability to predict the CCD, is vital for constraining carbon cycle and climate dynamics over multimillion-year time scales. We suggest that metabolic effects caused by elevated temperatures of the Paleogene must be accounted for when the long-term carbon and phosphorus cycles are considered. Otherwise, the predicted organic carbon burial is inconsistent with records of δ13C, temperature, pCO2, and carbonate accumulation rates. Another fundamental implication from our model outcome is that the observed long-term deepening of the CCD across the Cenozoic does not require an increase in silicate and carbonate weathering fluxes. Our results show that the CCD trends were largely decoupled from silicate and carbonate weathering fluxes during the Cenozoic. We propose that the decoupling developed partially because of the increasing proportion of carbonate buried in the open ocean relative to the shelf carbonate deposition, as the sea level regressed concomitant with global cooling and the formation of continental ice sheets. In addition, combination and interplay of various environmental and physical factors (e.g., temperature, sea-level rise, CaCO3 saturation, nutrient supply and concentration, upwelling, and circulation changes) affected the proliferation and distribution of CaCO3-producing organisms during the Cenozoic.


All simulations were initiated with modern steady-state fluxes and boundary conditions (see table S1) and with an initial, modern atmospheric CO2 concentration of 300 parts per million by volume (ppmv), as provided by GEOCARB. Temperature change in the model is prescribed on the basis of the δ18O variations rather than being predicted on the basis of the Earth system sensitivity and atmospheric pCO2 concentration. The rationale behind this approach is that the ocean carbonate chemistry is temperature dependent (obeys highly predictive thermodynamic laws), and we have a continuous and relatively reliable temperature proxy (δ18O) spanning the entire studied interval. In addition, our new modeling framework follows the Arrhenius pattern for chemical reactions, which states that metabolic rates increase with temperature [e.g., (39)]. This is implemented in the model as the temperature-dependent remineralization via dynamic particulate organic carbon attenuation factor, which is derived from the empirical data (40).

The attenuation of the organic carbon flux in the water column in our model is a function of water depth and is represented by the Martin curveF=Fzp×(zzp)b(1)where zp is the depth of production zones (−100 m) and Fzp is the amount of organic carbon exported from the production zone. z is the ocean depth, negative downward (see the Supplementary Materials). b is the coefficient that describes the strength of the decay of the export flux, and its default value is 0.858 (41). However, the more recent data question the validity of the Martin curve and propose a wider range of values for b (42), which, among other factors, may depend on the ambient water temperature (fig. S5). The Martin curve temperature dependency in this study is modeled after the following equations that describes the dependency of the attenuation factor on temperature (40)b(t)=0.062×T(t)+0.303(±2σ)(2)where b(t) is the Martin curve attenuation coefficient at time t, and T is the median water temperature of the upper 500 m of the water column at time t. σ (0.16, calculated in this study) is the SD of residuals (disagreement between the linear regression model and the dataset), and 2σ is the 95% confidence interval. The original formulation of Eq. 2 does not include the SD portion, but it is calculated in this study to produce the model error envelope. To be consistent with the default, modern-day LOSCAR setup, where remineralization rate in the intermediate water (the upper 1000 m) is 78%, the b coefficient at time t = 0 (modern) needs to be 0.6737, which requires T of about 6°C. The abovementioned b value is also necessary to produce the modern organic carbon burial rate of 5 × 1012 mol year−1 (17), which satisfies the preindustrial steady state.

The cumulative organic carbon burial over a depth interval [i.e., intermediate (m) and deep ocean (d)] is calculated similar to the approach described by Bjerrum et al. (43); however, our study assumes that organic carbon exported from the surface ocean attenuates according to the Martin curveFbg,j=Fzp(t)z1z2β(z)×dA*dz×(zzp)bdz(3)Fcrain,j=Fzp(t)z1z2dA*dz×(zzp)bdz(4)Fcrem,m=Fzp(t)×(1A*×(Dtzp)b)Fbc,m(5)Fcrem,d=Fzp(t)(Fcrem,m+Fbc,m+Fbc,d)(6)where j ∈ {m, d}, and m and d are intermediate and deep ocean boxes, respectively. Fbg, Fcrain, and Fcrem are organic carbon burial, organic carbon rain, and organic carbon remineralization, respectively. β is the burial efficiency, which is a function of depth and oxygen (fig. S6). dA*dz is the derivative of the normalized hypsographic curve (fig. S6). For organic carbon fluxes (rain, remineralization, and burial) in the intermediate water, z1 is −1000 m (Dt; table S3) and z2 is −100 m. z1 and z2 for the deep ocean box are −4500 and −1000 m, respectively.

Weatherability estimate

It has been proposed that the silicate weathering flux may have remained relatively constant during the Cenozoic and that the changes in weatherability (kw) drove the atmospheric CO2 concentration (21). On the basis of our reference scenario, we carried out a simple calculation estimating the change in weatherability necessary to explain the maximum drop in CO2 over the Cenozoic. The maximum atmospheric pCO2 in our reference scenario occurs at 53 Ma (∼1960 ppmv), while the lowest is the preindustrial (t = 0) CO2 set to 300 ppmv. Let silicate weathering depend on kw, a suite of dimensionless GEOCARB parameters (kGEO for simplicity; see the Supplementary Materials for all parameters), and on atmospheric CO2FSi=kw×kGEO×(CO2(t)CO2(0))n(7)where n is the weathering feedback constant set to 0.6 (15). As the weatherability scenario assumes no changes in silicate weathering, we set FSi = const. For simplicity, we also keep kGEO constant and ignore carbonate weathering, and then the kw at time t reduces tokw(t)kw(0)=1(CO2(t)CO2(0))n(8)

Therefore, weatherability at t = 53 Ma is 0.32, and weatherability at t = 0 is 1. It follows that to explain our reference scenario, global weatherability would have to increase by about three times from Early Eocene to present, which is within the range of potential variability (44, 45).


Supplementary material for this article is available at

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 the editor D. Lea, L. Kump, and one anonymous reviewer for comments and suggestions that improved the manuscript. Funding: This work was supported by NSF grant OCE16-58023 to R.E.Z. Author contributions: N.K. wrote the manuscript and expanded and modified the LOSCAR model originally developed by R.E.Z. R.E.Z. provided guidance during model development and contributed to refining the manuscript text. 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.

Stay Connected to Science Advances

Navigate This Article