The 405 kyr and 2.4 Myr eccentricity components in Cenozoic carbon isotope records

Cenozoic stable carbon (δ¹³C) and oxygen (δ¹⁸O) 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 δ¹⁸O are partly forced by ice-volume changes that have mostly occurred since the Oligocene. The cyclic δ¹³C–δ¹⁸O co-variation is found in both ice-free and glaciated climate states, however. Consequently, there should be a mechanism that forces the δ¹³C 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 δ¹³C mainly varies with a 405 kyr period in the model, the dynamics of dissolved inorganic carbon in the oceans and of atmospheric CO₂ 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) δ¹³C and δ¹⁸O 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.

1 Introduction

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 10⁴–10⁵ year timescales. This has been shown in numerous palaeoclimate proxy records, notably stable carbon isotope (δ¹³C) and oxygen isotope (δ¹⁸O) 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; Westerhold et al., 2011; Lauretano et al., 2015) and glaciated (Oligocene-Miocene; Pälike et al., 2006a; Holbourn et al., 2015a) climate regimes within the Cenozoic (Fig. 1). On this timescale, δ¹³C co-varies in phase with δ¹⁸O, 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 δ¹⁸O primarily represent changes in global ice volume and deep-sea temperatures. Changes in δ¹³C 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 δ¹³C and δ¹⁸O 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).

Figure 1 Long-term high-resolution benthic foraminiferal stable carbon and oxygen isotope records show 405 kyr cyclicity. (a) Palaeoceanographic reconstruction of 30 Ma with approximate site locations. Map generated with GPlates, using the rotation frame and tectonic reconstruction of (Matthews et al., 2016). Data from references in the legend show the 405 kyr periodicity in the δ¹⁸O (b) and δ¹³C (c) records throughout the Cenozoic at various sites in the Pacific. Darker lines going through the records represent 20-point moving averages while darker lines at arbitrary heights represent the 405 kyr filtered records (at 2.47±0.15 Myr⁻¹). The Westerhold et al. (2011) record was converted to the GTS2012 (Ogg, 2012) timescale by converting magnetic chron ages and placing the Paleocene–Eocene Thermal Maximum (PETM) onset at 56.0 Ma. Note the ∼16 Myr x-axis break in the Eocene.

On timescales >10⁶ years, the δ¹³C 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 (C_org) (Berner et al., 1983). Increases in ¹³C-depleted C_org burial are caused by increased delivery of clays through river runoff, as C_org 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 C_org burial would lead to δ¹³C maxima in the water column. However, we observe the opposite phase relation between eccentricity and δ¹³C 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 δ¹³C.

Various other forms of organic carbon sequestration may be the cause of the recorded phasing. For example, increased burial of C_org 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 ¹³C-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 C_org. This methane, which has δ¹³C 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 ¹³C-enriched relative to the warm phase. This would also result in the co-occurrence of 405 and 100 kyr eccentricity maxima and relatively negative δ¹³C and δ¹⁸O values. However, in Paleocene–Eocene sediment records, methane hydrate dynamics cannot explain the co-occurrence and cyclicity of δ¹³C and δ¹⁸O because they would result in a slope of Δδ¹³C−Δδ¹⁸O that is very different from the one observed (Zeebe et al., 2017). They may play a role in the Oligocene and Miocene, (Pälike et al., 2006a; Holbourn et al., 2015a; 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 δ¹³C and δ¹⁸O 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 δ¹³C 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.

Table 1 Possible mechanisms that can cause δ¹³C–δ¹⁸O to co-occur as a result of input eccentricity ϵ forcing maxima (↑) and minima (↓).

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 δ¹³C and δ¹⁸O. We decide to explore C_org burial in the oceans. In nature, continental margins provide the largest C_org storage potential over astronomical timescales (Berner, 1982). C_org burial is the product of sediment accumulation and organic carbon contents (corrected for porosity). C_org 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 C_org burial (e.g. oxygen) are secondary to sediment deposition (Hedges and Keil, 1995). Therefore, we focus our modelling efforts on linear astronomical forcing of C_org burial as a starting point. Since the model does not resolve shallow versus deep C_org burial, we force the long-term balance between burial and oxidation of C_org 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 δ¹³C and δ¹⁸O 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 C_org burial, coinciding with minima in global exogenic δ¹³C. Paillard (2017), in a study on Pleistocene glacial–interglacial variability, argued that silicate weathering alone could not simultaneously account for the large δ¹³C changes and the small pCO₂ changes, while nutrient or ecology-induced changes in C_org burial would not result in a change to these cycles during the onset of glaciations. Therefore Paillard (2017) focused on monsoon-driven changes in C_org 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 C_org burial. Zeebe et al. (2017) forced C_org 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 (Pälike et al., 2006a; Ma et al., 2011; Paillard, 2017; 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 (g₄−g₃) (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 C_org burial in marine sediment with astronomical forcing in this modelling study to assess the response of the carbon cycle.

2 Material and methods

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 (LOSCAR; Zeebe, 2012). In the model, C_org 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), pCO₂, pH, and the δ¹³C, 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 (similar to Berner et al., 1983). Long-term carbon cycling is implemented in the model as a flux of volcanic degassing and metamorphic processes ($F_{\mathrm{vc}}^{\mathrm{0}}$) and weathering of organic matter from the rock reservoir ($F_{\mathrm{OC}}^{\mathrm{in}}$), with atmospheric CO₂ being consumed by carbonate weathering (F_cc, Reaction R1) and silicate weathering (F_si, Reaction R2).

$$\begin{array}{}\text{(R1)}& {\displaystyle }& {\displaystyle }{\mathrm{CaCO}}_{\mathrm{3}}+{\mathrm{CO}}_{\mathrm{2}}+{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\rightleftharpoons {\mathrm{Ca}}^{\mathrm{2}+}+\mathrm{2}{\mathrm{HCO}}_{\mathrm{3}}^{-}\text{(R2)}& {\displaystyle }& {\displaystyle }{\mathrm{CaSiO}}_{\mathrm{3}}+\mathrm{2}{\mathrm{CO}}_{\mathrm{2}}+{\mathrm{H}}_{\mathrm{2}}\mathrm{O}\rightleftharpoons {\mathrm{Ca}}^{\mathrm{2}+}+{\mathrm{SiO}}_{\mathrm{2}}+\mathrm{2}{\mathrm{HCO}}_{\mathrm{3}}^{-}\end{array}$$

These weathering reactions respond as a negative feedback to variations in the CO₂ concentration as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }{F}_{\mathrm{cc}}=F_{\mathrm{cc}}^{\mathrm{0}}(p{\mathrm{CO}}_{\mathrm{2}}/p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{0}}{)}^{n_{\mathrm{cc}}},\text{(2)}& {\displaystyle }& {\displaystyle }{F}_{\mathrm{si}}=F_{\mathrm{si}}^{\mathrm{0}}(p{\mathrm{CO}}_{\mathrm{2}}/p{\mathrm{CO}}_{\mathrm{2}}^{\mathrm{0}}{)}^{n_{\mathrm{si}}},\end{array}$$

where F_cc and F_si are the carbonate and silicate weathering fluxes respectively and the superscript “0” refers to the initial (steady-state) value of the weathering flux and CO₂ respectively. The n_cc and n_si parameters control the strength of the weathering feedback. Steady state is reached when the silicate weathering flux balances the CO₂ degassing flux from volcanism; i.e. F_si=F_vc. These are all slow processes, where the steady-state balance is restored after a perturbation on timescales of the order of 10⁵–10⁶ 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 C_org that is produced in surface oceans is fully remineralized in the intermediate- and deep-ocean boxes (see Zeebe, 2012, for a full description). We note that Komar and Zeebe (2017) included a long-term, coupled C_org–O₂–PO₄ feedback in LOSCAR. However, a simplified approach is used here (see below), since our simple forcing does not invoke additional controls on C_org burial, such as productivity.

2.2 Orbital forcing

We spun up the modern (i.e. post-Paleocene–Eocene) set-up of LOSCAR with atmospheric CO₂ concentrations of 450 ppmv with default parameter values listed in Zeebe (2012) for 10 Myr. Initial CO₂ concentrations of 450 ppmv are roughly consistent with averaged late Paleogene to Neogene proxy estimates (e.g. Foster et al., 2012; Zhang et al., 2013). The equilibrium was then perturbed by linear scaling of the C_org burial rate with the orbital forcing signal at each time step (Eq. 3), with higher burial occurring with weaker forcing (F_inp, some form of eccentricity–tilt–precession (ETP), eccentricity, insolation, or clipped insolation).

$$\begin{array}{}\text{(3)}& {\displaystyle }& {\displaystyle }Z={\displaystyle \frac{F_{\mathrm{inp}}-min{F}_{\mathrm{inp}}}{max{F}_{\mathrm{inp}}-min{F}_{\mathrm{inp}}}}\text{(4)}& {\displaystyle }& {\displaystyle }{F}_{\mathrm{orb}}=\mathrm{2}(Z-\stackrel{\mathrm{‾}}{Z})\text{(5)}& {\displaystyle }& {\displaystyle }{F}_{\mathrm{gb}}=F_{\mathrm{gb}}^{\mathrm{0}}(\mathrm{1}-\mathit{\alpha }{F}_{\mathrm{orb}}),\end{array}$$

where Z is the input forcing, F_inp, scaled between 0 and 1. The scaled orbital forcing, F_orb, is Z minus its mean value, $\stackrel{\mathrm{‾}}{Z}$, multiplied by 2 so that it is scaled between −1 and 1. F_gb is the C_org flux resulting from the initial steady-state C_org flux, $F_{\mathrm{gb}}^{\mathrm{0}}$, being scaled with adjustable parameter α and F_orb. 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 C_org.

To force these changes in C_org burial, different combinations of the orbital parameters can be used. F_inp 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 C_org 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 CO₂ concentration, the impact of the forcing on the C_org 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 δ¹³C and δ¹⁸O (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; Tian et al., 2013; Holbourn et al., 2015a), 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 ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ isotopes, planktonic foraminiferal-, nannofossil-, radiolarian-, and/or magnetostratigraphic-datum events, and known ages of changes in wt % CaCO₃ 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 $\mathrm{1}:\mathrm{0.5}:-\mathrm{0.4}$) starting with the 405 kyr cycle and then moving on to the shorter astronomical frequencies, while Holbourn et al. (2015a) tuned their δ¹⁸O minima to ETP ($\mathrm{1}:\mathrm{1}:+\mathrm{0.2}$) maxima and Tian et al. (2013) tuned δ¹⁸O 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 δ¹³C and δ¹⁸O 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 C_org burial. This was achieved with the Blackman–Tukey method, using a Bartlett window, with 30 % of the series, and a time step of 10⁻⁵ cycles kyr⁻¹ 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⁻¹), performing a Hilbert transform which connects all the peaks in the record, and applying the MTM to the result (adapted from 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 δ¹³C using a time step of 50 kyr and a window size of 0.7 Myr (similar to Pälike et al., 2006a; Liebrand et al., 2016).

Figure 2 LOSCAR model output time series for orbitally forced carbon burial with orbital forcing strength α=0.5 and noise level 0.2. (a) 2.4 Myr (0.42±0.15 Myr⁻¹) and (b) 405 kyr (frequency of 2.46±0.15 Myr⁻¹) bandpass filters are shown (same colours, arbitrary y axis). (c) The ETP ($\mathrm{1}:\mathrm{0.5}:-\mathrm{0.4}$) construct used as orbital forcing in this scenario (gold). (d) Atmospheric CO₂ (purple, parts per million volume, ppm) and (e) the total exogenic carbon inventory (brown, Eg (10¹⁸ g)). (f) Pacific CCD (blue, km; note the cutoff y axis due to large outliers, explained in Sect. S1) (g) Deep Pacific δ¹³C (green, ‰ VPDB, reversed axis). Deep Pacific temperature (°C), and DIC of the deep Pacific model output are omitted from this plot because their pattern mirrors that of C_tot and pCO₂ (Fig. S10).

3 Results and discussion

3.1 Model output

All tracers in the model output show cyclic alternations as a result of implementing astronomical forcing of C_org 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 C_tot in the ocean basins closely matches changes in atmospheric CO₂, 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 CO₂ concentrations, and high C_tot co-occurring with a deepening of the CCD for the 405 kyr and 2.4 Myr cycles. DP δ¹³C co-varies approximately in anti-phase (i.e. ¹³C-depleted DIC co-occurring with high C_tot) 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 $\mathrm{1}:\mathrm{0.5}:-\mathrm{0.4}$ ETP but also for various other combinations of ETP (e.g. $\mathrm{1}:\mathrm{0.5}:-\mathrm{0.4}$, $\mathrm{0.1}:\mathrm{0}:\mathrm{0.1}$, $\mathrm{1}:\mathrm{1}:\mathrm{1}$), 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 C_org 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).

Figure 3 MTM spectral analysis of LOSCAR model output with orbital forcing strength α=0.5 and noise level 0.2. Shown are the model tracers and the ETP forcing used on absolute (a) (log axis) and normalized (b) (linear axis) spectral power. Relative spectral power was calculated by dividing each value by their maximum value. Insets (c) and (d) are the same as (a) and (b) but zoomed in on the frequency axis. Periods of interest (kyr) are labelled and indicated with dashed lines.

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 CO₂, deep Pacific temperature, and C_tot). The other strong eccentricity cycles consist of a 1 Myr cycle that is related to a Mercury–Jupiter (g₁−g₅) resonance and the 405 kyr cycle that is related to Venus–Jupiter (g₂−g₅).

Interestingly, the low-pass filtering through the carbon cycle model is different for the different tracers. The DIC, atmospheric CO₂, C_tot, and temperature show very high relative dominance of the long eccentricity cycles (2.4, 1, and 0.7 Myr), while δ¹³C – 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.

Table 2 Calculated residence times of the main carbon cycle processes.

The δ¹³C gradients of carbon in the oceans are modulated by changes in the (¹³C-depleted) C_org pump. The biological pump operates much more rapidly than silicate- and carbonate weathering but remains constant in the model. Therefore, changes in the δ¹³C gradients are dominated by the total C_org 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 CO₂ 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 CO₂ in the atmosphere (of the order of several years) – due to its relatively small reservoir size and fast cycling – causes atmospheric CO₂ 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 CO₂ 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 CO₂ (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, δ¹³C minima were compared to ETP maxima (i.e. −δ¹³C was compared to ETP), while the other tracers (CCD, C_tot, CO₂) were not adjusted. The high-frequency range shows similar lags to the DIC and C_tot. 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 −δ¹³C 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 −δ¹³C 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 δ¹³C leads over DIC. This is because the DIC response is determined by ocean chemistry and thus weathering, whereas the δ¹³C signal is mostly influenced by residence time and C_org cycling.

Figure 4 Phase lag (°, y axis) of the various modelled carbon cycle tracers (see legend; L stands for low-latitude surface, I for intermediate, and D for deep ocean box of the Pacific (P)) with respect to the ETP forcing (of strength 0.5 and noise 0.2) as a function of frequency (cycles kyr⁻¹, bottom x axis) and period (kyr, top x axis). Note minor break in x axes to better show the low-frequency range. Error bars represent 2σ of the variation of the phase over the specified frequency intervals (red shaded regions). Note that −δ¹³C is shown, such that a phase lag (positive values) represents a lag of δ¹³C minima with respect to ETP maxima.

Figure 5 (a) δ¹³C composite record (grey line) with a low-pass filter for frequencies below 3.1 Myr (dashed purple line). (b) The adjusted composite (purple) in the time domain. (c) The same records in the frequency domain. Dashed smooth lines represent the AR-1 fits for the δ¹³C and detrended δ¹³C records (grey and purple respectively) with associated 95 % confidence intervals (light and dark red respectively). Main orbital periods (kyr) are labelled and marked with dotted lines. (d) A zoomed-in version of (c) on the frequency axis, with the red shaded region indicating the range of the low-pass filter.

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 $F_N=\mathrm{1}/\left(\mathrm{2}\mathrm{\Delta }t\right)=\mathrm{0.125}{\mathrm{kyr}}^{-\mathrm{1}}$ and $F_R=\mathrm{1}/\left(N\mathrm{\Delta }t\right)=\mathrm{0.046}{\mathrm{Myr}}^{-\mathrm{1}}$ 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 (Pälike et al., 2006a; Tian et al., 2013; Holbourn et al., 2015a). Despite detrending, very high spectral power remains in the very low-frequency range (<0.5 Myr⁻¹) (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 δ¹³C 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⁻¹), 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.

Figure 6 Evolutive multitaper harmonics analysis (c) of the data composite δ¹³C (grey) after detrending (black) the <1 Myr signal (red line) in the time (a) and frequency (b) domains. Periods of interest are labelled (b) on the same axis as (c). A window size of 0.7 Myr was used, with a step size of 10 kyr (lines on the right side in c). Methods adapted from Pälike et al. (2006a).

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 CO₂ (Fig. S10) are larger than the temperature estimates based on the δ¹⁸O composite data (which would be ∼1 °C). This is because the carbon cycle model (not a climate model) relates temperature to CO₂ directly through an input climate sensitivity of 5 °C per CO₂ doubling, whereas in reality temperature is a function of a number of variables, including but not limited to pCO₂. 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 pCO₂, 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 (e.g. Rohling et al., 2012; von der Heydt et al., 2016). The lack of a strong 2.4 Myr cycle in δ¹⁸O records (and its associated temperature) therefore does not necessarily imply that the 2.4 Myr cycle cannot be dominant in pCO₂.

Unfortunately, even the available pCO₂ records for the Pliocene–Pleistocene (e.g. 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 pCO₂ proxy records are required to establish a possible 2.4 Myr cyclicity. Since the resulting model pCO₂ 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 CaCO₃ 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 δ¹³C response is overall slower compared to data records. For example, the lag of δ¹³C 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 C_org burial forcing using an artificial combination of eccentricity, precession, and tilt, we are able to generate a δ¹³C 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 C_org should be considered as a potential candidate driving δ¹³C and δ¹⁸O variability on eccentricity timescales.

3.4 Outlook

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 δ¹³C of C_org shows a less pronounced dampening of high-frequencies compared to other model tracers such as the CCD, total exogenic carbon, temperature, or atmospheric CO₂. 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.

4 Conclusions

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 δ¹³C.

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 C_org burial in the marine realm could be a mechanism that explains this link in both a glaciated and an ice-free world. As C_org burial is primarily controlled by sediment accumulation (Berner, 1982; Hedges and Keil, 1995; Berner, 2006), the direct top–down control for C_org burial would be river runoff, which corresponds to seasonal precipitation patterns (see also Paillard, 2017; Zeebe et al., 2017). Although the exact mechanistic link between eccentricity forcing and sediment deposition remains elusive, we demonstrate that marine C_org 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 CO₂. 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 pCO₂ spectrum was dominated by the 2.4 Myr eccentricity could be evaluated using long, high-resolution proxy records of pCO₂ or other carbonate chemistry parameters.