Deposition Models for Chronological Records

Christopher Bronk Ramsey
Research Laboratory for Archaeology, University of Oxford, Dyson Perrins
Building, South Parks Road, Oxford, OX1 3QY

Abstract
Long term sedimentary sequences provide a wealth of useful information for research
into the palaeo-environment, especially in relation to past climate change. Shorter
records provide similar information in many archaeological contexts. However if suchrecords are to be used to provide precise timing of events, and more critically the
relative timing between different records, methods are needed to provide accurate
and precise age-depth models for these sequences. Given the imprecision of individual calibrated radiocarbon determinations it is necessary to use the information
we have about the deposition process to refine our chronologies and also to provide interpolation between dated levels in any sequence. Even with layer counted
chronologies, the uncertainties are sometimes hard to quantify.
This paper outlines a range of Bayesian models for deposition which have been
implemented in the computer programme OxCal. These models can be used to
combine information from the sediments themselves with radiocarbon or other direct
dating information. Such models have the potential to integrate information between
different records and provide a coherent chronology on which to base environmental
or archaeological research.
Key words: Dating, Statistics, Radiocarbon, Calibration, Sedimentation, Lateglacial, Quaternary environments, Ice cores, Varves

1

Introduction

In order to integrate the time-scales of different palaeo-climate records it isnecessary that the chronologies for each record are determined with sufficientprecision and accuracy. This entails the incorporation of information from anumber of different sources including absolute age information, relative age
information and cross correlations between records. This requires the buildingof some sort of model whether ad hoc or formal.

Preprint submitted to Quaternary Science Reviews

13 September 2006

Absolute age information can come from a variety of sources. Sometimes wehave documentary information relating events to specific years; more often we
have to deal with ages from scientific dating methods with their inevitable
uncertainties. Relative age information likewise can be very precise as is thecase for annually laminated ice cores or it can be derived from assumptions
about deposition rates. When dealing with cross correlations between cores
we have to be careful not to introduce circularity, for example by assumingthat climatic signals are synchronous. That said, some signals (such as theprimary deposition of tephra particles) are synchronous in terms of a timescale measured in years.
If we are to achieve the required degree of precision and accuracy we needto treat the uncertainties associated with the various forms of information
correctly and specifically we need to be aware when these uncertainties arecorrelated in some way, and when they are truly independent.

2

Methodological background

In order to be able to put together information of the kind that is useful in thiscontext a mathematical framework is required. The one that is employed inthis paper is the Bayesian approach which allows us to combine information ofa probabilistic nature in a formal way and which has been applied in various
forms in palaeo-climate records recently (e.g. see Walker et al. (2003) andBlockley et al. (2004). Bayesian methods have been applied to constrained
radiocarbon dates for many years now (Buck et al. (1991), Buck et al. (1992))
and with greater frequency since the availability of software to perform suchanalyses (Bronk Ramsey (1995), Bronk Ramsey (2001), Jones and Nicholls(2002)). For details of Bayesian theory readers should look elsewhere (e.g. see
Gilks et al. (1996), D’Agostini (2003) and particularly Buck et al. (1996) andBuck and Millard (2001)).
There are some aspects of Bayesian nomenclature which it is useful to clarify inthis context. In particular, there are two terms frequently used which need to
be understood. The first of these is the prior model. This is essentially what weknow about any system before we start some set of measurements. The second
is the likelihood which describes the information that we get from the set ofmeasurements itself. Bayes theorem tells us how to combine these types of
information. In the context of a depositional model of the type described here
we will build a prior model which incorporates information on, for example, theorder of deposition and depth information. The actual dating information canthen be introduced in the form of probability distribution functions (PDFs)
which represent the likelihood that any one sample has a particular age.
2

What we are aiming to do mathematically is then to find a representative set of
possible ages for each depth point in a sedimentary sequence. In practice, as theages are highly correlated by depth, this is a task which is usually impossible
to perform analytically and so we use a method called Markov Chain Monte
Carlo (MCMC) sampling to build up a distribution of possible solutions (seeGilks et al. (1996)). Using Bayes theorem the algorithms employed essentiallysample over all possible solutions with a probability which is a product of theprior and likelihood probabilities. The resulting distributions are referred toas the posterior probability densities and take account of both the deposition
model and the actual age measurements made.
Implementation of models of this kind are too complex for most people totackle from scratch. The models discussed here are all implemented in thesoftware package OxCal (v4.0) which is freely available 1 . This program allows
the information necessary for analysis to be entered and then automates themodel specification and analysis stages. The full mathematical specification
for the model is also given for reference. The MCMC analysis in OxCal v4.0
exclusively uses the Metropolis-Hastings algorithm (see, for example Gilks
et al. (1996)) rather than the mix of Gibbs-sampling and Metropolis-Hastings
used in previous versions of the program. This paper does not go into detailsof the MCMC implementation itself, only the formulation of the deposition
models.

3

Deposition Models

In order to specify a deposition model completely we need to deal both withany direct dating information (referred to here using the Bayesian term likelihood) and the information that relates the different components of the sedimentary sequence (the prior model). In practice the exact dividing line between the prior and the likelihood can be arbitrary. For example it can bequestioned if the information giving the local reservoir offset for a region ispart of the likelihood or part of the prior. Fortunately it makes no difference
mathematically to the end result though it does affect some of the diagnostic
measures (such as the agreement index - see section 4.1.1).

3.1

Likelihoods

It is by means of the likelihoods that the absolute date information is entered.
This information can come in a variety of forms. One of the most common is1

programme available for use online at https://c14.arch.ox.ac.uk/oxcal

3radiocarbon date information. This is entered into the model in the form of acalibrated date distribution (see for example Bronk Ramsey, 2001). In OxCalsuch a distribution can easily be entered when setting up a model with the
R Date function. So for example:
R_Date("A",3450,28);
will indicate that an event A has a likelihood characterised by a radiocarbon
date of 3450 ± 28. For the calibration of the radiocarbon date a calibration
data-set is required (such as Reimer et al. (2004) or Hughen et al. (2004)) andin the case of marine samples a local ∆R offset. For applications like this it isimportant that the ∆R offset is treated as a correlated uncertainty in relation
to the depositional model on all of the relevant dates (as recommended byJones and Nicholls, 2001); this is the case in OxCal v4.0,where the ∆R offset
is treated as a model parameter with a normally distributed prior, but not inprevious versions of the software.
Radiocarbon is unusual in generating such complicated likelihood distributions. For most dating methods a uniform, or normal (Gaussian) distribution
can be used. The following are all allowed definitions of age likelihoods of thisform:
Age(N("B",3000,30));
Age(U("C",2970,3030));
Age(Top_Hat("D",3000,30));
the first being a normal distribution about 3000 and the next two are different
ways of describing uniform distributions. The form of these distributions isshown in figure 1. We will return to examples of their use later.

3.2

The deposition process

Before looking at possible mathematical formulations for prior models of thedeposition it is worth looking at the deposition process itself in some detailto see what features of real depositional systems we would like to include in
our models. There are three main aspects of any deposition process which weneed to take into account:• complexity in the underlying mechanisms• random elements to the processes• abrupt changes to different modes of deposition
The first two of these alone might lead us to give up on producing deposition models. The third is a further complication. However, if we look at real4

Fig. 1. Typical likelihood distributions for chronological information (see section
3.1.

sedimentary sequences a number of factors are apparent. Firstly the abruptchanges in modes of deposition are usually recognisable in terms of the sediments themselves. In other words the expertise of the sedimentologist can beused to specify these deposition boundaries and apply these as prior information in our model. In other cases such discontinuities are clearly apparent in
the radiocarbon data themselves. Furthermore we know that in many cases(for example peat growth or lake sedimentation), despite the complexity andrandomness of the processes, the deposition rate is approximately constant.
In the end what this means is that the randomness of the process is on a
sufficiently fine scale that the deposition appears continuous and uniform. Of
course this is only true for as long as any one mode of deposition persists.
In other cases deposition is very random and piece-wise (for example slumps,
deposition from major flood events or archaeological site formation).
In some instances the actual deposition rates are known exactly. This istheoretically the case in tree-ring sequences, varved lake sediments and ice
cores. Such situations can easily be dealt with a defined deposition rate (seeBronk Ramsey et al. (2001) and the D Sequence function of OxCal). An extension to this is where the deposition rate is known approximately (see theV Sequence function in OxCal Bronk Ramsey (1995)).
One more general model that has often been applied is one where the accumulation rate is unknown but assumed to be completely constant (see for example
Kilian et al. (2000) or Christen et al. (1995)); in this paper this will be referred to as the uniform deposition model (or implemented as U Sequence in
OxCal v4). Another general model is that which merely constrains the dated
5

events to be in the order that they are within the sediment (implemented as a
Sequence in Oxcal - Bronk Ramsey (1995)); this is applicable in cases where
we can make no direct use of the depth information at all because the random
nature of the deposition process is very extreme.
The uniform deposition models (constant accumulation rate) have been further elaborated by allowing a number of changes in deposition rate to be
automatically located within the sequence (see Aguilar et al. (2002), Blaauw
et al. (2003) and Blaauw and Christen (2005)). However this is still inherently
a uniform deposition model albeit broken up into a number of sections. Others have suggested using other specific non-linear curves (as in Christen et al.
(1995) and Walker et al. (2003)) to define the deposition process but thesestill do not take into account the random nature of the processes. Mixed effectregression has also been applied (Heegaard et al. (2005)) but this does nottake into account the complexity of the radiocarbon calibration information.
An alternative approach is suggested in this paper. If we consider the deposition process as being inherently random then it would be useful to see what
we expect as the functional form for deposition. In order to understand this we
will consider a very simple example of water collecting in a rain gauge during
a period of prolonged and fairly steady rain. Because of the finite size of therain drops the height of water in the gauge will not rise continuously but injumps as each drop enters the gauge. The time gap between drops arriving will
also vary: this process can be described in terms of a Poisson process (that is
one in which the events are essentially random). This type of age-depth model
will be referred to as a P Sequence.
Figure 2 shows how the P Sequence model would work in practice. The larger
the increments the more variable the actual deposition is likely to be, whereas
with very fine increments the deposition becomes almost constant and approximates to the constant deposition model (U Sequence). The increment
size for a particular model is given by a parameter k which gives the numberof accumulation events per unit depth.
Table 1 gives a summary of the types of deposition model discussed here andtheir OxCal functions. Appendix A gives the detailed mathematical formulation of these deposition models.

3.3

Applicability to real situations

It is necessary to consider whether or not the P Sequence methodology issuitable for real depositional scenarios. In addition we need to have some way
to work out what the correct choice of k (increments per unit length) is forany particular situation.
6z

P_Sequence-coarse
P_Sequence-fine
U_Sequence00tFig. 2. Examples of deposition scenarios under three different models. In the
U Sequence model the material builds up at a constant rate. In the P Sequencemodel there is random variation from an approximately constant rate; if the sediment is coarse (large increments, or low k parameter) the variation can be considerable; for finer sediments (small increments, or high k parameter) the accumulation
is much more constant and approximates to the U Sequence model.

OxCal key word
D Sequence
V Sequence
U Sequence

P Sequence

Sequence

Description
Age gaps between points known precisely (as in
the case of tree rings, varved sediments or ice
layersAge gaps between points known approximately
with normally distributed uncertainty
Deposition assumed to be a function of anotherparameter z (usually defined as depth)
Deposition assumed to be random giving approximate proportionality to z (usually definedas depth); more flexible than the U Sequence
but less flexible than the Sequence
Ages of specified events assumed to be in thecorrect order but no use made of depth information

Equation
A.4
A.6
A.24

A.27

A.22

Table 1
The full range of depositional models considered in this paper with their OxCal key
words and equations

7

There may be some cases where the actual theoretical framework of random
Poisson processes can really be applied. If the deposition is due to essentiallyrandom events and each deposition event is of similar size then the model canbe applied directly. The k parameter is defined by the size of the deposition
events. It is useful to consider such situations as they give us an idea of whata reasonable range of values might be for k. If we consider a coarse sedimentcontaining gravel the maximum suitable value for k is likely to be around
100 m−1 whereas a fine sediment might well have a value of up to 1000 m−1
or possibly even higher. This is because of the granularity of the sediments
themselves. We expect to have variability at least at the centimetre level in agravel with grain sizes of that order whereas a much finer sediment is likely
to have a more regular sedimentation, with noise on a finer scale.
In most circumstances, however, the deposition process is more complex than
this. In particular, although it may be random in nature, the size of deposition
events is likely to be variable and may also (as in the case of varved sediments)
have a cyclical nature to it. However, in these cases the actual functional
form of the P Sequence model may still be very appropriate. In section A.2.3
it can be seen that the probability distribution for an interpolated point isapproximately normally distributed. Since many complex processes do leadto populations of variation with normal characteristics this leads us to expectthat the functional form of the model is probably applicable in cases where thetheoretical framework is not right. The difficulty in such cases is in choosing
the correct parameter k.
There are cases where it might be possible to estimate an appropriate k value
from direct measurements. In particular where there is an exposed section
or multiple cores with stratigraphic marker layers you can use the variability
in distances between the layers (see figure 3 and equation A.20). If this is
not possible it may be necessary to resort to estimating k from the dating
information itself (see section 4.1.2, Blockley et al. (2006b) and Blockley et al.
(2006a)) but clearly this is not as satisfactory since there are dangers of circularreasoning.
As an example of another case where this model might be justified, it is also
worth considering the Poisson process in relation to deposition of regular layers(such as ice or lake laminations). In such cases the deposition is usually fairly
regular rather than being random. However, we can use a model like the
Poisson model outlined above to deal with deviations from such a regularprocess. We might well expect such deviations (such as missing or double
layers) to be random processes. Supposing we have a segment of sedimentwith m uncertain randomly disposed layers we would expect the fractionalq
uncertainty in number of uncertain layers at the midpoint to be 1/(2 (m))√(from equation A.17) which gives an absolute uncertainty of m/2. So for
example an ice core or varved segment with 100 uncertain layers but with8

Fig. 3. Method for estimating the k parameter for exposed sections or multiple
cores: three stratigraphic layers are identified and the distances ui and ui between
them measured for each section; Ui = ui /(ui + vi ) and Vi = vi /(ui + vi ) are thencalculated; from these the mean values µU , µV and the unbiased estimate of thestandard deviation σU , σV can be calculated; for each sequence the estimate of thek parameter is then given by µU µV /[σU σV (pi +qi )] from equation A.20; for example
if the three u values are measured to be 45, 50 and 48 cm and the three v values
55, 63 and 72 cm, we get estimates for k for the three sections of about 3 cm−1 or
300 m−1 .

a known start and end point would have an uncertainty of 5 years at thecentre assuming the disposition of the errors was random. Now supposing
those uncertain layers are spread over M total years. The absolute uncertainty
is still the√same but the fractional uncertainty is lower by a factor m/M and isequal to m/(2M ); if we use the P Sequence model defined above to model
this using the layer count as the parameter z we require a k = M/m.

3.4

Model definition

Appendix A covers the mathematical background to the types of depth model
discussed in this paper. Here we will look at how the models can be definedin practice using the Chronological Query Language (CQL) notation of OxCal (Bronk Ramsey, 1998). The definition of the three depositional models
described here is performed using three functions:• Sequence(name)• P Sequence(name,k[,s])9P_Sequence("",1000,50){
Boundary("Bottom"){};
R_Date("",1010,25){ z=0.65; };R_Date("",887,25) { z=0.61; };R_Date("",979,25) { z=0.57; };R_Date("",848,25) { z=0.53; };R_Date("",809,25) { z=0.49; };Boundary("Change"){ z=0.47; };R_Date("",743,25) { z=0.44; };R_Date("",595,25) { z=0.38; };R_Date("",613,25) { z=0.32; };R_Date("",485,25) { z=0.26; };R_Date("",395,25) { z=0.20; };Boundary("Top"){};};
Fig. 4. Typical model specification for a short sedimentary sequence• U Sequence(name[,s])where name is the label applied to the sequence in the model, k is the assumed
number of events per unit length (using the same units as for the depth z)
and s is the (optional) number of automatically generated depth model points
per unit length.
As an example assume that we have a sequence of 9 radiocarbon dated levels
in a depositional sequence which runs from 0.675m deep up to 0.225m deep.
The deposition is such that we expect the depositional unit granularity to bearound a mm (or 1000 m−1 ). We would like to have a depth model generated
at increments of about every 2cm (or 50 m−1 ). The sedimentology suggeststhat the deposition changes at a depth of 0.48m. The P Sequence model forthis can then be defined as in figure 4, where the radiocarbon date for thelowest dated sample is at a depth of 0.65m and has a radiocarbon date of
1010±25 BP and so on. The Sequence or U Sequence models can be definedby replacing the top line with one of:• Sequence("")• U Sequence("",50)
If the z values for the boundaries are not defined then the z is assumed tobe the same as that for the next dated item in the sequence except in thecase of the last boundary where it is assumed to be the same as the item
just preceding it. In practice this means that for a single segment sequence
the z values for the boundaries do not have to be specified unless you wish
to extrapolate beyond the dated points in the sequence and only those for10

OxCal v4beta2 Bronk Ramsey (2006); r:5 IntCal04 atmospheric curve (Reimer et al 2004)

OxCal v4beta2 Bronk Ramsey (2006); r:5 IntCal04 atmospheric curve (Reimer et al 2004)

0.2 R_Date(395,25)

0.2 R_Date(395,25)

R_Date(485,25)

R_Date(485,25)

0.3

0.3
R_Date(613,25)

R_Date(613,25)

R_Date(595,25)
Depth (m)Depth (m)R_Date(595,25)
0.4
R_Date(743,25)

0.4
R_Date(743,25)

Boundary Change
0.5

R_Date(809,25)

0.5

R_Date(848,25)

R_Date(848,25)

R_Date(979,25)
0.6

R_Date(979,25)
0.6

R_Date(887,25)
R_Date(1010,25)

800

R_Date(809,25)

R_Date(887,25)
R_Date(1010,25)

1000

1200
1400
Calendar Date (AD)

1600

800

1000

1200
1400
Calendar Date (AD)1600

Fig. 5. This shows a depth model output for the example given in figure 4; in lightgrey you can see the distributions for the single calibrated dates (the likelihoods); in
darker grey, the marginal posterior distributions are shown which take into accountthe depth model; the depth model curves are envelopes for the 95% and 68% HPD
ranges; the left hand plot shows the model based on a step change in deposition
rate exactly as defined in figure 4, and the right hand model that is based on
an exponential rise in deposition over time (as described in section A.4.3); in thisinstance both models give very similar resultschanges in deposition need to be specified.

In these models the underlying assumption is that the probability of deposition
is constant over time. We can assume that the deposition rate is linearly risingin any segment by replacing the first Boundary with a Zero Boundary or for
exponential rises with a Tau Boundary. See section A.4 for more details.

Having defined the model in this way, the OxCal program will then calculate
the mathematical functions to be applied to give the appropriate prior, as
outlined in Appendix A. The MCMC analysis is performed automatically togenerate the posterior probability distributions for each event. The highest
probability density (HPD) range is also calculated and presented in a table
and can be plotted as an age-depth model (as shown in figure 5 for the example
specified in figure 4).
114Applications

Having defined the way in which these deposition models can be specified it
is worth looking at their application in practice. In this context we will lookat two specific examples.
The first example to be considered is a the varved lake sediment from Soppensee (Hajdas et al. (1993) and Hajdas et al. (2000)). This is an interesting
example because it shows that the approaches outlined here can be appliedboth in terms of depth and in terms of varve count. The second example, isthat of ice-core layer counting which is clearly critical to the work of integration of climate records.

4.1

Lake sediment deposition

The example of Soppensee (Hajdas et al., 1993) is chosen here because it has
a number of interesting features:• there is a good sequence of radiocarbon dates• the sediment is varved through much of the sequence• there are clearly some discrepancies between the radiocarbon and varve
chronologies• the Laacher See Tephra (LST) is present in the record giving some independent check of the chronology• the chronology of the sequence has interesting implications for the timingof climatic events in the Last Glacial Interglacial Transition.
In this paper the formulation of the model for the sequence will be discussed
in detail. Another paper in this volume (Blockley et al., 2006b) covers theimplications of the derived chronology.

4.1.1

Identification of outliers

In dealing with a real example like this, the problems associated with any
sedimentary sequence become immediately apparent. The first of these is thefact that there are clearly some outliers in the sequence. There are different
approaches that can be used for identifying outliers. Some statisticians recommend the use of automatic outlier detection (Christen, 1994). This approachis statistically sound but it does have some difficulties. The method involves
the assignment of a prior probability for any sample being an outlier and theanalysis then allows samples to be down-weighted if need be. The results of
the analysis depend on the outlier probability - and this is clearly very difficult12

to determine in any real situation. The approach also automates an operationthat perhaps ought to have some expert judgement applied to it. In some sediments, for example intrusive material (younger than context) is less likely tobe present than residual material (older than context). The other approach is
to identify outliers either by eye (where they are obvious) or using statisticalmeasures but always mediated by a knowledge of the sedimentary processes.
The final model is then based on a specific subset of the results and so it
is clearer what information has been used to arrive at the final results. It is
often useful to see how robust the results are to the inclusion or exclusion ofparticular key measurements.
Another method is to look for cases where the posterior probability distribution has a very small overlap with the likelihood probability distribution.
In OxCal such an overlap is calculated in the form of the agreement index(Bronk Ramsey (1995)) which gives the relative likelihood of the posterior
from the model to that of a ‘null’ model (where there are no constraints orinterrelationships between the events). See chapter 9 of Gilks et al. (1996) fora discussion of these types of measure. The threshold for acceptable agreement index is 60%. This level has been chosen by empirical comparison to χ2
tests and is criticised by some Bayesian statisticians as being arbitrary. In theend the choice is between an arbitrary prior probability for a sample being
an outlier or an arbitrary cut-off in acceptable overlap between the posterior
and the likelihood; there is no totally satisfactory solution. In this paper theagreement index method has been employed as this is what is implementedin OxCal. The relationship between the practicalities of rejecting dates using
this threshold and the fitting of different deposition models is discussed in
Blockley et al. (2006b).
Further ideas of how outliers might be treated are discussed in section 5.
In this illustrative case we wish to look at a variety of different models. So in
deciding which samples to exclude from the analysis we selected those whichwere inconsistent with the other dates in most of the models (including thestraightforward Sequence). This meant that there were still measurements
included which were inconsistent with an entirely uniform deposition rate(using either varve count or depth).
In the end, in addition to the outliers already identified by Hajdas et al. (1993),
further data-points were rejected on the basis that they gave low agreementwith the selected models (Blockley et al., 2006b).

4.1.2

Construction of models

Because this study was drawing on already published material it was thought
most appropriate to use the different sedimentary units (I to VI) of the deposit
13defined by Hajdas et al. (1993).
Having defined the data points to be included in the model and the boundaries, all that remains is to select the mathematical model to be employed.
One possibility here would be to use the exact age differences from the varve
counting to define the relative age of the samples using the known age gaps
and either a D Sequence or a V Sequence, (see section A.1) however this presupposes that the uncertainties in the varve counting are either non-existentor very well characterised - and certainly earlier in the sequence this is clearly
not the case.
Instead five different models were looked at. The first of these was simply to
use the order of deposition to constrain the dates in a Sequence. This is themost conservative of the models and almost certainly gives much wider ranges
than are realistic as neither the depth information nor the varve counts wereused. However even this model gives results which are more precise than thecalibrated dates alone.
The next two models assume uniform deposition (U Sequence) either by depth
or in terms of varve count. The meaning of these is fairly clear: in the case of
depth we assume that within each segment the same depth of sediment is laid
down per unit time; in the case of varves we assume now, not that the varves
are necessarily annual but that the number counted will be proportional to thetime gap between any two points on the core. Such models are almost certainly
too prescriptive since they allow for no fluctuations in deposition from year
to year and they do not take into account that uncertainty in varves is likely
to be due to the occasional missing or double varve in an otherwise annual
sequence.
The final two models use the P Sequence methodology outlined above for either deposition against depth or against varve number. The problem in suchcases is the estimation of the k parameter to obtain a suitable degree of flexibility for the models. The k parameter for depth was chosen to be 300m−1 or
3cm−1 which implies a depositional unit of about 3mm (see fig 6). For varves
the k parameter was chosen to be 0.25 which implies a depositional unit of four
varves (or on average about 2mm). These were the highest k which gave a satisfactory agreement with the actual dating information (using the agreement
indices). This is not a very satisfactory method of choosing this parameter
since it relies on the dating information being able to discriminate between
different models and there is a danger of circular arguments. For example twodates will always be consistent with a completely uniform deposition - eventhough this might be far from the truth.
It is worth looking at the implied uncertainty in these prior P Sequence models to see what they mean in practice. For the depth sequence model, if we had
14

P_Sequence("",3){
Boundary("VI/VII"){};
R_Date("ETH-6809",12150,90){ z=633.5; };R_Date("ETH-6808",11930,90){ z=631.5; };R_Date("ETH-6807",12040,90){ z=629.5; };R_Date("ETH-6806",11385,90){ z=611; };R_Date("ETH-6805",11300,85){ z=610; };R_Date("ETH-5305", 11380, 105){ z=607; };R_Date("ETH-6933",11470,70){ z=606; };R_Date("ETH-6932",11160,60){ z=600.5; };Boundary("V/VI"){ z=597; };Event("LST"){ z=595; };R_Date("ETH-5290",10760,105){ z=594; };R_Date("ETH-7703",10440,100){ z=577; };Boundary("IV/V"){ z=571.5; };R_Date("ETH-6929",10400,70){ z=569; };R_Date("ETH-6803",9965,75){ z=550.5; };R_Date("ETH-7710",10135,100){ z=547; };R_Date("ETH-7701",9970,100){ z=542.5; };Boundary("III/IV"){ z=539.5; };Boundary("II/I"){ z=524.5; };R_Date("ETH-6623",9595,70){ z=524.5; };R_Date("ETH-6622",9625,65){ z=523.5; };R_Date("ETH-7700",9530,95){ z=519.5; };...
R_Date("ETH-7211",6325,50){ z=398.5; };R_Date("ETH-7353",6180,55){ z=391.5; };R_Date("ETH-7210/7352",6190,40){ z=390.5; };Boundary(){};};
Fig. 6. Main elements of the model specification for the analysis of Soppensee forthe case of a P Sequence model applied against depth; there is a long sequenceof radiocarbon dates between depths of 519.5 and 398.5 cm deep not shown here(marked with ...).

a 1m segment with the age known at both ends the√ fractional age uncertainty(from equation A.17) at the centre would be 1/(2 300) which is 2.9% or about
60 years (assuming an average deposition per year of about 0.5mm). For thevarve depth sequence, for 2000 varves (about√ 1m) the relative age uncertainty
at the centre of a segment would be 1/(2 0.25 × 2000) which is about 2.2%
or about 45 years. In practice a better approach for an ab initio study might
be to estimate these uncertainties directly by making a detailed analysis ofpart of the core or multiple cores as shown in figure 3. It is also hoped that,15

with experience of using such models, better methods for estimating k froma knowledge of the sediments will be developed. It is also interesting to notethat the value found empirically for this depth sequence is right in the middleof the expected range of 100-1000m−1 (see section 3.3).

4.1.3

Analysis results

A typical depth model output is shown in figure 7. Here you can see how
the analysis allows the depth model to be slightly non-linear even in the top
section (I) and allows changes in deposition rates at boundaries. The resultsfrom the analysis of the Soppensee core are discussed in much more detail in
Blockley et al. (2006b). However there are several useful lessons that can beextracted.
• the U Sequence models diverge slightly and show poor agreement with theactual dating information unless more outliers are rejected• the results of both P Sequence models are very close and allow much moreof the dates to be included in the model (as restricted by the agreement
index).
• the precision of the P Sequence models is close to that of the varve counting
precision and not much lower than the (artificially high) precision of the
U Sequence models• the Sequence model is much less precise but overlaps all of the others• as far as we can tell the P Sequence models are accurate in their estimationof the date of the Laacher See Tephra

4.2

Ice core layers

Finally we come to consider, in a more speculative way, the application of
the methodology outlined here to other sorts of records, which have little to
do with radiocarbon calibration. In particular we will look at the way the
P Sequence model might be applied to the estimation of interpolation andextrapolation of uncertainties in ice-core sequences.
In principle, or course, the annual layers laid down in ice cores allow the direct
counting of years back from some known date. However, in practice there are
errors in any such counting process. This is at least in part due to layers thatare difficult to distinguish and so have to be identified as uncertain layers
in the sequence. There are also other possible problems such as missing coresegments, ablation events or other re-workings of the primary deposition of
the ice but these are usually discussed and accounted for in initial countingof ice core layers. Here we will concentrate on the random errors associated16OxCal v4beta Bronk Ramsey (2006); r:5 IntCal04 atmospheric curve (Reimer et al 2004)

R_Date ETH-7210/7352
ETH-7353
ETH-7211
400 R_Date ETH-7387
R_Date ETH-7388
ETH-7212/7354
R_Date ETH-7213
R_Date ETH-7214
R_Date ETH-7216
ETH-7217
R_DateETH-7218
R_DateETH-6238
ETH-7389
ETH-7586
R_DateETH-7588
ETH-6236
R_DateETH-7589
R_Date ETH-7390
ETH-7590
ETH-7622
R_DateETH-7391
R_DateETH-7591
R_DateETH-7593
R_Date ETH-7592
ETH-7392
R_Date ETH-5296
ETH-6231
ETH-7393
R_Date450 R_Date ETH-6230
R_Date ETH-7623
ETH-6144/6152
R_DateETH-7395
ETH-6143
R_DateETH-7355
R_Date ETH-5291
ETH-6142

Depth (cm)

R_Date ETH-9641
Boundary I/II
R_Date ETH-6615
ETH-6614
R_Date500 R_Date ETH-6616
R_DateR_Date ETH-6936
ETH-6617
R_DateETH-6619
R_DateR_Date ETH-6620
ETH-7699
R_Date ETH-7700
R_DateETH-6622
Boundary
ETH-6623
II/I
Boundary III/IV
R_Date ETH-7701
R_Date ETH-7710
550 R_Date ETH-6803

R_DateETH-6929
Boundary
IV/V
R_Date ETH-7703
Laacher SeeTephra

R_Date ETH-5290
Boundary V/VI
600 R_Date ETH-6932
R_Date ETH-5305
R_Date ETH-6805
ETH-6806

R_DateETH-6807
R_DateR_Date ETH-6808
ETH-6809
650
16000

14000

12000

10000

8000

6000

Calendar Date (BP)Fig. 7. Age-depth model for Soppensee using depth as the variable and assuming
the deposition is a Poisson process (P Sequence); the input for this model is shownin figure 6

with uncertain layers which are harder to test for using physical analyses of
ice cores.

4.2.1

Binomial treatment

Appendix B gives a simple Binomial treatment of ice-core uncertainties. In
this we assume that each uncertain layer has a probability of 0.5 of being
real and 0.5 of being illusory. If, for example, there are a total of N = 1600
17

uncertain layers (as determined by the analysis of the original ice core layercount) the most likely number √of true layers amongst these is N/2 = 800 andthe 95% confidence range is ± N = ±40. From this it can be seen that fora core with a total of N = 1600 uncertain layers and 39200 definite layers
the 95% range would be 40000 ± 40 which is compared to the 100% range of40000 ± 800.
However, this simplified treatment has some problems associated with it. Most
critically, we do not know that the probability of an uncertain layer being realis exactly 0.5, it might in fact be 0.9 or 0.1, and clearly our lack of knowledgeof this probability far outweighs the uncertainty estimated from this type ofanalysis.

4.2.2

P Sequence treatment

Instead of the Binomial treatment we wish to consider the problem using a
P Sequence model against layer count. We will again assume, for the sake ofargument, that within a 40000 year sequence there are 1600 uncertain layers.
We can (other possible problems aside) then use these uncertainties to define
a maximum range of possible dates for each depth within the sequence. So
at 40000 layer counted ice core years, we know that the true age lies between
39200 and 40800 - and we will assume that any of these ages is equally likely.
Such a distribution can be defined in OxCal using the T op Hat function as
Age(Top Hat(40000,800)). We will then also use the result we derived in
section 3.3 that for m uncertain layers in a total of M , a suitable k for a
P Sequence is equal to M/m which in this case gives k = 25. We can thenbuild a model for the sequence as a whole.
Figure 8 shows two possible such models. The first model (I) assumes that weknow no more than this about the date of layers within the core. The top of thecore is fixed within a year. If we then look at the output from the deepest point
which is shown as A1 in figure 9 we can see that the posterior distribution isvery little altered by the model. This makes sense as although the distribution
is forced to be fairly uniform with respect to layer count the uncertain layers
can, in principle all be real or all false. However, model II shows that if, forexample, we have a well dated event for the 10000’th layer which is 10000 ±10 years old (with normally distributed errors), the model becomes muchtighter giving us a 95% HPD range of some 320 years at the earliest pointin the sequence. This is much more realistic than the 80 year range from the
Binomial treatment. Whether it is accurate or not in a particular instance ishard to ascertain. This is because there may by systematic effects that havenot been taken into account. We assume here that the error rate is independentof depth which is probably not reasonable. It may therefore, in a real case,be necessary to split the sequence into segments with similar characteristics.
18

P_Sequence("Model I",25){
Boundary(){};
A1=Age(Top_Hat(40000,800)){ z=40000; };B1=Age(Top_Hat(30000,600)){ z=30000; };C1=Age(Top_Hat(20000,400)){ z=20000; };D1=Age(Top_Hat(10000,200)){ z=10000; };E1=Age(N(0,1)){ z=0; };Boundary(){};};
P_Sequence("Model II",25){
Boundary(){};
A2=Age(Top_Hat(40000,800)){ z=40000; };B2=Age(Top_Hat(30000,600)){ z=30000; };C2=Age(Top_Hat(20000,400)){ z=20000; };D2=Age(N(10000,10)){ z=10000; };E2=Age(N(0,1)){ z=0; };Boundary(){};};
Fig. 8. Two model definitions for ice core sequences; in model I only the layer counting uncertainties are included; in model II a well dated event has been introduced
at 10000 years ago.

However, in cases where we have multiple cores and multiple tie points, all of
this information can, at least in principle, also be built into the model andwe have the potential to construct a truly integrated chronology although,in the case of Greenland, this would involve some assumptions about timetransgressions.

5

Further developments in modelling deposition processes

It is certainly possible to further refine these models and a few possibilities
will be discussed here. However, it should always be remembered that further
complexity comes at a cost - both in terms of analysis time and, more importantly, in terms of clarity. No model will fully encapsulate all of the details of
the real processes and it may be better to keep the model sufficiently simple
that it is relatively easy for the informed user to intuitively understand theassumptions underlying it.
The first area where we might wish to add extra aspects to the model is to deal
with dated samples that are not exactly the same age as the deposits in which
19

A1

A2

41000

40500
40000
Calendar Date (BP)39500

Fig. 9. Output from the models defined in figure 8; A1 gives the distribution for thestart of the sequence assuming only the ice layer information; in light grey is the
traditional maximum error estimate and in dark grey the results of the Bayesiananalysis showing very little change; A2 shows the result of including an age constraint at 10000 years ago on the estimate for the date of the start of the core; inboth cases the 95% HPD range is shown.

they are found. There are various different ways in which such samples canbe treated. Above (section 4.1.1) we looked at the treatment of some results
as outliers, either identified by the user, or by the use of automatic outlier
detection.
One possible elaboration of the models described here would be to add theautomatic outlier detection methods developed by Christen (1994). Howevera different approach has been taken in OxCal v4. This allows any date to beoffset by some unknown amount - we can then model the offsets for a wholesedimentary sequence in a number of different ways:• The most straight-forward is simply to assume that they are part of a normal
distribution of unknown variance; the analysis will then find that variancewhile constraining the mean to be zero (that is the dates are assumed to
be on average correct). Although this is simple it perhaps shares some ofthe same problems of the automatic outlier detection in that if the reasonsfor offsets are not understood, we cannot be sure that they average to zero.
However it might be applicable in cases of bi-directional bioturbation. Thisapproach was tested in Blockley et al. (2006a).
• The other two distributions of offsets that can be modelled are suited to thesituation where we know that the samples are suspected of being residual(older than their context) or intrusive (younger than their context). In thesecases we can assume that they are distributed as an exponential tail either tothe older or younger side of the sequence. Again the analysis can determine
the time constant of this exponential tail.
These methods of dealing with offsets would seem to offer a reasonable description of what is going in the depositional processes.
20Other information might also be brought into a model using additional parameters. For example, in the example outlined in Figure 4 we might haveinformation about the relative duration of both segments of the deposition.
In OxCal we can apply a prior probability distribution to the ratio of thesedurations. This is done by adding the following two lines to the model definition:
Ratio=(Top-Change)/(Change-Bottom);Ratio&=N(1.0,0.1);
The first of these lines defines another parameter ‘Ratio’ which is the ratio ofthe duration of the upper segment to that of the lower segment; the secondline then assigns a prior probability to this ratio which is a normal distribution
with a mean of 1.0 and a standard deviation of 0.1. In principle any functional
form can be given to this prior. The ability to add such information is potentially useful, however in practice it is rare that such information is know inquantitative form.
Two other developments are possible, although they have not yet been implemented. The first is to allow for the actual depth of a sample in a sediment
to be uncertain - effectively an alternative to dealing with any offsets in thetime dimension. The second, which is more difficult, is to make the k parameter of the P Sequence model a variable - so that the analysis can find themost appropriate value automatically. This is attractive given the difficulty inassessing a suitable value. That said, there are some disadvantages in such an
approach:• a suitable prior distribution for k would have to be defined anyway• the value of k would become strongly dependent on the treatment of offsets
or outliers since if the offsets are large enough any distribution could be
fitted to a straight line
• the resultant model output would be dependent on a series of assumptions
that it is hard to comprehend easily - and perhaps equally hard to justify.
So although further mathematical elaborations of the models are theoreticallypossible, these should not be expected to solve all of the problems in dealing
with depositional models. Ultimately if you have a number of offset dates
whose cause and distribution is hard to ascertain and you also do not know
how regular your deposition is, you do not have enough information to produce
a good age depth model and no statistical sleight of hand is going to help with
this.
As discussed in Blockley et al. (2006b) modelling needs to be combined withother stratigraphical information and laboratory analyses. However, when thisis achieved, robust high precision age models can be produced.
21

6

Conclusions

In this paper a number of different depth models have been discussed. These
come into two broad categories: the specific models for use in instances where
the deposition rate is known either precisely (D Sequence) or approximately(V Sequence) and the more general models where the actual deposition rateis unknown.
The generalised models cover a whole spectrum from those where the deposition is assumed to be exactly uniform (U Sequence) through to cases where all
we can say is that the deposition took place in a particular order (Sequence).
The P Sequence fills this spectrum with a k parameter defining the stiffnessof the model; if k = 0 this gives a model identical to a Sequence and if k = ∞(infinity) the model is the same as a U Sequence. From a theoretical perspective the P Sequence would seem to be the most realistic in all real cases sinceno deposition is truly uniform and neither is it such that the depth can tell us
nothing (as in the case of the Sequence).
In support of this, all of the deposition models discussed here have been assessed on synthetic sedimentary sequences and the P Sequence found to be
the most effective, in general terms, at reproducing the true sedimentation
rate in all cases considered (Blockley et al., 2006a).
The implementation of these models in OxCal, allows us to produce age-depth
models for a whole range of chronological records. It also, equally importantly,
allows us to study the sedimentation process to better understand issues suchas residuality, intrusion and changes in sedimentation rate.

7

Acknowledgements

Thanks are due to English Heritage for funding these developments in OxCal
though a research grant (Historic Environment Enabling Programme Grant
3164) and the Natural Environment Research Council for institutional support
through Services and Facilities. I am also very grateful to all of those who havecontributed to the ideas behind this paper through discussions and previous
publications; in particular I would like to thank Alex Bayliss, Maartin Blaauw,
Simon Blockley, Caitlin Buck, Zhiyu Guo, Andrew Millard, Geoff Nicholls andPeter Steier all of whom have been of particular help.

22

APPENDICESAMathematical formulation of deposition models

This section will deal with the mathematical formalism that underlies thedeposition models defined in this paper. The aim is to model the deposition
events that led the record as it is currently found. The events are labelledfrom the start (bottom) of the record i = a to the end (top) i = b. They are
assumed to have occurred at times ti (which are calendar dates and not to be
confused with the radiocarbon dates) and to correspond to specific depths zi
in the record. In some models the actual depths are not used but in all cases
we assume that the events occurred in a specific order - that is that:ti < ti+1a≤i<bfor all i such that

(A.1)

This can be expressed in terms of an unnormalised prior which is defined as:
pH (t) = pH (ta , ta+1 , · · · , tb ) =

1 ta < ta+1 < · · · < tb
0 otherwise(A.2)

We will also assume that z changes monotonically with time. For clarity weassume that z increases with increasing t:
zi < zi+1a≤i<bfor all i such that

(A.3)

although z could equally well be defined to be depth below the surface without
affecting the models (the OxCal implementation will accept z as height or
depth). The overall model can usefully be split into three main components:• the likelihood distributions p(yi |ti ) for individual parameters (see section
3.1) where, for example, direct dating information yi is included for a specificQdate ti . The product of all of these will be denoted by p(y|t) = i p(yi |ti ).
• the constraints on the parameters (in the case of sedimentary sequences
these are given in equation A.1 or, as a prior, in equation A.2)
• other informed prior probability functions relating the specific parameters(which are described in this section)
In all cases we will assume that the uninformed prior for any parameter p(ti )is constant. Their product p(t) will therefore also be constant.
23

A.1

Known age gaps

Before discussing generalised deposition models it is worth considering a coupleof special cases. The first of these is the case where the exact time interval
is known between two different points in the deposition sequence. This istheoretically the situation with a varved sediment or a tree ring sequence. In
such cases we have information which defines the time gap precisely (referredto here as a D Sequence. We have:ti+1 = ti + gi

for all i such that

a≤i<b

(A.4)

where gi is the gap between one point and the next. In practice such a model
has only one independent variable (which can be chosen to be the start of thesequence, ta ) since all of the relative ages are exactly defined. This type ofmodel is already widely used (see Bronk Ramsey et al. (2001) and Christen(2003)). If we take the prior for ta to be uniform the posterior is proportional
to the product of the likelihoods with t being a a function only of ta and g:p(t|y, g) ∝ p(y|t(ta , g))(A.5)

Some allowance for variation can be introduced where this gap is only approximately known with a mean gap of gi and an uncertainty in that figure of si .
Then we have:ti+1 = ti + gi ± si!
(ti+1 − ti − gi )2
1
p(ti , ti+1 |gi , si ) ∝ √ exp −
2s2i
si 2π
and the informed prior probability contains a factor proportional to:p(t|g, s) ∝b−1
Yi=a si1√(ti+1 − ti − gi )2
exp −
2s2i
2π

!

(A.6)

Using Bayes theorem, the posterior probability then becomes:p(t|y, g, s) ∝ p(y|t)p(t|g, s)pH (t)(A.7)

Such models have been implemented in OxCal (Bronk Ramsey, 1995) as a
V Sequence. These models do not use depth information at all explicitly, except in as much as they are subject to the constraint given in equation A.1 andcan only be used for age-depth models in very particular circumstances. Other24

functional forms might be more appropriate and this can also be catered for.
The usual way of doing this is to define parameters for the intervals between
events di = ti+1 − ti and then assign a likelihood function to the parameterdi . This additional information can be used to further constrain the Sequencemodel described in the following section.

A.2

Interpolation between points

In order to make use of depth information we need some method of interpolation between points of known age. Let us suppose that we know the (calendarage) time ti of depth zi and the time ti+1 of depth zi+1 . What we would like toknow is what the probability distribution function is for the time t at depth
z such that zi < z < zi+1 (and thus ti < t < ti+1 ).
A.2.1

Simple ordered deposition (Sequence)

Consider the case of the model in which all we assume is that the deposition
is monotonic (which in a sedimentary sequence simply means that there areno erosional episodes). This model will be referred to here as a Sequence. In
this case for any intermediate depth all we know is that the time of deposition
must lie somewhere between ti and ti+1 . In this case there is a uniform priorfor t over the range ti to ti+1 :p(t|z) =A.2.2

1
ti+1 −ti0ti < t < ti+1
otherwise(A.8)

Uniform deposition (U Sequence)

Next let us consider the model of uniform deposition that is often applied (seeexample 2 in Christen et al. (1995) or see Kilian et al. (2000)). In this model(U Sequence) the time at any particular depth is exactly defined throughlinear interpolation over the range zi ≤ z ≤ zi+1 :t = ti + (ti+1 − ti )A.2.3

(z − zi )(zi+1 − zi )(A.9)

Poisson model deposition (P Sequence)

Finally we turn to the Poisson mediated deposition model (P Sequence). In
this case to derive a probability distribution we need to know the size ofthe incremental jumps in depth. In the case of the rain analogy this is theheight gain in the rain gauge due to one rain drop, with every drop having
25

the same size. We will introduce this as a parameter k which is the numberof increments per unit length - so if it takes 1000 rain drops to raise the levelin the gauge by 0.01 m then k = 105 m−1 . The total number of increments
occurring in the range zi to zi+1 is then simply N = k(zi+1 − zi ) and thenumber arriving in the range zi to z is n = k(z − zi ). Strictly speaking N andn should be integers since z is only able to change by increments of 1/k. The
arrival times of the N increments are assumed to form a Poisson process, withexponentially distributed intervals between arrivals, and with an unknown but
constant arrival rate.
The arrival time of each increment has a uniform probability distributionfunction:f (t0 ) =1
(ti+1 −ti )0ti < t0 < ti+1
otherwise(A.10)and therefore the probability that any such increment takes place between ti
and t is:
q=

Zt
ti

0

0

f (t )dt =t − titi+1 − ti!

(A.11)

Thus the probability that n of them precede t and N − n succeed t is givenby:N!
p(n|t, N ) ∝n!(N − n)!

t − titi+1 − ti!nti+1 − tti+1 − ti!N −n= N Cn q n (1 − q)N −n!
1(n − N q)2≈q
exp −2N q(1 − q)
2πN q(1 − q)

(A.12)
for large N(A.13)

which is a likelihood function for n. We now use Bayes theorem, observing
that N Cn is independent of t and assuming that the uninformed prior p(t) is
uniform. From equation A.12 we find an informed prior for t:p(t|n, N ) ∝ p(n|t, N )p(t)∝ q n (1 − q)N −n(t − ti )k(z−zi ) (ti+1 − t)k(zi+1 −z)p(t|z, k) ∝(ti+1 − ti )k(zi+1 −zi )26(A.14)

By taking the derivative of equation A.14 (with regard to t) and equating tozero we see that this function reaches a maximum when:zi+1 − zz − zi=t − titi+1 − t(A.15)

which makes sense as this is equal to the interpolated point from the U Sequence(since equation A.15 can be rearranged to give equation A.9). If we now assume that N is large and that we are close to this maximum the distributionp(t|z, k) approximates to a normal distribution with respect to t with a mean:
ti+1 − tiµt = ti + (z − zi )zi+1 − zi!

and a standard deviation of:σt = (ti+1vuu (z − zi )(zi+1 − z)q− t )ti(zi+1 − zi )2

1
k(zi+1 − zi )(A.16)

This is useful as it allows us to see the effect of different k parameters on therigidity of the depth model. In particular we can see that for a midpoint thefractional uncertainty σt /(ti+1 − ti ) is:
1
1 1
σt
1
= q
= √
(ti+1 − ti )2 k(zi+1 − zi )2 N

(A.17)

More generally we can estimate k where we have some measure of the variability of the sedimentation rate in different sections. We are most likely toknow how variable z is for some specific t. We can use equation A.13, and theassumption that we are not too far from the maximum probability, to estimate
the expected standard deviation in z:σz =σn 1 q=N q(1 − q)
k
k v
√ u
Nu
t (t − ti )(ti+1 − t)=k(ti+1 − ti )2√ v
uNut (z − zi )(zi+1 − z)≈k(zi+1 − zi )2vuu (z − zi )(zi+1 − z) 1
σz√≈t
(zi+1 − zi )(zi+1 − zi )2N

27(A.18)

We can make measurements that allow us to estimate this if we have multiple
cores sharing chrono-stratigraphic markers or an exposure with identifiable
layers. Let us define:u = zi+1 − z
v = z − zi
U = u/(u + v)
V = v/(u + v)

(A.19)

By measuring u and v in several places (see figure 3 in the main text) we canfind mean and standard deviations for U and V (the means being defined asµU , µV and the standard deviations as σU , σV ). We then have from equation
A.18:
σz(zi+1 − zi )√1
= UV √N
1
σU σV ≈ µU µV
N
µU µV
N≈
σU σV

σU = σV ≈And then for each sedimentary sequence we arrive at an estimate for the kparameter, kest :kest =µU µV
1
σU σV (u + v)

(A.20)

To summarise the nature of interpolation using this P Sequence model: theprobability peaks at the same point as that derived from the U Sequencemodel and, like the Sequence model, falls to zero as t approaches either ti or
ti+1 . As k tends to zero the distribution becomes identical to the Sequencemodel and as it tends to infinity it effectively becomes the same as theU Sequence model. At intermediate values where k is high the distribution
becomes very similar to a normal (Gaussian) distribution.

A.3

Overall prior probabilities

Let us now set up the entire prior probability associated with the three generalised models under consideration here (Sequence, U Sequence and P Sequence)28

and see how they might be used to calculate posterior probability densities.
In all cases we will consider the sedimentary sequence in sections which aredefined by boundaries. These sections can essentially be treated independently
although there may be some constraints, such as the start of one section beingof the same age as the end of the previous one.
For each segment of interest we have the event a (at time ta ) which is the startof the segment, an event b (at time tb ) which marks the end of the segment, andall of the intermediate events of interest a + 1 · · · b − 1 (at times ta+1 · · · tb−1 ).
Within the segment we assume that the deposition process remains the sameand that we have the constraints ti < ti+1 (and thus ta < tb ).

A.3.1

Simple ordered deposition (Sequence)In the case of the simple Sequence, the prior probability density for any single
event ta+1 · · · tb−1 occurring at time t0 is given by:
0fu (t ) =1
(tb −ta )0

ta < t0 < tbotherwise(A.21)

And choosing the prior for ta and tb to be uniform, the overall informed priorprobability, subject to the constraints given in equation A.1, is proportionalto:pu (t) ∝

b−1
Y

1
1
=
(tb − ta )b−a−2
i=a+1 tb − ta

(A.22)

and as the total number of possible combinations of ti is proportional to (tb −ta )b−a−2 this prior is neutral with respect to the length of the segment. The
posterior is then calculated as:pu (t|y) ∝ p(y|t)pu (t)pH (t)(A.23)In OxCal such a model would be defined using the commands:
Sequence(){
Boundary("a");...
Boundary("b");};29

A.3.2

Uniform deposition (U Sequence)In the case of the uniform sequence, or U Sequence, the deposition rate is
assumed to be exactly defined by:dz∝ fu (t) = constant within the range ta < t < tb
dt
And so integrating with respect to time, and choosing ta and tb to be our
independent variables with uniform priors, intermediate times in the range
ta < t < tb are precisely determined by the times of the end segments:t = ta + (tb − ta )z − zazb − za(A.24)

This gives a posterior based only on the likelihoods and the constraints fromequations A.24 and A.1:
pu (t|y, z) ∝ p(y|t(ta , tb , z))pH (ta , tb )(A.25)In OxCal such a model would be defined using the commands:
U_Sequence(){
Boundary("a"){ z= za; };...
Boundary("b"){ z=zb; };};where za and zb are suitably defined.
In this case the prior is again neutral with respect to the difference between
ta and tb .

A.3.3

Poisson model deposition (P Sequence)

Finally we come to the case of the Poisson mediated P Sequence. Here thetotal number of events is assumed to be k(zb − za ) − 1 plus the start and end
events (a and b). This is made up of both modelled events ti with defined depth
zi and N other accumulation events (which do not have model parameters
associated with them) that are assumed to result in the accumulation of thesediment. The difference in z between each accumulation event is 1/k and so
between two modelled events occurring at times ti and ti+1 there are ni =Pk(zi+1 − zi ) − 1 accumulation events. Strictly speaking N = bi=a ni and ni
30should be integer, however in the OxCal implementation they are taken to bereal numbers with a minimum of 0. The modelled events are a subset of theaccumulation events. The accumulation events can occur at any time between
ta and tb and thus the probability distribution function, as for the Sequencemodel, is given by equation A.21. The chance that any one will occur between
time ti and time ti+1 is given by the integral of the function fu (t0 ) from equation
A.21:
tZi+1
ti

ti+1 − tifu (t )dt =tb − ta00(A.26)and so the probability that exactly ni = k(zi+1 − zi ) − 1 of them occur in thesame interval is given by:
pu (zi , zi+1 |ti , ti+1 , ta , tb , k) ∝ N Cni

ti+1 − titb − tani 

ti+1 − ti1−tb − taN −ni

Furthermore the probability of exactly ni events occurring in each of the intervals (ti to ti+1 ) is:pu (z|t, k) ∝ N !

b−1
Y1 ti+1 − titb − tai=a ni !
k(zi+1 −zi)−1

and we use Bayes theorem to give:pu (t|z, k) ∝ pu (z|t, k)pu (t)where the informed prior pu (t) is from equation A.22 (still valid for the modelled events), we find that the overall informed prior:
Qb−1

pu (t|z, k) ∝a(ti+1 − ti )k(zi+1 −zi)−1(tb − ta )k(zb −za )−1(A.27)

Thus for a single segment of this form the overall posterior would be:pu (t|y, z, k) ∝ p(y|t)pu (t|z, k)pH (t)In OxCal such a model would be defined using the commands:
P_Sequence(k){Boundary("a"){ z= za; };31(A.28)

...
Boundary("b"){ z=zb; };};where k, za and zb are suitably defined.
Once again this prior is neutral to the overall length of the segment tb − ta
because of the number of possible combinations of k(zb − za ) + 1 events (which
includes the two boundaries) is proportional to (tb − ta )k(zb −za )−1 .
In the case of the U Sequence and the P Sequence the depth dependency ofthe model is independent of the number of modelled events (assuming this isless than the total number of events k(zb −za )−1 postulated in the sequence; if
it is greater then the formulation of the model will break down). However, thisis not the case in the Sequence model where the more events are introduced
into the model, the more rigidly uniform the deposition will be. Another wayof looking at the P Sequence is that it is the same as a Sequence model where
one event is defined for each unit of deposition (ie each unit of depth 1/k). Andagain the U Sequence corresponds to the limit of an infinite event density.

A.3.4

Multiple segments

So far we have only considered one segment of a deposition sequence. In manycases it may make sense to break a sedimentary sequence into several segments.
In such cases each segment is treated independently although frequently theend of one segment will be coincident with the start of the next. Where theboundaries between segments are made is down to the judgement of the sedimentologist since this needs to be based on an assessment of the sediments
themselves. What they should look for is evidence that the mode of sedimentation has changed. Because this is a matter of judgement it may be necessary
to try more than one model and see how sensitive the end results are to theseboundaries.
Mathematically one property that we would like the models to have is that ourprior for the overall length of the depositional sequence should be independentof the number of intermediate boundaries. If the first boundary is at time ta
and the final boundary at tm , each of the internal boundaries has a uniform
prior which is defined as:p(t) =

1
tm −ta0ta < t < tmotherwiseand for n internal boundaries between the two outer boundaries ta and tm ,32

the overall prior is proportional to:p(ta , tm |n) =1(tm − ta )n(A.29)

which, since the total number of possible combinations of n + 2 boundary
dates is proportional to (tm − ta )n gives a prior for the overall span whichis independent of the number of internal boundary divisions defined for themodel. In OxCal this factor is applied as long as the ‘Uniform span prior’
option is set.
The final factor in the prior comes from the fact that in practice the overall
sedimentary sequence is sometimes constrained within outer limits. In the casethat there is a lower limit llima for ta and an upper limit ulimm for tm thenumber of possible combinations of solutions is proportional to:g(ta , tm , llima , ullimm ) = (ulimm − llima ) − (tm − ta )(A.30)and so the reciprocal of this can be added as a prior factor (these latter two factors were suggested in Nicholls and Jones (2001) and Jones and Nicholls (2002)
and are incorporated in OxCal prior to version 4 as described in Bronk Ramsey (2001)). If there is also an upper limit ulima for ta and a lower limit llimm
for tm the number of possible combinations becomes:g(ta , tm , llima , llimm , ullima , ullimm )= min( (ulimm − llima ) − (tm − ta ),(tm − ta ) − (llimm − ulima ),ulima − llima ,ulimm − llimm )(A.31)and so the reciprocal of this quantity:p(ta , tm |llima , llimm , ullima , ullimm )=1g(ta , tm , llima , llimm , ullima , ullimm )(A.32)

that is added as a prior factor in OxCal v4. When ulima = ulimm and llima =llimm the factors from equations A.30 and A.31 are identical.
33

A.4

Other underlying deposition scenarios

So far we have considered deposition processes which are on average uniform.
There are some situations in which this is clearly inappropriate. In particular where we have sedimentation in depressions (for example filling of ditchfeatures in an archaeological context), peat deposition (see example 1 in Christen et al. (1995)), accumulation restarting after a hiatus or sedimentation rate
varying with climatic instability (as in Blockley et al. (2004)) we would expect
some long-term trends. Likewise if we look at age depth models for ice cores,
these would need to take into account the compression of the ice at lower
levels.
The mathematical formulation for two further sets of deposition models will
be considered here: linear and exponential rises in accumulation rate (theseare easily extended to cover falls in accumulation rate).

A.4.1

Linear rise in accumulation rate (start point known)

If we postulate that the deposition rate starts from zero at time ta and depth
za and increases linearly until time tb and depth zb , we can replace the priorgiven in equation A.21 with the function:0

fl (t ) =

 

2
0

t0 −ta(tb −ta )2

ta < t0 < tbotherwise(A.33)and thus we find that the prior for a simple sequence based on this model isthe constraint prior pH (t) multiplied by:pl (t) ∝

b−1
Y
i=a+1

ti − ta(tb − ta )2!

(A.34)Such a model would be defined in OxCal using the code:
Sequence(){
Zero_Boundary("a");...
Boundary("b");};34For a U Sequence we define the accumulation rate within the time ta < t < tb
to be:dz∝ fl (t) ∝ (t − ta )dtwhich, if we once more choose ta and tb to be our independent variables, gives
an exact equation for t(z) in the range za < z < zb :st = ta + (tb − ta )z − zazb − za(A.35)

This model would be defined in OxCal using the code:
U_Sequence(){
Zero_Boundary("a"){ z= za; };...
Boundary("b"){ z=zb; };};where za and zb are suitably defined.
For a P Sequence based on this deposition scenario we need to to replace
equation A.26 with the integral of fl (t0 ) from equation A.33:
tZi+1
ti

t2 − t2i − 2ta ti+1 + 2ta ti
fl (t )dt = i+1(tb − ta )20!

0(ti+1 − ta )2 − (ti − ta )2=(tb − ta )2!

(A.36)

We can then follow through the same arguments that led to equation A.27
to find that the overall prior associated with the deposition model is the constraint prior pH (t) multiplied by:2
2 k(zi+1 −zi )−1 Qb−1
i=a ((ti+1 − ta ) − (ti − ta ) )i=a+1 (ti − ta )(A.37)(tb − ta )2(k(zb −za )−1)

Qb−1

pl (t|z, k) ∝Such a model would be defined in OxCal with the code:
P_Sequence(k){35Zero_Boundary("a"){ z= za; };...
Boundary("b"){ z=zb; };};where k, za and zb are suitably defined.

A.4.2

Linear rise in accumulation rate (start point not known)

If za (the depth at which the accumulation rate started from zero) is unknown
we can allocate a uniform prior for ta and assume that all of the deposition
events occur in the depth range, za+1 to zb and in the time interval ta+1 to tb .
The underlying prior probability for an accumulation event then becomes:0fl∗ (t ) =

 

2
0

t0 −ta(tb −ta )2 −(ta+1 −ta )2



ta+1 < t0 < tbotherwise(A.38)For a U Sequence we define the accumulation rate within the time ta+1 < t <tb to be:dz
∝ fl∗ (t) ∝ (t − ta )dtHowever in this case za is unknown. We therefore have three independent
variables which we will choose to be ta (the time at which deposition starts)
ta+1 (the first time for which we have depth information) and tb (the end of
the sequence). Solving the differential equation subject to these conditions weobtain an exact equation for t(z) in the range za+1 < z < zb :vuut = ta + (tb − ta )t1 −z − zbza+1 − zb!

ta − ta+1
1−
ta − tbThis model would be defined in OxCal using the code:
U_Sequence(){
Zero_Boundary("a");...
Boundary("b"){ z=zb; };};where zb is suitably defined.
36

2 !

(A.39)

For a P Sequence we need to substitute for equation A.26 the integral of thefunction fl∗ (t0 ):tZi+1fl∗ (t0 )dt0 =ti(ti+1 − ta )2 − (ti − ta )2(tb − ta )2 − (ta+1 − ta )2(A.40)

and using this the overall prior is the constraint prior pH (t) multiplied by:pl∗ (t | za+1 , · · · , zb , k)
Qb−1

∝

i=a+1



k(zi+1 −zi )−1((ti+1 − ta )2 − (ti − ta )2 )(ti − ta )((tb − ta )2 − (ta+1 − ta )2 )k(zb −za+1 )−1(A.41)Such a model would be defined in OxCal with the code:
P_Sequence(k){Zero_Boundary("a");...
Boundary("b"){ z=zb; };
};

A.4.3

Exponential rise in accumulation rate

If instead the deposition rate is expected to increase exponentially we canformulate the model in the following way. We define tb as the endpoint of the
distribution and ta to be one time constant earlier. In this case, the constraints
are different from those given in equation A.1 in that the only constraint on ta
is that ta < tb . The prior for a single deposition event, which can be substituted
for the function in equation A.21, is then:fe (t0 ) =
 e−(t0 −tb )/(ta −tb )tb −ta

0−∞ < t0 < tb

(A.42)

otherwise

Using this, the overall prior for a single-segment sequence, including the constraints is:pe (t) ∝ pH (ta , tb )pH (ta+1 , · · · , tb )

b−1
Y
i=a+1

37



b
exp − ttai −t
−tb







tb − ta(A.43)

Such a model would be defined in OxCal using the code:
Sequence(){
Tau_Boundary("a");...
Boundary("b");};For a U Sequence we define the accumulation rate within the time range−∞ < t < tb to be:dz∝ fe (t) ∝ e−(t−tb )/(ta −tb )dt

We again have three independent variables which we will choose to be ta(one time constant before tb ) ta+1 (the first time for which we have depth
information) and tb (the end of the sequence). Solving the differential equation
subject to these conditions we obtain an exact equation for t(z) in the range
za+1 < z < zb :z − zbt = tb − (ta − tb ) ln 1 −za+1 − zb!

ta+1 − tb
1 − exp −ta − tb!

(A.44)

This model would be defined in OxCal using the code:
U_Sequence(){
Tau_Boundary("a");...
Boundary("b"){ z=zb; };};where zb is suitably defined.
For a P Sequence based on this exponentially rising deposition model weconsider only those deposition events which occur between the depths za+1
and zb . For those the normalised probability is:fe∗ (t0 ) =












t0 −t
a
b

t

−t

exp − t −tbb
1−exp − a+1
t −t
a

b





(tb −ta )0

ta+1 < t0 < tbotherwise

38(A.45)

We must once again substitute for equation A.26 the integral of the functionfe∗ (t0 ):
tZi+1



fe∗ (t0 )dt0 =

ti





−tb
b
− exp − ttai −t
exp − ti+1
ta −tb−tb

−tb
1 − exp − ta+1
ta −tb

(A.46)

We can then follow through the same arguments that led to equation A.27 tofind that the overall prior associated with the deposition model including theconstraints is proportional to:pe (t|za+1 , · · · , zb , k) ∝pH (ta , tb )pH (ta+1 , · · · , tb ) ×
Qb−1

i=a+1

(A.47)

 

 
b


k(zi+1 −zi )−1 exp − tti −t


a −tb

 exp − ti+1 −tb − exp − ti −tb

ta −tb

ta −tb−tb
1 − exp − ta+1
ta −tbtb −tak(zb −za+1 )−1

Such a model would be defined in OxCal with the code:
P_Sequence(k){Tau_Boundary("a");...
Boundary("b"){ z=zb; };};A.5

Summary of mathematical methods

We are able to derive a mathematical prior probability for all of the models
described in table 1. In many cases this is actually only a factor in this probability (in other words it is unnormalised) but fortunately for MCMC analysis
this is all that is required.
Table 1 refers to the equations specifying the prior for each of the main models(together with factors from equations A.29 and A.32 if required). When taken
together with the likelihood distributions (see section 3.1) and the constraints(as given in equation A.1) these define the probability function for the model
as a whole.
The mathematical formulation for Sequence, U Sequence and P Sequencemodels can also be derived for linear and exponential increase or decrease as
39

shown here. In all cases they can be derived so that the prior for the overall
scale of the sequence is neutral. These models are all implemented in OxCal
v4.

B

Binomial treatment of Ice-Core uncertainties

For each uncertain layer which might correspond to either zero or one year
we can treat this as a Bernoulli trial which is either false or true. We assume
the probability that the layer is real is p and the probability that it is not is
q = 1 − p. For a whole segment of core with N uncertain layers it then follows
that the probability that there are exactly n real extra layers in the segment
is given by:
Pp (n|N ) =N!
pn (1 − p)N −nn!(N − n)!

The mean of this distribution is µ = N p as you would expect. It also turns outthat if N is high enough this approximates to a normal or Gaussian distribution
with a variance of:
σ 2 = N p(1 − p)that is we get:(n − N p)2
exp −Pp (n|N ) ≈ q
2N p(1 − p)
2πN p(1 − p)1!

In practice if we assume that p (and therefore q) are 0.5 the situation further
simplifies and for N uncertain layers we get the following results:• The mean expected number of layers is N/2• The probability√ distribution is approximately normal with a standard uncertainty of N /2√
• The 95% confidence range is ± N

40

References
Aguilar, D., Litton, C., O’hagan, A., 2002. Novel statistical model for a piecewise linear radiocarbon calibration curve. Radiocarbon 44 (1), 195 – 212.
Blaauw, M., Christen, J., 2005. Radiocarbon peat chronologies and environmental change. Journal Of The Royal Statistical Society Series C-Applied
Statistics 54, 805 – 816.
Blaauw, M., Heuvelink, G., Mauquoy, D., van der Plicht, J., van Geel, B., 2003.
A numerical approach to c-14 wiggle-match dating of organic deposits: best
fits and confidence intervals. Quaternary Science Reviews 22 (14), 1485 –1500.
Blockley, S., Blaauw, M., Bronk Ramsey, C., van der Plicht, J., 2006a. Assessing uncertainties in age modelling sedimentary records in the Lateglacialand Early Holocene. Quaternary Science Reviews in press.
Blockley, S., Bronk Ramsey, C., Lane, C., 2006b. Depth modelling of theSoppensee core. Quaternary Science Reviews ibid, xx – xxx.
Blockley, S., Lowe, J., Walker, M., Asioli, A., Trincardi, F., Coope, G., Pollard, A., Donahue, R., 2004. Bayesian analysis of radiocarbon chronologies:
examples from the European Lateglacial. Journal of Quaternary Science 19,
159–175.
Bronk Ramsey, C., 1995. Radiocarbon calibration and analysis of stratigraphy:
The OxCal program. Radiocarbon 37 (2), 425 – 430.
Bronk Ramsey, C., 1998. Probability and dating. Radiocarbon 40 (1), 461 –474.
Bronk Ramsey, C., 2001. Development of the radiocarbon calibration program
OxCal. Radiocarbon 43 (2), 355 – 363.
Bronk Ramsey, C., van der Plicht, J., Weninger, B., 2001. ‘Wiggle-Matching’radiocarbon dates. Radiocarbon 43 (2), 381 – 389.
Buck, C., Cavanagh, W., Litton, C. (Eds.), 1996. Bayesian Approach to Interpreting Archaeological Data. John Wiley & Sons.
Buck, C., Kenworthy, J., Litton, C., Smith, A., 1991. Combining archaeological and radiocarbon information - a Bayesian-approach to calibration.
Antiquity 65 (249), 808 – 821.
Buck, C., Litton, C., Smith, A., 1992. Calibration of radiocarbon results pertaining to related archaeological events. Journal of Archaeological Science
19 (5), 497 – 512.
Buck, C., Millard, A. (Eds.), 2001. Tools for Constructing Chronologies, Lecture Notes in Statistics. Vol. 177. Springer.
Christen, J., 1994. Summarizing a set of radiocarbon determinations - a robust
approach. Applied Statistics-Journal of The Royal Statistical Society Series
C 43 (3), 489 – 503.
Christen, J., 2003. Bwigg: an internet facility for bayesianradiocarbon
wiggle-matching.
Internet
Archaeology
7,
http://intarch.ac.uk/journal/issue13/christen index.html.
Christen, J., Clymo, R., Litton, C., 1995. A Bayesian approach to the use of41

C-14 dates in the estimation of the age of peat. Radiocarbon 37 (2), 431 –441.
D’Agostini, G., 2003. Bayesian inference in processing experimental data: principles and basic applications. Reports on Progress in Physics 66 (9), 1383 –1419.
Gilks, W., Richardson, S., Spiegelhalter, D. (Eds.), 1996. Markov Chain Monte
Carlo in Practice. Chapman and Hall.
Hajdas, I., Bonani, G., Zolitschka, B., 2000. Radiocarbon dating of varve
chronologies: Soppensee and Holzmaar Lakes after ten years. Radiocarbon
42 (3), 349 – 353.
Hajdas, I., Ivy, S., Beer, J., Bonani, G., Imboden, D., Lotter, A., Sturm,
M., Suter, M., 1993. AMS radiocarbon dating and varve chronology of lake
Soppensee - 6000 to 12000 C-14 years BP. Climate Dynamics 9 (3), 107 –116.
Heegaard, E., Birks, H., Telford, R., 2005. Relationships between calibrated
ages and depth in stratigraphical sequences: an estimation procedure by
mixed-effect regression. Holocene 15 (4), 612 – 618.
Hughen, K., Baillie, M., Bard, E., Beck, J., Bertrand, C., Blackwell, P., Buck,
C., Burr, G., Cutler, K., Damon, P., Edwards, R., Fairbanks, R., Friedrich,M., Guilderson, T., Kromer, B., McCormac, G., Manning, S., Bronk Ramsey, C., Reimer, P., Reimer, R., Remmele, S., Southon, J., Stuiver, M., Talamo, S., Taylor, F., van der Plicht, J., Weyhenmeyer, C., 2004. Marine04
marine radiocarbon age calibration, 0-26 cal kyr BP. Radiocarbon 46 (3),1059 – 1086.
Jones, M., Nicholls, G., 2001. Reservoir offset models for radiocarbon calibration. Radiocarbon 43 (1), 119 – 124.
Jones, M., Nicholls, G., 2002. New radiocarbon calibration software. Radiocarbon 44 (3), 663 – 674.
Kilian, M., van geel, B., van der plicht, J., 2000. C-14 AMS wiggle matching ofraised bog deposits and models of peat accumulation. Quaternary Science
Reviews 19 (10), 1011 – 1033.
Nicholls, G., Jones, M., 2001. Radiocarbon dating with temporal order constraints. Journal Of The Royal Statistical Society Series C-Applied Statistics
50, 503 – 521.
Reimer, P., Baillie, M., Bard, E., Bayliss, A., Beck, J., Bertrand, C., Blackwell,P., Buck, C., Burr, G., Cutler, K., Damon, P., Edwards, R., Fairbanks, R.,
Friedrich, M., Guilderson, T., Hogg, A., Hughen, K., Kromer, B., Mccormac,
G., Manning, S., Bronk Ramsey, C., Reimer, R., Remmele, S., Southon,J., Stuiver, M., Talamo, S., Taylor, F., van der plicht, J., Weyhenmeyer,C., 2004. Intcal04 terrestrial radiocarbon age calibration, 0-26 cal kyr BP.
Radiocarbon 46 (3), 1029 – 1058.
Walker, M., Coope, G., Sheldrick, C., Turney, C., amd Blockley SPE, L. J.,
Harkness, D., 2003. Devensian Lateglacial environmental changes in Britain:
a multi-proxy record from Llanilid, South Wales, UK. Quaternary Science
Reviews 22, 475–520.
42