The 405 kyr and 2 . 4 Myr eccentricity components in Cenozoic carbon isotope records

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 covariations are 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 5 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 10 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 time scale 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 to 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 is present as 15 an amplitude modulator of the 405 and 100 kyr eccentricity cycles.


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 4 -10 5 year timescales.This has been shown in numerous palaeoclimate proxy records, notably stable carbon isotope (δ 13 C) and oxygen isotope (δ 18 O) ratios of deepsea 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, δ 13 C co-varies in phase with δ 18 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 δ 18 O primarily represent changes in global ice volume and deep-sea tem-peratures.Changes in δ 13 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 δ 13 C and δ 18 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).
On timescales > 10 6 years, the δ 13 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 13 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 δ 13 C maxima in the water column.However, we observe the opposite phase relation between eccentricity and δ 13 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 δ 13 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 13 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 δ 13 C values of down to −70 ‰ VPDB (Vienna Pee Dee Belemnite), can be temporarily stored in methane hydrates (Dickens, 2003(Dickens, , 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, result-ing in an exogenic carbon pool that is 13 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 δ 13 C and δ 18 O values.However, in Paleocene-Eocene sediment records, methane hydrate dynamics cannot explain the co-occurrence and cyclicity of δ 13 C and δ 18 O because they would result in a slope of δ 13 C− δ 18 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 δ 13 C and δ 18 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 δ 13 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.
It is unknown which mechanism is responsible for the cooccurrence of the 405 and 100 kyr eccentricity cycles in climate records, but we do know that they must link astronomical forcing to changes in δ 13 C and δ 18 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 δ 13 C and δ 18 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  Requires presence of ice sheets 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 δ 13 C. Paillard (2017), in a study on Pleistocene glacial-interglacial variability, argued that silicate weathering alone could not simultaneously account for the large δ 13 C changes and the small pCO 2 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 4 − g 3 ) (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.

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.

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 2 , pH, and the δ 13 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 0 vc ) and weathering of organic matter from the rock reservoir (F in OC ), with atmospheric CO 2 being consumed by carbonate weathering (F cc , Reaction R1) and silicate weathering (F si , Reaction R2).
These weathering reactions respond as a negative feedback to variations in the CO 2 concentration as follows: 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 2 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 2 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 5 -10 6 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 deepocean boxes (see Zeebe, 2012, for a full description).We note that Komar and Zeebe (2017) included a long-term, coupled C org -O 2 -PO 4 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.

Orbital forcing
We spun up the modern (i.e.post-Paleocene-Eocene) setup of LOSCAR with atmospheric CO 2 concentrations of 450 ppmv with default parameter values listed in Zeebe (2012) for 10 Myr.Initial CO 2 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-tiltprecession (ETP), eccentricity, insolation, or clipped insolation). (3) 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, 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 0 gb , 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 2 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.

Data compilation
A high-resolution (22 Myr long, 34-12 Ma, median sampling interval of 4 kyr) composite record of benthic foraminiferal δ 13 C and δ 18 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.009N, 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 40 Ar/ 39 Ar isotopes, planktonic foraminiferal-, nannofossil-, radiolarian-, and/or magnetostratigraphic-datum events, and known ages of changes in wt % CaCO 3 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 1 : 0.5 : −0.4) starting with the 405 kyr cycle and then moving on to the shorter astronomical frequencies, while Holbourn et al. (2015a) tuned their δ 18 O minima to ETP (1 : 1 : +0.2) maxima and Tian et al. (2013) tuned δ 18 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 δ 13 C and δ 18 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.

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 compari-son 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 readding 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 −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 (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 δ 13 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).

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 2 , 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 2 concentrations, and high C tot co-occurring with a deepening of the CCD for the 405 kyr and 2.4 Myr cycles.DP δ 13 C co-varies approximately in anti-phase (i.e. 13 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 1 : 0.5 : −0.4 ETP but also for various other combinations of ETP (e.g. 1 : 0.5 : −0.4,0.1 : 0 : 0.1, 1 : 1 : 1), eccentricity alone, and unclipped/clipped 65  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 2 , 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 1 − g 5 ) resonance and the 405 kyr cycle that is related to Venus-Jupiter (g 2 − g 5 ).
Interestingly, the low-pass filtering through the carbon cycle model is different for the different tracers.The DIC, atmospheric CO 2 , C tot , and temperature show very high relative dominance of the long eccentricity cycles (2.4, 1, and 0.7 Myr), while δ 13 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.
The δ 13 C gradients of carbon in the oceans are modulated by changes in the ( 13 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 δ 13 C gradients are dominated by the total C org inventory, on which we directly impose our forcing.

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 2 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 2 in the atmosphere (of the order of several years) -due to its relatively small reservoir size and fast cycling -causes atmospheric CO 2 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 2 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 2 (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, δ 13 C minima were compared to ETP maxima (i.e.−δ 13 C was compared to ETP), while the other tracers (CCD, C tot , CO 2 ) were not adjusted.The highfrequency 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 −δ 13 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 −δ 13 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 δ 13 C leads over DIC.This is because the DIC response is determined by ocean chemistry and thus weathering, whereas the δ 13 C signal is mostly influenced by residence time and C org cycling.

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 = 1/(2 t) = 0.125 kyr −1 and F R = 1/(N t) = 0.046 Myr −1 respectively.
The dominant 2.4 Myr cyclicity observed in our model output seems absent in the δ 13 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 −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)(26)(27)(28)(29)(30)(31)(32)(33)(34).This could be related to the ∼ 1 Myr obliquity cycle.

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 2 (Fig. S10) are larger than the temperature estimates based on the δ 18 O composite data (which would be ∼ 1 • C).This is because the carbon cycle model (not a climate model) relates temperature to CO 2 directly through an input climate sensitivity of 5 • C per CO 2 doubling, whereas in reality temperature is a function of a number of variables, including but not limited to pCO 2 .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 2 , 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 δ 18 O records (and its associated temperature) therefore does not necessarily imply that the 2.4 Myr cycle cannot be dominant in pCO 2 .
Unfortunately, even the available pCO 2 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 2 proxy records are required to establish a possible 2.4 Myr cyclicity.Since the resulting model pCO 2 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 3 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 δ 13 C response is overall slower compared to data records.For example, the lag of δ 13 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 δ 13 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 δ 13 C and δ 18 O variability on eccentricity timescales.

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 palaeorecord, as modelled δ 13 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 2 .Ideally this would require highresolution 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.

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 δ 13 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 is primarily controlled by sediment accumulation (Berner, 1982;Hedges and Keil, 1995;Berner, 2006), the direct topdown 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 2 .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 2 spectrum was dominated by the 2.4 Myr eccentricity could be evaluated using long, high-resolution proxy records of pCO 2 or other carbonate chemistry parameters.
Code and data availability.The LOSCAR model source code is available upon request to loscar.model@gmail.com.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 (Pälike et al., 2006b;Tian et al., 2013;Holbourn et al., 2015b).The R code that was used to assemble the composite is available in the above-mentioned GitHub repository.
Author contributions.Modelling work was conducted by MJC and IJK.Data compilation was generated by IJK.All authors provided feedback on the paper.
Competing interests.The authors declare that they have no conflict of interest.

Figure 1 .
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 δ 18 O (b) and δ 13 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 −1 ).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.

Figure 3 .
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.

Figure 4 .
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 −1 , 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 −δ 13 C is shown, such that a phase lag (positive values) represents a lag of δ 13 C minima with respect to ETP maxima.

Figure 5 .
Figure 5. (a) δ 13 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 δ 13 C and detrended δ 13 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.

Figure 6 .
Figure 6.Evolutive multitaper harmonics analysis (c) of the data composite δ 13 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).

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

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