In radiocarbon applications, it is often necessary to bind atmospheric radiocarbon curves from different sources. Here, I show an example on how to bind a pre- and a post-bomb curve as well as a forecasted curve for which atmospheric radiocarbon data have not been released yet.

As an example, we will consider the case in which we need an atmospheric radiocarbon curve from the years 500 to 2020 AD. This implies binding first the IntCal09 curve with one of the post-bomb curves, and then doing a forecast from the year 2009 to 2020.

Binding pre- and post-bomb curves

SoilR includes already three prebomb datasets: IntCal09, IntCal13, and IntCal20, and curves for the post-bomb period stored as a list in the object Hua2021. SoilR also includes a handy function called bind.C14curves that allows you to bind a pre- and a post-bomb curve. This function has three arguments: prebomb, postbomb, and time.scale. The prebomb argument accepts either IntCal09, IntCal13, or IntCal20, and we recommend to use IntCal20 because it is the curve currently recognized as official by the radiocarbon community. The postbomb argument requires the specification of the hemispheric zone as described in Hua et al. (2021). The time.scale argument determines whether the results are reported in a time scale in years before present (BP) or anno domini (AD).

In this example, we will bind the IntCal20 and the post-bomb curve for the northern hemisphere zone 2, which corresponds to most of the temperate zone in this hemisphere.

preandpost=bind.C14curves(prebomb=IntCal20,postbomb=Hua2021$NHZone2,time.scale="AD")

The combined dataset goes from about 53,000 years before present (-53050 AD) to the year 2019.375. The dataset is therefore very large:

dim(preandpost)

## [1] 10440     3

For example, if we are only interested in values of Delta14C starting from the year 500, we can select only the values of interest and plot the results

bombcurve=preandpost[preandpost[,1]>=500,1:2]
plot(bombcurve[,1:2],type="l")

plot of chunk unnamed-chunk-3

Forecast of the bomb curve

We use the R package forecast to predict the values of atmospheric radiocarbon until the year 2030.

library(forecast)

The forecast algorithm requires equally spaced values for the radiocarbon time series, however the original time series is not equally spaced, for this reason we need to create a new series based on an interpolation procedure using a spline function

y0=1975 # Initial year used for detecting trend in time series
q=4 #Quaterly time-scale
y=seq(y0,2019.375,by=1/q)
lambda = 1.21 * (10^(-4))

# Function to transform dataseries to F' from F
FpfromF=function(X){ # X must be a data.frame with first column time and second column Fraction Modern
  y=X[,1]
  fM=X[,2]
  Fp=fM*exp((1950-y)*lambda)
  return(data.frame(year=y,Fp=Fp))
}

FpNZ2=FpfromF(Hua2021$NHZone2[,c(1,4)]) #Absolute fraction modern F'
qNZ2=spline(FpNZ2,xout=y) #Spline interpolation of the NH_Zone 2 dataset at a quaterly basis
NZ2=ts(qNZ2$y,start=y0,freq=q) #Transformation into a time-series object

Now, we use the ets function to fit an exponential smoothing state space model, and then perform a forecast based on this model

mNZ2<-ets(NZ2) #Fits an exponential smoothing state space model to the time series
foryrs=10 # Number of years to forecast
fNZ2<-forecast(mNZ2, h=foryrs*q, level=c(69,90)) #Uses the fitted model to forecast 10 years into the future

Now we are ready to combine both curves

bc=data.frame(Year=c(bombcurve[-dim(bombcurve)[1],1],
                     seq(tsp(fNZ2$mean)[1],tsp(fNZ2$mean)[2], by=1/tsp(fNZ2$mean)[3])),
              Delta14C=c(bombcurve[-dim(bombcurve)[1],2],as.numeric((fNZ2$mean)-1)*1000))  # As the forecast is given in fraction modern, we need to convert it to D14C.

This data.frame is now ready to use in SoilR. We can now plot this final dataset.

plot(bc,type="l")

plot of chunk unnamed-chunk-8

And if we want to see the fraction modern curve along with the forecast, we can take the forecasted fraction modern values and add it to the fraction modern from the NHZone2 whose last record is for 2019.375.

Fmodern=data.frame(Year=c(Hua2021$NHZone2$Year,
                     seq(tsp(fNZ2$mean)[1],tsp(fNZ2$mean)[2], by=1/tsp(fNZ2$mean)[3])),
              Fmodern=c(Hua2021$NHZone2$mean.F14C,as.numeric(fNZ2$mean)))


plot(Fmodern, type="l", xlim=c(1980, 2030), ylim=c(0.9,1.4))

plot of chunk unnamed-chunk-9