Orbital tuning as an inverse problem - Wiley Online Library

Click Here

PALEOCEANOGRAPHY, VOL. 25, PA2203, doi:10.1029/2009PA001769, 2010

for

Full Article

Orbital tuning as an inverse problem: Chronology of the early Aptian oceanic anoxic event 1a (Selli Level) in the Cismon APTICORE A. Malinverno,1 E. Erba,2 and T. D. Herbert3 Received 23 March 2009; revised 25 November 2009; accepted 15 December 2009; published 4 May 2010.

[1] Orbital tuning, the process of fitting sedimentary cycles to orbital periodicities, can estimate with high resolution the timing and duration of key events in the geological record. We formulate here orbital tuning as the inverse problem of finding the variation in sedimentation rate that matches sediment cycles with orbital periodicities. Instead of obtaining a single best estimate, we apply a Bayesian formulation and define a probability distribution of sedimentation rate variations that result in powerful orbital periodicities. By sampling this distribution with a Monte Carlo method, we quantify the uncertainty in the inferred sedimentation rates due to uncertainties in the tuning periods and to components of the sedimentary signal that are unrelated to orbital cycles. The method is applied to the chronology of a 30 m interval in the Cismon APTICORE borehole (Southern Alps, Italy) that includes the early Aptian oceanic anoxic event 1a (Selli Level), estimated to last 1.11 ± 0.11 Ma (95% interval). The d 13C record shows a sudden negative shift of about −1‰ at the base of the Selli Level (22–47 ka) followed by a recovery to preshift values in ∼240 ka, consistent with a scenario where light carbon is quickly added and then flushed out of the ocean‐ atmosphere system. Citation: Malinverno, A., E. Erba, and T. D. Herbert (2010), Orbital tuning as an inverse problem: Chronology of the early Aptian oceanic anoxic event 1a (Selli Level) in the Cismon APTICORE, Paleoceanography, 25, PA2203, doi:10.1029/2009PA001769.

1. Introduction [2] Over the last few decades, time series analyses of Quaternary sediment sequences have conclusively shown periodicities linked to Milankovitch orbital cycles [e.g., Hays et al., 1976; Martinson et al., 1987]. This discovery sparked considerable interest in using orbital cycles to measure time in the geological record. Numerous studies found cycles in Mesozoic sediments that could be interpreted as the climate response to orbital cycles, and cyclostratigraphy started its development [e.g., Schwarzacher and Fischer, 1982; Herbert and Fischer, 1986; Fischer et al., 1990; D’Argenio et al., 2004; Hinnov and Ogg, 2007]. A key technique in cyclostratigraphy is orbital tuning, which is the process of matching orbital periodicities with sedimentary cycles in sequences that have variable sedimentation rates [e.g., Martinson et al., 1987; Muller and MacDonald, 2000; Herbert, 2001; Meyers and Sageman, 2007]. Tuning converts stratigraphic depth to time, and provides crucial information on the chronology of events recorded in the sequence.

1

Lamont-Doherty Earth Observatory, Columbia University, Palisades, New York, USA. 2 Dipartimento di Scienze della Terra, Università degli Studi di Milano, Milan, Italy. 3 Department of Geological Sciences, Brown University, Providence, Rhode Island, USA. Copyright 2010 by the American Geophysical Union. 0883‐8305/10/2009PA001769

[3] We present here an inverse method for orbital tuning that infers sedimentation rates by fitting orbital cycles to a climate‐related sediment property. A key feature of the method is that it accounts for uncertainties in the resulting chronology due to nonorbital components in the data and to uncertainty in the tuning frequencies. Uncertainty quantification is essential to test correlations between different sequences and to properly integrate data from different sources of chronological information (e.g., from a geomagnetic polarity timescale and from cyclostratigraphy). An effective integration can be achieved by weighing each piece of chronological information on the basis of the respective uncertainty. A similar approach to invert for sedimentation rate, named “average spectral misfit,” has been proposed by Meyers and Sageman [2007] and Meyers [2008]; the two methods will be compared later. [4] The method is applied to establish the chronology of a 30 m thick interval in the Cismon APTICORE borehole, drilled in a Lower Cretaceous (Barremian‐Aptian) sequence of the Southern Alps, northern Italy [Erba and Larson, 1998; Erba et al., 1999]. The interval studied includes the early Aptian oceanic anoxic event (OAE) 1a, which is one of the most extreme examples of global anoxia in the mid‐ Cretaceous greenhouse world [Schlanger and Jenkyns, 1976; Jenkyns, 1980, 1999; Arthur and Premoli Silva, 1982; Sliter, 1989; Bralower et al., 1994, 2002; Erba, 2004]. The lithological expression of the global OAE 1a in Italy is the “Livello Selli” [Coccioni et al., 1987], a low‐ CaCO3 interval containing black shales and radiolarian layers. We date with our inverse method the duration of the Selli Level, of the associated carbon isotopic anomaly, and

PA2203

1 of 16

PA2203

MALINVERNO ET AL.: ORBITAL TUNING AS AN INVERSE PROBLEM

PA2203

Figure 1. Well log and other data in the 8–38 m stratigraphic depth interval of the Cismon APTICORE borehole. The Schlumberger Formation Microimager (FMI) image shows low resistivities in dark colors and high resistivities in light colors. The FMI uncalibrated resistivity curve is computed from the horizontally averaged image value and compares closely with the lower‐resolution spherically focused resistivity curve (SFLU). Low‐resistivity intervals such as the Selli Level correspond to intervals of low bulk density and high porosity. The CaCO3 content is from infrared reflectance measurements above and below the Selli Level [Li et al., 2008] and from carbonate content measurements in the Selli Level. Here d13Ccarb after Erba et al. [1999] and isotopic stages C2–C7 according to the definition of Menegatti et al. [1998]. of biotic changes preceding, contemporaneous to, and following global oceanic anoxia. The dates we obtain are crucial for future modeling of the carbon cycle, of the carbonate system, and of biota evolution. [5] In addition, we obtain an estimate of the period of obliquity in the Early Cretaceous. The periods of obliquity and precession in the geological past are expected to be less than their present value due to a shorter Earth‐Moon distance and a shorter length of day [Berger et al., 1992; Néron de Surgy and Laskar, 1997; Laskar et al., 2004], and cyclostratigraphic analyses provide useful constraints to orbital calculations [Park and Herbert, 1987; Pälike and Shackleton, 2000; Lourens et al., 2001].

2. Background: OAE 1a [6] The organic C‐enriched sediments of OAE 1a are associated with a complex d13C anomaly that has been

recognized worldwide in pelagic [Weissert, 1989; Menegatti et al., 1998; Bralower et al., 1999; Erba et al., 1999; Bellanca et al., 2002; Price, 2003; Ando et al., 2008], neritic [Jenkyns, 1995], and terrestrial records [Gröcke et al., 1999; Jahren et al., 2001; Ando et al., 2002]. Major changes in marine biota further characterize the OAE 1a episode: high‐ resolution data indicate a speciation of calcareous plankton [Erba, 1994, 2004; Premoli Silva et al., 1999; Leckie et al., 2002] preceding and accompanying a considerable decrease in calcification of calcareous nannoplankton culminating in the nannoconid crisis (in the sense of Erba [1994]) as well as a drastic reduction in planktonic foraminifers [Premoli Silva et al., 1999; Leckie et al., 2002]. Evolutionary pulses, abrupt changes in abundance and diversification, and calcification crises of planktonic communities have been interpreted as the result of oceanic fertilization and excess pCO2, both presumably related to the emplacement of the

2 of 16

PA2203


PA2203

Figure 2. The weight fraction of CaCO3 [Li et al., 2008] and the logarithm of the FMI resistivity curve are highly correlated because carbonate‐rich intervals have lower porosity and therefore higher electrical resistivity than carbonate‐poor intervals. The CaCO3 and log resistivity curve have been standardized to have the same variance for comparison.

Ontong Java Plateau [Erba, 1994, 2004; Larson and Erba, 1999; Leckie et al., 2002; Weissert and Erba, 2004; Méhay et al., 2009]. This event could have triggered a large‐ scale dissociation of methane hydrates, causing a sudden negative shift in the d13C record near the base of OAE 1a [Jahren et al., 2001; Bellanca et al., 2002; Jenkyns, 2003].

3. Selli Level in the Cismon Borehole [7] We concentrate here on the interval 8–38 m in stratigraphic depth of the Cismon borehole (see Erba et al. [1999] for the conversion of well depth to stratigraphic depth). The continuous sequence in this interval includes a 5 m thick Selli Level, defined on a lithostratigraphic basis as the interval containing black shales with high organic carbon content (Figure 1). The Selli Level at Cismon is composed of grey to dark grey marly limestones alternating with radiolarian‐rich calcarenite horizons and black shale layers with low CaCO3 content. Below and above the Selli Level, the lithology consists of light grey to dark grey limestones. The CaCO3 content curve in Figure 1 is a composite of infrared reflectance data from Li et al. [2008] in the limestones above and below the Selli Level and of sparser actual measurements of CaCO3 content in the Selli Level. We plot the measured CaCO3 in the Selli Level because it displays clearly a carbonate‐poor interval (“CaCO 3 dissolution phase” in Figure 1) that is not evident in the infrared reflectance data. [8] Similarly to other Cretaceous pelagic limestone sequences in Italy [Fischer and Schwarzacher, 1984; Herbert and Fischer, 1986], the Cismon section displays distinct decimeter‐scale cycles that have been related to orbital periodicities [Herbert, 1992; Li et al., 2008]. These lithological cycles are clearly visible on electrical resistivity images obtained by the Schlumberger Formation Microimager (FMI) downhole tool (Figure 1). The FMI tool generates a detailed image with a horizontal and vertical resolution of 0.5–1 cm using 192 button electrodes mounted

on pads kept in contact with the borehole wall [Luthi, 2001]; the Cismon FMI data are uniformly sampled every ∼1.2 cm in stratigraphic depth. The color in the FMI image of Figure 1, where dark shades are low resistivities, closely matches the color variation of the rocks, where darker intervals contain higher fractions of organic‐rich clays. In the downhole log data of the Cismon borehole, the Selli Level is an interval of low electrical resistivity, low bulk density, and high porosity compared to the limestone above and below (Figure 1). [9] For our tuning exercise, we use the logarithm of the FMI resistivity curve computed from the horizontal average of the image data (Figure 1). The logarithms of the measured FMI resistivities track closely the measured CaCO3 content in the Cismon borehole (Figure 2). Locklair and Sageman [2008] also used FMI data for cyclostratigraphic analysis and observed a similar relationship between resistivity and CaCO3 content in a chalk‐marl Cretaceous sequence (Western Interior basin, U.S.A.). The likely explanation for this correlation is that carbonate‐rich intervals tend to be more cemented and have lower porosity compared to intervals that contain less CaCO3 and more clay minerals. The bulk electrical resistivity of a sedimentary rock is dominated by its content of pore water, which is far more conductive than sediment grains. Experiments show that the electrical resistivity of water‐saturated sedimentary rocks is proportional to a negative power of their porosity [Archie, 1942; Ellis and Singer, 2007]. Therefore, the logarithm of the resistivity is inversely proportional to porosity, and resistivity closely tracks sedimentary cycles of varying CaCO3 content. [10] An overall chronology of the Cismon section has been established by Channell et al. [1979, 2000], who identified and correlated polarity chrons M0 to M9 to the M sequence geomagnetic polarity timescale [Channell et al., 1995]. The Selli Level, however, is in the long normal Cretaceous interval above the top of M0. Channell et al. [2000] estimated sedimentation rates of 9–12 m/Ma in the carbonate‐rich intervals immediately below the Selli Level, but the low CaCO3 content of the Selli Level suggests lower

3 of 16

PA2203


Figure 3. (a) A sedimentation rate model with constant rates in each layer corresponds to (b) an age model consisting of linear segments. If the data only constrain sedimentation rates, the age model has no absolute age information and only gives a “floating” timescale, as indicated by the gray arrows. sedimentation rates. The chronology of the Selli Level can be determined by cyclostratigraphy, and Herbert [1992] obtained a duration of about 1.2 Ma for the Selli Level in the Cismon outcrop by counting 100 kyr bed bundles. In a recent study, Li et al. [2008] identified precession cycles in a spectral analysis of cm‐scale CaCO3 data from the Cismon core, and estimated a similar duration of 1.23 Ma for the Selli Level. These durations imply an overall sedimentation rate of about 4 m/Ma for the 5 m thick Selli Level.

4. Orbital Tuning as an Inverse Problem 4.1. Background [11] Orbital tuning consists of estimating the variation in sedimentation rates by matching orbital signals. In principle, tuning can be carried out by letting sedimentation rates vary freely from point to point to maximize the amplitude of cycles at orbital frequencies. The problem with this kind of tuning is that fluctuations in sediment properties that are unrelated to orbital cycles can be mapped to orbital frequencies by arbitrary variations in sedimentation rates. Reliable tuning should be minimal, in the sense that local fluctuations in sedimentation rate should be kept as small as possible [Muller and MacDonald, 2000]. Minimal tuning can be accomplished using a simple sedimentation rate model with a few layers of constant sedimentation rates, which corresponds to an age model composed of linear segments (Figure 3). We choose here the sedimentation rate model of Figure 3 because it has a few parameters and because it allows for step changes in sedimentation rates (as expected at the boundaries of the Selli Level). The method we describe, however, can be used with any parameterized sedimentation rate model. [12] There are two main ways to infer a sedimentation rate/age model from sequences of climate‐related sediment properties. The first approach is to obtain an absolute age model by correlating sedimentary sequences with orbital

PA2203

cycle curves (e.g., insolation). This approach has been widely applied in the Cenozoic, where the absolute timing of orbital cycles can be predicted reliably [Laskar et al., 2004]. The most recent version of the Neogene geological timescale [Lourens et al., 2004] is a compilation of numerous such individual cyclostratigraphic studies. Inverse methods to estimate age models by correlation are described by Martinson et al. [1982], Brüggemann [1992], Yu and Ding [1998], and Lisiecki and Lisiecki [2002]. [13] In the geological past beyond ∼50 Ma, the absolute timing of calculated orbital cycles is not expected to be accurate because of the chaotic evolution of the solar system [Laskar, 1999; Laskar et al., 2004]. The second approach to orbital tuning obtains sedimentation rate models by matching cycles in sedimentary records with characteristic orbital periodicities. This approach has been applied to the Mesozoic and beyond and results in “floating” age models [Hinnov, 2004] that need independent information to obtain absolute ages (Figure 3). The determination of sedimentation rates is typically made on the basis of a spectral analysis of the sedimentary sequence [e.g., Herbert and Fischer, 1986; Park and Herbert, 1987; Olsen and Kent, 1999; Meyers and Sageman, 2007]. The inverse method we describe is applicable to this spectrally based approach. [14] Orbital cycles used for spectrally based tuning include long eccentricity (with a present‐day period of ∼400 ka), short eccentricity (periods of ∼120 ka and ∼95 kyr), obliquity (∼41 ka), and precession (∼19–23 ka) [e.g., Hinnov, 2004]. In the geological past, the periods of some of these cycles are expected to differ from the present‐day value. Tidal dissipation decreases the rotation rate of the Earth and increases the Earth‐Moon distance, causing the obliquity and precession periods to be lower in the past [Berger et al., 1992; Néron de Surgy and Laskar, 1997; Pälike and Shackleton, 2000; Lourens et al., 2001; Laskar et al., 2004]. On the other hand, the periods of eccentricity should be stable in the Mesozoic [Berger et al., 1992], with the long eccentricity being the most reliable [Laskar et al., 2004]. 4.2. Bayesian Inversion [15] Solving an inverse problem consists of inferring values of model parameters on the basis of measured data. The parameters of our sedimentation rate model can be written as a vector m = (u, z), which is composed of a vector u of sedimentation rates and a vector z of depths of the layer interfaces. For example, the model of Figure 3 consists of seven parameters: four sedimentation rates and three layer interface depths. The measurements are a vector d of FMI log resistivities (Figure 1). A common way to implement this kind of inversion is to define an objective function that measures how well the data in d match orbital cycles for a given choice of the model parameters m. For example, the objective function may measure how much spectral power is concentrated in the frequency bands that should contain orbital signals [e.g., Meyers and Sageman, 2007]. The solution of the inverse problem is the sedimentation rate model that maximizes this objective function. [16] We also wish to quantify the uncertainty that is inherent to the solution. An already mentioned source of

4 of 16

PA2203


uncertainty in tuning is that the frequencies of orbital cycles in the geological past may not be accurately known. In addition, power spectra of sediment property sequences typically consist of peaks corresponding to orbital cycles superimposed on a background variability that has a red spectrum with more power at lower frequencies [Mann and Lees, 1996; Wunsch, 2004; Meyers et al., 2008]. Because of these background fluctuations, different sedimentation rate models may result in similar spectral peaks at orbital frequencies, similar values of the objective function, and a corresponding uncertainty in the solution. [17] To quantify the uncertainty of the inferred sedimentation rate model, we apply the Bayesian formulation that has been widely used in geophysical inverse problems [Tarantola and Valette, 1982; Jackson and Matsu’ura, 1985; Duijndam, 1988; Mosegaard and Tarantola, 1995]. In this formulation, the solution of the inverse problem is a posterior probability distribution that can be written from Bayes’s rule as pðmjdÞ ¼ const:pðmÞpðdjmÞ;

ð1Þ

where p(m∣d) is the posterior distribution of the model parameters in m given the data in d, p(m) is the prior distribution and p(d∣m) the likelihood function. The constant in equation (1) normalizes the posterior distribution and is not important in this context. 4.3. Outline of the Inversion Procedure [18] In the following, we first define the likelihood p(d∣m), which measures how well a given sedimentation rate model m matches sediment cycles with orbital tuning frequencies (section 4.4). The likelihood is initially determined for a single layer with a uniform sedimentation rate (section 4.5) and then for a model with several layers as in Figure 3 (section 4.6). The inversion procedure consists of two steps: (1) establish an initial sedimentation rate model (section 4.6) and (2) perturb the initial model with a Monte Carlo algorithm to quantify the uncertainty of sedimentation rates and layer interface depths (section 4.7). 4.4. Definition of the Likelihood Function [19] In the Bayesian formulation, the prior distribution accounts for what is known about the model parameters a priori and the likelihood function accounts for the information provided by the data. We use here simple uniform distributions for the prior (i.e., p(m) is constant between prescribed minimum and maximum values), so that the posterior distribution is fully described by the likelihood function. The likelihood is formally the probability of measuring the observed data in d for a given value of the model parameters m, and for Nf tuning frequencies it is pðdjmÞ ¼

Nf Z Y

’j ð f Þpðdjm; f Þdf ;

ð2Þ

j¼1

where ’j(f) is a probability distribution describing the uncertainty of the jth tuning frequency and p(d∣m, f) is the

PA2203

likelihood of a sinusoidal signal of temporal frequency f in a red noise background: Cd ð f Þ ; pðdjm; f Þ ¼ const: exp Cb ð f Þ

ð3Þ

where Cd(f) is the periodogram (the squared Fourier transform) of the data and Cb(f) the periodogram of the red noise background. Equation (3) is a simple extension of the formulas of Bretthorst [1988, 2001] for the likelihood of a sinusoidal signal in a background of Gaussian white noise. Definitions of the periodogram for evenly and unevenly sampled data, of the fitting procedure for the periodogram of the red noise background, and of the distributions of the tuning frequencies ’j(f) are in the Appendices. [20] The likelihood in equation (2) depends on the sedimentation rate model because the periodogram of the data is in terms of temporal frequency and needs to be calculated for data converted from depth to age using the parameters in m. The likelihood will be greater when the periodogram ratio in equation (3) is greater, i.e., the more the periodogram of the data rises over the background in the frequency bands of orbital cycles. 4.5. Likelihood of a Uniform Sedimentation Rate [21] To illustrate the calculation of the likelihood, consider the simplest possible sedimentation rate model: a uniform sedimentation rate in a given depth interval. Figure 4 shows the data and background periodograms of a 5 m interval in the FMI log resistivity data. The periodograms in Figures 4b and 4c are a function of spatial frequency fz, which for a given sedimentation rate u corresponds to a temporal frequency f ¼ ufz :

ð4Þ

[22] The periodogram ratios Cd(f)/Cb(f) computed for two possible values of sedimentation rate are in Figures 4d and 4e. A sedimentation rate that aligns peaks in the periodogram ratio with orbital frequencies will result in a large value of the likelihood. The orbital tuning frequencies are three eccentricities (periods of 404, 124, and 94.9 ka), obliquity (37 ka), and a broad interval for precession (centered on 20.3 ka). These tuning frequencies are obtained from the orbital solution of Laskar et al. [2004] at 120 Ma; for details, see Appendix C. The likelihood function computed as a function of sedimentation rate has a clear peak around u = 7.1 m/Ma. [23] The likelihood in Figure 4f is analogous to the “average spectral misfit” of Meyers and Sageman [2007] and Meyers [2008]. These two measures of fit are not identical; in particular, the average spectral misfit depends on the difference between observed and predicted peak orbital frequencies, while the likelihood is a function of how much above‐background spectral power is near orbital frequencies. In general, however, the higher the likelihood of a sedimentation rate, the lower the average spectral misfit will be (a low spectral misfit signifies good fit, and thus will be a

5 of 16

PA2203


PA2203

Figure 4. (a) Calculation of the likelihood of the sedimentation rate u for log resistivity data in a 5 m interval. (b) The periodogram of the data Cd(fz) (black line) and a fitted background periodogram Cb(fz) (gray line) and (c) the ratio Cd(fz)/Cb(fz). Spatial frequencies fz can be converted to temporal frequencies for different values of sedimentation rate u. (d) Here u = 3.5 m/Ma results in peaks of the periodogram ratio that are not aligned with the tuning frequencies, and the likelihood of u computed from equation (2) is small. Tuning periods are at the top of the plot, and the gray scale background shows the probability distributions of the tuning frequencies. (e) Here u = 7.1 m/Ma results in peaks of the periodogram ratio that match the 94.9 ka and 37 ka tuning frequencies, resulting in (f) a high likelihood. mirror image of the likelihood), and these two approaches are complementary. 4.6. Likelihood of a Sedimentation Rate/Age Model [24] We now extend the calculation of the likelihood to a sedimentation rate model such as that of Figure 3. To split the Cismon borehole interval into a stack of layers, we start from three layers that should have different sedimentation rates: the Selli Level and the limestones above and below. This initial subdivision, however, results in multimodal likelihoods, which suggest that the sedimentation rate is not constant in each layer. By splitting these three layers and adjusting the depths of the interfaces, we obtain the six‐layer model of Figure 5. The likelihood in each layer (Figure 5c) generally shows a single prominent maximum, which is attained when the 94.9 ka short eccentricity and the 37 ka obliquity cycles are both matched (Figure 5b). The maximum likelihood sedimentation rates in the Selli Level are 4–5 m/Ma and are lower than those in the surrounding limestones (up to

9 m/Ma), in agreement with previous studies [Herbert, 1992; Channell et al., 2000; Li et al., 2008]. [25] While the results in Figure 5 give the likelihood of sedimentation rates in each layer, we need to evaluate the likelihood of the whole sedimentation rate model. This calculation is illustrated in Figure 6. The model in Figure 6a transforms depth in the FMI log resistivity data to time, and the resulting periodogram of the data in temporal frequency is used to compute the likelihood of the sedimentation rate model from equation (2). The orbital tuning frequencies are short eccentricity (period of 94.9 ka) and obliquity (37 ka), which mostly affect the likelihood in each layer (Figure 5b); details are in Appendix C. [26] It should be noted that the calculation of the likelihood assumes that sedimentation was continuous throughout the depth interval spanned by the model. The 8–38 m interval in the Cismon borehole contains no major unconformities [Erba et al., 1999], but small gaps in sedimentation cannot be completely ruled out. Hiatuses cause a

6 of 16

PA2203


PA2203

Figure 5. Likelihood of the sedimentation rate u in six depth intervals of the Cismon APTICORE. (a) FMI log resistivity curve with data in the Selli Level rescaled to have the same mean and variance as in the rest of the section. (b) Image of the periodogram ratio Cd(fz)/Cb(fz) in each depth interval and (c) image of the likelihood of sedimentation rates. The dotted lines in Figure 5c show maximum likelihood sedimentation rates, and the dotted lines in Figure 5b mark the corresponding spatial frequencies for tuning periods of 404, 124, 94.9, 37, and 20.3 ka.

split in orbital spectral peaks [Meyers and Sageman, 2004], and in principle can be detected; we do not examine this problem here. [27] Although the sedimentation rate model of Figures 5c and 6a results in sizable spectral peaks at orbital frequencies, there may be alternative models that give similar results. The solution of the inverse problem is a posterior distribution of sedimentation rate models (equation (1)), and to quantify the posterior uncertainty we need to explore the range of models that give prominent orbital signals. 4.7. Markov Chain Monte Carlo Sampling [28] Rather than striving to find a best value for the parameters of the sedimentation rate model, we follow the inverse recipe advocated by Tarantola [2006]: consider many possible models consistent with a priori information and keep those that fit the data (in our case, result in powerful components at the orbital tuning frequencies). In practice, this exploration of the space of possible models can be efficiently

carried out by Markov chain Monte Carlo sampling [Gilks et al., 1996]. Applications of Markov chain Monte Carlo to geophysical inverse problems are described by Mosegaard and Tarantola [1995], Sen and Stoffa [1995], Sambridge and Mosegaard [2002], Malinverno [2002], Malinverno and Leaney [2005], and Sambridge et al. [2006]. [29] We apply the Markov chain Monte Carlo Metropolis algorithm, where a candidate model is obtained first by small modifications of the last model sampled, and the candidate is then accepted or rejected [Metropolis et al., 1953; Chib and Greenberg, 1995]. The sampling chain starts from the initial model of Figure 6a and proceeds by proposing a candidate model mcand obtained from the current model m by changing slightly the sedimentation rates and the depth of one of the layer interfaces. (The interfaces at the top and base of the Selli Level are kept fixed because they correspond to a distinct lithological change). We use the procedure described by Mosegaard and Tarantola [1995], where the generation of candidate models follows a random walk

7 of 16

PA2203


PA2203

Figure 6. (a) Calculation of the likelihood of a sedimentation rate model. The model transforms the depth coordinate in (b) the FMI log resistivity data to (c) time. (d) The periodogram Cd(f) of the data (black line) and the background periodogram Cb(f) (gray line) in temporal frequency. (e) The periodogram ratio is used to compute the likelihood of the sedimentation rate model (equation (2)); tuning periods are at the top of the plot, and the gray scale background shows the probability distributions of the tuning frequencies. that samples the prior distribution. The chain samples the candidate model with a probability of acceptance a defined as ¼ min

pðdjmcand Þ ;1 : pðdjmÞ

ð5Þ

If the candidate model has a likelihood greater than or equal to that of the current model, it is always accepted; otherwise, it is accepted with a probability equal to the ratio of the likelihoods. If the candidate is rejected, the chain samples again the current model. This simple algorithm will return a sample distributed as in the posterior distribution [Mosegaard and Tarantola, 1995]. [30] Besides the parameters of the sedimentation rate model, there is also uncertainty about the period of obliquity in the Early Cretaceous. For example, while at 120 Ma the orbital model of Laskar et al. [2004] gives an obliquity period of 37 ka, Berger et al. [1992] estimate a longer period of 38.6 ka. These differences are due to different assumptions about tidal dissipation in the geological past.

To account for this uncertainty, we allow the period of obliquity to vary in an interval that spans a range of possible Early Cretaceous values (34–39 ka). In practice, this is done by changing the period of obliquity in addition to the parameters of the sedimentation rate model during the Monte Carlo sampling.

5. Monte Carlo Sampling Results and Discussion 5.1. Sedimentation Rate and Age Model [31] The results of sampling sedimentation rate and age models are in the images of Figure 7, which are two‐ dimensional histograms of 50,000 samples. As noted earlier, the age models obtained by spectrally based tuning give a “floating” timescale that needs independent ties to absolute age. All the sampled age models in Figure 7b have been constrained to have an age of 121 Ma at the base of magnetic polarity chron M0. This age is one of the constraints used by Channell et al. [1995] in their CENT94 geomagnetic polarity timescale, and is supported by recent radio-

8 of 16

PA2203


PA2203

Figure 7. Results of Monte Carlo sampling: (a) image showing a two‐dimensional histogram of the sampled sedimentation rate models and (b) image showing the corresponding age models, which all have an age of 121 Ma at the base of M0. The continuous line in Figure 7b is the mean of the sampled age models. Independently determined ages shown for comparison are for the base and top of M0 (triangles; CENT94 timescale of Channell et al. [1995]) and for selected nannofossil events (circles; Table 1). metric dating [He et al., 2008]. To match the GTS2004 timescale [Gradstein et al., 2004; Ogg et al., 2008], which has an age of 125 Ma for the base of M0, 4 Ma should be added to all absolute ages estimated here. The solid line in Figure 7b is the mean of the sampled age models, and is available for download in the auxiliary material.1 Ages of first occurrence of selected nannofossils [Bralower et al., 1995; Channell et al., 1995, 2000] and of the “nannoconid crisis” [Erba, 1994, 1996, 2004; Channell et al., 2000] that were obtained independently from our tuning (Table 1) are consistent with the sampled age models. [32] It should be stressed that the uncertainties in the age model of Figure 7 and in the rest of this paper assume an exact age of 121 Ma at the base of M0. If an uncertainty were assigned to this age tie, all age uncertainties would correspondingly increase. The total uncertainty can be easily computed by summing the variances due to uncertainties in sedimentation rates and in the age tie, which were obtained independently. For example, suppose that the tuning results gave an age uncertainty of ±300 ka for a given horizon when the age of the base M0 was assumed to be exact. If the actual uncertainty for the age of the base M0 was ±400 ka, 1 Auxiliary materials are available at ftp://ftp.agu.org/apend/pa/ 2009pa001769.

the age uncertainty of the horizon would increase to (3002 + 4002)1/2 = 500 ka. On the other hand, the uncertainty in the duration of an interval (e.g., the duration of polarity chron M0) is not affected by the age tie, because it only depends on uncertainties in the sedimentation rate model. [33] Major changes in sedimentation rates occur at base and top of the Selli Level as well as in the underlying Maiolica and overlying Scaglia Variegata (Figure 7). Specifically, a considerable decrease in sedimentation rate from 9 to 6.5 m/Ma correlates with the start of the nannoconid decline [Channell et al., 2000; Erba and Tremolada, 2004] and an increase in the amplitude of resistivity fluctuations (Figure 1). The base of the Selli Level is marked by a drop in sedimentation rate to 4.8 m/Ma and minimum values (down to 4 m/Ma) are recorded in its upper part. The lower part of the Selli Level, however, contains a carbonate dissolution interval and would be expected to have the lowest sedimentation rates if the other sedimentary inputs were the same. A possible explanation is that the drop in carbonate was compensated by an increased burial of organic carbon‐ rich black shale [Erba et al., 1999]. Moreover, Tejada et al. [2009] documented a gradual increase in the 187Os/188Os ratio in the coeval portion of the Selli Level at Gorgo a Cerbara, that they interpret as due to an increase in global weathering rates. Consequently, an increase in the fine‐

9 of 16

PA2203


PA2203

Table 1. Estimated Ages of Nannofossil Events in the Literature Event

Stratigraphic Depth in Cismon APTICORE (m)

Age (Ma)

14.28 15.27 18.12 23.9 30.36 32.45 33.48

118.2 118.5 119 120.34 121 121.325 121.725

FO Braarudosphaera africanaa FO Corollithion achylosuma FO Eprolithus floralisa Start nannoconid crisisb FO Nannoconus truittiia FO Rucinolithus irregularisa FO Flabellites oblongusa a

After Bralower et al. [1995] and Channell et al. [1995, 2000]. After Erba [1994, 1996, 2004] and Channell et al. [2000].

b

grained terrigenous input at the beginning of OAE1a may have contributed to attenuate the drop in sedimentation rate in the lowermost part of the Selli Level. Sedimentation rates increase again above the Selli Level and reach values of 8 m/Ma where marly limestones prevail (from 12 m upwards). 5.2. Chronology of Selli Level [34] Ages and durations (and their uncertainties) for the Selli Level in the Cismon borehole listed in Table 2 are obtained from a histogram of the sampled sedimentation rate and age models. For example, Figure 8 shows the histogram of sampled durations of the Selli Level. The duration reported in Table 2 is the average of the sampled durations, and the 95% uncertainty bounds are computed from the interval containing the central 95% of the samples. The Selli Level in the Cismon borehole is estimated to last 1.11 ± 0.11 Ma, which is consistent with the results of Herbert [1992] and Li et al. [2008]. The base and top of the Selli Level are dated to 120.21 and 119.11 Ma, respectively (Table 2). The uncertainties in these dates are small (±40– 120 ka), but as noted earlier they assume no uncertainty in the age of 121 Ma assigned to the base of polarity chron M0. For comparison to other sections, it may be more useful to note that the time interval between the top of M0 and the base of the Selli Level is 298 ± 15 ka, somewhat less than

the 500 ka estimated by Herbert [1992] by counting 100 ka bed bundles in the Cismon outcrop. 5.3. Chronology of C Isotope Intervals [35] Figure 9a plots as a function of age the d13C measured on the carbonate fraction in the Cismon borehole [Erba et al., 1999] and the carbon isotope intervals C2–C7 originally defined by Menegatti et al. [1998]. Starting near the base of the Selli Level, d13C suddenly decreases by about −1‰ in interval C3. d13C then gradually increases back to its pre‐Selli Level value in interval C4, remains approximately constant during interval C5, and increases again by about 1‰ in interval C6. The 1.59 Ma long C7 interval (Table 2) is a positive excursion postdating OAE 1a worldwide. It denotes a major perturbation of the carbon cycle, presumably due to excess burial of organic matter although such deposits have not been documented yet. [36] The durations of the C isotope intervals C3–C7 obtained from Monte Carlo sampling are listed in Table 2. These durations are similar to those estimated by Li et al. [2008]; the main differences are that we obtain a shorter C4 (239 ka instead of 330 ka) and a slightly longer C3 (46.7 instead of 41 ka). Li et al. [2008] note that in the Cismon borehole the C3 interval in the d 13C of organic carbon [van Bruegel et al., 2007] covers a shorter time span, that they

Figure 8. Histogram of durations of the Selli Level sampled by the Monte Carlo algorithm. The mean sampled duration is 1.11 Ma. The dashed lines encompass the central 95% of the sampled durations. 10 of 16

PA2203


PA2203

Table 2. Chronology of Selli Level and of Selected Intervals Horizon

Stratigraphic Depth (m)

Age ± 95% Bounds

18.64 23.67

119.11 ± 0.12 Ma 120.21 ± 0.04 Ma

Depth Span (m)

Duration ± 95% Bounds

23.67–18.64 23.30–22.60 23.90–17.00 31.50–23.90 23.90–23.67 18.64–17.00 25.75–23.67 29.17–25.75 23.75–23.30 18.49–8.32 19.95–18.49 22.27–19.95 23.49–22.27 23.75–23.49 23.60–23.49

1.11 ± 0.11 Ma 137 ± 22 ka 1.44 ± 0.11 Ma 1.08 ± 0.06 Ma 33 ± 1.7 ka 302 ± 23 ka 298 ± 15 ka 490 ± 25 ka 83.9 ± 17 kaa 1.59 ± 0.07 Ma 349 ± 47 ka 510 ± 70 ka 239 ± 39 ka 46.7 ± 13.7 kaa 21.5 ± 3.5 ka

Top of Selli Level Base of Selli Level Interval Selli Level CaCO3 dissolution phase Nannoconid crisis Start nannoconid decline–start nannoconid crisis Start nannoconid crisis–base Selli Level Top Selli Level–end nannoconid crisis Top M0–base Selli Level Polarity Chron M0 Base C3 (d 13Ccarb)–start CaCO3 dissolution phase C7 (d 13Ccarb) C6 (d 13Ccarb) C5 (d 13Ccarb) C4 (d 13Ccarb) C3 (d 13Ccarb) C3 (d 13Corg)

a Duration uncertainty accounts for an uncertainty in the depth to the base of C3, which can be 17–34 cm thick depending on the sample chosen to define the base.

estimate as 27 ka. Our estimate for the duration of the C3 interval in the d13C of organic carbon is 21.5 ± 3.5 ka (Table 2). [37] The most striking feature of the d13C curve (Figure 9a) is the negative excursion at the base of the Selli Level that comprises isotopic intervals C3 and C4 and lasts a total of about 290 ka. The negative d13C shift that takes place over a few tens of ka in C3 can be explained by a sudden injection of light carbon in the ocean‐atmosphere reservoir. The mechanism typically invoked for this sharp negative carbon isotope excursion in OAE 1a is the dissociation of preexisting marine methane hydrates, whose carbon has a d 13C ≈ −60‰ [Jahren et al., 2001; Beerling et al., 2002; Jenkyns, 2003]. Alternative origins for the isotopically light carbon (e.g., from volatilization of C in terrestrial organic matter or from volcanic CO2) are viewed as unlikely because of the unrealistically large quantities of carbon required [Jahren et al., 2001]. On the other hand, volcanic CO2 emitted during the formation of the Ontong‐ Java large igneous province [Larson and Erba, 1999; Price, 2003; Weissert and Erba, 2004] could have warmed the ocean and triggered methane release from gas hydrates [Jahren et al., 2001; Bellanca et al., 2002; Jenkyns, 2003]. After a rapid introduction of isotopically light carbon, the negative d13C anomaly is expected to gradually disappear over a timescale that is 2–3 times the 100 ka residence time of carbon in the ocean and atmosphere [Berger and Vincent, 1986; Dickens, 2000; Beerling et al., 2002]. The 240 ka duration of the C4 interval in the Cismon d 13C record is consistent with this interpretation. 5.4. Chronology of Biotic Changes [38] The most dramatic change in marine biota associated with OAE 1a is a biocalcification crash, expressed by the nannoconid crisis in pelagic settings [Bralower et al., 1994, 1999; Erba, 1994, 2004; Erba and Tremolada, 2004] coeval

with a global shallow‐water carbonate demise [Föllmi et al., 1994, 2006; Wissler et al., 2003; Skelton et al., 2003; Weissert and Erba, 2004]. The nannoconid crisis lasted 1.44 ± 0.11 Ma; its onset preceded OAE 1a by 33 ± 1.7 ka and calcification resumed 302 ± 23 ka after the end of the Selli Level deposition. The first sign of calcification difficulty is dated at ∼121.3 Ma and shortly follows the onset of a major speciation episode among calcareous nannoplankton, which coincides with the base of nannofossil zone NC 6 (Figure 1) and is dated at ∼121.4 Ma (Figure 9b). The diversification of new coccolith taxa and the biocalcification crisis of the heavily calcified nannoconids are possibly linked to rapidly rising CO2 and might record the Ontong Java Plateau emplacement [Larson and Erba, 1999; Premoli Silva et al., 1999; Erba, 2004; Walczak et al., 2006]. [39] The decrease in nannoconid abundance occurred in successive steps in several sections [Erba, 1994; Erba and Tremolada, 2004], seemingly marking major phases of the Ontong Java Plateau construction. Repeated pulses of volcanogenic CO2 are thought to be responsible for the progressive biocalcification deficit. In this early phase, calcareous nannoplankton experienced conditions increasingly hampering calcite secretion, as shown by the origination of tiny coccolith species and a decrease of the heavily calcified nannoliths. Only in the core of the d13C negative spike, a carbonate dissolution phase is detected as a consequence of a major CCD shoaling. This acidification climax lasted 137 ± 22 ka and is evidenced by a minimum in calcareous nannoplankton and planktonic foraminiferal abundance and species richness [Premoli Silva et al., 1999; Leckie et al., 2002]. The lysocline then rapidly deepened during the last phase of C4, as suggested by the increase in CaCO3 and recovery of calcareous plankton. Rapid and massive carbon input results in such an oscillation of the lysocline in carbon cycle models that include a weathering feedback [e.g., Dickens, 2000].

11 of 16

PA2203


PA2203

the 38.6 ka in the orbital model of Berger et al. [1992] and the 39.2 ka estimated by Park and Herbert [1987] in a spectral analysis of the Aptian‐Albian sequence at Piobbico (central Italy). It should be stressed, however, that our estimate rests on the assumed short eccentricity period. If the short eccentricity period differed from 94.9 ka, the estimated obliquity period would change accordingly.

6. Conclusions

Figure 9. (a) Chronology of stable carbon isotope record and (b) stratigraphic events in the Cismon APTICORE borehole. The solid curve in Figure 9a uses the mean age model of Figure 7, which assumes an age of 121 Ma at the base of M0. The gray scale image of the d 13C record in Figure 9a is a two‐dimensional histogram obtained using all the sedimentation rate models sampled by the Monte Carlo algorithm. The d13Ccarb data are after Erba et al. [1999] and isotopic stages C2–C7 are according to the definition of Menegatti et al. [1998]. 5.5. Obliquity Period at 120 Ma [40] Figure 10 plots the likelihood of sedimentation rate models as a function of the obliquity periods sampled by the Monte Carlo algorithm assuming a short eccentricity period of 94.9 ka. If the data did not constrain the period of obliquity, the likelihood would not display a clear maximum in Figure 10. The data are informative about the obliquity period, and the sampled obliquity period with the maximum likelihood is 37.5 ka. The envelope of the highest likelihood values (dashed line in Figure 10) approximates the likelihood of the obliquity period, and shows that the 95% uncertainty of the estimated obliquity is at least ±1 ka. The 37.5 ka obliquity period estimated here is therefore consistent with the 37 ka obliquity period of the Laskar et al. [2004] orbital solution at 120 Ma, which is shorter than

[41] Orbital tuning is the process of estimating a chronology by matching sediment cycles with orbital periodicities in a sequence that has variable sedimentation rates. Orbital tuning can be viewed as an inverse problem: infer sedimentation rate models that give rise to prominent orbital cycles in the data. A sedimentation rate model is a parameterized variation of sedimentation rate with depth in the section, and this study uses a model consisting of a stack of layers with constant sedimentation rates. Any sedimentation rate model has a corresponding age model that is “floating,” in the sense that it needs independent age information to provide an absolute chronology. [42] This paper describes an inverse method to tune a sedimentary sequence that is based on a Bayesian formulation and uses a Markov chain Monte Carlo algorithm to sample the posterior distribution of sedimentation rate models that result in substantial orbital cycles in the data. The variation in the sampled sedimentation rate models measures the uncertainty in the tuning due to fluctuations in the sedimentary record that are not related to astronomical cycles and to uncertainties in the orbital periods. [43] The method is applied to establish the chronology of the “Selli Level” (OAE 1a) in the Cismon APTICORE borehole, drilled in an Early Cretaceous sequence of the southern Alps, Italy. The Monte Carlo sampling results give a duration for the anoxic sedimentation in the Selli Level of 1.11 ± 0.11 Ma. When the sampled age models are tied to an age of 121 Ma for the base of magnetic polarity chron M0 [Channell et al., 1995], the most likely age interval for the Selli Level deposition is 120.21 to 119.11 Ma. The d 13C data measured on the carbonate fraction in the Cismon borehole show a sudden negative shift of about −1‰ at the base of the Selli Level that takes place over an estimated 46.7 ± 13.7 ka. The duration of this negative shift in the d13C data for organic carbon is shorter, 21.5 ± 3.5 ka. After the negative shift, d13C returns to its preshift value in 239 ± 39 ka. These durations are consistent with a scenario where a large quantity of isotopically light carbon is added quickly to the ocean‐atmosphere system (e.g., by the dissociation of marine gas hydrates), followed by a slower recovery on a timescale comparable to the ∼100 ka residence time of carbon in the ocean‐atmosphere reservoir. [44] Within the negative d13C spike a carbonate dissolution phase of 137 ± 22 ka was presumably triggered by excess CO2, ocean acidification, and CCD shoaling. The biotic response to the OAE 1a environmental perturbation is documented by multiple changes in nannoplankton abundance and composition during a ∼2.6 Ma interval encompassing the Selli Level. While the nannoconid decline and

12 of 16

PA2203

PA2203


where w = 2p f is angular frequency and Rð!Þ ¼

N X

di cosð!ti Þ;

ðA2Þ

di sinð!ti Þ:

ðA3Þ

i¼1

I ð!Þ ¼

N X i¼1

Figure 10. Likelihood of sedimentation rate models sampled by the Monte Carlo algorithm as a function of the obliquity period. The sampled obliquity period with maximum likelihood (circled) is at 37.5 ka. The envelope of the highest likelihood values (dashed line) approximates the likelihood of the obliquity period.

The coefficients in equations (A2) and (A3) can be efficiently computed with fast Fourier transform algorithms [e.g., Press et al., 1992]. [48] If the data were acquired with a nonuniform sampling in depth or if the sedimentation rate model contains different sedimentation rates, the ages t1, t2, …, tN will not be uniformly sampled. In this case, the periodogram of the data in equation (3) is the Lomb‐Scargle periodogram [Lomb, 1976; Scargle, 1982] Cd ð f Þ ¼

2 1 R2LS ð!Þ ILS ð!Þ þ ; 2 CLS ð!Þ SLS ð!Þ

ðA4Þ

where crisis imply a hampered calcification, the coeval origination of coccoliths possibly testify to a biomineralization strategy of calcareous nannoplankton under excess CO2. Nannoconid abundance resumed 302 ± 23 ka after termination of global anoxia and might indicate CO2 drowning below threshold values hindering calcification of large nannoliths. [45] The period of obliquity in the geological past is expected to be less than the present value of 41 ka due to the effects of tidal dissipation, which progressively decreases the Earth’s rotation rate and increases the Earth‐Moon distance. We estimate the obliquity period in the early Cretaceous (∼120 Ma) by letting it be one of the variable parameters sampled by the Monte Carlo algorithm. For an assumed short eccentricity period of 94.9 ka, the sampling results give a best estimate of the obliquity period of 37.5 ka, which is close to the 37 ka in the orbital solution of Laskar et al. [2004].

RLS ð!Þ ¼

Cd ð f Þ ¼

1 2 R ð!Þ þ I 2 ð!Þ ; N

di cosð!ti Þ;

ðA5Þ

di sinð!ti Þ;

ðA6Þ

i¼1

ILS ð!Þ ¼

N X i¼1

CLS ð!Þ ¼

N X

cos2 ð!ti Þ;

ðA7Þ

sin2 ð!ti Þ;

ðA8Þ

i¼1

SLS ð!Þ ¼

N X i¼1

and t is a time delay chosen to make the sine and cosine basis orthogonal on the times ti at each frequency w as "P # N 1 1 i¼1 sinð2!ti Þ ¼ tan : PN 2 i¼1 cosð2!ti Þ

Appendix A: Periodogram of the Data [46] The periodogram of the data Cd(f) in equation (3) is a function of temporal frequency, and thus must be computed for data that are a function of age. Sediment property data (in our case, log electrical resistivities) in the vector d = d1, d2, …,dN are measured at depths z1, z2, …, zN. For the calculation of the periodogram, these depths must be converted to ages t1, t2, …, tN with a sedimentation rate model. [47] If the data were acquired with a uniform sampling in depth and if the sedimentation rate model is assumed constant in the whole depth interval, the ages ti will also be uniformly sampled. In this case, the periodogram of the data in equation (3) is the Schuster periodogram [e.g., Bretthorst, 1988; Muller and MacDonald, 2000]

N X

ðA9Þ

If the ages in the data are evenly sampled, t = 0, CLS(w) ≈ N/2, SLS(w) ≈ N/2, and the Lomb‐Scargle periodogram in equation (A4) becomes the same as the Schuster periodogram in equation (A1). Algorithms to compute the Lomb‐Scargle periodogram are given by Press et al. [1992] and Muller and MacDonald [2000].

Appendix B: Background Periodogram [49] The background periodogram Cb(f) in equation (3) is defined as

ðA1Þ

13 of 16

Cb ð f Þ ¼

a f02 þ f 2

b

ðB1Þ

PA2203

PA2203


Table C1. Orbital Frequencies at 120 Ma Used in Tuning Orbital Cycle Long eccentricity Short eccentricity Short eccentricity Obliquity Precession

Mean Period (ka)

Mean Frequency (cycles/Ma)

Standard Deviation (cycles/Ma)

404 124 94.9 37 20.3

2.475 8.065 10.54 27.03 49.26

0.11 0.28 0.14 0.18 5

and is similar to the power spectrum of an AR(1) process (autoregressive process of order 1), which has often been used to represent the red noise component of climate‐related records [e.g., Mann and Lees, 1996; Meyers et al., 2008]. The parameter a is related to the total variance of the red noise, and at high frequencies f ≫ f0 the background periodogram is proportional to f−2b. The exponent b always equals 1 in an AR(1) process, and we found that values appreciably different from 1 were needed to match the background in the Cismon log resistivity data. [50] The three parameters of the background periodogram (a, b, and f0) are obtained by fitting the periodogram of the data. Mann and Lees [1996] note that a statistically robust fit is required to avoid strong periodic signals affecting the estimate of the background. In our case, the orbital signals are not very powerful, and we obtained good results by fitting parameters to the logarithm of the data periodogram with nonlinear least squares optimization (Levenberg‐ Marquardt [see Press et al., 1992]) and by adjusting a so that the integral of the fitted periodogram matched the integral of the data periodogram.

Appendix C: Orbital Tuning Frequencies [51] The distributions ’j(f) (equation (2)) that quantify uncertainty in the tuning frequencies are taken to be normal distributions, and Table C1 lists the mean orbital frequencies and their standard deviations. The mean periods and frequencies in Table C1 are computed from the orbital solutions of Laskar et al. [2004] in a 5 Ma interval centered on 120 Ma.

In practice, we computed periodograms of the orbital solutions in 1 Ma windows, determined the frequencies of periodogram peaks corresponding to the orbital periods in each window, and computed the average peak frequencies in the full 5 Ma interval. The standard deviations of the eccentricity and obliquity frequencies in Table C1 are the residual standard deviations on the averages. Precession has three major spectral peaks whose relative amplitude varies with time, and the standard deviation of the precession frequency in Table C1 spans the range of these distinct spectral peaks. [52] In the whole 30 m interval studied, we used for tuning only the 94.9 ka short eccentricity and the obliquity because these two are the most powerful orbital cycles in the data, as also observed in other Cretaceous sequences in Italy [e.g., Park and Herbert, 1987]. The calculation of the likelihood of a sedimentation rate model requires computing the Lomb‐Scargle periodogram (Appendix A), which is more intensive the more tuning frequencies are used. Using only two tuning frequencies shortens considerably computing times for Monte Carlo sampling. To account for uncertainties of tidal dissipation in the geological past, the obliquity period is allowed to vary during sampling (see section 4.7). [53] Acknowledgments. The late Roger Larson had the foresight of acquiring the high‐resolution downhole log data that were vital to this study. Thanks to Jim Channell for helpful discussions and to Jill Weinberger, Trevor Williams, Sean Higgins, and Jim Murray for help with the downhole log data. Jerry Dickens, Stephen Meyers, and an anonymous reviewer gave many productive suggestions. This is LDEO contribution 7318.

References Ando, A., T. Kakegawa, R. Takashima, and T. Saito (2002), New perspective on Aptian carbon isotope stratigraphy: Data from d 13C records of terrestrial organic matter, Geology, 30, 227–230. Ando, A., H. Kaiho, H. Kawahata, and T. Kakegawa (2008), Timing and magnitude of early Aptian extreme warming: Unraveling primary d 18 O variation in indurated pelagic carbonates at Deep Sea Drilling Project Site 463, central Pacific Ocean, Palaeogeogr. Palaeoclimatol. Palaeoecol., 260, 463–476. Archie, G. E. (1942), The electrical resistivity log as an aid in determining some reservoir characteristics, Trans. Am. Inst. Min. Metall. Pet. Eng., 146, 54–62. Arthur, M. A., and I. Premoli Silva (1982), Development of widespread organic carbon‐rich strata in the Mediterranean Tethys, in Nature and Origin of Cretaceous Carbon‐Rich Facies, edited by S. O. Schlanger and M. B. Cita, pp. 8– 54, Academic, London. Beerling, D. J., M. R. Lomas, and D. R. Gröcke (2002), On the nature of methane gas‐hydrate

dissociation during the Toarcian and Aptian oceanic anoxic events, Am. J. Sci., 302, 28–49. Bellanca, A., E. Erba, R. Neri, I. Premoli Silva, M. Sprovieri, F. Tremolaca, and D. Verga (2002), Palaeoceanographic significance of the Tethyan ‘Livello Selli’ (early Aptian) from the Hybla Formation, northwestern Sicily: Biostratigraphy and high‐resolution chemostratigraphic records, Palaeogeogr. Palaeoclimatol. Palaeoecol., 185, 175–196. Berger, A., M. F. Loutre, and J. Laskar (1992), Stability of the astronomical frequencies over the Earth’s history for paleoclimate studies, Science, 255, 560–566. Berger, W. H., and E. Vincent (1986), Deep‐sea carbonates: Reading the carbon‐isotope signal, Geol. Rundsch., 75, 249–269. Bralower, T. J., M. A. Arthur, R. M. Leckie, W. V. Sliter, D. J. Allard, and S. O. Schlanger (1994), Timing and paleoceanography of oceanic dysoxia/anoxia in the late Barremian to early Aptian, Palaios, 9, 335–369. Bralower, T. J., R. M. Leckie, W. V. Sliter, and H. R. Thierstein (1995), An integrated Creta-

14 of 16

ceous microfossil biostratigraphy, in Geochronology, Timescales, and Stratigraphic Correlation, edited by W. A. Berggren et al., Spec. Publ. SEPM Soc. Sediment. Geol., 54, 65–79. Bralower, T. J., E. Cobabe, B. Clement, W. V. Sliter, C. L. Osburn, and J. Longoria (1999), The record of global change in mid‐Cretaceous (Barremian‐Albian) sections from the Sierra Madre, northeastern Mexico, J. Foraminiferal Res., 29, 418–437. Bralower, T. J., et al. (2002), Leg 198 summary, Proc. Ocean Drill. Program Initial Rep., 198, 1–148. Bretthorst, G. L. (1988), Lecture Notes in Statistics, vol. 48, Bayesian Spectrum Analysis and Parameter Estimation, Springer, New York. Bretthorst, G. L. (2001), Nonuniform sampling: Bandwidth and aliasing, in Bayesian Inference and Maximum Entropy Methods in Science and Engineering, edited by J. T. Rychert, G. J. Erickson, and C. R. Smith, AIP Conf. Proc., 567, 1–28.

PA2203


Brüggemann, W. (1992), A minimal cost function method for optimizing the age‐depth relation of deep‐sea sediment cores, Paleoceanography, 7, 467–487. Channell, J. E. T., W. Lowrie, and F. Medizza (1979), Middle and Early Cretaceous magnetic stratigraphy from the Cismon section, northern Italy, Earth Planet. Sci. Lett., 42, 133–166. Channell, J. E. T., E. Erba, M. Nakanishi, and K. Tamaki (1995), Late Jurassic–Early Cretaceous time scales and oceanic magnetic anomaly block models, in Geochronology, Timescales, and Stratigraphic Correlation, edited by W. A. Berggren et al., Spec. Publ. SEPM Soc. Sediment. Geol., 54, 51–63. Channell, J. E. T., E. Erba, G. Muttoni, and F. Tremolada (2000), Early Cretaceous magnetic stratigraphy in the APTICORE drill core and adjacent outcrop at Cismon (Southern Alps, Italy), and correlation to the proposed Barremian‐Aptian stratotype, Geol. Soc. Am. Bull., 112, 1430–1443. Chib, S., and E. Greenberg (1995), Understanding the Metropolis‐Hastings algorithm, Am. Stat., 49, 327–335. Coccioni, R., O. Nesci, M. Tramontana, F. C. Wezel, and E. Moretti (1987), Descrizione di un livello‐guida “radiolaritico‐bituminoso‐ ittiolitico” alla base delle Marne a Fucoidi nell’Appennino Umbro‐Marchigiano, Boll. Soc. Geol. Ital., 106, 183–192. D’Argenio, B., A. G. Fischer, I. Premoli Silva, H. Weissert, and V. Ferreri (Eds.) (2004), Cyclostratigraphy: Approaches and Case Histories, Spec. Publ. SEPM Soc. Sediment. Geol., 81, 311 pp. Dickens, G. R. (2000), Methane oxidation during the late Paleocene thermal maximum, Bull. Soc. Geol. Fr., 171, 37–49. Duijndam, A. J. W. (1988), Bayesian estimation in seismic inversion, part I: Principles, Geophys. Prospect., 36, 878–898. Ellis, D. V., and J. M. Singer (2007), Well Logging for Earth Scientists, Springer, Dordrecht, Netherlands. Erba, E. (1994), Nannofossils and superplumes: The early Aptian “nannoconid crisis,” Paleoceanography, 9, 483–501. Erba, E. (1996), The Aptian stage, Bull. Inst. R. Sci. Nat. Belg., 66, 31–43. Erba, E. (2004), Calcareous nannofossils and Mesozoic oceanic anoxic events, Mar. Micropaleontol., 52, 85–106. Erba, E., and R. L. Larson (1998), The Cismon APTICORE (Southern Alps, Italy): A “reference section” for the Lower Cretaceous at low latitudes, Riv. Ital. Paleontol. Stratigr., 104, 181–192. Erba, E., and F. Tremolada (2004), Nannofossil carbonate fluxes during the Early Cretaceous: Phytoplankton response to nutrification episodes, atmospheric CO2, and anoxia, Paleoceanography, 19, PA1008, doi:10.1029/ 2003PA000884. Erba, E., J. E. T. Channell, M. Claps, C. Jones, R. Larson, B. Opdyke, I. Premoli Silva, A. Riva, G. Salvini, and S. Torricelli (1999), Integrated stratigraphy of the Cismon APTICORE (Southern Alps, Italy): A “reference section” for the Barremian‐Aptian interval at low latitudes, J. Foraminiferal Res., 29, 371–391. Fischer, A. G., and W. Schwarzacher (1984), Cretaceous bedding rhythms under orbital control?, in Milankovitch and Climate: Understanding the Response to Astronomical Forcing, edited by A. Berger et al., pp. 163– 175, D. Reidel, Dordrecht, Netherlands.

Fischer, A. G., P. L. de Boer, and I. Premoli Silva (1990), Cyclostratigraphy, in Cretaceous Resources, Events and Rhythms, edited by R. N. Ginsburg and B. Beaudoin, pp. 139–172, Kluwer Acad., Dordrecht, Netherlands. Föllmi, K. B., H. Weissert, M. Bisping, and H. Funk (1994), Phosphogenesis, carbon‐ isotope stratigraphy and carbonate platform evolution along the northern Tethyan margin, Geol. Soc. Am. Bull., 106, 729–746. Föllmi, K. B., A. Godet, S. Bodin, and P. Linder (2006), Interactions between environmental change and shallow water carbonate buildup along the northern Tethyan margin and their impact on the Early Cretaceous carbon isotope record, Paleoceanography, 21, PA4211, doi:10.1029/2006PA001313. Gilks, W. R., S. Richardson, and D. J. Spiegelhalter (Eds.) (1996), Markov Chain Monte Carlo in Practice, Chapman and Hall, London. Gradstein, F. M., J. G. Ogg, and A. G. Smith (Eds.) (2004), A Geologic Time Scale 2004, Cambridge Univ. Press, Cambridge, U. K. Gröcke, D. R., S. P. Hesselbo, and H. C. Jenkyns (1999), Carbon‐isotope composition of Lower Cretaceous fossil wood: Ocean‐atmosphere chemistry and relation to sea‐level change, Geology, 27, 155–158. Hays, J. D., J. Imbrie, and N. J. Shackleton (1976), Variations in the Earth’s orbit: Pacemaker of the ice ages, Science, 194, 1121–1132. He, H., Y. Pan, L. Tauxe, H. Qin, and R. Zhu (2008), Toward age determination of the M0r (Barremian–Aptian boundary) of the Early Cretaceous, Phys. Earth Planet. Inter., 169, 41–48, doi:10.1016/j.pepi.2008.07.014. Herbert, T. D. (1992), Paleomagnetic calibration of Milankovitch cyclicity in Lower Cretaceous sediments, Earth Planet. Sci. Lett., 112, 15–28. Herbert, T. D. (2001), Orbitally tuned time scales, in Encyclopedia of Ocean Sciences, edited by J. Steele, S. Thorpe, and K. Turekian, pp. 2048–2054, Academic, San Diego, Calif. Herbert, T. D., and A. G. Fischer (1986), Milankovitch climatic origin of mid‐Cretaceous black shale rhythms in central Italy, Nature, 321, 739–743. Hinnov, L. A. (2004), Earth’s orbital parameters and cycle stratigraphy, in A Geologic Time Scale 2004, edited by F. M. Gradstein, J. G. Ogg, and A. G. Smith, pp. 55–62, Cambridge Univ. Press, Cambridge, U. K. Hinnov, L. A., and J. G. Ogg (2007), Cyclostratigraphy and the astronomical time scale, Stratigraphy, 4, 239–251. Jackson, D. D., and M. Matsu’ura (1985), A Bayesian approach to nonlinear inversion, J. Geophys. Res., 90, 581–591. Jahren, A. H., N. C. Arens, G. Sarmiento, J. Guerrero, and R. Amundson (2001), Terrestrial record of methane hydrate dissociation in the Early Cretaceous, Geology, 29, 159–162. Jenkyns, H. C. (1980), Cretaceous anoxic events: From continents to oceans, J. Geol. Soc. London, 137, 171–188. Jenkyns, H. C. (1995), Carbon isotope stratigraphy and paleoceanographic significance of the Lower Cretaceous shallow water carbonates of Resolution Guyot, mid‐Pacific Mountains, in Proc. Ocean Drill. Program Sci. Results, 143, 89–97. Jenkyns, H. C. (1999), Mesozoic anoxic events and palaeoclimate, Zentralbl. Geol. Palaeontol., Teil 1, 7–9, 943–949. Jenkyns, H. C. (2003), Evidence for rapid climate change in the Mesozoic‐Paleogene

15 of 16

PA2203

greenhouse world, Philos. Trans. R. Soc. London, Ser. A, 361, 1885–1916. Larson, R. L., and E. Erba (1999), Onset of the mid‐Cretaceous greenhouse in the Barremian‐ Aptian: Igneous events and the biological, sedimentary, and geochemical response, Paleoceanography, 14, 663–678. Laskar, J. (1999), The limits of Earth orbital calculations for geological time‐scale use, Philos. Trans. R. Soc. London, Ser. A, 357, 1735–1759. Laskar, J., P. Robutel, F. Joutel, M. Gastineau, A. C. M. Correia, and B. Levrard (2004), A long‐term numerical solution for the insolation quantities of the Earth, Astron. Astrophys., 428, 261–285. Leckie, R. M., T. J. Bralower, and R. Cashman (2002), Oceanic anoxic events and plankton evolution: Biotic response to tectonic forcing during the mid‐Cretaceous, Paleoceanography, 17(3), 1041, doi:10.1029/2001PA000623. Li, Y.‐X., T. J. Bralower, I. P. Montañez, D. A. Osleger, M. A. Arthur, D. M. Bice, T. D. Herbert, E. Erba, and I. Premoli Silva (2008), Toward an orbital chronology for the early Aptian Oceanic Anoxic Event (OAE1a, ∼120 Ma), Earth Planet. Sci. Lett., 271, 88–100. Lisiecki, L. E., and P. A. Lisiecki (2002), Application of dynamic programming to the correlation of paleoclimate records, Paleoceanography, 17(4), 1049, doi:10.1029/2001PA000733. Locklair, R. E., and B. B. Sageman (2008), Cyclostratigraphy of the Upper Cretaceous Niobrara formation, Western Interior, U.S.A.: A Coniacian–Santonian orbital timescale, Earth Planet. Sci. Lett., 269, 540–553. Lomb, N. R. (1976), Least squares frequency analysis of unequally spaced data, Astrophys. Space Sci., 39, 447–462. Lourens, L. J., R. Wehausen, and H.‐J. Brumsack (2001), Geological constraints on tidal dissipation and dynamical ellipticity of the Earth over the past three million years, Nature, 409, 1029–1033. Lourens, L., F. J. Hilgen, N. J. Shackleton, J. Laskar, and D. S. Wilson (2004), The Neogene period, in A Geologic Time Scale 2004, edited by F. M. Gradstein, J. G. Ogg, and A. G. Smith, pp. 409–440, Cambridge Univ. Press, Cambridge, U. K. Luthi, S. M. (2001), Geological Well Logs: Their Use in Reservoir Modeling, Springer, New York. Malinverno, A. (2002), Parsimonious Bayesian Markov chain Monte Carlo inversion in a nonlinear geophysical problem, Geophys. J. Int., 151, 675–688. Malinverno, A., and W. S. Leaney (2005), Monte Carlo Bayesian look‐ahead inversion of walkaway vertical seismic profiles, Geophys. Prospect., 53, 689–703. Mann, M. E., and J. M. Lees (1996), Robust estimation of background noise and signal detection in climatic time series, Clim. Change, 33, 409–445. Martinson, D. G., W. Menke, and P. Stoffa (1982), An inverse approach to signal correlation, J. Geophys. Res., 87, 4807–4818. Martinson, D. G., N. G. Pisias, J. D. Hays, J. Imbrie, T. C. Moore Jr., and N. J. Shackleton (1987), Age dating and the orbital theory of the ice ages: Development of a high‐resolution 0 to 300,000‐ year chronostratigraphy, Quat. Res., 27, 1–29. Méhay, S., C. E. Keller, S. M. Bernasconi, H. Weissert, E. Erba, C. Bottini, and P. A. Hochuli (2009), A volcanic CO2 pulse triggered the Cretaceous Oceanic Anoxic Event

PA2203


1a and a biocalcification crisis, Geology, 37, 819–822, doi:10.1130/G30100A.1. Menegatti, A. P., H. Weissert, R. S. Brown, R. V. Tyson, P. Farrimond, A. Strasser, and M. Caron (1998), High‐resolution d 13C stratigraphy through the early Aptian “Livello Selli” of the Alpine Tethys, Paleoceanography, 13, 530–545. Metropolis, N., A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller (1953), Equation of state calculations by fast computing machines, J. Chem. Phys., 21, 1087–1092. Meyers, S. R. (2008), Resolving Milankovitchian controversies: The Triassic Latemar Limestone and Eocene Green River Formation, Geology, 36, 319–322, doi:10.1130/G24423A.1. Meyers, S. R., and B. B. Sageman (2004), Detection, quantification, and significance of hiatuses in pelagic and hemipelagic strata, Earth Planet. Sci. Lett., 224, 55–72, doi:10.1016/j. epsl.2004.05.003. Meyers, S. R., and B. B. Sageman (2007), Quantification of deep‐time orbital forcing by average spectral misfit, Am. J. Sci., 307, 773–792. Meyers, S. R., B. B. Sageman, and M. Pagani (2008), Resolving Milankovitch: Consideration of signal and noise, Am. J. Sci., 308, 770–786. Mosegaard, K., and A. Tarantola (1995), Monte Carlo sampling of solutions to inverse problems, J. Geophys. Res., 100, 12,431–12,447. Muller, R. A., and G. J. MacDonald (2000), Ice Ages and Astronomical Causes: Data, Spectral Analysis and Mechanisms, Springer, London. Néron de Surgy, O., and J. Laskar (1997), On the long term evolution of the spin of the Earth, Astron. Astrophys., 318, 975–989. Ogg, J. G., G. Ogg, and F. M. Gradstein (2008), The Concise Geologic Time Scale, Cambridge Univ. Press, Cambridge, U. K. Olsen, P. E., and D. V. Kent (1999), Long‐period Milankovitch cycles from the Late Triassic and Early Jurassic of eastern North America and their implications for the calibration of the Early Mesozoic time scale and the long‐term behavior of the planets, Philos. Trans. R. Soc. London, Ser. A, 357, 1761–1786. Pälike, H., and N. J. Shackleton (2000), Constraints on astronomical parameters from the geological record for the last 25 Myr, Earth Planet. Sci. Lett., 182, 1–14.

Park, J., and T. D. Herbert (1987), Hunting for paleoclimatic periodicities in a geologic time series with an uncertain time scale, J. Geophys. Res., 92, 14,027–14,040. Premoli Silva, I., E. Erba, G. Salvini, D. Verga, and C. Locatelli (1999), Biotic changes in Cretaceous anoxic events, J. Foraminiferal Res., 29, 352–370. Press, W. H., B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling (1992), Numerical Recipes in C: The Art of Scientific Computing, 2nd ed., Cambridge Univ. Press, New York. Price, G. D. (2003), New constraints upon isotope variation during the Early Cretaceous (Barremian–Cenomanian) from the Pacific Ocean, Geol. Mag., 140, 513–522. Sambridge, M., and K. Mosegaard (2002), Monte Carlo methods in geophysical inverse problems, Rev. Geophys., 40(3), 1009, doi:10.1029/ 2000RG000089. Sambridge, M., K. Gallagher, A. Jackson, and P. Rickwood (2006), Trans‐dimensional inverse problems, model comparison and the evidence, Geophys. J. Int., 167, 528–542, doi:10.1111/j.1365-246X.2006.03155.x. Scargle, J. (1982), Studies in astronomical time series analysis, II: Statistical aspects of spectral analysis of unevenly spaced data, Astrophys. J., 263, 835–853. Schlanger, S. O., and H. C. Jenkyns (1976), Cretaceous oceanic anoxic events: Causes and consequences, Geol. Mijnbouw, 55, 179–184. Schwarzacher, W., and A. G. Fischer (1982), Limestone‐shale bedding and perturbations of the Earth’s orbit, in Cyclic and Event Stratification, edited by G. Einsele and A. Seilacher, pp. 72–95, Springer, Berlin. Sen, M. K., and P. L. Stoffa (1995), Advances in Exploration Geophysics, vol. 4, Global Optimization Methods in Geophysical Inversion, Elsevier, Amsterdam. Skelton, P. W., R. A. Spicer, S. P. Kelley, and I. Gilmour (2003), The Cretaceous World, Cambridge Univ. Press, Cambridge, U. K. Sliter, W. V. (1989), Aptian anoxia in the Pacific basin, Geology, 17, 909–912. Tarantola, A. (2006), Popper, Bayes and the inverse problem, Nat. Phys., 2, 492–494. Tarantola, A., and B. Valette (1982), Inverse problems = quest for information, J. Geophys., 50, 159–170. Tejada, M. L. G., K. Suzuki, J. Kuroda, R. Coccioni, J. J. Mahoney, N. Ohkouchi, T. Sakamoto, and

16 of 16

PA2203

Y. Tatsumi (2009), Ontong Java Plateau eruption as a trigger for the early Aptian oceanic anoxic event, Geology, 37, 855–858, doi:10.1130/G25763A.1. van Bruegel, Y., S. Schouten, H. Tsikos, E. Erba, G. D. Price, and J. S. Sinninghe Damsté (2007), Synchronous negative carbon isotope shifts in marine and terrestrial biomarkers at the onset of the early Aptian oceanic anoxic event 1a: Evidence for the release of 13 C‐ depleted carbon into the atmosphere, Paleoceanography, 22, PA1210, doi:10.1029/ 2006PA001341. Walczak, P. S., R. A. Duncan, L. J. Clarke, and E. Erba (2006), Cretaceous anoxic event 1a linked with submarine plateau volcanism: Geochemical evidence from marine sedimentary sections, Eos Trans. AGU, 87(52), Fall Meet. Suppl., Abstract PP41B‐1191. Weissert, H. (1989), C‐isotope stratigraphy, a monitor of paleoenvironmental change: A case study from the Early Cretaceous, Surv. Geophys., 10, 1–61. Weissert, H., and E. Erba (2004), Volcanism, CO2 and palaeoclimate: A Late Jurassic–Early Cretaceous carbon and oxygen isotope record, J. Geol. Soc. London, 161, 695–702. Wissler, L., H. Funk, and H. Weissert (2003), Response of Early Cretaceous carbonate platforms to changes in atmospheric carbon dioxide levels, Palaeogeogr. Palaeoclimatol. Palaeoecol., 200, 187–205. Wunsch, C. (2004), Quantitative estimate of the Milankovitch‐forced contribution to observed Quaternary climate change, Quat. Sci. Rev., 23, 1001–1012. Yu, Z. W., and Z. L. Ding (1998), An automatic orbital tuning method for paleoclimate records, Geophys. Res. Lett., 25, 4525–4528. E. Erba, Dipartimento di Scienze della Terra, Università degli Studi di Milano, Via M a ng ia ga ll i 3 4, I‐ 2 01 33 M i l a n o, It al y. ([email protected]) T. D. Herbert, Department of Geological Sciences, Brown University, Box 1846, Providence, RI 02912, USA. (timothy_herbert@ brown.edu) A. Malinverno, Lamont‐Doherty Earth Observatory, Columbia University, 61 Rt. 9W, Palisades, NY 10964, USA. (alberto@ldeo. columbia.edu)