The Plio-Pleistocene climatic evolution as a consequence of orbital forcing on the carbon cycle

Since the discovery of ice ages in the 19th century, a central question of climate science has been to understand the respective role of the astronomical forcing and of greenhouse gases, in particular changes in the atmospheric concentration of carbon dioxide. Glacial–interglacial cycles have been shown to be paced by the astronomy with a dominant periodicity of 100 ka over the last million years, and a periodicity of 41 ka between roughly 1 and 3 million years before present (Myr BP). But the role and dynamics of the carbon cycle over the last 4 million years remain poorly understood. In particular, the transition into the Pleistocene about 2.8 Myr BP or the transition towards larger glaciations about 0.8 Myr BP (sometimes referred to as the mid-Pleistocene transition, or MPT) are not easily explained as direct consequences of the astronomical forcing. Some recent atmospheric CO2 reconstructions suggest slightly higher pCO2 levels before 1 Myr BP and a slow decrease over the last few million years (Bartoli et al., 2011; Seki et al., 2010). But the dynamics and the climatic role of the carbon cycle during the Plio-Pleistocene period remain unclear. Interestingly, the δ13C marine records provide some critical information on the evolution of sources and sinks of carbon. In particular, a clear 400 kyr oscillation has been found at many different time periods and appears to be a robust feature of the carbon cycle throughout at least the last 100 Myr (e.g. Paillard and Donnadieu, 2014). This oscillation is also visible over the last 4 Myr but its relationship with the eccentricity appears less obvious, with the occurrence of longer cycles at the end of the record, and a periodicity which therefore appears shifted towards 500 kyr (see Wang et al., 2004). In the following we present a simple dynamical model that provides an explanation for these carbon cycle variations, and how they relate to the climatic evolution over the last 4 Myr. It also gives an explanation for the lowest pCO2 values observed in the Antarctic ice core around 600–700 kyr BP. More generally, the model predicts a two-step decrease in pCO2 levels associated with the 2.4 Myr modulation of the eccentricity forcing. These two steps occur respectively at the PlioPleistocene transition and at the MPT, which strongly suggests that these transitions are astronomically forced through the dynamics of the carbon cycle.


Introduction
The idea that the orbital parameters of the Earth may influence climate has a long history, linked mostly to the development of theories of ice ages (e.g.Paillard, 2015).But it is clear from geological records that astronomical climatic variations have occurred throughout the Earth's history, with or without ice being present on Earth.It is therefore certain that the astronomical parameters are influencing climate through other mechanisms than the growth and decay of ice sheets.This is well known concerning records of monsoons or records of low-latitude precipitations, which are strongly influenced by precession.A very illustrative example is given by the Mediterranean sapropels (Lourens et al., 1996;Hilgen et al., 1999) which are used to calibrate the 40 Ar / 39 Ar radiochronometers (Kuiper et al., 2008).Similarly, a 400 kyr oscillation is observed in the δ 13 C of the foraminifera recovered from marine records, throughout most of the Cenozoic (Pälike et al., 2006;Cramer et al., 2003;Sexton et al., 2011;Billups et al., 2004;Wang et al., 2010).This oscillation is present in the benthic records, but also in many plank-Published by Copernicus Publications on behalf of the European Geosciences Union.tic ones, suggesting that these δ 13 C variations are linked to global ocean δ 13 C changes.This persistent oscillation was recently used to reconstruct the evolution of the Earth's carbon over the last 100 Myr (Paillard and Donnadieu, 2014).A key difficulty is to understand the dynamics of this cycle.In particular, during the last million years these oscillations appear to stretch and the relationship with eccentricity becomes less clear (e.g.Wang et al., 2004Wang et al., , 2010)), as illustrated in Fig. 1.
Before 1 Myr BP when ice sheets remained medium sized, the cyclicity appears locked to eccentricity, with high eccentricity values associated with decreasing or low values in δ 13 C.This phase relationship appears consistent with earlier time periods, with the chronology of Cenozoic marine cores being sometimes based on the association of high eccentricity and low δ 13 C values (e.g.Pälike et al., 2006;Cramer et al., 2003).A simple deduction is that, most probably, the dynamics behind this oscillation are essentially stable and linked to the astronomical forcing before 1 Myr BP, but it is strongly disturbed by the large Quaternary glaciations afterwards.This observation has major implications on the possible mechanisms, as we will see further on.
There is no consensus on the cause of these δ 13 C oscillations, but monsoons or the associated low-latitude precipitations are known to respond to precessional forcing, and therefore to be modulated by the 400 kyr eccentricity cycles.Still many factors may contribute to the evolution of the carbon cycle on these timescales, like erosion, vegetation dynamics, ocean biogeochemical or dynamical changes.It was therefore suggested that the δ 13 C cycles could be caused by the modulation of weathering in monsoonal regions (Pälike et al., 2006) or by ecological shifts in marine organisms, possibly linked to nutrient availability (Wang et al., 2004;Rickaby et al., 2007).It is worth emphasizing that during the last million years, if the link with eccentricity is less obvious, there are clear indications that these δ 13 C shifts are associated with major changes in the Earth carbon cycle.For instance, carbonate deposition exhibits major changes, well correlated with these global δ 13 C changes (Bassinot et al., 1994;Wang et al., 2004), and the record of atmospheric pCO 2 from Antarctic ice cores also shows a 10 to 20 ppm long-term modulation with a minimum level around 0.6-0.7 Myr BP and a maximum around 0.3-0.4Myr BP (Lüthi et al., 2008) in phase with the long-term carbonate preservation cycle.A mechanistic modelling of these 400 to 500 kyr cycles is therefore a critical missing element in our understanding of climate-carbon evolution over the Plio-Pleistocene period.
With a simple ocean box model (Russon et al., 2010), it was shown that silicate weathering alone could not account for the simultaneously observed rather large δ 13 C changes (> 0.4 ‰) and rather small pCO 2 variations (< 20 ppm) in this frequency band during the last million years.Furthermore, with silicate weathering only, the model-predicted phase relationships were also inconsistent with observations of δ 13 C, carbonate deposition and pCO 2 .Changes in organic matter fluxes are therefore a necessary ingredient in order to account for the observed rather large δ 13 C changes.A possible mechanism could therefore be linked to ocean organic matter burial, associated with changes in nutrient supply or ecological shifts (Rickaby et al., 2007).But it is then very difficult to explain why this mechanism would change drastically with the occurrence of major glaciations, as suggested by Fig. 1.We will therefore build our model on a different perspective, involving a more direct link between monsoons and organic matter burial, that should be strongly affected by sea level changes.
Organic matter burial takes place mostly on the continental shelves.Recent reassessments of riverine carbon fluxes to the ocean have emphasized the role of the erosion of continental organic carbon in the overall balance (e.g.Galy et al., 2007;Hilton et al., 2015).When investigating the influence of monsoons on the carbon cycle, it is natural to have a closer look at river discharges in monsoonal areas.Carbon budgets on major present-day erosional systems have provided some contrasted results, with riverine organic matter being either a net carbon source for the ocean (Burdige, 2005) or a net sink through organic carbon burial in sedimentary fans (Galy et al., 2007).The first study was mostly based on the Amazon basin, while the second estimation is from the Himalayan system.The differences are likely linked to different river basin configurations and different sedimentary deposition dynamics.This dramatically highlights the impact of geomorphology on terrestrial organic carbon burial, and suggests that the long-term global balance might be different in a context of large glacial-interglacial sea level variations like the last million years, when compared to earlier periods with much smaller sea level changes.Our conceptual model is therefore built on the impact of monsoon-driven terrestrial organic matter burial on the global carbon cycle.

Conceptual model
We are interested in the evolution of the global Earth carbon, that is, the carbon content of the atmosphere, the ocean and the biosphere, which amounts today approximately to C = 40 000 PgC (petagrams of carbon, i.e. 10 15 gC).This evolution results from possible imbalances between the volcanic inputs V , the oceanic carbonate deposition flux D associated with silicate weathering and its alkalinity flux to the ocean W , and the organic carbon burial B.  Bartoli et al., 2011) and alkenone isotopes (pink and blue lines for the min and max envelope, from Seki et al., 2010).δ 13 C in cores 1143 (red: Wang et al., 2004); 849 (blue: Mix et al., 1995); 846 (green: Shackleton et al., 1995).The same δ 13 C records filtered at 400 kyr (bandpass = 2.5 Myr −1 ), eccentricity (grey, from Laskar et al., 2004) and filtered eccentricity (black).
and runoff, and transforms it into a HCO − 3 that finally reaches the ocean.When considering the "global" Earth surface budget C which includes the ocean and atmosphere, W has therefore no direct effect on C and does not appear in Eq. (1a) for dC/dt, but only as a source of alkalinity in Eq. (1b).On timescales larger than several millennia, if we assume that the oceanic calcium concentration does not change significantly over the last few millions of years, carbonate compensation will restore the oceanic carbonate content.Therefore, to first order, we can write Solving for D, this leads to the long-term evolution equation for carbon: For simplicity, we will assume that the main stabilizer of the carbon system is the silicate weathering, with a fixed relaxation time τ C , i.e.W = C/τ C .Solving the present-day equilibrium with δ 13 Eq = 0 ‰ as a typical value for carbonates, we easily deduce typical equilibrium values for the fluxes: B 0 = V /5; C Eq = (8/5)τ C V = 40 000 PgC.If we assume a relaxation time τ C of 200 kyr (Archer et al., 1997), we obtain V = (5/8)C Eq /τ C = 125 TgC yr −1 and B 0 = 25 TgC yr −1 .For a larger value τ C = 400 kyr (Archer, 2005), we would get V = 62 TgC yr −1 .There is no consensus on the actual total carbon emissions from volcanism (including all aerial and submarine sources), but these values for V (or τ C ) span more or less the range of current estimates from about 40 to 175 TgC yr −1 (Burton et al., 2013).
It must be stressed that B stands for all organic carbon fluxes, whether they correspond to organic carbon burial (positive contributions to B) or to organic matter oxidation (negative contributions to B).While the long-term average equilibrium value B 0 needs to be positive to account for the isotopic balance as shown above, this is not necessarily always the case for the instantaneous values of B, as we will illustrate in what follows with the astronomical forcing.Indeed, B represents a sum of positive and negative terms whose individual absolute magnitudes are much larger than the long-term net value B 0 .For instance, the oxidation of petrogenic organic carbon alone will contribute negatively to B, with a magnitude that may be as large as 40 TgC yr −1 (Blair et al., 2003).
The isotopic 13 C budget can be written as where δ 13 C is the isotopic composition of ocean carbon, δ 13 V the isotopic composition of the volcanic carbon input, δ 13 B the isotopic composition of organic matter and δ 13 D the isotopic composition of marine carbonates.This can be rewritten as If we neglect isotopic fractionation during carbonate precipitation (in other words, δ 13 D = δ 13 C) and more generally during carbonate compensation, we finally obtain In the following we will assume a constant −5 ‰ volcanic source δ 13 V , as well as a constant −25 ‰ organic matter value δ 13 B (e.g.Porcelli and Turekian, 2010).
In order to translate the total carbon content C into an equivalent pCO 2 level, we will use a simple scaling.Indeed, if we assume, to first order, that C may represent the carbon content of a well-mixed ocean, then from chemical equilibrium pCO 2 should be proportional to HCO − 3 2 / CO 2− 3 .After carbonate compensation (i.e.assuming that CO 2− 3 remains constant) and considering that C is dominated by bicarbonates HCO − 3 under standard pH conditions, we end up with the approximate scaling that pCO 2 varies roughly as C 2 , or pCO 2 = 280 (C/40 000) 2 (in ppm).To reproduce a multi-million year trend, we need to add one explicitly in the weathering relaxation: W = C/τ C = ( C + C Eq − γ t)/τ C , with the coefficient γ set to 1.2 TgC yr −1 to obtain the desired pCO 2 levels at the start of the simulation, i.e. about 350 ppm at 4 Myr BP, according to current estimates (Bartoli et al., 2011;Seki et al., 2010).The model is integrated from an arbitrary initial condition at 5 Myr BP and the first 1 Myr is discarded.
In the following, we describe how carbon burial B should vary with monsoons, and what consequences these variations have on the total carbon content C as well as on carbonate isotopes δ 13 C.In order to represent the monsoon's response to astronomical forcing, we introduce a simple truncation of the precessional forcing: where e is the eccentricity and ω the climatic precession.
Indeed, soil erosion or sediment transport are dominated by intense events, not by the average climate.Such a nonlinear response can be mimicked in a simple way by the above expression that accounts only for positive monsoonal forcing, not for a negative one.Consequently, the model will be influenced by the amplitude modulation of the precessional forcing, i.e. the eccentricity.To avoid useless parameters, we furthermore introduce the normalization which results in a precessional forcing F (t) with amplitude 1 and 0 mean.
We implicitly account for a slow terrestrial organic carbon reservoir (soil) as "buried organic carbon".It is reasonable to assume that monsoon or enhanced precipitation will favour primary production and soil formation.But this recent soil together with older soils and with petrogenic organic carbon (Galy et al., 2008) will be eroded and transported to the ocean through enhanced river discharges.If the corresponding carbon is remineralized in the ocean without too much burial in the alluvial fan, the net perturbation of the burial flux is likely to be negative (i.e.net "old" soil erosion and remineralization).We will refer to this case as the "Amazon-like" situation, with the perturbation F (t) being subtracted to the baseline burial B 0 by writing B = B 0 − aF (t).In contrast, if most of the organic carbon is buried and preserved in the sediment, then the perturbation is likely to be positive, since it induces a net "recent" soil formation and burial.We call this the "Himalayan-like" situation, with now B = B 0 + aF (t).Before 1 Myr BP and the associated major sea level changes, the river fans and continental shelves should evolve mostly in a progradational way (see scheme in Fig. 3), a situation which a priori favours organic carbon remineralization, while aggradational situations are likely to be more frequent in the late Pleistocene, with therefore a possible temporary reversal of the organic carbon burial.

Results
Our first simulations, with B = B 0 − aF (t), correspond to a perpetual "Amazon-like" situation.They correspond to experiment a (black lines) with no trend in the total carbon, and experiment b (blue lines), with an explicit linear trend in carbon.The value of the parameter a is chosen in order to obtain approximately the correct amplitude for these simulated 400 kyr oscillations (a = 50 TgC yr −1 ).Still, as can be seen in Fig. 2, we obtain a surprisingly good match between the simulated and observed δ 13 C, with overall very similar cycles.More specifically, the δ 13 C black and blue simulated curves are superimposed and almost undistinguishable, since the linear trend added to the carbon cycle has almost no impact on the δ 13 C.They are both most of the time within the range of observed values (grey curves).The two main exceptions occur at about 0.3 and 2.3 Myr BP, with the simulated δ 13 C being significantly too high.In experiment a (black lines), pCO 2 is oscillating around its equilibrium value, with two significant negative excursions occurring near 2.5 Myr BP and near 0.5 Myr BP.These lower values are directly linked to the ∼ 2.4 Myr modulation of eccentricity.Obviously, with fixed or periodic parameters, this model cannot simulate a long-term decreasing trend in carbon.When explicitly adding such a linear decreasing trend (experiment b, blue lines), the two minima described above become two decreasing steps.The first one, occurring around  (2004).Sea level curve LR04 (black line) from Lisiecki and Raymo (2005) used to compute the river incision z MIN defined as the previous sea level minima (blue line).The geomorphological variable s used from experiment c (red lines) relaxed to its prescribed maximum value s MAX ∼ z 3 MIN (black line).The orange shaded areas correspond to the aggradation regimes (i.e.s < 0.85s MAX ).Total carbon C rescaled as pCO 2 for experiments a (black, precessional forcing only), b (blue, similar experiment but with a linear trend in carbon), and c (red, using the geomorphological dynamics from Eq. 3).Carbon isotopic composition δ 13 C for experiments a (black), b (blue), and c (red).In grey, the min and max values of the 13 C records from Fig. 1.The 400 kyr filtered values of δ 13 C results (blue and red) together with the range of filtered records (grey).The 400 kyr filtered eccentricity as in Fig. 1.In order to obtain these results, we chose τ C = 200 kyr (Archer et al., 1997)  2.8 Myr BP, is coincident with the Plio-Pleistocene transition and the development of Northern Hemisphere glaciations.The second one near 0.8 Myr BP is coincident with the mid-Pleistocene transition (MPT) and the significant amplification of glaciations.Note that the timing of these two steps is directly linked to the astronomical forcing: it does not depend at all on the specifics of the trend that we used here.Two similar pCO 2 decreasing episodes are also seen in the data (Fig. 1), though it is difficult to associate them with a precise timing, due to the difficulties in accurately reconstructing pCO 2 from indirect proxies.
In order to account for the observed departure of the δ 13 C oscillations from a simple eccentricity forcing, we need to introduce a retroaction of Quaternary sea level changes onto the sedimentary dynamics of alluvial fans and continental shelves, and consequently onto organic carbon burial.As explained above, we will reverse the sign of our burial flux perturbation, and change it into B = B 0 + aF (t) when some conditions are met on the geomorphology of river outputs.In particular, it is necessary to account for a changing reservoir size that can be filled with sediments in an aggradational way.Indeed, at the first major sea level drop, rivers are incising though the river and fan bedrock, thus providing room for the accumulation of sediments loaded with organic carbon.This volume should be filled progressively with sedimentary organic carbon up to a point when further river incision, and consequent aggradation of sediment, no longer affect the global organic carbon but only move sedimentary carbon from one place to another.In other words, we will assume that the global "Himalayan-like" situation (i.e.net organic carbon burial) is only a transient situation, linked to the first occurrence of a sea level minimum.In order to illustrate this mechanism, we add a new equation for the slow geomorphological reservoir S for organic carbon in river beds or river fans.We define its maximal size S MAX from the observed sea level changes using the reference stack LR04 (Lisiecki  top) With small sea level changes, we assume that the dominant sedimentary regime is progradation, with rather small organic carbon burial in coastal areas.The net effect of precessional forcing is (old) soil erosion, therefore a net transfer of carbon to the ocean-atmosphere.(bottom) With large sea level changes during the late Quaternary, the dominant sedimentary regime can switch temporarily to aggradation just after major sea level drops and river incisions.During these transitory phases, the net effect of precessional forcing is reversed, with net organic carbon burial in river beds and river fans.and Raymo, 2005) In other words, the sedimentary organic carbon reservoir S grows at the pace of the above-mentioned astronomical perturbation F 0 (t) up to its maximal size S MAX .After a short transient period, this reservoir remains therefore equal to this maximum value S MAX in the absence of major sea level drops, as during the pre-Quaternary period.In contrast, for each significant sea level drop, S MAX increases abruptly and we start a new transient phase whose duration is linked to parameter b.When S is small compared to the maximal reservoir size S MAX , then the aggradational scheme is favoured, with river beds and deltaic net organic carbon burial.But when S is close to its maximum value, we switch back to a mostly sedimentation scheme, meaning that potential sea level changes will no longer affect net global organic carbon burial: Using this simple crude criterion, we obtain the results shown in Fig. 2 (experiment c, red lines).As expected, this simple model does switch from the background "Amazon-like" or progradational burial mode to a "Himalayan-like" or aggradational mode, after each significant sea level drop, and most notably at two time periods, the first one between 2.4 and 2.5 Myr BP (as a consequence of the Plio-Pleistocene transition) and the second and largest one between 350 and 650 kyr BP (as a consequence of the MPT).The start of these transient periods is directly linked to sea level drops, according to the LR04 forcing, while the duration of these transients is linked both to the 0.85 S MAX threshold and the b parameter, whose values are chosen to qualitatively better match the δ 13 C data.For the results shown in Fig. 2, b = (160 kyr) −1 .Indeed, this sedimentary switch mechanism allows for a much better agreement with measured δ 13 C around 0.3 and 2.3 Myr BP, while the first simulations were systematically too high at this time, as illustrated by the difference between the blue and red curves in Fig. 2. We also simulate correctly the δ 13 C maximum around 500 kyr BP and the occurrence of two broad "500 kyr" cycles over the last million years.With this burial mode switching mechanism, we are also able to predict an absolute minimum in carbon content, or long-term pCO 2 , around 600 kyr BP, in rather good agreement with the long-term trend of pCO 2 measured in Antarctic ice cores.Indeed, pCO 2 from the Dome C record is about 5 to 10 ppm lower before the MPT (between 400 and 800 kyr BP), which is also what we obtain in our experiment c.

Discussion
When variations in B, as determined by parameter a, are smaller than the baseline value B 0 , the model cannot reproduce the oceanic amplitude of δ 13 C observed in marine benthic records.The observed 400 kyr signal in δ 13 C records therefore requires major changes in the organic carbon burial, with almost no global net burial, but net oxidation episodes, during maxima of precessional forcing.This strong forcing therefore implies significant oscillations in the Earth carbon cycle for this time frequency, up to 4 or 5 % in total carbon content.This is translated here into 10 to 20 ppm variations of pCO 2 using our simple scaling, but it is very likely that these changes would be much larger, when accounting for interactions between pCO 2 and climate.Indeed, colder climates are more favourable to oceanic carbon storage, as observed during the last glacial cycles.According to this mechanism, in the ordinary sedimentary regime (progradation), we obtain changes in the carbon cycle with pCO 2 maxima and δ 13 C minima associated directly with eccentricity maxima.This is indeed consistent with long Cenozoic records (e.g.Pälike et al., 2006).
When we allow for changes in the sedimentary regime triggered by sea level changes, the model can also reproduce more peculiar features.Indeed, up to now it has been difficult to explain the last two long-term cycles observed in the marine δ 13 C, each being approximately 500 kyr long, with a maximum now (δ 13 Cmax-I), a well-marked maximum at about 500 kyr BP (δ 13 Cmax-II) and a previous one around 1000 or 1100 kyr BP (δ 13 Cmax-III).In the model described here, these two long oscillations are generated from the eccentricity forcing, but with an abrupt switch to aggradation mode at about 620 kyr BP caused by the sea level drop at MIS 16.This switch reverses the phase of the 400 kyr carbon oscillation during a few hundred thousands of years.Interestingly, this also induces a slight minimum in the carbon (or pCO 2 ) results, consistent with the observed low pCO 2 values observed in the Antarctic ice core around 600-700 kyr BP.This mechanism also allows for simulated marine δ 13 C in better agreement with data at about 2.4 Myr BP.
It has already been noted (Wang et al., 2004) that the climatic evolution over the last million years, in particular the MPT (about 0.8 Myr BP) and the Mid-Brunhes Event (MBE, about 0.4 Myr BP) is associated with the carbon isotopic maxima δ 13 Cmax-II and δ 13 Cmax-III.This is a strong indication of a possible causal link between the long-term well-recognized eccentricity forcing on the carbon cycle and the Plio-Pleistocene climatic evolution.There is therefore a strong incentive to build a mechanistic astronomical theory of the carbon cycle.But a prerequisite towards understanding this long-term precessionally forced carbon cycle and its climatic consequences is to explain the observed changes during the Quaternary, in terms of δ 13 C, and simultaneously in the atmospheric CO 2 levels (Lüthi et al., 2008).The model results outlined above are a first step in this direction.
As detailed above, the fact that the 400 kyr carbon isotope cycle is perturbed during the Pleistocene strongly points towards a major role for organic matter burial over continental shelf areas being affected by sea level changes.Obviously, this model is far too simple to represent faithfully the complexities of sedimentary dynamics in coastal areas, its consequences on organic matter preservation, on carbon cycle and ultimately on climate.Furthermore, we provide here no explanation for the prescribed multi-million-year decreasing trend in pCO 2 .There is unfortunately no clear consensus on the actual mechanisms involved, though this trend has been often attributed to long-term changes in continental weathering linked either to mountain uplift (Raymo and Ruddiman, 1992), to continental drift or mantle degassing rate (Lefebvre et al., 2013).Furthermore, we considered only sea level changes as a potential feedback on organic matter burial in coastal areas.Obviously, many other important climatic feedbacks would also play a role.For instance, increased temperature would probably reduce net primary production as a consequence of increased stratification, and therefore reduce organic carbon deposition in coastal sediments, but it would also decrease oxygen concentrations and consequently would favour organic matter preservation.Similarly, stronger monsoon events would enhance the delivery of nutrients to the continental shelves, and therefore biological productivity.This would in addition deliver more fine-grained clay minerals that are necessary to seal and preserve organic matter from oxidation.This would work opposite to our continental soil-carbon mechanism for which enhanced monsoons lead to more organic carbon oxidation in agreement with the isotopic records.But, as a proof of concept, our model is chosen as minimalistic as possible.It does not attempt to include all potentially important mechanisms.

Conclusion
Our basic assumptions are primarily based on recent reassessments of riverine organic carbon inputs to the ocean.With the above conceptual model, we demonstrate that simple mechanistic assumptions can account for the major patterns of the observed global evolution of carbon and carbon isotopes over this time period.First, enhanced precessional forcing linked to high eccentricity leads to more continental organic carbon being washed out and remineralized, and therefore a net decrease in overall organic carbon burial.Second, this mechanism is temporarily reversed following major sea level drops associated with glaciations.This model was built on the premises that changes in organic matter or petrogenic organic carbon fluxes are responsible for the 400 kyr oscillations observed in Cenozoic 13 C records, and that the large sea level variations occurring during the Quaternary strongly affect this process.Continental margins and sedimentary fans are a very likely key component, as illustrated by our simple conceptual model.But obviously, many complex processes are involved in the interactions between organic matter burial or oxidation, monsoons and sea level changes.The geomorphological mechanism described here is one possibility which allows us, for the first time, to account both for the persistent 400 kyr oscillation observed in 13 C records during the Cenozoic and its change during the last million years.It also suggests the occurrence of possibly significant CO 2 drops at about 0.8 Myr BP (mid-Pleistocene transition) and at about 2.8 Myr BP (Plio-Pleistocene transition), which would ultimately link the timing of these transitions to the astronomical forcing.Our model also provides a possible explanation for the puzzling shifted level in the CO 2 records associated with the MBE.

Figure 2 .
Figure2.From top to bottom: precessional forcing F 0 (t) = Max(0, −e sin ω) (black line) fromLaskar et al. (2004).Sea level curve LR04 (black line) fromLisiecki and Raymo (2005) used to compute the river incision z MIN defined as the previous sea level minima (blue line).The geomorphological variable s used from experiment c (red lines) relaxed to its prescribed maximum value s MAX ∼ z 3 MIN (black line).The orange shaded areas correspond to the aggradation regimes (i.e.s < 0.85s MAX ).Total carbon C rescaled as pCO 2 for experiments a (black, precessional forcing only), b (blue, similar experiment but with a linear trend in carbon), and c (red, using the geomorphological dynamics from Eq. 3).Carbon isotopic composition δ 13 C for experiments a (black), b (blue), and c (red).In grey, the min and max values of the 13 C records from Fig.1.The 400 kyr filtered values of δ 13 C results (blue and red) together with the range of filtered records (grey).The 400 kyr filtered eccentricity as in Fig.1.In order to obtain these results, we chose τ C = 200 kyr(Archer et al., 1997) or equivalently V = 125 TgC yr −1 .The trend (experiments b and c) is set to γ = 1.2 TgC yr −1 to induce a drift from about 350 to about 280 ppm.The amplitude of the organic matter burial perturbation (experiments a, b and c) is set to a = 50 TgC yr −1 .The filling rate of the sedimentary reservoir (experiment c) is set to b = (160 kyr) −1 .The model is integrated from an arbitrary initial condition at 5 Myr BP and the first 1 Myr is discarded.
Figure2.From top to bottom: precessional forcing F 0 (t) = Max(0, −e sin ω) (black line) fromLaskar et al. (2004).Sea level curve LR04 (black line) fromLisiecki and Raymo (2005) used to compute the river incision z MIN defined as the previous sea level minima (blue line).The geomorphological variable s used from experiment c (red lines) relaxed to its prescribed maximum value s MAX ∼ z 3 MIN (black line).The orange shaded areas correspond to the aggradation regimes (i.e.s < 0.85s MAX ).Total carbon C rescaled as pCO 2 for experiments a (black, precessional forcing only), b (blue, similar experiment but with a linear trend in carbon), and c (red, using the geomorphological dynamics from Eq. 3).Carbon isotopic composition δ 13 C for experiments a (black), b (blue), and c (red).In grey, the min and max values of the 13 C records from Fig.1.The 400 kyr filtered values of δ 13 C results (blue and red) together with the range of filtered records (grey).The 400 kyr filtered eccentricity as in Fig.1.In order to obtain these results, we chose τ C = 200 kyr(Archer et al., 1997) or equivalently V = 125 TgC yr −1 .The trend (experiments b and c) is set to γ = 1.2 TgC yr −1 to induce a drift from about 350 to about 280 ppm.The amplitude of the organic matter burial perturbation (experiments a, b and c) is set to a = 50 TgC yr −1 .The filling rate of the sedimentary reservoir (experiment c) is set to b = (160 kyr) −1 .The model is integrated from an arbitrary initial condition at 5 Myr BP and the first 1 Myr is discarded.

Figure 3 .
Figure 3. Simple scheme of the two different geomorphological dynamics considered here.(top)With small sea level changes, we assume that the dominant sedimentary regime is progradation, with rather small organic carbon burial in coastal areas.The net effect of precessional forcing is (old) soil erosion, therefore a net transfer of carbon to the ocean-atmosphere.(bottom) With large sea level changes during the late Quaternary, the dominant sedimentary regime can switch temporarily to aggradation just after major sea level drops and river incisions.During these transitory phases, the net effect of precessional forcing is reversed, with net organic carbon burial in river beds and river fans.
by finding the previous sea level minimum z MIN (i.e. the lower envelope) with the scaling S MAX ∼ z 3 MIN since it represents a volume of sediment (see Fig. 2): if S < S MAX : dS/dt = bF 0 (t) otherwise : S = S MAX .