MatCal: Open Source Bayesian 14C Age Calibration in Matlab

Bryan C. Lougheed; Stephen P. Obrochta

(1) Overview

Introduction

Radiocarbon (¹⁴C) age calibration is a necessary process for estimating the true age range of a ¹⁴C determination. Atmospheric ¹⁴C activity (Δ¹⁴C) has not been constant throughout history, due to changes in the rate of ¹⁴C production in the upper atmosphere, as well as changes in the Earth’s carbon cycle. Simply knowing the ¹⁴C activity of a given sample and applying the half-life of ¹⁴C is, therefore, in itself not sufficient to determine the true age range; one must also have knowledge of past Δ¹⁴C.

Reconstructions of past Δ¹⁴C have been carried out using independently dated terrestrial archives, including tree-ring, speleothem and lake records. Additionally, independently dated marine archives have also been used to indirectly reconstruct past Δ¹⁴C. The aforementioned Δ¹⁴C records can be combined to produce what is known as a calibration curve, a continuous record of past atmospheric ¹⁴C activity. The most commonly used calibration curves are produced by the IntCal group of scientists, the most recent versions of which were published in 2013 [, ]. The latest calibration curves are IntCal13 (a global atmospheric calibration curve), SHCal13 (a southern hemisphere atmospheric calibration curve) and Marine13 (a calibration curve for the surface mixed layer of the global oceans). All calibration curves include a 1 sigma uncertainty range of possible ¹⁴C ages (¹⁴C yr BP) corresponding to a certain true age, or calibrated age (Cal yr BP).

¹⁴C age determinations produced by radiocarbon laboratories are given in the form of a normal (Gaussian) distribution, whereby mean and 1 sigma values for the ¹⁴C age are reported. The age calibration process involves using the shape of the calibration curve and its associated 1 sigma value to produce a probable calibrated age range for the ¹⁴C age determination. A number of software tools offering ¹⁴C calibration functionality have been authored and are widely used by researchers. Two of the most widely used such tools are the standalone calibration programs Calib [] and the age modelling software OxCal []. Additionally, ¹⁴C calibration routines are built into the R-based, open source age modelling software packages Clam [] and Bacon []. More recently, a Python-based, open source ¹⁴C age calibration script (IOSACal) has also been written [].

Many researchers within geosciences and archaeology make use of Matlab when analysing data. Some of these researchers also carry out ¹⁴C determinations as part of their work. MatCal provides a ¹⁴C calibration function that can easily be implemented into a Matlab-based workflow.

Implementation and architecture

The first methods used for ¹⁴C age calibration, including in the earliest versions of Calib, involved using the so called ‘intercept method’, where one essentially uses the calibration curve as a lookup table to find the calibration curve values directly associated with median, 1 sigma and 2 sigma values of the normal distribution of the ¹⁴C age determination. Later, OxCal [] and BCal [] pioneered ¹⁴C age calibration by using the normal distribution of the ¹⁴C age determination and the 1 sigma range of the calibration curve to calculate a probability density function (PDF) of calibrated ages, which is subsequently analysed using the Bayesian highest posterior density (HPD) method, resulting in a superior calibrated age range determination []. MatCal also employs HPD analysis of a PDF.

The matcal function loads a chosen calibration curve into the Matlab workspace and then linearly interpolates the full range of calibration curve calibrated years to annual resolution. Following Bronk Ramsey [], an age probability (P) is calculated (Eq. 1) for each annual calibrated year (t), where the normal distribution of the ¹⁴C age determination is represented by D ± σ_D (input variables c14age and c14err), and the calibration curve ¹⁴C age for each calibrated year (t) is represented as f(t) ± σ_f(t). Both the calibration curve ¹⁴C data and the user inputted ¹⁴C data are converted to F14C space when calculating P(t), which is a more suitable approach, especially for older ages [].

(Equation 1)

P (t) = e x p (− (D − f (t)) 2 2 (σ D 2 + σ f 2 (t)) σ D 2 + σ f 2 (t)

M1 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ P\left( t \right) = \frac{{exp\left( {\frac{{ - {{\left( {D - f\left( t \right)} \right)}^2}}}{{2(\sigma _D^2 + \sigma _f^2\left( t \right)}}} \right)}}{{\sqrt {\sigma _D^2 + \sigma _f^2\left( t \right)} }} \] \end{document}

A resulting n by 2 array containing calibrated years and associated probabilities (normalised to between zero and one) is subsequently analysed using the Bayesian highest posterior density (HPD) method in order to calculate the calibrated age probability range associated with the 1 sigma (68.27%) and 2 sigma (95.45%) confidence levels (CL). Specifically, the HPD method is implemented as follows: the aforementioned n by 2 array is sorted by probability, after which a cumulative P is calculated. All calibrated years with a cumulative probability falling between 1-CL and 1 are subsequently selected and sorted by calibrated year, after which all continuous annual intervals of calibrated years represent the discrete HPD confidence intervals. The probability associated with each HPD interval is determined by taking the sum of the probability for all calibrated years within the interval itself. All relevant information is also plotted in a Matlab figure with handle ‘14’ (Figure 1).

Figure 1

An example of a ¹⁴C calibration output plot from MatCal. The red shaded area indicates the ¹⁴C age normal distribution. The darker and lighter blue shading respectively indicates the 1 sigma and 2 sigma credible intervals of the calibrated age PDF. The darker and lighter grey shaded areas indicate the calibration curve 1 sigma and 2 sigma confidence intervals, respectively. The green error bars (only shown in the case of IntCal13) indicate the raw ¹⁴C data used to construct the calibration curve.

It is also possible to include a reservoir age (R(t) or ΔR) and reservoir age uncertainty (σ_R) when calibrating a ¹⁴C determination, using the input variables resage and reserr. Reservoir age is applied before the calibration process commences, by altering the ¹⁴C determination (D ± σ_D) to D_r ± σ_r as follows:

D r = D − R and σ r = (σ D 2 + σ R 2) ½

M2 \documentclass[10pt]{article} \usepackage{wasysym} \usepackage[substack]{amsmath} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage[mathscr]{eucal} \usepackage{mathrsfs} \usepackage{pmc} \usepackage[Euler]{upgreek} \pagestyle{empty} \oddsidemargin -1.0in \begin{document} \[ {\rm{D}_{\rm{r}} = {\rm{D}} - {\rm{R}} \ {\rm{and}} \ {\sigma _r} = {(\sigma _D^2 +}\sigma _R^2)}^{1/2} \] \end{document}

Quality control

MatCal has been tested against the two most commonly used ¹⁴C calibration tools, the aforementioned Calib and Oxcal (Table 1). MatCal is found to agree very well with both tools. The very minor differences are most likely due to the slightly different methods used to interpolate the calibration curve. OxCal, for example, interpolates the calibration curve to sub-annual resolution, whereas MatCal interpolates to annual resolution. Seeing as the pre-bomb calibration curves are published in five and ten year resolution, an equal case can be made to apply either annual or sub-annual interpolation. It can be noted that both MatCal and OxCal provide superior calibration for the very oldest dates and/or dates with very large uncertainty, due to the calculation of calibrated age probability in F14C space.

Table 1

¹⁴C age calibration in MatCal 2.0, Calib 7.0.0 and OxCal 4.2 and using IntCal13. Calibration results shown are the 1σ calibrated age ranges in Cal yr BP.

¹⁴C age	MatCal 2.0	Calib 7.0.0	OxCal 4.2

500 ± 30	536 – 514	536 – 514	537 – 513
10000 ± 60	11610 – 11520 11507 – 11326	11610 – 11519 11508 – 11326	11610 – 11519 11508 – 11326
15000 ± 80	18346 – 18109	18345 – 18109	18346 – 18110
25000 ± 100	29177 – 28861	29175 – 28859	29178 – 28860
35000 ± 120	39758 – 39342	39757 – 39339	39760 – 39341
45000 ± 250	48801 – 47953	48792 – 47944	48801 – 47953

Further quality control has been carried out for Matlab in the form of user testing, which has led to the implementation of a number of error and warning messages. For example, if the user attempts to calibrate a ¹⁴C age which produces a calibrated age PDF partially exceeding the range of the chosen calibration curve (e.g. <0 cal yr BP or >50000 cal yr BP in the case of IntCal13), a warning message is returned to the command window and also displayed on any subsequent output plots. Furthermore, user interpretable error messages are returned when, for example, an invalid calibration curve name is selected or if invalid and/or not enough input parameters are set.

(2) Availability

Operating system

Any system capable of running Matlab R2012a or higher.

Programming language

Tested as working in Matlab R2012a, likely to work in earlier versions also.

Additional system requirements

Unknown.

Dependencies

The MatLab plotting toolbox is required if the user decides to make use of the plotting functions.

List of contributors

Bryan C. Lougheed (wrote ¹⁴C calibration routine, plotting module and manuscript)

Stephen P. Obrochta (wrote input parser and rewrote parts of code to increase compatibility across multiple Matlab versions)

Software location

Archive (e.g. institutional repository, general repository) (required – please see instructions on journal website for depositing archive copy of software in a suitable repository)

Name: Figshare

Persistent identifier: doi: 10.6084/m9.figshare.3380710.v1

Licence: GNU General Public License Version 2

Publisher: Bryan C. Lougheed

Version published: 2.0

Date published: 2016-10-04

Code repository

Name: GitHub

Identifier:https://github.com/bryanlougheed/MatCal/

Licence: GNU General Public License Version 2

Date published: 2016-10-04

Language

English

(3) Reuse potential

The matcal function is very suitable for implementation in a Matlab-based geosciences project. For example, users might want to experiment with rapidly calibrating ¹⁴C ages using multiple calibration curves or multiple reservoir ages. Previously, researchers may have had to implement calibration in an external software package, but now it is possible to implement the ¹⁴C calibration into a Matlab workflow. For example, one may wish to batch calibrate multiple ¹⁴C ages, which can be easily done by calling the matcal function within a loop.

The calibrated age probability density function (calprob) and highest posterior density credible intervals (p68_2 and p95_4) outputs are also ideal for implementation within a geochronological age modelling routine, e.g. Monte Carlo techniques.

[B1] Reimer, P J Bard, E Bayliss, A Beck, J W Blackwell, P G Bronk Ramsey, C et al. (2013). IntCal13 and Marine13 radiocarbon age calibration curves 0–50,000 years cal BP Radiocarbon 55(4): 1869–1887, DOI: https://doi.org/10.2458/azu_js_rc.55.16947

[B2] Hogg, A G Hua, Q Blackwell, P G Niu, M Buck, C E Guilderson, T P et al. (2013). SHCal13 Southern Hemisphere Calibration, 0–50,000 Years cal BP Radiocarbon 55(4): 1889–903, DOI: https://doi.org/10.2458/azu_js_rc.55.16783

[B3] Stuiver, M and Reimer, P J (1993). Extended 14C data base and revised CALIB 3.0 14C age calibration program Radiocarbon 35: 215–30, DOI: https://doi.org/10.1017/S0033822200013904

[B4] Bronk Ramsey, C (1995). Radiocarbon calibration and analysis of stratigraphy; the OxCal program Radiocarbon 37(2): 425–30, DOI: https://doi.org/10.1017/S0033822200030903

[B5] Blaauw, M (2010). Methods and code for “classical” age-modelling of radiocarbon sequences Quat Geochronol 5(5): 512–8, DOI: https://doi.org/10.1016/j.quageo.2010.01.002

[B6] Blaauw, M Christen, J A et al. (2011). Flexible paleoclimate age-depth models using an autoregressive gamma process Bayesian Anal 6(3): 457–74, DOI: https://doi.org/10.1214/ba/1339616472

[B7] Costa, S and Gutiérrez-Roig, M (2016). IOSACal: v0.3. Zenodo DOI: https://doi.org/10.5281/zenodo.49858

[B8] Buck, C E, Christen, J A and James, G N (1999). BCal: an on-line Bayesian radiocarbon calibration tool Internet Archaeology 7 Available at: http://intarch.ac.uk/journal/issue7/buck.

[B9] Telford, R J, Heegaard, E and Birks, H J B (2004). The intercept is a poor estimate of a calibrated radiocarbon age The Holocene 14(2): 296–8, DOI: https://doi.org/10.1191/0959683604hl707fa

[B10] Bronk Ramsey, C (2008). Radiocarbon dating: Revolutions in Understanding Archaeometry 50(2): 249–75, DOI: https://doi.org/10.1111/j.1475-4754.2008.00394.x

Journal of Open Research Software

Software Metapapers

MatCal: Open Source Bayesian ¹⁴C Age Calibration in Matlab

Abstract