1 Data and Models
Our approach uses forward-models of each aspect of the system that results in the paleoclimate data, and then learns all unknown quantities in a joint inferential analysis. The components we include are an astronomical forcing model (a) that drives a climate model (b); an archive model (c) that describes which quantities are stored and how they are stored; and a measurement model (d) relating observations (e) to the true values. This is summarised in Figure 1. The models used here are relatively simple, but even for these the inference is computationally challenging. However, in principle each component could be replaced by a more sophisticated choice, and, notwithstanding computational challenges, the inference methodology described in this article could still be used.
1.1 (a) Forcing
The prevailing theory is that the glacial–interglacial cycle is primarily driven by the seasonal and spatial variation of incoming solar radiation, termed “insolation”, due to variations in the Earth’s orbit around the Sun. The orbit is characterised by precession, obliquity, and eccentricity. Precession refers to the angle, , made between the point of perihelion (the point of the orbit when the Earth is closest to the Sun) and the vernal point marking the spring equinox, and as such determines when in the seasonal cycle the Earth is closest to the Sun. Obliquity is the angle between the equator and the orbital plane, and determines the insolation contrast between summer and winter. Eccentricity measures how much the Earth’s orbit deviates from a perfect circle (indicated by zero eccentricity), and hence modulates the effect of precession. It is often convenient to refer to climatic precession, , which combines the effects of eccentricity and precession in order to indicate the effect on the Northern Hemisphere summer insolation. Climatic precession can be complemented with coprecession, , to effectively compute insolation at any time of year at any latitude . In this article we compute these quantities using the algorithm in , which is suitably accurate over the past 1 Myr.
Phenomenological models are typically forced by an orbital forcing function that summarises the effect of the seasonal and spatial distribution of insolation . Here we use a forcing function of the form
where , , and , are the normalised climatic precession, coprecession, and obliquity respectively. The parameter weights the linear combination, and the model is unforced when is set to zero.
1.2 (b) Climate model
The astronomical forcing alone does not explain all of the features of the glacial-interglacial cycle, and so the internal dynamics of the climate system must also be considered [19, 20]. The approach we follow is to model the Earth’s climate as a dynamical system forced by the variation in the insolation . The complexity of the climate model is necessarily limited by the computational cost of generating model simulations. Many simple phenomenological models have been proposed, which typically comprise of a small number of differential equations representing hypothesised relationships between different aspects of the climate. Here we take the CR14 model (a modified version of ), which treats the climate model as a forced oscillator, i.e., the Earth’s climate would fluctuate between hot and cold periods in the absence of forcing, but the oscillation is paced by the astronomical forcing. In addition, we represent atmospheric variability as a stochastic process, resulting in the following SDE (suppressing dependence on time)
in which is taken to be ice volume, and is a non-physical variable acting to switch between glacial and interglacial states. The parameter scales the ice accumulation/ablation process depending on the sign, scales a linear feedback process, stabilises the system when positive, is an inverse time scale, controls the ratio of time scales between and , and and scale the stochastic fluctuations.
1.3 (c) Archive model
There are no direct measurements of the Earth’s climate or ice sheet extent over the time scale of interest. Information about the past state of the climate is stored in climate archives. For example, the quantity of in the ocean varies as a function of global temperature and sea ice extent. This information can be extracted from sediment cores, for example, in order to study the climate history. In general, larger values of are indicative of a cold climate with a large amount of sea ice. Here we assume a simple linear model between ice volume (), and (denoted ), i.e.,
The records we use in this article are taken from benthic sediment cores. Sediment accumulation is inherently stochastic, so as a starting point we model sediment accumulation by
where is the amount of sediment (in m), a standard Brownian motion, and the parameters and
represent the mean and variance of the sediment accumulation process. Under this model, sediment can accumulate and dissipate, but whenthe trend is for linear accumulation over time.
Equation 3 can be inverted to give a model for the time, , corresponding to some core depth , where identifies core slices. Notationally we take larger values of to be nearer the top of the core (), as the present, as the present sediment level, and as the top of the core. When a core is sampled at depth , the climate information recorded corresponds to the most recent time at which . This gives a first passage time problem under the time reversal of Equation 3
, the solution of which is an inverse Gaussian distribution, so that
We can obtain
using Bayes theorem by noting thatfollows Equation 4 conditioned on present values. Note that even though our sediment accumulation model is non-monotonic, age monotonically increases with core depth.
Following the sediment accumulation process the sediment is subject to post-depositional effects such as core compaction, which we can include by extending the model. In order to model compaction, we utilise the linear porosity model of , which should be suitable for the core depths of interest . This compaction adjustment introduces two additional parameters that need to be inferred, the gradient, , and the intercept, . The model transforms depth measurements into a non-compacted equivalent, using
which can then be substituted into Equation 4.
1.4 (d) Measurement model
The final component in the forward model relates the output of the climate model to measurements taken from paleoclimate archives. For
we use Gaussian white noise measurement error,
where denotes the observation at depth , , and scales the amount of measurement error. In addition, there are often features contained in cores that have been independently dated, providing valuable information for age estimation. Geomagnetic reversals, for instance, have been independently dated to high accuracy, and are frequently observable in benthic cores. However, these are rare events, with the most recent being the Brunhes-Matuyama (BM) reversal, which occurred approximately 780 kyr ago. Where the BM reversal can be observed within a core, we use 780 kyr as the age estimate of the associated
observation, and assume Gaussian error with a standard deviation of 2 kyr.
1.5 (e) Data
We use the ODP677  and ODP846  benthic sediment cores. Both cores have been dated both as part of an orbitally tuned scheme , and a non-orbitally tuned scheme , giving two different dated records for comparison. The BM reversal is identifiable in both ODP677 (at 30.4 m) and ODP846 (at 28.7 m), and so we take the measurements at these depths as the starting values. The number of observations following the BM reversal are and , for ODP677 and ODP846 respectively.
1.6 (f) Inverse problem
We employ a Bayesian approach that simultaneously estimates observation ages, model parameters, climate states, and chooses between models by estimating Bayes factors. Formally, we target the posterior distribution
where is the user defined prior distribution (provided in the supplementary information), and , , and are densities induced by the forward models described in (a)–(d). The normalizing constant, , is termed the model evidence. The ratio of model evidence terms between two models, i.e. where and are model identifiers, is the Bayes factor of model over , and indicates the relative explanatory power between the two models. Bayes factors are a commonly used tool for performing model selection in a Bayesian framework. Standard interpretations of the Bayes factor are described in . The Bayes factors provide a principled way to undertake model selection, such as comparing two different phenomenological models, or different astronomical forcings. Since the posterior distribution has no analytical solution, we use a Monte Carlo approach which characterizes the posterior distribution with a large number of random samples, each one of which consists of a set of parameter values, climate reconstructions, and observation ages. Further details are given in the Methods section.
2 Results and Discussion
2.1 Joint inferential analysis reliably infers unobserved components from synthetic data
We begin with a simulation study to demonstrate the ability of our joint inferential analysis for age estimation, state estimation, parameter estimation, and model selection. We simulate a set of observations from the forward models, and then attempt to recover the true parameters, state, and observation ages. We also estimate the Bayes factors between the forced and unforced model to ensure that we can determine the importance of the astronomical forcing. Specifically, observation times were drawn from an imagined core of length 32m sampled at 0.1 m intervals, giving observations in total (giving 980 quantities to infer). The first observation is a noisy measurement of the true age, where the noise is sampled from a Gaussian distribution with mean zero and a standard deviation of 2 kyr, as if we had observed the BM reversal.
shows the simulated observations, and compares them with the 95% highest density region (HDR) intervals for the estimated observation times and climate states, as well as the marginal posterior distributions of the parameters and their true values. The linear trend has been removed from the age versus depth plot so that the variation is more clearly visible. Note that the age-depth relationship is not a simple linear relationship (which would be a horizontal line), and in particular, there is a large period of time in which little sediment is deposited in the middle of the record. Despite this hiatus, the true observation times are in regions of high posterior probability density throughout the dataset, showing that we are able to recover the observation ages. Likewise, the majority of the true values for both the observable and unobservable state variables lie within the 95% HDR intervals throughout the core. As would be expected, the posterior variance for the unobservable state is larger than for the observable state, particularly when the system switches between glacial and interglacial periods. Finally, the true parameter values lie in regions of high posterior probability density.
We repeat the analysis for these data but now using an unforced version of the CR14 model (with the forcing parameters set to 0). We focus on the Bayes factor to determine whether we can infer the importance of the astronomical forcing. The Bayes factor is in favour of the forced model, suggesting that even with the age uncertainty, the data strongly support the forced model [15, 23]. In other words, even though we have a relatively small number of noisy age-depth measurements from a single core, we are able to correctly infer the importance of the astronomical forcing in the climate record, even after accounting for the uncertainty that arises from estimating 17 model parameters, the depth-time relationship, and the true climate signal. To ensure that this is a reasonable result, we can perform the model selection experiment on simulated data that has been generated with the unforced model. In this case, we find the Bayes factor is in favour of the unforced model, demonstrating that we are inferring the forcing parameters, and not simply assuming that the forcing plays a crucial role. This performance is consistent across different realisations of the synthetic data.
2.2 Reconstructions from ODP677 and ODP846 are consistent with LR04 but not H07
We now analyse ODP677 and ODP846, which are shown in Figure 3 with the estimated sequence of 95% HDR intervals for the observation ages. We include the age estimates from the LR04  and H07  stacks for comparison. The age uncertainties are larger than in the simulation study, likely as a result of model discrepancies. Between the two cores, the age uncertainties are typically larger in ODP846 than ODP677, with the mean standard deviation of the age estimates being 6.5 kyr in ODP846 and 3.5 kyr in ODP677. Both are smaller than the previous uncertainty estimates, which were up to 11 kyr [7, 11, 8]. Additionally, the most uncertain estimates are not necessarily at the mid-point between age control points (such as the present, or geomagnetic reversals), which has previously been assumed [11, 8]
. Rather, the age control points only seem to constrain the ages within a few metres of the core. Our age estimates for ODP677 are consistent with the LR04 age estimates, which lie in credible intervals throughout the sediment core. On the other hand, the H07 estimates deviate greatly from our estimates between 11 m and 16 m. Our age estimates for ODP846 are consistent with both LR04 and H07, primarily due to the larger variance in the age estimates. However, the LR04 estimates are notably closer to the posterior mean. For both datasets it can be seen that using a linear age-depth relationship will lead to poor estimates of observation times.
The sequence of 95% HDR intervals of the normalized ice volume vs time are shown in Figure 4. We include the LR04 and H07 stacks for comparison. It is reassuring that our ice-volume reconstructions from ODP677 and ODP846 are remarkably similar despite being obtained independently. The similarities with LR04 are again very striking, whereas H07 is out of agreement between 200 and 400 kyr ago. The likely reason is that the ages in the H07 stack are purely depth derived, and since this period is distant from the age-control points provided by the BM reversal and the present (core-top), the H07 reconstruction has low accuracy here.
Our approach has the advantage of using information on climate forcing, while preserving the possibility to test the alternative hypothesis that the astronomical forcing has no influence. Repeating the experiment using the unforced CR14 model yields Bayes factors in favour of the forced model against the unforced model of approximately in ODP677, and in ODP846. In other words, the forced model is strongly supported by ODP677, while evidence about orbital forcing from ODP846 is weaker, but not contradictory.
2.3 Joint inferential analysis, as opposed to a multi-stage analysis, is essential
In  it was demonstrated that model selection experiments are sensitive to the choice of age model by showing that two sets of age estimates consistent with the reported uncertainty lead to different conclusions, and it was thus argued that the age uncertainty must be incorporated into any analysis. In order to further illustrate this statement we propose the following experiment. In performing a joint inferential analysis we have randomly generated realizations of the observation ages of each sediment core. By sampling from these realizations we can obtain plausible age estimates that differ from each other consistently with the age uncertainty estimates. We can then repeat the analysis with the fixed age estimates in order to test whether our conclusions significantly differ when the age uncertainty is ignored. Taking four realizations of the ages for ODP677, and then obtaining the Bayes factor in favour of the forced over the unforced model, yields Bayes factors of approximately , , , and . Likewise for ODP846 we obtain Bayes factors of , , , and ; values less than one are evidence in favour of the unforced model. In ODP677 we always favour the forced model, but with varying degrees of confidence. In ODP846 our conclusions can change drastically, likely because the age uncertainty is greater. In either case it is demonstrated that ignoring the age uncertainty can significantly alter the results of such analyses. Due to strong couplings between the different components of the system, accurately characterizing and propagating uncertainties requires a joint inferential analysis.
We have investigated an approach to calibrating dynamical climate models and testing between competing hypotheses via model selection, whilst jointly fitting an age model to a sediment core. Performing a joint inferential analysis in this manner is highly challenging, requiring state of the art statistical methods and intensive computation (the analyses presented here each took six days on a standard desktop computer). Nevertheless there are notable advantages in undertaking a joint inferential analysis over splitting the analysis between multiple stages. Firstly, the joint inferential analysis both estimates climate dynamics and forcings, and uses this information to constrain the ages, without the risk of circular reasoning. Secondly, we are able to characterize a range of uncertainties, and push these uncertainties forward into our conclusions, making those conclusions more reliable.
There are several ways in which the approach presented here could potentially be extended. Firstly, in each experiment we only utilize data from a single sediment core, and a natural extension is to combine observations from multiple cores. However, this is a non-trivial extension; with multiple cores the order in time of the observations is unknown, and so developing effective proposal distributions in a Monte Carlo approach is significantly more difficult. Secondly, numerous sources of uncertainty are not accounted for here. Examples include identification error for the BM reversal, and bioturbation in the sediment core. Such sources of uncertainty could be incorporated into the analysis by extending the models. Finally, the models considered here are all relatively simple, and could be replaced with more complex models. A primary obstacle in these extensions is the computational cost involved, which is already high when using simple models. Fortunately, sequential Monte Carlo methods such as those used here are amenable to parallelization, giving the possibility of dramatically improving computation times.
Our focus for model selection was testing between a forced and an unforced model for the glacial–interglacial cycle, but the joint inferential approach is relevant to a much wider range of investigations. Dynamical models can be used to investigate, for example, changes in oscillation regimes such as the mid-Pleistocene transition, abrupt changes during glacial periods such as Dansgaard–Oeschger events, and relationships between different climate variables such as ice volume and . In each of these cases the calibration of the dynamical model will be sensitive to the inferred ages of the observations. Testing between different hypotheses therefore requires effective joint quantification of age and model uncertainties, which can be achieved by performing a joint inferential analysis. Ultimately the work presented here is a first step towards more comprehensive analyses of paleoclimates that we hope will be built upon in the future.
We employ the sequential Monte Carlo squared (SMC) algorithm , which was previously implemented in  to test between competing phenomenological models of the glacial-interglacial cycle in the absence of any age uncertainty. Here we have extended the target distribution to also estimate the observation ages. SMC is advantageous as it requires little user input in selecting tuning parameters, and so can be applied to multiple models and data sets with relative ease. The full algorithm, details, and code can be found in the supplementary information.
JC was funded by an EPSRC PhD studenthsip. RDW was supported by the EPSRC-funded Past Earth Network (Grant number EP/M008363/1). MC is supported by the Belgian National Fund of Scientific Research.
-  Blaauw M (2012) Out of tune: the dangers of aligning proxy archives. Quaternary Science Reviews 36:38–49.
-  Shackleton NJ, et. al. (1984) Oxygen isotope calibration of the onset of ice-rafting and history of glaciation in the North Atlantic region. Nature 307:620–623.
-  Emiliani C (1995) Pleistocene temperatures. Journal of Geology 63(6):538-578.
-  Shackleton NJ (1967) Oxygen isotope analyses and Pleistocene temperatures re-assessed. Nature 215:15–17.
-  Tingley M, Craigmile P, Haran M, Li B, Mannshardt E, Rajaratnam B (2012) Piecing together the past: statistical insights into paleoclimatic reconstructions. Quaternary Science Reviews 35:1–22.
-  Imbrie J, et. al. (1984) The orbital theory of Pleistocene climate: Support from a revised chronology of the marine record. Milankovitch and Climate (Part 1), eds Berger AL, Imbrie J, Hays J, Kukla G, Saltzman B (Springer, New York), pp 269–305.
-  Lisiecki L, Raymo M (2004) A Pliocene–Pleistocene stack of 57 globally distributed benthic records. Paleoceanography 20:PA1003.
-  Huybers P (2007) Glacial variability over the last two million years: an extended depth-derived age- model, continuous obliquity pacing, and the pleistocene progression. Quaternary Science Reviews 26:37–55.
-  Milankovitch M (1941) Kanon der erdbestrahlung und seine anwendung auf das eiszeitenproblem (Belgrade, Königlich Serbische Akademie).
-  Milankovitch M (1998) Canon of insolation and the ice-age problem (Beograd, Narodna Biblioteka Srbije).
-  Huybers P, Wunsch C (2004) A depth-derived Pleistocene age model: Uncertainty estimates, sedimentation variability, and nonlinear climate change. Paleoceanography 19:PA1028.
-  Crucifix M (2012) Oscillators and relaxation phenomena in Pleistocene climate theory. Philosophical Transactions of the Royal Society A 370:1140–1165.
-  Crucifix M (2013) Why could ice ages be unpredictable? Climate of the Past, 9:2253–2267.
-  Feng F, Bailer-Jones CAL (2015) Obliquity and precession as pacemakers of Pleistocene deglaciations. Quaternary Science Reviews 122:166–179.
-  Carson J, Crucifix M, Preston S, Wilkinson RD (2018) Bayesian model selection for the glacial-interglacial cycle. Journal of the Royal Statistical Society: Series C 67(1):25–54.
-  Shackleton NJ, Berger A, Peltier WR (1990) An alternative astronomical calibration of the lower Pleistocene timescale based on ODP site 677. Transactions of the Royal Society of Edinburgh: Earth Sciences 81(4):251–261.
-  Mix A, Le J, Shackleton NJ (1995) Benthic foraminiferal stable isotope stratigraphy of site 846: 0-1.8 ma. Proceedings of the Ocean Drilling Program, Scientific Results 138:839–854.
-  Berger AL (1978) Long term variations of daily insolation and Quaternary climate changes. Journal of Atmospheric Sciences 35:2362–2367.
-  Imbrie J, Imbrie JZ (1980) Modelling the climatic response to orbital variations. Science 207(4434):943–953.
-  Raymo ME (1997) The timing of major climate terminations. Paleoceanography 12(4):577–585.
-  Bahr DB, et. al. (2001) Exponential approximations to compacted sediment porosity profiles. Computers and Geosciences 27(6):691–700.
-  Singer BS, Pringle MS (1996) Age and duration of the Matuyama-Brunhes geomagnetic polarity reversal from incremental analysis of lavas. Earth and Planetary Science Letters 139(1996):47–61.
-  Kass RE, Raftery AE (1995) Bayes factors. Journal of the American Statistical Association 90(430):773–795.
-  Chopin N, Jacob PE, Papaspiliopoulos O (2013) SMC: an efficient algorithm for sequential analysis of state-space models. Journal of the Royal Statistical Society: Series B 75(3):397–426.