rice

1 Introduction

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

2 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)

3 Calibration curves

3.1 Plots

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)

To find the IntCal20 radiocarbon year belonging to a certain calendar year, e.g., 900 cal BP:

calBP.14C(900)
## [1] 936  10

4 Calculations

This package also provides functions related to radiocarbon ‘realms’. First there are two functions to calculate radiocarbon ages from pMC values (in this case of a postbomb date):

pMC.age(150, 1)
## [1] -3257    53

and the other way round:

age.pMC(-2300, 40)
##           y    sdev
## [1,] 133.15 0.66138

The same for calculations in the F14C realm:

F14C.age(.150, .01)
##          y   sdev
## [1,] 15240 518.44

and the other way round:

age.F14C(-2300, 40)
##           y      sdev
## [1,] 1.3315 0.0066138

To transfer \(\Delta^{14}C\) (a proxy for atmospheric 14C concentration at t cal BP) to F14C, and the other way around:

F14C.D14C(0.71, t=4000)
## [1] 151.8406
D14C.F14C(152, 4000)
## [1] 0.7100983

These functions can be used to investigate \(\Delta^{14}C\) over time:

cc <- rintcal::ccurve()
cc.Fmin <- age.F14C(cc[,2]+cc[,3])
cc.Fmax <- age.F14C(cc[,2]-cc[,3])
cc.D14Cmin <- F14C.D14C(cc.Fmin, cc[,1])
cc.D14Cmax <- F14C.D14C(cc.Fmax, cc[,1])
par(mar=c(4,3,1,3), bty="l")
plot(cc[,1]/1e3, cc.D14Cmax, type="l", xlab="kcal BP", ylab="")
mtext(expression(paste(Delta, ""^{14}, "C")), 2, 1.7)
lines(cc[,1]/1e3, cc.D14Cmin)
par(new=TRUE)
plot(cc[,1]/1e3, (cc[,2]+cc[,3])/1e3, type="l", xaxt="n", yaxt="n", col=4, xlab="", ylab="")
lines(cc[,1]/1e3, (cc[,2]-cc[,3])/1e3, col=4)
axis(4, col=4, col.axis=4)
mtext(expression(paste(""^{14}, "C kBP")), 4, 2, col=4)

4.1 Contamination

The above functions can be used 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 14C BP would be contaminated with 1% of modern carbon (F14C=1)?

contaminate(5000, 20, .01, 1)
##           y   sdev
## [1,] 4930.9 19.779

The effect of different levels of contamination can also 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()

4.2 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,] 1685 1710 12.9
## [2,] 1719 1732  8.1
## [3,] 1758 1758  0.1
## [4,] 1804 1823  8.2
## [5,] 1832 1892 50.8
## [6,] 1906 1927 14.8

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(130,10)

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

4.3 Marine offsets

Dates on marine material will often have to be calibrated with the Marine20 calibration curve5, and many coastal locations will have an additional regional reservoir offset (deltaR) (Reimer and Reimer 2006)6. 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

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)

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

4.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 Burnell7, Joy Division8 and the Hallstatt Plateau9):

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)


  1. Stuiver, R., Polach, H.A., 1977. Discussion: reporting of 14C data. Radiocarbon 19, 355-363 http://dx.doi.org/10.1017/S0033822200003672↩︎

  2. Reimer, P.J., Brown, T.A., Reimer, R.W., 2004. Discussion: reporting and calibration of post-bomb 14C Data. Radiocarbon 46, 1299-1304 http://dx.doi.org/10.1017/S0033822200033154↩︎

  3. Millard, R., 2014. Conventions for reporting radiocarbon determinations. Radiocarbon 56, 555-559 http://dx.doi.org/10.2458/56.17455↩︎

  4. Reimer, P.J., et al., 2020. The IntCal20 Northern Hemisphere radiocarbon age calibration curve (0-55 cal kBP). Radiocarbon 62, 725-757 http://dx.doi.org/10.1017/S0033822200032999↩︎

  5. Heaton, T.J., et al., 2020. Marine20-the marine radiocarbon age calibration curve (0-55,000 cal BP). Radiocarbon 62, 779-820 http://dx.doi.org/10.1017/RDC.2020.68↩︎

  6. Reimer, P.J., Reimer, R.W., 2006. A marine reservoir correction database and on-line interface. Radiocarbon 43, 461-463. http://dx.doi.org/10.1017/S0033822200038339↩︎

  7. https://www.cam.ac.uk/stories/journeysofdiscovery-pulsars↩︎

  8. https://www.radiox.co.uk/artists/joy-division/cover-joy-division-unknown-pleasures-meaning/↩︎

  9. https://en.wikipedia.org/wiki/Hallstatt_plateau↩︎