The 405 kyr and 2.4 Myr eccentricity components in Cenozoic carbon isotope records
- 1Department of Earth Sciences, Faculty of Geosciences, Utrecht University, Princetonlaan 8a, 3584 CB, Utrecht, the Netherlands
- 2Department of Oceanography, University of Hawai`i at Mānoa, 1000 Pope Road, HI 96822, Honolulu, USA
Correspondence: Ilja J. Kocken (email@example.com)
Cenozoic stable carbon (δ13C) and oxygen (δ18O) isotope ratios of deep-sea foraminiferal calcite co-vary with the 405 kyr eccentricity cycle, suggesting a link between orbital forcing, the climate system, and the carbon cycle. Variations in δ18O are partly forced by ice-volume changes that have mostly occurred since the Oligocene. The cyclic δ13C–δ18O co-variation is found in both ice-free and glaciated climate states, however. Consequently, there should be a mechanism that forces the δ13C cycles independently of ice dynamics. In search of this mechanism, we simulate the response of several key components of the carbon cycle to orbital forcing in the Long-term Ocean-atmosphere-Sediment CArbon cycle Reservoir model (LOSCAR). We force the model by changing the burial of organic carbon in the ocean with various astronomical solutions and noise and study the response of the main carbon cycle tracers. Consistent with previous work, the simulations reveal that low-frequency oscillations in the forcing are preferentially amplified relative to higher frequencies. However, while oceanic δ13C mainly varies with a 405 kyr period in the model, the dynamics of dissolved inorganic carbon in the oceans and of atmospheric CO2 are dominated by the 2.4 Myr cycle of eccentricity. This implies that the total ocean and atmosphere carbon inventory is strongly influenced by carbon cycle variability that exceeds the timescale of the 405 kyr period (such as silicate weathering). To test the applicability of the model results, we assemble a long (∼22 Myr) δ13C and δ18O composite record spanning the Eocene to Miocene (34–12 Ma) and perform spectral analysis to assess the presence of the 2.4 Myr cycle. We find that, while the 2.4 Myr cycle appears to be overshadowed by long-term changes in the composite record, it is present as an amplitude modulator of the 405 and 100 kyr eccentricity cycles.
During the Cenozoic (past 66 Myr) the Earth's climate system gradually shifted from an ice-free greenhouse world to an ice-house world with massive ice sheets at both poles (Zachos et al., 2001). Orbital forcing – quasi-cyclic changes in the climate system as a result of changes in the orbital configuration of the Earth (Milankovitch cyclicity) – caused superimposed variability in the climate system and carbon cycle on 104–105 year timescales. This has been shown in numerous palaeoclimate proxy records, notably stable carbon isotope (δ13C) and oxygen isotope (δ18O) ratios of deep-sea foraminifera. Particularly distinctive is the imprint of the 405 kyr eccentricity term (Westerhold et al., 2011) on climate and the carbon cycle, under both ice-free (Paleocene–Eocene; Lauretano et al., 2015; Westerhold et al., 2011) and glaciated (Oligocene-Miocene; Holbourn et al., 2015a; Pälike et al., 2006a) climate regimes within the Cenozoic (Fig. 1). On this timescale, δ13C co-varies in phase with δ18O, although the latter also responds strongly to the short eccentricity cycles (95 and 126 kyr). Similar to the shorter eccentricity cycles, the 405 kyr eccentricity cycle has a minor effect on total annual insolation (Laskar et al., 2004) but dominantly acts on the climate system as an amplitude modulator of precession, causing a greater seasonal contrast during eccentricity maxima. Changes in foraminiferal δ18O primarily represent changes in global ice volume and deep-sea temperatures. Changes in δ13C represent ocean-wide changes in the partitioning of carbon pools with different isotopic signatures among ocean, atmosphere, and rock reservoirs. The co-variance of oxygen and carbon isotopes therefore implies a coupling of climate and the carbon cycle, where δ13C and δ18O minima correspond with eccentricity maxima, or periods of increased seasonal contrast. It is important to note, however, that orbital-scale variability is a function of both cyclic forcing and other, probably non-cyclic variability on different timescales (e.g. from tectonics).
On timescales >106 years, the δ13C of the global exogenic carbon pool is largely controlled by carbon burial on the continental margins, where the isotopic composition of the buried carbon is mainly determined by the ratio between inorganic carbon and organic carbon (Corg) (Berner et al., 1983). Increases in 13C-depleted Corg burial are caused by increased delivery of clays through river runoff, as Corg is primarily buried in association with clay particles (Hedges and Keil, 1995). One might suspect that clay supply and burial increase during periods of increased seasonal contrast (eccentricity maxima) because this would optimize weathering, erosion, and transport. In this scenario, the resulting enhanced Corg burial would lead to δ13C maxima in the water column. However, we observe the opposite phase relation between eccentricity and δ13C in palaeo-records. Alternatively, clay formation may be optimized during phases of year-round relatively wet conditions in soils to facilitate silicate weathering and clay formation. If the system is not limited by transport rates, this scenario is more likely to occur during periods of low seasonal contrast (eccentricity minima), resulting in the observed phase relation between eccentricity and δ13C.
Various other forms of organic carbon sequestration may be the cause of the recorded phasing. For example, increased burial of Corg may occur on land during eccentricity minima because constant precipitation above peatlands provides conditions for the development of year-round soil anoxia (Kurtz et al., 2003; Zachos et al., 2010). The resulting 13C-enriched global exogenic carbon pool would be consistent with the phase relation we see in sediment records. This mechanism may be problematic, however, because it does not explain the consistent patterns throughout the Cenozoic, as terrestrial basins were unlikely stable over tens of millions of years. Another factor that might influence the total mass and isotopic composition of inorganic carbon in the oceans involves submarine methane hydrate dynamics. Microbes in continental margin sediments produce biogenic methane from Corg. This methane, which has δ13C values of down to −70 ‰ VPDB (Vienna Pee Dee Belemnite), can be temporarily stored in methane hydrates (Dickens, 2003, 2011). The total mass of methane stored in hydrates is affected by deep-sea temperature and could thus be controlled by orbital cycles. Assuming a steady sediment methane production, more biogenic methane is directly mixed into the exogenic carbon pool under relatively warm climates (Dean et al., 2018). Under colder climates, more methane is trapped in hydrates, resulting in an exogenic carbon pool that is 13C-enriched relative to the warm phase. This would also result in the co-occurrence of 405 and 100 kyr eccentricity maxima and relatively negative δ13C and δ18O values. However, in Paleocene–Eocene sediment records, methane hydrate dynamics cannot explain the co-occurrence and cyclicity of δ13C and δ18O because they would result in a slope of Δδ13C−Δδ18O that is very different from the one observed (Zeebe et al., 2017). They may play a role in the Oligocene and Miocene, (Holbourn et al., 2015a; Pälike et al., 2006a; Tian et al., 2013) but could also be related to changes in ice volume (see below). Given the intermediate ocean-depth temperature response to orbital variability and seasonal contrast in a fully coupled climate model (Lunt et al., 2011), methane hydrate dynamics could be a plausible explanation for coupled δ13C and δ18O variability.
Alternatively, ice-sheet dynamics have been proposed as the origin of the 405 kyr cycle for the past 35 Ma (de Boer et al., 2014). In a 3-D ice sheet and climate model, increased insolation with higher surface-air temperature triggers a reduction in ice volume through sub-shelf melting, which can be linked to a decrease in δ13C values through changes in the carbon cycle (de Boer et al., 2014). While the 405 kyr cycle in stable isotopes can be simulated adequately for the past 5 Myr (de Boer et al., 2014), astronomical pacing of the carbon cycle in the early Paleogene – with the absence of permanent ice sheets – means that this mechanism alone does not explain the link between the carbon cycle, astronomical forcing, and climate variability throughout the Cenozoic.
It is unknown which mechanism is responsible for the co-occurrence of the 405 and 100 kyr eccentricity cycles in climate records, but we do know that they must link astronomical forcing to changes in δ13C and δ18O. We decide to explore Corg burial in the oceans. In nature, continental margins provide the largest Corg storage potential over astronomical timescales (Berner, 1982). Corg burial is the product of sediment accumulation and organic carbon contents (corrected for porosity). Corg burial is thus linearly related to sediment accumulation, which in turn is approximately linearly related to sediment delivery to the ocean in the absence of major changes in sea level (Berner, 1982). Other factors governing Corg burial (e.g. oxygen) are secondary to sediment deposition (Hedges and Keil, 1995). Therefore, we focus our modelling efforts on linear astronomical forcing of Corg burial as a starting point. Since the model does not resolve shallow versus deep Corg burial, we force the long-term balance between burial and oxidation of Corg directly, without implying which particular mechanism is responsible.
Previous modelling studies have attempted to unravel the origin of the 405 kyr cyclic patterns found in Cenozoic δ13C and δ18O records. Pälike et al. (2006a) hypothesized that the cycles are driven by variations in marine biosphere productivity in response to changes in insolation. Ma et al. (2011), however, focused their simulations on the effects of increased weathering with eccentricity maxima, resulting in increased inorganic carbon input into the oceans, which, in combination with higher nutrients, leads to increased productivity and ultimately increased Corg burial, coinciding with minima in global exogenic δ13C. Paillard (2017), in a study on Pleistocene glacial–interglacial variability, argued that silicate weathering alone could not simultaneously account for the large δ13C changes and the small pCO2 changes, while nutrient or ecology-induced changes in Corg burial would not result in a change to these cycles during the onset of glaciations. Therefore Paillard (2017) focused on monsoon-driven changes in Corg burial, which may be similar to our approach – as the monsoonal changes lead to changes in precipitation, which cause an increased riverine clay flux, resulting in increased Corg burial. Zeebe et al. (2017) forced Corg burial in marine sediments directly with various orbital forcing scenarios and noise to understand a sediment record from the Paleocene and Eocene. These and previous modelling studies demonstrate how high-frequency spectral power in input forcing is dampened, resulting in increased power in the lower frequencies (Ma et al., 2011; Paillard, 2017; Pälike et al., 2006a; Zeebe et al., 2017). The long residence time of carbon in the oceans (Broecker and Peng, 1982) effectively filters the input signal, resulting in a strongly muted response on short timescales, which allows only the lowest frequencies in the input forcing to become apparent in the record.
The question remains if this low-pass filtering also affects the much longer 2.4 Myr eccentricity cycle, for which indications have been found in geological records for the Eocene (Lourens et al., 2005), Paleocene (Westerhold et al., 2011), Oligocene (Pälike et al., 2006a; Valero et al., 2014), and Miocene (Holbourn et al., 2007; Liebrand et al., 2016). The 2.4 Myr cycle is present in the astronomical solution as both an amplitude modulator (AM) of the 124 and 95 kyr eccentricity cycles as well as a true eccentricity cycle and is caused by orbit–orbit resonance between Mars and Earth (g4−g3) (Laskar et al., 2011). This cycle is never unambiguous in palaeoclimate records because long time series are required and the long period of this cycle overlaps with the timescale of slow tectonically driven climate change (>1 Myr), and it can therefore be obscured by the latter. The signal can also be lost when data are detrended prior to spectral analysis.
Here we aim to investigate the general role of eccentricity forcing on the carbon cycle across the Cenozoic and specifically provide a model–data comparison with Oligocene and Miocene carbon isotope records. We linearly force Corg burial in marine sediment with astronomical forcing in this modelling study to assess the response of the carbon cycle.
We investigate the response of the carbon cycle to variations in orbital parameters by conducting a modelling study with simple linear forcing of organic carbon burial, using the Long-term Ocean-atmosphere-Sediment CArbon cycle Reservoir model (Zeebe, 2012). In the model, Corg burial and weathering were forced with various aspects and combinations of the orbital solution by Laskar et al. (2004). We use spectral analysis techniques to study cyclic patterns in the model output. Finally, we compare modelling results to a long-term high-resolution data composite of published benthic foraminiferal isotope records from the Pacific Ocean.
2.1 LOSCAR model
The LOSCAR model simulates and tracks changes in the partitioning of carbon throughout the ocean basins, sediments, and atmosphere over both short (centuries) and long (millions of years) timescales (Zeebe, 2012). The model consists of a number of easily adjustable ocean boxes: in the set-up used here, it comprises an atmosphere, the Atlantic, Indian, and Pacific oceans (each with a surface, intermediate, and deep box), and one general high-latitude surface ocean box. All oceans are coupled to a prescribed number of sediment boxes (13 for each ocean here). A thermohaline circulation is imposed in the model by way of flow rates between specific ocean boxes. Various tracers (e.g. dissolved inorganic carbon (DIC), pCO2, pH, and the δ13C, as well as the carbonate compensation depth (CCD)) are tracked for the different ocean boxes in the model through time.
Carbon cycling in the model operates on various timescales that are relevant here, including volcanism and the long-term weathering feedbacks and the short-term gas exchange (Berner et al., 1983). Long-term carbon cycling is implemented in the model as a flux of volcanic degassing and metamorphic processes () and weathering of organic matter from the rock reservoir (), with atmospheric CO2 being consumed by carbonate weathering (Fcc, Reaction 0) and silicate weathering (Fsi, Reaction 0).
These weathering reactions respond as a negative feedback to variations in the CO2 concentration as follows:
where Fcc and Fsi are the carbonate and silicate weathering fluxes respectively and the superscript “0” refers to the initial (steady-state) value of the weathering flux and CO2 respectively. The ncc and nsi parameters control the strength of the weathering feedback. Steady state is reached when the silicate weathering flux balances the CO2 degassing flux from volcanism; i.e. Fsi=Fvc. These are all slow processes, where the steady-state balance is restored after a perturbation on timescales of the order of 105–106 years (Zeebe, 2012).
Carbon cycling also occurs on a much more rapid timescale in the oceans and atmosphere through the biological pump, air–sea exchange, and remineralization of carbon. In the original LOSCAR model (version 2.0.4), the biological pump is decoupled from the long carbon cycle described above, meaning that Corg that is produced in surface oceans is fully remineralized in the intermediate- and deep-ocean boxes (Zeebe, 2012). We note that Komar and Zeebe (2017) included a long-term, coupled Corg–O2–PO4 feedback in LOSCAR. However, a simplified approach is used here (see below), since our simple forcing does not invoke additional controls on Corg burial, such as productivity.
2.2 Orbital forcing
We spun up the modern (i.e. post-Paleocene–Eocene) set-up of LOSCAR with atmospheric CO2 concentrations of 450 ppmv with default parameter values listed in Zeebe (2012) for 10 Myr. Initial CO2 concentrations of 450 ppmv are roughly consistent with averaged late Paleogene to Neogene proxy estimates (Foster et al., 2012; Zhang et al., 2013). The equilibrium was then perturbed by linear scaling of the Corg burial rate with the orbital forcing signal at each time step (Eq. 3), with higher burial occurring with weaker forcing (Finp, some form of eccentricity–tilt–precession (ETP), eccentricity, insolation, or clipped insolation).
where Z is the input forcing, Finp, scaled between 0 and 1. The scaled orbital forcing, Forb, is Z minus its mean value, , multiplied by 2 so that it is scaled between −1 and 1. Fgb is the Corg flux resulting from the initial steady-state Corg flux, , being scaled with adjustable parameter α and Forb. Thus, weaker input forcing is associated with increased burial, which in nature could – if weathering is not limiting (Sect. 1) – be explained by increased clay formation, delivery, and deposition and associated burial of Corg.
To force these changes in Corg burial, different combinations of the orbital parameters can be used. Finp was implemented as eccentricity or precession separately or as artificial combinations of eccentricity, obliquity (tilt), and precession (together ETP), which is created by normalizing (subtracting the mean, then dividing by the standard deviation) and adding the ETP components after scaling them with different factors. Forcing the model with individual components and ETP may seem rather artificial but allows us to study the impact of individual orbital components and infers no mechanistic understanding of the underlying processes – such as the location of deep water formation. We also force with various insolation schemes with different latitudes and seasons.
To test the robustness of the carbon cycle response to changes in orbitally forced Corg burial, a parameterized amount of white noise (scaled between −1 and 1) was added to the burial rate. Noise was generated and applied when a time interval larger than 1 kyr since the previous noise addition had passed. As additional sensitivity analyses, several (newly introduced) model parameters were changed to study their effect on model output. The atmospheric CO2 concentration, the impact of the forcing on the Corg burial (α), the amount of noise, as well as the time step after which noise was added were adjusted to find how these affect model output. When using ETP as forcing, the relative make-up of the ETP curve, was varied, while for summer insolation different ways of artificial clipping were applied so as to simulate a non-linear response resulting in increased spectral power in the eccentricity components.
2.3 Data compilation
A high-resolution (22 Myr long, 34–12 Ma, median sampling interval of 4 kyr) composite record of benthic foraminiferal δ13C and δ18O (Fig. 1) was compiled using data from Ocean Drilling Program (ODP) Leg 199 Site 1218 (8∘53.378′ N, 135∘22.00′ W; 4.8 km water depth; Pälike et al., 2006a) and Integrated ODP (IODP) Site U1337 (3∘50.009′ N, 123∘12.352′ W; 4.2 km water depth; Holbourn et al., 2015a; Tian et al., 2013), located in the equatorial East Pacific (Fig. 1a). We use records from the Pacific because it represents the largest ocean reservoir by volume and compare these to output from LOSCAR's Pacific box. Initial age models for the records were established from 40Ar∕39Ar isotopes, planktonic foraminiferal-, nannofossil-, radiolarian-, and/or magnetostratigraphic-datum events, and known ages of changes in wt % CaCO3 contents. The age models of these records are based on similar dating and tuning strategies and we therefore assume that they are compatible for the purposes of this long-term study. Initial age models were subsequently improved through astronomical tuning to the orbital solution (Laskar et al., 2004). Pälike et al. (2006a) tuned their record to an ETP model (with the ratio ) starting with the 405 kyr cycle and then moving on to the shorter astronomical frequencies, while Holbourn et al. (2015a) tuned their δ18O minima to ETP () maxima and Tian et al. (2013) tuned δ18O minima to obliquity maxima (which represents the same phase relation).
Areas of overlap in the various records comprising the data composite were checked for consistency in phase relationship, absolute value (Sect. S3), and amplitude variation. The δ13C and δ18O composite data from Pälike et al. (2006a) were limited to the higher-resolution range between 20 and 34 Ma to keep the resolution comparable to the other records. Due to minor inconsistencies in the two age models for site U1337, data from Tian et al. (2013) were excluded in the overlapping depth region because of the lower resolution and shorter record length compared to the Holbourn et al. (2015a) record.
2.4 Spectral analysis
All analyses were performed using the R programming language (R Core Team, 2018), including various functions of the Astrochron package (Meyers, 2014), unless indicated differently. To identify periodicity in both the model output and the data composite, the multitaper method (MTM) by Thomson (1982) was applied. First, data and model output were linearly re-sampled to the median sampling interval of the data composite (4 kyr). For the MTM analysis, two tapers were used and the data or model output were zero-padded with a length of 5 times the input length. The data composite was also filtered for periods larger than 3.1 Myr to remove longer-scale climate trends that likely have causes other than orbital forcing to allow for an easier model–data comparison of records and spectra. The value of 3.1 Myr was chosen based on a small trough in spectral output of the raw data composite around this value (Sect. 3). This was achieved by applying a low-pass filter (for periods below 3.1 Myr), subtracting this signal from the original composite and re-adding the mean of the low-passed signal. Cross-correlation was performed between model output and input forcing to estimate phase lags and provide insight into the dominant mechanisms determining model response to orbital forcing of Corg burial. This was achieved with the Blackman–Tukey method, using a Bartlett window, with 30 % of the series, and a time step of 10−5 cycles kyr−1 using Analyseries 2.0.8 (Paillard et al., 1996). Amplitude modulation of shorter cycles was investigated by filtering for the short eccentricity cycles (0.01±0.003 kyr−1), performing a Hilbert transform which connects all the peaks in the record, and applying the MTM to the result (Zeebe et al., 2017). To evaluate how the frequency distribution and amplitude modulation changes through time in the composite record, evolutive multitaper harmonics analysis (EHA) was performed on the 1 Myr filtered, negative composite of δ13C using a time step of 50 kyr and a window size of 0.7 Myr (Liebrand et al., 2016; Pälike et al., 2006a).
3.1 Model output
All tracers in the model output show cyclic alternations as a result of implementing astronomical forcing of Corg burial in the model (Figs. 2 and S10 for all model output). Model tracers do not show much inter-basin variability; thus, we will consider the deep Pacific box (DP) to assess the generalized oceanic behaviour of tracers. As expected, the total exogenic carbon inventory Ctot in the ocean basins closely matches changes in atmospheric CO2, dissolved inorganic carbon (DIC), and temperature. While the carbonate compensation depth (CCD) of the DP shows cyclic extreme outliers due to model design (Sect. S1 in the Supplement), the underlying pattern of change is also similar, with high temperatures, high CO2 concentrations, and high Ctot co-occurring with a deepening of the CCD for the 405 kyr and 2.4 Myr cycles. DP δ13C co-varies approximately in anti-phase (i.e. 13C-depleted DIC co-occurring with high Ctot) with these tracers but shows more higher-frequency fluctuations.
The MTM analysis of model output shows significant spectral power at orbital periods (Fig. 3). Most notably, decreased spectral power occurs in the higher frequencies (precession, obliquity) in the forcing signal, while the lower frequencies (long-term eccentricity cycles), which are sometimes only marginally present in the input forcing, gain spectral power. This low-pass filtering of spectral power, occurs for our main scenario of ETP but also for various other combinations of ETP (e.g. , , ), eccentricity alone, and unclipped/clipped 65 and 30∘ N summer insolation (Sect. S4). Results are qualitatively similar for different contributions of E, T, and P in the ETP curve as forcing. Furthermore, the results are largely unaffected by the introduction of relatively strong white noise, which is transformed into red noise, which affects Corg burial by up to the same levels as the astronomical forcing, and different initial modelling conditions. Together, these sensitivity studies indicate that the model results in terms of spectral response are very robust (Sect. S4).
By comparing the dominant frequencies in our output to the characteristic frequencies of orbital solution from Laskar et al. (2004), the corresponding orbit–orbit resonances between the planets were found (Table S1 in the Supplement). Surprisingly, the 2.4 Myr cycle in carbon and oxygen isotopes is strongly present in the model output (including DIC, atmospheric CO2, deep Pacific temperature, and Ctot). The other strong eccentricity cycles consist of a 1 Myr cycle that is related to a Mercury–Jupiter (g1−g5) resonance and the 405 kyr cycle that is related to Venus–Jupiter (g2−g5).
Interestingly, the low-pass filtering through the carbon cycle model is different for the different tracers. The DIC, atmospheric CO2, Ctot, and temperature show very high relative dominance of the long eccentricity cycles (2.4, 1, and 0.7 Myr), while δ13C – and to a lesser extent CCD – show more spectral power in the shorter eccentricity range (405, 127, and 95 kyr) (Fig. 3). To investigate why spectral power is filtered differently for these tracers, residence times were calculated for the key carbon cycle feedbacks in the model by dividing the reservoir sizes by the respective fluxes (See Table 2). Silicate weathering, with a residence time of 655 kyr, is the slowest feedback system in our model. Carbonate weathering operates on timescales of the same order of magnitude (280 kyr), while oceanic carbonate compensation is an order of magnitude faster (10 kyr). This means that silicate weathering most strongly low-pass filters the spectral power, resulting in a dominant 2.4 Myr eccentricity cycle.
The δ13C gradients of carbon in the oceans are modulated by changes in the (13C-depleted) Corg pump. The biological pump operates much more rapidly than silicate- and carbonate weathering but remains constant in the model. Therefore, changes in the δ13C gradients are dominated by the total Corg inventory, on which we directly impose our forcing.
3.1.1 Phase relations of model tracers
Cross-spectral analysis between the key carbon cycle tracers in the model output and the ETP forcing reveals different phasing for the various tracers, as well as for the different timescales under consideration (Fig. 4). A phase of 0∘ implies that the tracer is changing exactly in phase with the ETP forcing, while at 180 or −180∘ they are in anti-phase. Positive (negative) values indicate that the variable is lagging (leading) with respect to the ETP forcing.
Broadly, we observe that phasing remains rather constant for each tracer for periods <405 kyr, while absolute lag time decreases as a function of frequency (Fig. S7). The atmospheric CO2 concentration and ocean temperature show the smallest lags – or fastest response – compared to the oceanic tracers. This is because the much shorter residence time of CO2 in the atmosphere (of the order of several years) – due to its relatively small reservoir size and fast cycling – causes atmospheric CO2 concentrations to respond rapidly to changes in sea surface carbon. The fast temperature response occurs because of the simple climate sensitivity parameterization used in the model. Atmospheric CO2 and ocean temperature therefore appear to be in quasi-steady state with respect to the other tracers in the model. Furthermore, cycling in the deep (D) and intermediate (I) Pacific (P) boxes is slower than the low-latitude surface (L) boxes for DIC. This is also a result of relative reservoir size and air–sea exchange. The CCD has a negative phase for the higher frequencies (precession, obliquity, and 100 kyr eccentricity), indicating anti-phasing relative to the other tracers. With ETP maxima, the CCD shoals at these shorter timescales due to relatively rapid carbonate compensation as a short-term response to increased atmospheric CO2 (Table 2). On longer timescales (>100 kyr), however, the CCD shifts to in-phase cycling – deepening of the CCD with ETP maxima – as a result of increased weathering and increased riverine input of DIC and TA, which increases the oceanic DIC inventory but also restores alkalinity to the ocean.
For the phase analysis, δ13C minima were compared to ETP maxima (i.e. −δ13C was compared to ETP), while the other tracers (CCD, Ctot, CO2) were not adjusted. The high-frequency range shows similar lags to the DIC and Ctot. Surprisingly, the DP box shows a slightly faster response than the LP and IP boxes, in fact more akin to the atmospheric reservoir for the precession frequencies. With decreasing frequency (i.e. the 405 kyr to 2.4 Myr eccentricity components), this relation between the deep and shallower boxes switches and the −δ13C starts to lag more and more, reaching an overshoot of more than half a period for the 2.4 Myr cycle in the DP box, causing the 2.4 Myr-cycle in −δ13C to appear to be leading the signal in the ETP forcing (Fig. 2). This can be demonstrated to be an overshoot response (Sect. S5). It is also interesting to note that δ13C leads over DIC. This is because the DIC response is determined by ocean chemistry and thus weathering, whereas the δ13C signal is mostly influenced by residence time and Corg cycling.
3.2 Data composite
The composite record consists of a carbon isotope record of relatively high resolution (median sampling interval of ∼4 kyr), spanning the Oligocene and most of the Miocene (34–12 Ma). However, there is a lower-resolution stretch of data between 20.0 and 21.5 Ma (Fig. 1). The highest and lowest (Nyquist and Rayleigh) detectible frequencies are thus and respectively.
The dominant frequencies in the data composite correspond to the main orbital periods of eccentricity (405, 124, and 95 kyr), obliquity (41 kyr), and precession (23, 22, and 19 kyr) (Fig. 5), similar to the results of the individual records (Holbourn et al., 2015a; Pälike et al., 2006a; Tian et al., 2013). Despite detrending, very high spectral power remains in the very low-frequency range (<0.5 Myr−1) (Fig. 5c). This could be caused by non-periodic long-term climatic variations or possibly by bundles of 2.4 Myr cycles (Boulila et al., 2012).
The dominant 2.4 Myr cyclicity observed in our model output seems absent in the δ13C composite record when performing MTM spectral analysis. The EHA records a strong influence of the 405, 127, and 95 kyr cycles throughout the composite record, which vary in intensity over 2.4 Myr intervals (Fig. 6). Furthermore, the envelope of the 100 and 405 kyr filtered record also changes at 2.4 Myr intervals (Fig. S8). This indicates that the 2.4 Myr cycle acts purely as an AM on the other eccentricity cycles in the composite record and not as a true cycle.
As a true cycle, the 2.4 Myr signal may be overshadowed by spectral power in the lower frequencies (Fig. 5). If these are filtered out (periods <3.1 Myr, red shaded region in Fig. 5c), the AR-1 fit and 95 % confidence levels are lowered on the low-frequency range (<20 Myr−1), but spectra are unaffected. The relatively high spectral peak in the 2.4 Myr range of the low-pass filtered data composite is therefore interpreted to be a combined result of a weak (overshadowed) signal in the data and an artefact produced by cutting off (low-pass filtering) the even lower frequencies, leaving only this frequency range in the spectral outcome. The obliquity and precession components also change through time but contribute less power to the EHA spectrum. Qualitatively, a weak 1 Myr bundling can be inferred for these frequencies, especially for the older range (25–34 Ma). This could be related to the ∼1 Myr obliquity cycle.
3.3 Model–data comparison
The model generates quite comparable time series to the high-passed data composite in terms of amplitude, phasing, and associated spectra. The distribution of spectral power in astronomical components is similar and can be matched to fit the data composite more closely by adjusting the strength of the orbital forcing, the make-up of the ETP, or the noise level (Sect. S4). However, while the 2.4 Myr eccentricity cycle is predominantly present as a true cycle in the model output, the composite record instead contains a strong presence of this cycle as amplitude modulator of the shorter eccentricity cycles. Upon further inspection, we also find ∼2.4 Myr bundling of the short eccentricity cycles in the EHA of model output (Fig. S9). This is probably the result of the eccentricity component in the ETP, as this bundling is absent when forcing with, e.g., 65∘N summer insolation. The changes in temperature (∼3 ∘C) that were calculated from model output CO2 (Fig. S10) are larger than the temperature estimates based on the δ18O composite data (which would be ∼1 ∘C). This is because the carbon cycle model (not a climate model) relates temperature to CO2 directly through an input climate sensitivity of 5 ∘C per CO2 doubling, whereas in reality temperature is a function of a number of variables, including but not limited to pCO2. For example, insolation directly affects the climate system components (e.g. atmosphere), which in turn affect temperature. These insolation-driven changes to the climate system are not part of the temperature response to pCO2, however. This behaviour cannot be captured by the carbon cycle model. In addition, this climate sensitivity strongly depends on the strengths of various feedbacks at various timescales, which likely also vary through time (Rohling et al., 2012; von der Heydt et al., 2016). The lack of a strong 2.4 Myr cycle in δ18O records (and its associated temperature) therefore does not necessarily imply that the 2.4 Myr cycle cannot be dominant in pCO2.
Unfortunately, even the available pCO2 records for the Pliocene–Pleistocene (Bartoli et al., 2011; Seki et al., 2010) are of insufficient resolution (∼200 kyr, potential for aliasing) and duration (∼5 Myr, at best two cycles) to test for a 2.4 Myr cycle. Furthermore, potentially cyclic changes could be overshadowed due to proxy uncertainty. New, long pCO2 proxy records are required to establish a possible 2.4 Myr cyclicity. Since the resulting model pCO2 can be adjusted by changing input parameters of forcing, noise, etc. (Sect. S4), comparison in terms of absolute values is not meaningful.
On long timescales, the model's CCD output can be compared to reconstructions of the Pacific CCD based on CaCO3 accumulation rates (CAR) from various sites (Pälike et al., 2012). These records indicate a deep and stable CCD during the Oligocene–Miocene up to about 18.5 Ma, but our modelled CCD is highly variable on 405 kyr and 2.4 Myr timescales – ranging from ∼10 to 100 m if we exclude model-set-up-induced outliers (Sect. S1). The CCD reconstructions for the Eocene, however, are highly variable on these timescales (changing by up to ∼500 m) and potentially show indications of long-term eccentricity-forced cycles. We performed spectral analyses on the CAR for the different sites as well as the inferred Pacific CCD (Figs. S1 and S2). In these records, long-term trends and step changes are controlled by long-term non-cyclic forcings such as volcanism, sea level, or long-term changes in the climate (Pälike et al., 2012). Therefore, we remove the long-term trend using a loess fit prior to spectral analysis. These analyses provide some evidence for 2.4 Myr cyclicity in the longest continuous CAR records for Site 1333 and Site 1334 (Fig. S1) as well as for the global CCD reconstruction (Fig. S2). The evidence is far from strong, however, because also these records are of low resolution – insufficient to resolve the short (100 and 405 kyr) eccentricity cycles and certainly precession. Furthermore, there are indications that the 2.4 Myr eccentricity cycle that is present in the astronomical solution is expressed as a 1.2 Myr AM of the short eccentricity cycle during the Eocene, due to either a resonance with or a non-linear response of the residence time of gas hydrates (Galeotti et al., 2010), further complicating the matter. The fact that we see some evidence for the 2.4 Myr cycle, given these caveats, however, is noteworthy.
The phasing of the model and data records differ in terms of absolute lag (kyr): the modelled δ13C response is overall slower compared to data records. For example, the lag of δ13C to the 405 kyr eccentricity cycle is ∼190 kyr, which is much larger than the 40–60 kyr estimates for the Paleocene (Westerhold et al., 2011) and 20–30 kyr during the Oligocene (Pälike et al., 2006a). Of course the lag times of data records are calculated after the tuning of a record and are therefore likely underestimated. However, the large difference between modelled lag times and lag times calculated from tuned records may indicate a lack of understanding regarding the underlying mechanisms.
With simple linear Corg burial forcing using an artificial combination of eccentricity, precession, and tilt, we are able to generate a δ13C profile for the Eocene–Miocene that is very comparable to a long-term high-passed data composite (Fig. 5) in both the time (Fig. 2) and frequency domains (Fig. 3). This indicates that cyclic burial of Corg should be considered as a potential candidate driving δ13C and δ18O variability on eccentricity timescales.
Future work could focus on evaluating the different magnitudes of the spectral power changes towards lower frequencies between different carbon cycle variables in the palaeo-record, as modelled δ13C of Corg shows a less pronounced dampening of high-frequencies compared to other model tracers such as the CCD, total exogenic carbon, temperature, or atmospheric CO2. Ideally this would require high-resolution records (at least 5 times the Nyquist frequency for the period of interest, i.e. every 2.3 kyr in order to resolve precession) and long-term (ideally more than 5 2.4 Myr cycles, spanning 12 Myr) records for various climate variables and could prove quite a challenge.
The carbon cycle box model simulations presented here reveal that simple orbital forcing of organic carbon burial in the ocean can reproduce the dominant 405 kyr cycle that is observed in Cenozoic carbon isotope records. The good correspondence between modelled carbon isotopes and the palaeo-record confirms orbitally forced burial of organic carbon as a possible candidate for causing the observed 405 kyr cyclicity in δ13C.
The link between astronomical forcing, carbon cycle dynamics, and the climate system that may give rise to this consistent eccentricity-forced periodicity is still poorly understood. Here we show that astronomical forcing of Corg burial in the marine realm could be a mechanism that explains this link in both a glaciated and an ice-free world. As Corg burial is primarily controlled by sediment accumulation (Berner, 1982, 2006; Hedges and Keil, 1995), the direct top–down control for Corg burial would be river runoff, which corresponds to seasonal precipitation patterns (Paillard, 2017; Zeebe et al., 2017). Although the exact mechanistic link between eccentricity forcing and sediment deposition remains elusive, we demonstrate that marine Corg burial provides an alternative to other proposed mechanisms that link the carbon cycle to astronomical forcing, such as carbon storage in peatland soils (Kurtz et al., 2003), methane hydrate reservoirs, or ice-sheet dynamics (de Boer et al., 2014).
Importantly, this modelling study implies a significant role for the 2.4 Myr eccentricity cycle in modulating carbon cycle dynamics, especially in DIC and atmospheric CO2. The 2.4 Myr eccentricity cycle is also present in the data composite but acts dominantly as an amplitude modulator of the 127 and 95 kyr eccentricity cycles and to a lesser extent as a true cycle, possibly due to overshadowing by slower cycles. The hypothesis that the Cenozoic pCO2 spectrum was dominated by the 2.4 Myr eccentricity could be evaluated using long, high-resolution proxy records of pCO2 or other carbonate chemistry parameters.
The LOSCAR model source code is available upon request to firstname.lastname@example.org. The implementation of astronomical forcing and noise addition to the model are available as a patch file on the corresponding author's GitHub (Kocken, 2019). The data that were used to create the data composite are available from the respective references (Holbourn et al., 2015b; Pälike et al., 2006b; Tian et al., 2013). The R code that was used to assemble the composite is available in the above-mentioned GitHub repository.
The supplement related to this article is available online at: https://doi.org/10.5194/cp-15-91-2019-supplement.
Modelling work was conducted by MJC and IJK. Data compilation was generated by IJK. All authors provided feedback on the paper.
The authors declare that they have no conflict of interest.
Ilja J. Kocken thanks the kfHein fund for financial support. Margot J. Cramwinckel thanks the Molengraaff
Fund for financial support. Richard E. Zeebe was supported by US NSF grants OCE12-20615
and OCE16-58023. This work used data generated on samples provided by the
International Ocean Discovery Program (IODP) and was carried out under the
program of the Netherlands Earth System Science Centre (NESSC). We also thank
Frits Hilgen for extensive discussions and two anonymous reviewers for their
input on the paper.
Edited by: Yves Godderis
Reviewed by: two anonymous referees
Bartoli, G., Hönisch, B., and Zeebe, R. E.: Atmospheric CO2 decline during the Pliocene intensification of Northern Hemisphere glaciations, Paleoceanography, 26, 1–14, https://doi.org/10.1029/2010PA002055, 2011. a
Berner, R. A., Lasaga, A. C., and Garrels, R. M.: The carbonate-silicate geochemical cycle and its effect on atmospheric carbon dioxide over the past 100 million years, 283, 641–683, https://doi.org/10.2475/ajs.283.7.641, 1983. a, b
Boulila, S., Galbrun, B., Laskar, J., and Pälike, H.: A ~9 myr cycle in Cenozoic δ13C record and long-term orbital eccentricity modulation: Is there a link?, Earth Planet. Sc. Lett., 317/318, 273–281, https://doi.org/10.1016/j.epsl.2011.11.017, 2012. a
Broecker, W. and Peng, T.-H.: Tracers in the sea, Lamont-Doherty Geological Observatory, Columbia University, 1982. a
Dean, J. F., Middelburg, J. J., Röckmann, T., Aerts, R., Blauw, L. G., Egger, M., Jetten, M. S. M., de Jong, A. E. E., Meisel, O. H., Rasigraf, O., Slomp, C. P., in 't Zandt, M. H., and Dolman, A. J.: Methane feedbacks to the global climate system in a warmer world, Rev. Geophys., 56, 1–44, https://doi.org/10.1002/2017RG000559, 2018. a
de Boer, B., Lourens, L. J., and van de Wal, R. S. W.: Persistent 400,000-year variability of Antarctic ice volume and the carbon cycle is revealed throughout the Plio-Pleistocene, Nat. Commun., 5, 2999, https://doi.org/10.1038/ncomms3999, 2014. a, b, c, d
Dickens, G. R.: Rethinking the global carbon cycle with a large, dynamic and microbially mediated gas hydrate capacitor, Earth Planet. Sc. Lett., 213, 169–183, https://doi.org/10.1016/S0012-821X(03)00325-X, 2003. a
Dickens, G. R.: Down the Rabbit Hole: Toward appropriate discussion of methane release from gas hydrate systems during the Paleocene-Eocene thermal maximum and other past hyperthermal events, Clim. Past, 7, 831–846, https://doi.org/10.5194/cp-7-831-2011, 2011. a
Foster, G. L., Lear, C. H., and Rae, J. W. B.: The evolution of pCO2, ice volume and climate during the middle Miocene, Earth Planet. Sc. Lett., 341–344, 243–254, https://doi.org/10.1016/j.epsl.2012.06.007, 2012. a
Galeotti, S., Krishnan, S., Pagani, M., Lanci, L., Gaudio, A., Zachos, J. C., Monechi, S., Morelli, G., and Lourens, L.: Orbital chronology of Early Eocene hyperthermals from the Contessa Road section, central Italy, Earth Planet. Sc. Lett., 290, 192–200, https://doi.org/10.1016/j.epsl.2009.12.021, 2010. a
Ogg, J. G.: Chapter 5 – Geomagnetic Polarity Time Scale, in: The Geologic Time Scale 2012, edited by: Gradstein, F. M., Ogg, J. G., Schmitz, M., and Ogg, G., Elsevier, 88–93, https://doi.org/10.1016/B978-0-444-59425-9.00005-6, 2012. a
Holbourn, A., Kuhnt, W., Schulz, M., Flores, J.-A., and Andersen, N.: Orbitally-paced climate evolution during the middle Miocene “Monterey” carbon-isotope excursion, Earth Planet. Sc. Lett., 261, 534–550, https://doi.org/10.1016/j.epsl.2007.07.026, 2007. a
Holbourn, A., Kuhnt, W., Kochhann, K. G. D., Andersen, N., and Sebastian Meier, K. J.: Global perturbation of the carbon cycle at the onset of the Miocene Climatic Optimum, Geology, 43, 123–126, https://doi.org/10.1130/G36317.1, 2015a. a, b, c, d, e, f
Holbourn, A., Kuhnt, W., Kochhann, K. G. D., Andersen, N., and Meier, K. J. S.: Stable isotope record, carbonate and organic carbon content of the Miocene section of IODP Site 321-U1337, PANGAEA, https://doi.org/10.1594/PANGAEA.839743, 2015b. a
Komar, N. and Zeebe, R. E.: Redox-controlled carbon and phosphorus burial: A mechanism for enhanced organic carbon sequestration during the PETM, Earth Planet. Sc. Lett., 479, 71–82, https://doi.org/10.1016/j.epsl.2017.09.011, 2017. a
Kurtz, A. C., Kump, L. R., Arthur, M. A., Zachos, J. C., and Paytan, A.: Early Cenozoic decoupling of the global carbon and sulfur cycles, Paleoceanography, 18, 1090, https://doi.org/10.1029/2003PA000908, 2003. a, b
Laskar, J., Fienga, A., Gastineau, M., and Manche, H.: La2010: a new orbital solution for the long-term motion of the Earth, Astron. Astrophys., 532, 1–15, https://doi.org/10.1051/0004-6361/201116836, 2011. a
Lauretano, V., Littler, K., Polling, M., Zachos, J. C., and Lourens, L. J.: Frequency, magnitude and character of hyperthermal events at the onset of the Early Eocene Climatic Optimum, Clim. Past, 11, 1313–1324, https://doi.org/10.5194/cp-11-1313-2015,2015. a
Liebrand, D., Beddow, H. M., Lourens, L. J., Pälike, H., Raffi, I., Bohaty, S. M., Hilgen, F. J., Saes, M. J. M., Wilson, P. A., Dijk, A. E. V., Hodell, D. A., Kroon, D., Huck, C. E., and Batenburg, S. J.: Cyclostratigraphy and eccentricity tuning of the early Oligocene oxygen and carbon isotope records from Walvis Ridge Site 1264, Earth Planet. Sc. Lett., 1, 1–14, https://doi.org/10.1016/j.epsl.2016.06.007, 2016. a, b
Lourens, L., Sluijs, A., Kroon, D., Zachos, J. C., Thomas, E., Röhl, U., Bowles, J., and Raffi, I.: Astronomical pacing of late Palaeocene to early Eocene global warming events, Nature, 435, 1083–1087, https://doi.org/10.1038/nature03814, 2005. a
Lunt, D. J., Ridgwell, A., Sluijs, A., Zachos, J., Hunter, S., and Haywood, A.: A model for orbital pacing of methane hydrate destabilization during the Palaeogene, Nat. Geosci., 4, 775–778, https://doi.org/10.1038/ngeo1266, 2011. a
Ma, W., Tian, J., Li, Q., and Wang, P.: Simulation of long eccentricity (400-kyr) cycle in ocean carbon reservoir during Miocene Climate Optimum: Weathering and nutrient response to orbital change, Geophys. Res. Lett., 38, 1–5, https://doi.org/10.1029/2011GL047680, 2011. a, b
Matthews, K. J., Maloney, K. T., Zahirovic, S., Williams, S. E., Seton, M., and Müller, R. D.: Global plate boundary evolution and kinematics since the late Paleozoic, Global Planet. Change, 146, 226–250, https://doi.org/10.1016/j.gloplacha.2016.10.002, 2016. a
Pälike, H., Norris, R. D., Herrle, J. O., Wilson, P. A., Coxall, H. K., Lear, C. H., Shackleton, N. J., Tripati, A. K., and Wade, B. S.: The heartbeat of the Oligocene climate system, Science, 314, 1894–1898, https://doi.org/10.1126/science.1133822, 2006a. a, b, c, d, e, f, g, h, i, j, k, l
Pälike, H., Norris, R. D., Herrle, J. O., Wilson, P. A., Coxall, H., Lear, C. H., Shackleton, N. J., Tripati, A. K., and Wade, B. S.: Age models and stable isotope analysis on sediment core Site 199–1218 from the equatorial Pacific, PANGAEA, https://doi.org/10.1594/PANGAEA.547942, 2006b. a
Pälike, H., Lyle, M. W., Nishi, H., Raffi, I., Ridgwell, A., Gamage, K., Klaus, A., Acton, G., Anderson, L., Backman, J., Baldauf, J., Beltran, C., Bohaty, S. M., Bown, P., Busch, W., Channell, J. E., Chun, C. O., Delaney, M., Dewangan, P., Dunkley Jones, T., Edgar, K. M., Evans, H., Fitch, P., Foster, G. L., Gussone, N., Hasegawa, H., Hathorne, E. C., Hayashi, H., Herrle, J. O., Holbourn, A., Hovan, S., Hyeong, K., Iijima, K., Ito, T., Kamikuri, S. I., Kimoto, K., Kuroda, J., Leon-Rodriguez, L., Malinverno, A., Moore, T. C., Murphy, B. H., Murphy, D. P., Nakamura, H., Ogane, K., Ohneiser, C., Richter, C., Robinson, R., Rohling, E. J., Romero, O., Sawada, K., Scher, H., Schneider, L., Sluijs, A., Takata, H., Tian, J., Tsujimoto, A., Wade, B. S., Westerhold, T., Wilkens, R., Williams, T., Wilson, P. A., Yamamoto, Y., Yamamoto, S., Yamazaki, T., and Zeebe, R. E.: A Cenozoic record of the equatorial Pacific carbonate compensation depth, Nature, 488, 609–614, https://doi.org/10.1038/nature11360, 2012. a, b
Rohling, E. J., Sluijs, A., Dijkstra, H. A., Köhler, P., Van De Wal, R. S., Von Der Heydt, A. S., Beerling, D. J., Berger, A., Bijl, P. K., Crucifix, M., Deconto, R., Drijfhout, S. S., Fedorov, A., Foster, G. L., Ganopolski, A., Hansen, J., Hönisch, B., Hooghiemstra, H., Huber, M., Huybers, P., Knutti, R., Lea, D. W., Lourens, L. J., Lunt, D., Masson-Demotte, V., Medina-Elizalde, M., Otto-Bliesner, B., Pagani, M., Pälike, H., Renssen, H., Royer, D. L., Siddall, M., Valdes, P., Zachos, J. C., and Zeebe, R. E.: Making sense of palaeoclimate sensitivity, Nature, 491, 683–691, https://doi.org/10.1038/nature11574, 2012. a
Seki, O., Foster, G. L., Schmidt, D. N., Mackensen, A., Kawamura, K., and Pancost, R. D.: Alkenone and boron-based Pliocene pCO2 records, Earth Planet. Sc. Lett., 292, 201–211, https://doi.org/10.1016/j.epsl.2010.01.037, 2010. a
Tian, J., Yang, M., Lyle, M. W., Wilkens, R., and Shackford, J. K.: Obliquity and long eccentricity pacing of the Middle Miocene climate transition, Geochem. Geophy. Geosy., 14, 1740–1755, https://doi.org/10.1002/ggge.20108, 2013. a, b, c, d, e, f
Valero, L., Garcés, M., Cabrera, L., Costa, E., and Sáez, A.: 20 Myr of eccentricity paced lacustrine cycles in the Cenozoic Ebro Basin, Earth Planet. Sc. Lett., 408, 183–193, https://doi.org/10.1016/j.epsl.2014.10.007, 2014. a
von der Heydt, A. S., Dijkstra, H. A., van de Wal, R. S., Caballero, R., Crucifix, M., Foster, G. L., Huber, M., Köhler, P., Rohling, E., Valdes, P. J., Ashwin, P., Bathiany, S., Berends, T., van Bree, L. G., Ditlevsen, P., Ghil, M., Haywood, A. M., Katzav, J., Lohmann, G., Lohmann, J., Lucarini, V., Marzocchi, A., Pälike, H., Baroni, I. R., Simon, D., Sluijs, A., Stap, L. B., Tantet, A., Viebahn, J., and Ziegler, M.: Lessons on Climate Sensitivity From Past Climate Changes, Curr. Clim. Chang. Reports, 2, 148–158, https://doi.org/10.1007/s40641-016-0049-3, 2016. a
Westerhold, T., Röhl, U., Donner, B., McCarren, H. K., and Zachos, J. C.: A complete high-resolution Paleocene benthic stable isotope record for the central Pacific (ODP Site 1209), Global Biogeochem. Cy., 25, 1–13, https://doi.org/10.1029/2011WR010620, 2011. a, b, c, d, e
Zachos, J., Pagani, M., Sloan, L., Thomas, E., and Billups, K.: Trends, rhythms, and aberrations in global climate 65 Ma to present, Science, 292, 686–693, 2001. a
Zachos, J. C., McCarren, H., Murphy, B., Röhl, U., and Westerhold, T.: Tempo and scale of late Paleocene and early Eocene carbon isotope cycles: Implications for the origin of hyperthermals, Earth Planet. Sc. Lett., 299, 242–249, https://doi.org/10.1016/j.epsl.2010.09.004, 2010. a
Zeebe, R. E., Westerhold, T., Littler, K., and Zachos, J. C.: Orbital forcing of the Paleocene and Eocene carbon cycle, Paleoceanography, 32, 440–465, https://doi.org/10.1002/2016PA003054, 2017. a, b, c, d, e