Introduction

Radiocarbon dating requires a range of calculations, for example calibration¹²³⁴, translations between pMC, F14C, C14 age and D14C, and assessing the impacts of contamination. This package provides functions to do so in R.

Installation

On first usage of the package, it has to be installed:

install.packages('rice')

The companion data package ‘rintcal’ which has the radiocarbon calibration curves will be installed if it isn’t already. New versions of R packages appear regularly, so please re-issue the above command regularly to remain up-to-date, or use:

update.packages()

To obtain access to the calibration curves and radiocarbon functions, first the package has to be loaded:

library(rice)

Calibration curves

Calibration curves can be plotted:

draw.ccurve()

Or, comparing two calibration curves:

draw.ccurve(1000, 2020, BCAD=TRUE, cc2='marine20', add.yaxis=TRUE)

Or zooming in to between AD 1600 and 2000 (using the BCAD scale):

draw.ccurve(1600, 1950, BCAD=TRUE)

Interesting things happened after 1950, as can be seen by adding a postbomb curve:

draw.ccurve(1600, 2020, BCAD=TRUE, cc2='nh1')

The postbomb curve dwarfs the IntCal20 curve, so we could also plot both on separate vertical axes:

draw.ccurve(1600, 2020, BCAD=TRUE, cc2='nh1', add.yaxis=TRUE)

The calibration curves can also be plotted in the ‘realms’ of F14C, pMC or D14C (see below), e.g.:

draw.ccurve(50000, 35000, realm="D")

Realms

This package also provides functions to translate values between the radiocarbon-relevant ‘realms’ of cal BP (calendar years before AD 1950), BC/AD, C14, F14C, pMC and D14C. The following Table lists the available functions:

As an example of the above functions, the IntCal20 C14 age and error belonging to one or more cal BP or BCAD ages can be found (interpolating linearly where necessary):

calBPtoC14(10.5)

##      [,1] [,2]
## [1,]  168   11

BCADtoC14(1940:1950)

##       [,1] [,2]
##  [1,]  170   11
##  [2,]  174   11
##  [3,]  178   11
##  [4,]  181   11
##  [5,]  185   11
##  [6,]  188   11
##  [7,]  190   11
##  [8,]  193   11
##  [9,]  195   11
## [10,]  197   11
## [11,]  199   11

To translate between BC/AD and cal BP ages, we can use (the last example avoids 0 BC/AD):

BCADtocalBP(2024)

## [1] -74

BCADtocalBP(-1, zero=TRUE)

## [1] 1951

BCADtocalBP(-1, zero=FALSE)

## [1] 1950

D14C (Δ¹⁴C, a proxy for atmospheric ¹⁴C concentration at t cal BP) can be transferred to F¹⁴C, and the other way around:

D14CtoF14C(152, t=4000)

## [1] 0.7100983

F14CtoD14C(0.71, t=4000)

## [1] 151.8406

This can also be done with C14 ages:

C14toD14C(0.71, t=4000)

## [1] 622.1646

D14CtoC14(152, t=4000)

## [1] 2750.1

These functions can be used to investigate Δ¹⁴C over time:

x <- seq(0, 55e3, length=1e3)
cc <- calBPtoC14(x)
Dcc <- calBPtoD14C(x)

par(mar=c(4,3,1,3), bty="l")
plot(x/1e3, Dcc[,1]+Dcc[,2], type="l", xlab="kcal BP", ylab="")
mtext(expression(paste(Delta, ""^{14}, "C")), 2, 1.7)
lines(x/1e3, Dcc[,1]-Dcc[,2])

par(new=TRUE)
plot(x/1e3, (cc[,1]-cc[,2])/1e3, type="l", xaxt="n", yaxt="n", col=4, xlab="", ylab="")
lines(x/1e3, (cc[,1]+cc[,2])/1e3, col=4)
mtext(expression(paste(""^{14}, "C kBP")), 4, 2, col=4)
axis(4, col=4, col.axis=4, col.ticks=4)

Pooling dates

Sometimes, the same material is measured using multiple radiocarbon dates. If we can be sure that the dated material stems from one single age in time (e.g., multiple dates on the same single bone, or perhaps the cereal grains within one bowl which could be assumed to all stem from the same season), then we can check to which degree the dates agree using a Chi2-test (Ward & Wilson 1978)⁵. If they do agree, then a pooled mean can be calculated. For example, take the Shroud of Turin, which was dated multiple times in three different labs:

data(shroud)
shroud

##            ID   y er
## 1   AA-3367.1 591 30
## 2   AA-3367.2 690 35
## 3   AA-3367.3 606 41
## 4   AA-3367.4 701 33
## 5   Ox-2575.1 795 65
## 6   Ox-2575.2 730 45
## 7   Ox-2575.3 745 55
## 8  ETH-3883.1 733 61
## 9  ETH-3883.2 722 56
## 10 ETH-3883.3 635 57
## 11 ETH-3883.4 639 45
## 12 ETH-3883.5 679 51

pool(shroud$y,shroud$er) 

## ! Scatter too large to calculate the pooled mean
## Chisq value: 20.697, p-value 0.037 < 0.05

Zu <- grep("ETH", shroud$ID) # Zurich lab only
pool(shroud$y[Zu],shroud$er[Zu])

## pooled mean: 676.1 +- 23.7
## Chi2: 2.74476177594226, p-value 0.601 (>0.05, OK)

## [1] 676.1406  23.7443

Contamination

To calculate the effect of contamination on radiocarbon ages, e.g. what age would be observed if material with a “true” radiocarbon age of 5000 +- 20 ¹⁴C BP would be contaminated with 1% of modern carbon (F¹⁴C=1)?

contaminate(5000, 20, .01, 1)

##        [,1]   [,2]
## [1,] 4930.9 19.779

Or imagine that you measured a dinosaur bone, dating to far beyond the limit of radiocarbon dating, and the sample is very clean as it contains only 0.5% modern contamination:

contaminate(66e6, 1e6, .005, 1, F14C.er=0.01)

##       [,1]   [,2]
## [1,] 42561 8825.2

The effect of different levels of contamination can be visualised:

real.14C <- seq(0, 50e3, length=200)
contam <- seq(0, .1, length=101) # 0 to 10% contamination
contam.col <- rainbow(length(contam))
plot(0, type="n", xlim=c(0, 55e3), xlab="real 14C age", ylim=range(real.14C), ylab="observed 14C age")
for(i in 1:length(contam))
  lines(real.14C, contaminate(real.14C, c(), contam[i], 1, decimals=5), col=contam.col[i])
contam.legend <- seq(0, .1, length=6)
contam.col <- rainbow(length(contam.legend)-1)
text(50e3, contaminate(50e3, c(), contam.legend, 1), 
  labels=contam.legend, col=contam.col, cex=.7, offset=0, adj=c(0,.8))

If that is too much code for you, try this function instead:

draw.contamination()

We can also calculate the amount of contamination required to ‘explain away’ certain ages. For example, one of the measurements of the Shroud of Turin dates to 591 ± 30 C14 BP. How much contamination would have to be inferred for the material to really date to, say, AD 40?

decontaminate(591, BCADtoC14(40)[1], 1)

## To observe an age of 591 C14 BP, a sample with a true age of 1970 C14 BP would have to be contaminated with 67.39% of carbon with F14C=1

So we’d require the sample to have been contaminated by 35% modern carbon to still date to around AD 40. But what if the sample had been repaired (or burned, somehow ‘adding carbon’) in, say, AD 1400, thus adding material of an not-entirely-modern F14C value (i.e., taking into account both the atmospheric C14 concentrations in AD 1400 and the fact that some of the C14 will have decayed since then)?

decontaminate(591, BCADtoC14(40)[1], BCADtoF14C(1400)[1])

## To observe an age of 591 C14 BP, a sample with a true age of 1970 C14 BP would have to be contaminated with 98.22% of carbon with F14C=0.93173

This means that the dated sample would have to consist almost entirely of Medieval age material - which is exactly what was found by Damon et al. (1989)⁶.

Fractions

Sometimes, one needs to estimate a missing radiocarbon age from a sample which has C14 dates on both the entire sample and on fractions, but where one of the samples was too small to be dated. This can be used in for example soils separated into size fractions, where one of the samples turns out to be too small to be dated. Requires to have the bulk age, the ages of the dated fractions, and the carbon contents and weights of all fractions.

Cs <- c(.02, .05, .03, .04) # carbon contents of each fraction
wghts <- c(5, 4, 2, .5) # weights for all fractions, e.g., in mg
ages <- c(130, 130, 130, NA) # ages of all fractions. The unmeasured one is NA
errors <- c(10, 12, 10, NA) # errors, unmeasured is NA
fractions(150, 20, Cs, wghts, ages, errors) # assuming a bulk age of 150 +- 20 C14 BP 

## estimated C14 age of fraction 4: 519.3 +- 437.3

Calibration

Now on to calibration of radiocarbon dates. We can obtain the calibrated probability distributions from radiocarbon dates, e.g., one of 130 ± 10 C14 BP:

calib.130 <- caldist(130, 10, BCAD=TRUE)
plot(calib.130, type="l")

It is also possible to find the likelihood of a single calendar year for our radiocarbon age, e.g., 145 cal BP:

l.calib(145, 130, 10)

## [1] 0.0052052

For reporting purposes, calibrated dates are often reduced to their 95% highest posterior density (hpd) ranges (please report all, not just your favourite one!):

hpd(calib.130)

##      from   to perc
## [1,] 1927 1906 14.8
## [2,] 1892 1832 50.8
## [3,] 1823 1804  8.2
## [4,] 1758 1758  0.1
## [5,] 1732 1719  8.1
## [6,] 1710 1685 12.9

Additionally, calibrated dates are often reduced to single point estimates. Note however how poor representations they are of the entire calibrated distribution!

calib.2450 <- caldist(2450, 20)
plot(calib.2450, type="l")
points.2450 <- point.estimates(calib.2450)
points.2450

## weighted mean        median          mode      midpoint 
##        2539.9        2512.3        2666.0        2531.5

abline(v=points.2450, col=1:4, lty=2)

Want a plot of the radiocarbon and calibrated distributions, together with their hpd ranges?

calibrate(2450, 20)

Sometimes one would want to smooth a calibration curve to take into account the fact that material has accumulated over a certain time (e.g., a tree, or peat). To do so, a calibration curve can be smoothed to produce a tailor-made calibration curve, after which this one is used to calibrate the date:

mycurve <- smooth.ccurve(smooth=50)
calibrate(2450, 20, thiscurve=mycurve)

Calibrating ‘young’ radiocarbon dates (close to 0 C14 BP) can cause an error, because a bomb curve might be required to capture the youngest ages. Do not worry, there is an option to avoid that error:

try(calibrate(130,30))

## Error in calibrate(130, 30) : 
##   This appears to be a postbomb age (or is close to being one). Please provide a postbomb curve

calibrate(130, 30, bombalert=FALSE)

## Date falls partly beyond calibration curve and will be truncated!

It is also possible to analyse the calibrated probability distributions, e.g. what is the probability (between 0 and 1) that the date stems from material that is of the age of 150 cal BP or younger? Or that it is older than that age?

younger(150, 130, 10)

## [1] 0.7669044

older(150, 130, 10)

## [1] 0.2330956

Marine offsets

Dates on marine material will often have to be calibrated with the Marine20 calibration curve⁷, and many coastal locations will have an additional regional reservoir offset (deltaR) (Reimer and Reimer 2006)⁸. The on-line database at http://calib.org/marine/ is very useful for this; it features the radiocarbon ages and deltaR of many shells of known collection date. The data from this database were downloaded (in August 2024) and can be queried. For example, a map can be drawn with all shell data within certain coordinates:

myshells <- map.shells(S=54, W=-8, N=61, E=0) # the northern part of the UK

## using a medium-scale map. If you prefer a higher-resolution map, please download rnaturalearthhires using the command
## remotes::install_github('ropensci/rnaturalearthhires')

The output can also be queried:

head(myshells)

##       lon   lat  no taxonN   dR dSTD collected res res.error C14 er       lab
## 430 -3.15 55.97 524    146 -141   57      1851 386        58 511 57  SRR-1818
## 433 -5.00 54.17 527    125 -180   47      1890 334        47 444 47 SRR-1821a
## 434 -5.00 54.17 528    125 -129   63      1890 385        64 495 63 SRR-1821b
## 435 -6.00 55.83 529     48 -224   46      1890 290        46 400 46 SRR-1822a
## 436 -6.00 55.83 530     48 -233   52      1890 281        52 391 52 SRR-1822b
## 437 -5.33 57.83 531     90 -171   29      1900 346        30 442 29  SRR-357a
##                                                                                                                                                                                                                                                        ref
## 430 Harkness, D D,, 1983. The extent of the natural 14C deficiency in the coastal environment of the United Kingdom,. Journal of the European Study Group on Physical, Chemical and Mathematical Techniques Applied to Archaeology PACT 8 (IV.9):351-364. 
## 433 Harkness, D D,, 1983. The extent of the natural 14C deficiency in the coastal environment of the United Kingdom,. Journal of the European Study Group on Physical, Chemical and Mathematical Techniques Applied to Archaeology PACT 8 (IV.9):351-364. 
## 434 Harkness, D D,, 1983. The extent of the natural 14C deficiency in the coastal environment of the United Kingdom,. Journal of the European Study Group on Physical, Chemical and Mathematical Techniques Applied to Archaeology PACT 8 (IV.9):351-364. 
## 435 Harkness, D D,, 1983. The extent of the natural 14C deficiency in the coastal environment of the United Kingdom,. Journal of the European Study Group on Physical, Chemical and Mathematical Techniques Applied to Archaeology PACT 8 (IV.9):351-364. 
## 436 Harkness, D D,, 1983. The extent of the natural 14C deficiency in the coastal environment of the United Kingdom,. Journal of the European Study Group on Physical, Chemical and Mathematical Techniques Applied to Archaeology PACT 8 (IV.9):351-364. 
## 437 Harkness, D D,, 1983. The extent of the natural 14C deficiency in the coastal environment of the United Kingdom,. Journal of the European Study Group on Physical, Chemical and Mathematical Techniques Applied to Archaeology PACT 8 (IV.9):351-364. 
##                     taxon    feeding
## 430 Crassostrea virginica suspension
## 433 Crassostrea virginica suspension
## 434 Crassostrea virginica suspension
## 435 Crassostrea virginica suspension
## 436 Crassostrea virginica suspension
## 437 Crassostrea virginica suspension

shells.mean(myshells)

## wmean -143.4, error is max of weighted uncertainty (5.3) & sdev (58.3): 58.3

## [1] -143.4   58.3

You can also extract say the 20 shells closest to a coordinate, e.g., 120 East and 10 North:

myshells <- find.shells(120, 10, 20)

## using a medium-scale map. If you prefer a higher-resolution map, please download rnaturalearthhires using the command
## remotes::install_github('ropensci/rnaturalearthhires')

shells.mean(myshells, distance=TRUE)

## wmean -120.9, error is max of weighted uncertainty (8.4) & sdev (57.4): 57.4

## [1] -120.9   57.4

Multiple calibrations

You can also draw one or more calibrated distributions:

set.seed(123)
dates <- sort(sample(500:2500,5))
errors <- .05*dates
depths <- 1:length(dates)
my.labels <- c("my", "very", "own", "simulated", "dates")
draw.dates(dates, errors, depths, BCAD=TRUE, labels=my.labels, age.lim=c(0, 1800))

or add them to an existing plot:

plot(300*1:5, 5:1, xlim=c(0, 1800), ylim=c(5,0), xlab="AD", ylab="dates")
draw.dates(dates, errors, depths, BCAD=TRUE, add=TRUE, labels=my.labels, mirror=FALSE)

or get creative (inspired by Jocelyn Bell Burnell⁹, Joy Division¹⁰ and the Hallstatt Plateau¹¹):

par(bg="black", mar=rep(1, 4))
n <- 50; set.seed(1)
draw.dates(rnorm(n, 2450, 30), rep(25, n), n:1,
  mirror=FALSE, draw.base=FALSE, draw.hpd=FALSE, col="white",
  threshold=1e-28, age.lim=c(2250, 2800), ex=.8)