Illustration

We provide in Figure 3 an illustrative example of the need to incorporate calendar age uncertainty. For simplicity, in this example, we fit and model our spline in the same domain. We consider a straightforward underlying function which we wish to reconstruct from 100 noisy observations $\left\{ {({Y_i},{T_i})} \right\}_{i = 1}^{100}$

$$f(\theta ) = {e^{ - \theta }}\sin (4\pi \theta ) + \cos (2\pi \theta )),$$

Figure 3 The importance of recognizing calendar age uncertainty (and correctly representing covariance within that uncertainty) when constructing a curve based on data arising from records with different observed timescales. Panel (a) shows the true, underlying, calendar ages of the data in two different records; while panel (b) shows a joint shift in the observed calendar ages within record 2. In such a case that observed timescales in two records are offset from one another, if we ignore this calendar age uncertainty then our spline estimate will introduce spurious variation as it flips between the records as in panel (c). Conversely if we incorporate such calendar age uncertainty and accurately represent it then we can still recover the underlying function accurately, as shown in (d).

where our underlying ${\theta _i} \sim Beta(1.1,1.1)$ and observed ${Y_i} \sim N(\,f({\theta _i}{),0.05^2})$ for $i = 1, \ldots ,100$. Further, let us assume that these observations arise from two different sediment cores, 50 from core 1 (shown as black triangles) and 50 from core 2 (shown as red dots) as seen in panel (a). Within core 1, the calendar ages are known absolutely. However, within core 2 the observed timescale is somewhat shifted/biased so that when we observe the calendar ages $T$ in this core they all share the same joint shift from their true values, i.e.

$${T_i} = {\theta _i}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,{\rm{if}}\;i\;{\rm{is\,in\,core\,1,}}$$

$${T_i} = {\theta _i} + \psi\,\,\,\,\,\,{\rm{if}}\;i\;{\rm{is\,in\,core\,2,}}$$

where $\psi \sim N{(0,0.1^2})$. The observed pairs $\left\{ {({Y_i},{T_i})} \right\}_{i = 1}^{100}$ are presented in panel (b) showing the shift in timescale within core 2. If we attempt to reconstruct the function using splines based upon our observed $\left\{ {({Y_i},{T_i})} \right\}_{i = 1}^{100}$ without recognizing the calendar age uncertainty, i.e. assuming both cores are on the same timescale, then we obtain the estimate in the panel (c). The estimate is spuriously variable as the spline will try to pass near all of the data on their non-comparable timescales. Conversely, if we recognize that the timescale in core 2 may need to be shifted onto the timescale of core 1, we obtain the estimate shown in panel (d). If we inform our Bayesian method that the calendar age observations in core 2 are subject to uncertainty, it will estimate the size of the joint shift needed to register core 2 onto the true timescale of core 1 and simultaneously reconstruct the true underlying function.

Implications for Visually Assessing Curve Fit

While it is important to incorporate errors-in-variables in curve construction, and accurately represent any dependence in the calendar age uncertainties, this does cause some difficulty in visually assessing the quality of fit between the raw constituent data and the final IntCal curve. Since the proposed method can shift the calendar ages of the observed data left and right within their uncertainties as described, one cannot estimate the final curve’s goodness-of-fit by eye using only the radiocarbon age axis if the raw data are only plotted at their initial, observed calendar ages. The method will have tried to align shared features even if they occur at somewhat different observed times within the different sets. This should be taken into consideration when viewing the final curve against the raw data. We provide the posterior calendar age estimates for the true calendar ages ${\theta _i}$ for all the data. The difficulty of assessing the fit of the curve by eye is further compounded by the offsets in $^{\rm14}{\rm{C}}$ described in Section 4.4, in particular the marine reservoir ages which vary significantly over time and so will change as the MCMC updates the calendar ages of the data.

Blocking and Additive Errors

Our model for the tree-ring determinations, taking account of blocking, considers the observed ${F_i}$ as

$${F_i} = {1 \over {{m_i}}}\left\{ {\sum\limits_{j=0}^{{m_i} - 1} f({\theta _i} + j)} \right\} + {\varepsilon _i} + {\eta _i}\quad {\rm{for}}\ i=1, \ldots ,n$$

where ${m_i}$ is the number of annual rings in the block for that determination; ${\theta _i}$ is the most recent year (start ring) of the block; and ${\varepsilon _i} \sim N(0,\sigma _i^2)$ where ${\sigma _i}$ is the uncertainty reported by the laboratory. This implicitly assumes that the tree rings over the section of ${m_i}$ rings are approximately equal in width, so that approximately equal annual amounts of wood were deposited and have ended up in the final measured sample. The ${\eta _i} \sim N(0,\tau _i^2)$ are independent random effects, with ${\tau _i}$ unknown, that represent potential over-dispersion in the $^{\rm14}{\rm{C}}$ determinations.

Additive Errors Model

As described in Section 5.1 we model ${\tau _i}$, the standard deviation of the additive over-dispersion on the $i$th determination, as proportional to the square-root of the underlying ${F^{14}}{\rm{C}}$, i.e.

$$\tau _i^2 = {\tau ^2}\ {1 \over {{m_i}}}\sum\limits_{j = 0}^{{m_i} - 1} f({\theta _i} + j).$$

The individual values of ${\eta _i}$ are not of interest and so we integrate them out, updating only the overall constant of proportionality $\tau$. Formally, a structure where ${\tau _i}$ depends upon the underlying curve means this term should be incorporated in the update to $\bi{\beta}$ within our sampler—a change to the curve (via its spline coefficients $\bi{\beta}$) would mean a change to the additive uncertainty ${\tau _i}$ on all observations. However, this would make the $\bi{\beta}$ update no longer a Gibbs step. Since, given the overall multiplier $\tau$, the change in ${\tau _i}$ between two potential curves is minimal, when updating $\bi{\beta}$ we consider ${\tau _i}$ as fixed; see Section 4.3.3 for more details.

Floating Tree-Ring Sequences

Several late-glacial tree sequences measured by Aix, ETH and Heidelberg (Reinig et al. Reference Reinig, Nievergelt, Esper, Friedrich, Helle, Hellmann, Kromer, Morganti, Pauly, Sookdeo, Tegel, Treydte, Verstege, Wacker and Büntgen2018; Capano et al. Reference Capano, Miramont, Shindo, Guibal, Marschal, Kromer, Tuna and Bard2020 in this issue) were dated by wiggle-matching so, while their internal relative chronologies were known precisely, their absolute calendar ages were somewhat uncertain. For such a floating tree-ring sequence, the true calendar ages are:

$$\theta _j^{{\rm{float}}} = \xi + {\delta _j}\ {\rm{for}}\ j\ {\rm{in\ floating\ sequence}}.$$

where $\xi$ is the unknown start date (i.e. most recent year) of the floating sequence and ${\delta _j}$ is the precisely-known internal ring-count chronology. We place a discretized (integer) normal prior on $\xi \sim N(M,{N^2})$. Here, the prior mean $M$ and variance ${N^2}$ for the start date of each such tree-ring sequence were provided. All of the other dendrodated determinations were assumed to have their calendar age ${\theta _i}$ known absolutely.

Finding Values of ${\bi{F^{14}}{\bi{C}}}$ at Each Year Represented in a Block

We begin by describing how we can relate each blocked determination to the spline coefficients $\bi{\beta}$. Let ${\bi{\theta}}^A = ({\theta _1},{\theta _1} + 1, \ldots ,{\theta _1} + {m_1} - 1,{\theta _2}, \ldots ,{\theta _2} + {m_2} - 1, \ldots ,{\theta _n} + {m_n} - {1)^T}$ be a vector concatenating every calendar year, including any potential blocking, represented in all $n$ determinations. Also, define vectors ${\bi{c}}}_{\bi{\theta}}^A$ and ${\bi{e}}}_{\bi{\theta}}^A$ so that

$$c_{\theta ,l}^A = {1 \over {1000}}{e^{ - \theta _l^A/8267}}\quad {\rm{and}}\ e_{\theta ,l}^A = {e^{ - \theta _l^A/8267}}.$$

Then, create the matrix ${\cal B}_{\bi{\theta}}^A$ containing the spline bases evaluated at each value in ${\bi{\theta}}^A$:

$${\cal B}_{\bi{\theta}}^A = \left( {\matrix{ {{B_1}({\theta _1})} \ldots {{B_K}({\theta _1})} \cr {{B_1}({\theta _1} + 1)} \ldots {{B_K}({\theta _1} + 1)} \cr \vdots {} \vdots \cr {{B_1}({\theta _n} + {m_n} - 1)} \ldots {{B_K}({\theta _n} + {m_n} - 1)} \cr } } \right),$$

so that, at the calendar years ${\bi{\theta}} ^A$, the modeled value of ${\Delta ^{14}}{\rm{C}}$ is ${\bi{g}}}^A = {\cal B}_{\bi{\theta}} ^A\bi{\beta}$ where ${\bi{\beta}}$ represents the spline coefficients. The modeled value of ${F^{14}}{\rm{C}}$ at an individual calendar age $\theta _l^A$ then becomes $f_l^A = c_{\theta ,l}^Ag_l^A + e_{\theta ,l}^A$ and at the full ${\bi{\theta}} ^A$, ${{\bf{f}}}^A = {\cal B}_{\bi{\theta}} ^\dagger {\bi{\beta}} + {\bi{e}}}_{\bi{\theta}} $ where ${\cal B}_{\bi{\theta}}^\dagger$ is formed by multiplying each row in ${\cal B}_{\bi{\theta}} ^A$ by the corresponding element in ${\bi{c}}}_{\bi{\theta}} ^A$.

A Blocking Matrix

For each blocked determination, we now wish to average the values of F ¹⁴C according to the multiple years represented in that block. Consider an averaging matrix ${\bf{M}}$ corresponding to averaging the individual values in ${{\bf{f}}^A}$:

$${\bf{M}} = \left( {\matrix{ {1/{m_1}} \ldots {1/{m_1}} 0 \ldots 0 \ldots \ldots \cr 0 0 0 {1/{m_2}} \ldots {1/{m_2}} 0 \ldots \cr \vdots {} {} {} {} {} {} {} \cr } } \right).$$

Our observed vector of ${F^{14}}{\rm{C}}$ measurements ${\bf{F}}$ thus becomes

$$\eqalign{{\bf{F}} = {\bf{M}}{\bi{\cal B}}_{\bi{\theta}} ^{\dagger} {\bi{\beta}} + {\bf{M}}{\bi{e}}_{\bi{\theta}} + {\bi{\varepsilon}} + {\bi{\eta}}, \cr = {\bi{{\cal B}}}_{\bi{\theta}} ^ \star {\bi{\beta}} + {\bf{e}}_{\bi{\theta}} ^ \star + {\bi{\varepsilon}} + {\bi{\eta}}.}$$

Note that both the matrix ${\cal B}_{\bi{\theta}} ^ \star$ and ${\bf{e}}_{\bi{\theta}} ^ \star$ depend upon the true calendar ages. However, since for the vast majority of the tree-rings, the calendar ages are known exactly, one only has to calculate ${\cal B}_{\bi{\theta}} ^ \star$ and ${\bf{e}}_{\bi{\theta}} ^ \star$ once as it is then fixed. The final ${\bi{{\cal B}}}_{\bi{\theta}} ^ \star$ has dimension $K \times n$ and hence the Gibbs updating of ${\bi{\beta}}$ is the same speed as if one had no blocking. To incorporate the floating tree-ring sequences, for any potential sequence start date $\xi$ we only need update a small submatrix of ${\bi{{\cal B}}}_{\bi{\theta}} ^ \star$ and ${\bf{e}}_{\bi{\theta}} ^ \star$ corresponding to these observations. The required updates for all sequence start dates $\xi$ covering the prior range can be calculated and stored before commencing curve estimation.

Posterior

Using the above hierarchical structure, the joint posterior is proportional to

$$\scale95%\eqalign{[\xi ,{\bi{\beta}} ,\lambda ,\tau |{\bf{F}},{\bi{\theta}} ] \propto [{\bf{F}}|{\bi{\beta}},{\bi{\theta}},\xi ,\tau ][{\bi{\beta}}|\lambda ][\lambda ][\tau ][\xi ] \cr \propto \det {({{\bf{W}}_\tau })^{{1 \over 2}}}\ \exp \left\{ { - {1 \over 2}\left[ {\left\{ {{{({\bf{F}} - {\bf{e}}_{\bi{\theta}}^ \star - {\bi{{\cal B}}}_{\bi{\theta}}^ \star {\bi{\beta}})}^T}{{\bf{W}}_\tau }({\bf{F}} - {\bf{e}}_{\bi{\theta}}^ \star - {\bi{{\cal B}}}_{\bi{\theta}}^ \star {\bi{\beta}})} \right\} + {{{{(\xi - M)}^2}} \over {{N^2}}} + \lambda {{\bi{\beta}}^T}{\bf{D}}{\bi{\beta}}} \right]} \right\} \cr \quad \times {\lambda ^{A + {\rm{rank}}({\bf{D}})/2 - 1}} \times \exp \left( { - {\lambda \over B}} \right),}$$

where ${{\bf{W}}_\tau } = {\rm{diag}}{(\sigma _1^2 + \tau _1^2, \ldots ,\sigma _n^2 + \tau _n^2)^{ - 1}}$ includes the additive error representing the over-dispersion. We can sample from this using Metropolis-within-Gibbs but where most of the updates can be performed directly via Gibbs and only a few require Metropolis-Hastings (MH).

Gibbs Updating ${\bi{\beta}}|{\bf{F}},{\bi{\lambda}},{\bi{\xi}},{\bi{\tau}}$

$$\pi ({\bi{\beta}}|{\bf{F}},\lambda ,\tau ) \propto \exp \left\{ { - {1 \over 2}\left[ {\left\{ {{{({\bf{F}} - {\bf{e}}_{\bi{\theta}}^ \star - {\bi{{\cal B}}}_{\bi{\theta}}^ \star {\bi{\beta}})}^T}{{\bf{W}}_\tau }({\bf{F}} - {\bf{e}}_{\bi{\theta}}^ \star - {\bi{{\cal B}}}_{\bi{\theta}}^ \star {\bi{\beta}})} \right\} + \lambda {{\bi{\beta}}^T}{\bf{D}}{\bi{\beta}}} \right]} \right\}$$

For the purposes of maintaining a Gibbs update here, we treat ${{\bf{W}}_\tau } = {\rm{diag}}{(\sigma _1^2 + \tau _1^2, \ldots ,\sigma _n^2 + \tau _n^2)^{ - 1}}$ as fixed and independent of ${\bi{\beta}}$, i.e. we ignore the dependence of the individual over-dispersion ${\tau _i}$ on the curve. With this minor approximation, algebra gives a direct complete conditional

$${\bi{\beta}}|{\bf{F}},\lambda \sim MVN\left( {{\bf{Q}}{\bi{{\cal B}}}_{\bi{\theta}}^{ \star T}{{\bf{W}}_\tau }({\bf{F}} - {\bf{e}}_{\bi{\theta}}^ \star ),{\bf{Q}}} \right),$$

where

$${\bf{Q}} = ({\bi{{\cal B}}}_{\bi{\theta}}^{ \star T}{{\bf{W}}_\tau }{{\bi{{\cal B}}}^ \star } + \lambda {\bf{D}}{)^{ - 1}}.$$

Here large computational savings can be made if we calculate ${\bf{Q}}({\cal B}_{\bi{\theta}}^{ \star T}{{\bf{W}}_\tau }{\bf{F}})$ rather than $({\bf{Q}}{\cal B}_{\bi{\theta}}^{ \star T}){{\bf{W}}_\tau }{\bf{F}}$ since ${\bf{Q}}$ is $K \times K$ while ${\cal B}_{\bi{\theta}}^{ \star T}$ is $K \times n$ (and $n \gg K$).

Gibbs Updating ${\bi{\lambda}} |{\bi {\beta}},{\bi{A}},{\bi{B}}$

Due to conjugacy of the $\lambda$ prior:

$$\pi (\lambda |{\bi{\beta}}) \propto {\lambda ^{A + {\rm{rank}}({\bf{D}})/2 - 1}}\exp \left\{ { - \lambda \left( {{{{{\bi{\beta}}^T}{\bf{D}}{\bi{\beta}}} \over 2} + {1 \over B}} \right)} \right\},$$

i.e. $\lambda |{\bi{\beta}} \sim Ga(A',B')$ where $A' = A + {\rm{rank}}({\bf{D}})/2$ and $B' = {\left( {{{{{\bi{\beta}}^T}{\bf{D}}{\bi{\beta}}} \over 2} + {1 \over B}} \right)^{ - 1}}$.

MH Updating ${\bi{\xi}} |{\bf{F}},{\bi{\lambda}},{\bi {\beta}},{\bi{\tau}}$

To update the start of our floating tree-ring sequences, we use an MH step:

• Propose $\xi ' \sim MVN(\xi ,{\sigma _{{\rm{prop}}}})$;

• Calculate ${\bi{{\cal B}}}_{\bi{\theta}}^ \star$ and ${\bi{{\cal B}}}_{\bi{\theta}'}^ \star$ where ${\bi{\theta}}'$ are the calendar ages of the floating sequence with the proposed new start date $\xi' $. Similarly calculate ${\bf{e}}_{\bi{\theta}}^ \star$ and ${\bf{e}}_{\bi{{\theta'}}}^ \star$. Accept $\xi'$ according to Hastings Ratio:

$$HR = \min \left\{ {1,{{\pi (\xi' |{\bf{F}},{\bi{\beta}},\tau )} \over {\pi (\xi |{\bf{F}},{\bi{\beta}},\tau )}}} \right\},$$

where

$$\pi (\xi |{\bf{F}},{\bi{\beta}},\tau ) = \exp \left\{ - {1 \over 2}\left[ \left\{ ({\bf{F}} - {\bf{e}}_{\bi{\theta}}^ \star - {\bi{{\cal B}}}_{\bi{\theta}}^ \star {\bi{\beta}})^T{\bf{W}}_\tau ({\bf{F}} - {\bf{e}}_{\bi{\theta}}^ \star - {\bi{{\cal B}}}_{\bi{\theta}}^ \star {\bi{\beta}}) \right\} + {(\xi - M)^2 \over N^2} \right] \right\}.$$

MH Updating ${\bi{\tau}} |{\bf{F}},{\bi{\lambda}} ,{\bi{\beta}}$

This also requires an MH update step. Let the residuals between the observed ${\bf{F}}$ and fitted values ${\bf{\hat F}} = {\bi{{\cal B}}}_{\bi{\theta}}^ \star {\bi{\beta}} + {\bf{e}}_{\bi{\theta}}^ \star$ be

$${\bf{R}} = {\bf{F}} - {\bi{{\cal B}}}_{\bi{\theta}}^ \star {\bi{\beta}} - {\bf{e}}_{\bi{\theta}}^ \star .$$

Then, we have independently ${R_i} \sim N(0,{\sigma _i^2} + \tau _i^2)$ where ${\tau _i} = \tau \sqrt {{{\hat F}_i}}$. With our prior $\tau \sim \pi (\tau )$ we update:

• Sample $\tau ' \sim {N_ + }(\tau ,{\eta ^2})$;

• Accept according to Hastings Ratio:

$$HR = \min \left\{ {1,{{f({\bf{R}};\tau ')\pi (\tau ')\Phi (\tau /\eta )} \over {f({\bf{R}};\tau )\pi (\tau )\Phi (\tau '/\eta )}}} \right\},$$

where $f({\bf{R}};\tau ) = \prod\nolimits_{i = 1}^n \phi ({R_i},0,\sigma _i^2 + \tau _i^2)$ and ${\tau _i} = \tau \sqrt {{{\hat F}_i}}$. Note the asymmetric proposal adjustment.

Choice of Splines

For the portion of the curve from 0 cal kBP back to approximately 14 cal kBP, we placed 2000 knots at the unique calendar age quantiles of the constituent determinations (each blocked determination was represented by its midpoint) to permit variable smoothing dependent upon data density. In regions where there were more data this allowed us to pick out finer scale variation. Such a selection enabled close to annual resolution in the detail of the final curve while still maintaining computational feasibility in curve construction. We also placed additional knots in the vicinity of the known Miyake-type events at 774–5 AD, 993–4 AD and $ \sim$660 BC (Miyake et al. Reference Miyake, Nagaya, Masuda and Nakamura2012, Reference Miyake, Masuda and Nakamura2013; O’Hare et al. Reference O’Hare, Mekhaldi, Adolphi, Raisbeck, Aldahan, Anderberg, Beer, Christl, Fahrni, Synal, Park, Possnert, Southon, Bard and Muscheler2019) to enable us to capture more rapid variation at these times. Around the calendar age of each such event, we added an extra 4 jittered knots, i.e. after addition of small amounts of random noise—for the $ \sim$660 BC event we used a slightly larger jitter due to its slightly less certain timing. If there is no event at these times then the curve will not introduce one but, if there is, then these additional knots will allow it to be picked up more clearly. All knot locations then remained fixed throughout the sampler.

Outlier Screening

In order to identify potential outliers, after fitting a preliminary curve, scaled residuals (incorporating potential blocking) were calculated for each datum, i.e. for an annual measurement

$${r_i} = {{{F_i} - \hat f({\theta _i})} \over {\sqrt {\sigma _i^2 + \upsilon {{({\theta _i})}^2}} }},$$

where $\upsilon ({\theta _i})$ is the posterior standard deviation on the calibration curve at ${\theta _i}$ cal BP, the calendar age of that datum. These residuals were then combined for each dataset ${\cal K}$ into a scaled deviance, ${Z_{\cal K}} = \sum\nolimits_{i \in {\cal K}} r_i^2$, and compared with a $\chi _{{n_{\cal K}}}^2$ where ${n_{\cal K}}$ are the number of determinations in that set. The mean offset of each dataset from the preliminary curve was also calculated ${\nu _{\cal K}} = {1 \over {{n_{\cal K}}}}\sum\nolimits_{i \in {\cal K}} \left\{ {{F_i} - \hat f({\theta _i})} \right\}$. Sets which had low $p$-values for their deviance and high mean offsets were then discussed with the data set providers as to whether they should be included in the final IntCal or not (see Bayliss et al. Reference Bayliss, Marshall, Dee, Friedrich, Heaton and Wacker2020 in this issue, for more details).

Run Length

The MCMC was run for 50,000 iterations. The first 25,000 of these iterations were discarded as burn-in. The remaining 25,000 were used to create the final curve (thinned to every 10th) and passed to the older part of the curve to create a seamless merging between the curve sections as discussed in Section 4.4.1. Since the main update steps (i.e. $\bi{\beta}$ and $\lambda$) are Gibbs this was felt to be sufficient although convergence was further assessed by initializing the sampler at different values and comparing the resultant curve estimates visually.

Uncertain Calendar Ages

All our data in the older part of the curve have estimates (either in the form of noisy observations or priors) for their calendar ages rather than absolutely known values. For the corals and speleothems these estimates are obtained via U-Th dating and are considered independent; for the Cariaco Basin they are either provided by varve counting or elastic palaeoclimate tie-pointing (Heaton et al. Reference Heaton, Bard and Hughen2013); similar elastic palaeoclimate tie-pointing provides the ages for the Pakistan and Iberian margins; for Lake Suigetsu they are provided by wiggle matching; while the floating tree-ring sequences have internally known chronologies but no absolute age estimates. We assume all our calendar age estimates are (potentially multivariate) normal. For a dataset ${\bi{{\cal K}}}$ on which we have noisy observations ${{\bf{T}}_{\cal K}}$ of its calendar ages e.g. U-Th dated corals and the varve counted section of Cariaco, we model:

$${{\bf{T}}_{\cal K}} \sim MVN({{\bi{\theta}}_{\cal K}},{{\bi{\Psi}}_{{\cal K},T}})\ {\rm{so\ that}}\ [{{\bf{T}}_{\cal K}}|{{\bi{\theta}}_{\cal K}}] \propto \exp \left\{ { - {1 \over 2}\left\{ {{{({{\bf{T}}_{\cal K}} - {{\bi{\theta}}_{\cal K}})}^T}{\bi{\Psi}}_{{\cal K},T}^{ - 1}({{\bf{T}}_{\cal K}} - {{\bi{\theta}}_{\cal K}})} \right\}} \right\}.$$

Here ${{\bi{{\Psi}}}_{{\cal K},T}}$ is the pre-specified and fixed covariance matrix that encodes the dependence within the calendar age observations for that dataset. For such datasets we place an uninformative prior on ${\bi{{\theta}}_{\cal K}}$. For Lake Suigetsu and the datasets where calendar age estimates are obtained via palaeoclimate tie-pointing, no actual calendar ages were observed but rather we have a prior on their true values:

$$[{{\bi{\theta}}_{\cal K}}] \propto \exp \left\{ { - {1 \over 2}\left\{ {{{({{\bf{T}}_{\cal K}} - {{\bi{\theta}}_{\cal K}})}^T}{\bi{\Psi}}_{{\cal K},T}^{ - 1}({{\bf{T}}_{\cal K}} - {{\bi{\theta}}_{\cal K}})} \right\}} \right\},$$

and the ${{\bf{T}}_{\cal K}}$ and ${\Psi_{{\cal K},T}}$ now represent our prior mean and covariance for the calendar ages. As we have no datasets on which we have both priors and noisy calendar age observations, these two types of calendar age estimate become equivalent for the purposes of updating. Details of how to construct the various covariance matrices for the case of varve counting and wiggle matching can be found in Niu et al. (Reference Niu, Heaton, Blackwell and Buck2013); for the datasets with priors obtained by tie-pointing in Heaton et al. (Reference Heaton, Bard and Hughen2013); and for Lake Suigetsu in Bronk Ramsey et al. (Reference Bronk Ramsey, Heaton, Schlolaut, Staff, Bryant, Brauer, Lamb, Marshall and Nakagawa2020 in this issue).

Offsets: Reservoir Ages and Dead Carbon Fractions

Several of our datasets do not provide direct measurements of the atmosphere but are instead offset. The speleothem records contain carbon obtained from dripwater which has passed through old limestone and so is a mixture of atmospheric CO$_2$ and dissolved old (dead) carbon that has no $^{\rm14}{\rm{C}}$ activity. The offset in radiocarbon age this creates is called a dead carbon fraction (dcf) and is specific to the speleothem. Similarly the marine records have an atmospheric offset called a marine reservoir age (MRA) that arises due to both the limited gas exchange with the atmosphere and ocean circulation drawing up deep old carbon. These MRAs are location specific and vary over time. For both of these types of data we can incorporate the offset (either dcf or MRA) in the radiocarbon age domain as:

$${h_{{\rm{offset}}}}(\theta ) = {h_{{\rm{atmos}}}}(\theta ) + {r_{\cal K}}(\theta ),$$

where ${h_{{\rm{offset}}}}(\theta )$ is the radiocarbon age at time $\theta$ cal BP in the offset environment; ${h_{{\rm{atmos}}}}(\theta )$ our atmospheric calibration curve in the radiocarbon age domain; and ${r_{\cal K}}(\theta )$ our offset. This alters our observational model so that, in the ${F^{14}}{\rm{C}}$ domain, for determinations from dataset ${\cal K}$ and with ${F^{14}}{\rm{C}}$ domain calibration curve $f(\theta )$, the offset becomes a multiplier,

$${F_i} = {e^{ - {r_{\cal K}}({\theta _i})/8033}}f({\theta _i}) + {\varepsilon _i}.$$

These offsets must be estimated within curve construction to adaptively synchronize the different records. For IntCal20, speleothem dcfs were considered to be approximately constant over time but with an unknown level. MRAs were modeled as time-varying, providing a step forward from IntCal09 and IntCal13. Initial MRA estimates for each of our locations were obtained by creating a preliminary atmospheric calibration curve using the same Bayesian spline methodology but based only on the Hulu Cave record (Southon et al. Reference Southon, Noronha, Cheng, Edwards and Wang2012; Cheng et al. Reference Cheng, Edwards, Southon, Matsumoto, Feinberg, Sinha, Zhou, Li, Li, Xu, Chen, Tan, Wang, Wang and Ning2018). This Hulu-based curve was then used as a forcing for an enhanced Hamburg Large Scale Geostrophic (LSG) ocean general circulation model to provide estimates ${\delta _{\cal K}}(\theta )$ of MRA at each given location (see Butzin et al. Reference Butzin, Heaton, Köhler and Lohmann2020 in this issue, for details). Due to coastal effects, these LSG estimates were considered to define the shape of the MRA, i.e. relative changes over time, but subject to a constant potential coastal shift. Ideally we might cycle the process of building the curve and re-estimating MRAs several times but the LSG model had a run time of several weeks and so this was not feasible. Notwithstanding the above attempts to accurately model dcf and MRA, we recognized that there was further variation in the offsets over time we were not able to fully capture. This was incorporated through the addition of further independent variability in the offsets from one calendar age to the next. Consequently our model for the offsets is:

$$\eqalign{{\rm{Speleothems(dcfs)}}: \ {r_{\cal K}}({\theta _i}) \sim N({\nu _{\cal K}},\zeta _{\cal K}^2),\cr {\rm{Marine\ records(MRAs)}}:\ {r_{\cal K}}({\theta _i}) \sim N({\nu _{\cal K}} + {\delta _{\cal K}}({\theta _i}),\zeta _{\cal K}^2),}$$

where ${\zeta _{\cal K}}$ is the further independent variability in offsets from year-to-year around our LSG/constant dcf model and added to data before curve construction. We place priors on the constant dcf mean/coastal shifts ${\nu _{\cal K}}$:

$$\pi ({\nu _{\cal K}}) \sim N({\rho _{\cal K}},\omega _{\cal K}^2).$$

We do not aim to estimate the precise values of ${r_{\cal K}}({\theta _i})$ during curve construction. However, the parameter ν, containing the values of ${\nu _{\cal K}}$ which determine the mean dcf/MRA offset for each dataset, is updated within our MCMC sampler.

To fully specify our offset model, values of the prior mean dcf/coastal shift, ${\rho _{\cal K}}$, for each dataset; and ${\omega _{\cal K}}$, our uncertainty on the level of this shift, were estimated from the overlap between additional data available for each speleothem/marine location and the dendrodated curve from 0–14 cal kBP. Similarly, these overlaps provided an estimate for each ${\zeta _{\cal K}}$. See Section 5.2 for more details. Finally, for the two floating Southern Hemisphere kauri trees we assumed an offset of $43 \pm 23$$^{\rm14}{\rm{C}}$ yrs ($1\sigma$) based upon the North-South hemispheric offset estimated in SHCal13 (Hogg et al. Reference Hogg, Hua, Blackwell, Niu, Buck, Guilderson, Heaton, Palmer, Reimer, Reimer, Turney and Zimmerman2013). Being direct measurements of the NH atmosphere, both Suigetsu and the three Bølling-Allerød floating tree-ring sequences are not offset, i.e. for all $\theta$, ${r_{\cal K}}(\theta ) = 0$.

Offsets: Cariaco Basin

The LSG model was not able to adequately resolve the MRA within the geographically unique Cariaco Basin which has a shallow sill that potentially limits exchange with the wider ocean. A further model for the MRAs in this location was therefore needed. As it covered a short time period, the MRA for the Cariaco varved record (Hughen et al. Reference Hughen, Southon, Bertrand, Frantz and Zermeño2004) was modeled as for the speleothems, i.e. independently varying around a constant level. For the Cariaco unvarved record (Hughen and Heaton Reference Hughen and Heaton2020 in this issue), we modeled the ${F^{14}}{\rm{C}}$ domain multiplier corresponding to any offset as

$${\gamma _{\cal K}}(\theta ) = {e^{ - {r_{\cal K}}(\theta )/8033}} = \sum\limits_{k = 1}^{30} {\beta _{{\rm{C}},k}}{B_{{\rm{C}},k}}(\theta ),$$

a further Bayesian spline with 30 knots placed at jittered quantiles of the Cariaco prior calendar ages. The value of ${{\bi {\beta}}_{\rm{C}}} = ({\beta _{{\rm{C}},1}}, \ldots ,{\beta _{{\rm{C}},30}}{)^T}$ determining this spline was also estimated during curve construction but with a fixed and large smoothing parameter ${\lambda _{\rm{C}}}$. This approach meant rapid changes in the ¹⁴C determinations within Cariaco were considered as atmospheric signal while smoother, longer term drifts away from the other data that might otherwise introduce spurious features into the curve, would be resolved as time-varying reservoir effects.

Heavier Tailed Errors

Due to the diversity of the datasets in this older time period, as well as the difficulty in incorporating all their potential complexities, we wished to create a curve that was not overly influenced by single determinations. We therefore permit the possibility of heavier tailed errors in the F ¹⁴C determinations for the data in this older section of curve extending beyond 13,913 cal BP. This is incorporated by introducing an error multiplier where for each observation we maintain

$${F_i} = {e^{ - {r_{\cal K}}({\theta _i})/8033}}f({\theta _i}) + {\varepsilon _i},$$

but with each determination’s reported uncertainty ${\sigma _i}$ scaled according to a further parameter ${\kappa _i}$

$${\varepsilon _i}|{\kappa _i} \sim N\left( {0,{{\sigma _i^2} \over {{\kappa _i}}}} \right),\qquad \qquad {\rm{where}}\ {\kappa _i} \sim Ga({\rm{shape}} = {{{\varrho _{{\cal K}(i)}}} \over 2},{\rm{rate = }}{{{\varrho _{{\cal K}(i)}}} \over 2}),$$

and ${\cal K}(i)$ is the dataset of observation $i$. All observations $i$ belonging to the same dataset ${\cal K}$ therefore have the same shape and rate in the Gamma prior for their individual ${\kappa _i}$. This model, integrating out ${\kappa _i}$, equates to a Student’s $t$-distribution for an individual observation i,

$${\varepsilon _i} \sim t({\rm{location}} = 0,{\rm{scale}} = {\sigma _i},{\rm{df}} = {\varrho _{{\cal K}(i)}}).$$

Further, for each dataset ${\cal K}$, we place an independent hyperprior on the value of its particular ${\varrho _{\cal K}}$ parameter so that the observations belonging to that dataset are grouped,

$${\varrho _{\cal K}} \sim Ga({\rm{shape}} = {A_\varrho },{\rm{rate}} = {B_\varrho }),$$

to capture potential differences between the various records. Consequently, if a dataset is seen to contain several outliers, the other determinations in that same set will also be treated with more caution. We choose a subjective hyperprior encapsulating an expectation of a low level of heavy-tailed behavior by selecting, for each dataset, ${A_\varrho } = 100$ and ${B_\varrho } = {1 \over 2}$. We update both ${\bi {\kappa}} = ({\kappa _1}, \ldots ,{\kappa _n}{)^T}$ and ${\bi{\varrho}}$, the values of the various ${\varrho _{\cal K}}$, within our sampler.

Parallel Tempering

A significant concern with the previous random walk approach was how well the MCMC sampler mixed. With this random walk approach, we only updated the calibration curve one calendar year at a time, conditional on the value at all other times, which restricted the ability of the curve estimate to move significantly. Furthermore, this update was performed by MH meaning we would frequently reject proposed updates. By switching to Bayesian splines, we are instead able to update, conditional on the current calendar ages, the entire curve through its spline coefficients simultaneously via Gibbs. These dual changes, to update the entire curve at once and to do so via Gibbs sampling as opposed to MH, hopefully begin to address mixing concerns. However, due to the uncertain calendar ages it is likely the posterior for the curve (and the true calendar ages) will remain multi-modal. In an attempt to overcome any remaining concerns over mixing, we also incorporate parallel tempering. In tempering, we run multiple chains concurrently. One of these MCMC chains has as its target our posterior of interest while the others have modified, higher temperature, targets. These higher temperature targets are typically flattened versions of the posterior of interest designed so that, as the temperature increases, the corresponding MCMC chain mixes more easily. We then propose swaps between the states of the multiple chains so that the chains that mix well (i.e. those that run at higher temperatures) improve the mixing of the chains which do not (i.e. our posterior of interest). We are free to choose the elements of the likelihood we wish to temper and by how much. We only temper the likelihood of the observed data, i.e. $[{\bf{F}}|{\bi{\theta}},{\bi{\beta}},{\bi{\nu}},{{\bi{\beta}}_{\rm{C}}},{\bi{\kappa}}]$; and $[{\bf{T}}|{\bi{\theta}}]$ (or equivalent prior $[{\bi{\theta}}]$) since these were thought to be the main elements restricting mixing and doing so keeps the MCMC updates straightforward. Our coupled chains run with a pair of temperatures ${\alpha _1},{\alpha _2} > 1$ with each chain sampling from a modified posterior,

$$\eqalign{{[{\bi{\theta}},{\bi{\beta}},\lambda ,{\bi{\nu}},{{\bi{\beta}}_{\rm{C}}}, {\bi{\kappa}},{\bi{{\varrho}}}|{\bf{F}},{\bf{T}}]_{{\alpha _1},{\alpha _2}}} \propto {L_{{\alpha _1},{\alpha _2}}}({\bi{\theta}},{\bi{\beta}},\lambda ,{\bi{\nu}},{{\bi{\beta}}_{\rm{C}}},{\bi{\kappa}},{\bi{{\varrho}}}|{\bf{F}},{\bf{T}})\cr = \ {[{\bf{T}}|{\bi{\theta}}]^{1/{\alpha _1}}}{[{\bf{F}}|{\bi{\theta}},{\bi{\beta}},{\bi{\nu}},{{\bi{\beta}}_{\rm{C}}},{\bi{\kappa}}]^{1/{\alpha _2}}}[{\bi{\beta}}|\lambda ][{\bi{\nu}}][{{\bi{\beta}}_{\rm{C}}}][\lambda ][{\bi{\kappa}}|{\bi{{\varrho}}}][{\bi{{\varrho}}}][{\bi{\theta}}]^{1/{\alpha _1}},}$$

so higher temperatures give flatter posteriors. We show in Section 4.4.3 how these modifications are simply equivalent to increasing our observational noise and so we can maintain straightforward and fast updates for each chain. We run the chain at four temperatures, including the unmodified posterior of actual interest where ${\alpha _1} = {\alpha _2} = 1$, in parallel. We discuss the effect of tempering, and how it aids mixing, in more detail within the Supplementary Information.

Merging with Already Created Tree-Ring-Only Section of Curve

We need to ensure we smoothly merge our older section of calibration curve with the independently created tree-ring-only curve described in Section 4.3. This is achieved by selection of an appropriate overlapping set of knots for the two sections. We first identify, in the tree-ring-only curve, the calendar age at which to create the join. This is determined to be the calendar age of the knot at which, due to the ending of the tree-ring determinations, the curve uncertainty on the tree-ring-only curve begins to increase significantly. For IntCal20, this was selected to be at 13,913 cal BP. This knot and the two knots either side that were used to create the tree-ring-only curve (five tree-ring-only knots in total) are then copied into the knot basis for the older portion of curve. On each update of the spline coefficients for the older section of curve, we pass the three spline coefficients that relate to the value of a particular realization of the tree-ring-only section of curve at the join. The Gibbs update for the older spline coefficients is then performed conditional on these three tree-ring-only spline coefficients. This ensures that when the posterior realizations of the two curves are connected together for the final curve, each realization is itself a cubic spline with not only continuity but also smoothness over the transition between the two sections.

Required Simplifications—Blocking and Over-Dispersion

To make the creation of the older section of the curve computationally practical, we do not incorporate exact modeling of blocking or over-dispersion. The large number of uncertain calendar ages mean that including blocking is not feasible. Each update to ${\bi{\theta}}$ requires recalculation of the spline basis matrix ${{\bi{{\cal B}}}_{\bi{\theta}}}$. Doing this at every year in ${{\bi{\theta}}^A}$, i.e. every year represented in all the blocks, would be extremely slow. Instead, we consider each determination to represent the single year in its centre. This significantly reduces the number of calendar ages ${\bi{\theta}}$ at which we have to recalculate our spline basis matrix. This is likely to have little effect on the overall curve since the density of data in this older section is not sufficient to identify very fine-scaled features. Over-dispersion was also not included in modeling the older curve as we lack the required data density to reliably estimate it. Further, our previous estimate of over-dispersion related to tree-ring determinations only, as opposed to the broader range of $^{\rm14}{\rm{C}}$ sources used within the older section of curve. Instead our modification to permit heavier tailed uncertainties aims to deal with potential over-dispersion in curve estimation. Note however, that when calibrating atmospheric ¹⁴C determinations against the curve in this older time period, we would still expect them to have similar potential additional sources of variation (beyond the laboratory quoted uncertainty) as seen in the predominantly dendrodated section. We therefore maintain predictive intervals for our final, older curve estimate by adding back in the over-dispersion estimated within the 0–14 cal kBP tree rings to the posterior curve estimate.

Posterior

Let us define, where subscripts denote the unknown (and updated) variables on which the elements depend:

$$\eqalign{ \hskip-4{\gamma _{\cal K}}(\theta ) = {e^{ - {r_{\cal K}}(\theta )/8033}}\quad \quad \quad\; {\Gamma _{{\bi{\nu}},{{\bi{\beta}}_{\rm{C}}},{\bi{\theta}}}} = {\rm{diag}}({\gamma _{{\cal K}(1)}}({\theta _1}), \ldots ,{\gamma _{{\cal K}(n)}}({\theta _n})), \cr C(\theta ) = {1 \over {1000}}{e^{ - \theta /8267}}\quad \quad \;\;{{\bi{{\cal C}}}_{\bi{\theta}}} = {\rm{diag}}(C({\theta _1}), \ldots ,C({\theta _n})), \cr E(\theta ) = {e^{ - \theta /8267}}\quad\quad \quad \quad \;\; {{\bf{e}}_{\bi{\theta}}} = (E({\theta _1}), \ldots ,E({\theta _n}{))^T}.}$$

Note that our atmospheric-offset multipliers ${\gamma _{\cal K}}(\theta )$ are determined by either ${\nu _{\cal K}}$, the values of the corresponding mean dcf or coastal MRA shift, or ${{\bi{\beta}}_{\rm{C}}}$ in the case of the Cariaco unvarved record. Then, considering each datum as representing a single year and letting ${\bf{F}} = ({F_1}, \ldots ,{F_n}{)^T}$, we can calculate for each dataset ${\cal K}$ the offset-adjusted calibration curve:

$${f_{\cal K}}(\theta ) = {\gamma _{\cal K}}(\theta )\left[ {C(\theta ){\bi{B}}{{(\theta )}^T}{\bi{\beta}} + E(\theta )} \right],\ {\rm{and\ so}}$$

$${\bf{F}}|{\bi{\theta}},{\bi{\beta}},{\bi{\nu}},{{\bi{\beta}}_{\rm{C}}},{\bi{\kappa}} \sim MVN\left( {{\Gamma _{{\bi{\nu}},{{\bi{\beta}}_{\rm{C}}},{\bi{\theta}}}}\left[ {{{\bi{{\cal C}}}_{\bi{\theta}}}{{\bi{{\cal B}}}_{\bi{\theta}}}{\bi{\beta}} + {{\bf{e}}_{\bi{\theta}}}} \right],{\rm{diag}}\left( {{{\sigma _i^2} \over {{\kappa _i}}}} \right)} \right).$$

Our joint posterior then becomes

$$\scale97%\eqalign{[{\bi{\theta}},{\bi{\beta}},{\lambda},{\bi{\nu}},{{\bi{\beta}}_{\rm{C}}},{\bi{\kappa}},{{\bi{\varrho}}}|{\bf{F}},{\bf{T}}] \propto [{\bf{T}}|{\bi{\theta}}][{\bf{F}}|{\bi{\theta}},{\bi{\beta}},{\bi{\nu}},{{\bi{\beta}}_{\rm{C}}},{\bi{\kappa}}][{\bi{\beta}}|\lambda ][{\bi{\nu}}][{{\bi{\beta}}_{\rm{C}}}][\lambda ][{\bi{\kappa}}|{\bi{\varrho}}][{\bi{\varrho}}][{\bi{\theta}}] \cr \propto \prod\limits_{\bi{\cal K}} \left( {\prod\limits_{j \in {\bi{\cal K}}} {{\sqrt {{\kappa _j}} } \over {{\sigma _j}}}} \right) \times \exp \left\{ { - {1 \over 2}[\left\{ {{{({\bf{T}} - {\bi{\theta}})}^T}{\bi{\Psi}}_T^{ - 1}({\bf{T}} - {\bi{\theta}})} \right\} + } \right. \cr \quad \left\{ {{{({\bf{F}} - {\Gamma _{{\bi{\nu}},{{\bi{\beta}}_{\rm{C}}},{\bi{\theta}}}}\left[ {{{\cal C}_{\bi{\theta}}}{{\bi{\cal B}}_{\bi{\theta}}}{\bi{\beta}} + {{\bf{e}}_{\bi{\theta}}}} \right])}^T}{{\bf{W}}_{\rm{k}}}({\bf{F}} - {\Gamma _{{\bi{\nu}},{{\bi{\beta}}_{\rm{C}}},{\bi{\theta}}}}\left[ {{{\cal C}_{\bi{\theta}}}{{\bi{\cal B}}_{\bi{\theta}}}{\bi{\beta}} + {{\bf{e}}_{\bi{\theta}}}} \right])} \right\} + \sum\limits_{\bi{\cal K}} {{{{({\nu _{\bi{\cal K}}} - {\rho _{\bi{\cal K}}})}^2}} \over {\omega _{\bi{\cal K}}^2}} \cr \quad + \lambda {{\bi{\beta}}^T}{\bf{D}}{\bi{\beta}} + {\lambda _{\rm{C}}}{\bi{\beta}}_{\rm{C}}^T{{\bf{D}}_{\rm{C}}}{{\bi{\beta}}_{\rm{C}}}]\} \times {\lambda ^{A + {\rm{rank}}({\bf{D}})/2 - 1}} \times \exp \left( { - {\lambda \over B}} \right) \cr \quad \times \prod\limits_{\bi{\cal K}} \left( {\prod\limits_{j \in {\bi{\cal K}}} {{{{{{{\varrho _{\bi{\cal K}}}} \over 2}}^{{{{\varrho _{\bi{\cal K}}}} \over 2}}}} \over {\Gamma ({{{\varrho _{\bi{\cal K}}}} \over 2})}}\kappa _j^{{{{\varrho _{\bi{\cal K}}}} \over 2} - 1}{e^{ - {{{\varrho _{\bi{\cal K}}}} \over 2}{\kappa _j}}}} \right) \times \prod\limits_{\bi{\cal K}} {{B_\varrho ^{{A_\varrho }}} \over {\Gamma ({A_\varrho })}}\varrho _{\bi{\cal K}}^{{A_\varrho } - 1}{e^{ - {B_\varrho }{\varrho _{\bi{\cal K}}}}},}$$

where ${\bf{W}}_{{\bi{\kappa}}} = {\rm{diag}}{\Big({{\sigma _1^2} \over {{\kappa _1}}}, \ldots ,{{\sigma _n^2} \over {{\kappa _n}}}\Big)^{ - 1}}$ includes the precision multiplier representing the heavy tails, and ${\bi{\Psi}}_T$ is the combined covariance matrix for all the calendar ages. As for the predominantly dendrodated section of curve, we can sample from this using a Metropolis-within-Gibbs algorithm, and again much of the updating can be performed directly via Gibbs sampling.

Gibbs updating the calibration curve ${\bi \beta}}|{\bf{F}},{\bi{\theta}},{\bi{\lambda}},{\bi{\nu}},{{\bi{\beta}}_{\rm{C}}},{\bi{\kappa}}$

This update is very similar to the equivalent update in the tree-ring-based section of curve. Letting ${{\bf{F}}^ \star } = {\bf{F}} - {\Gamma _{{\bi{\nu}},{{\bi{\beta}}_{\rm{C}}},{\bi{\theta}}}{{\bf{e}}_{\bi{\theta}}}}$ then, were it not for our need to create a smooth transition from the trees, the posterior

$${\bi{\beta}}|{{\bf{F}}^ \star },{\lambda} ,{\bi{\theta}},{\bi{\nu}},{\bi{\kappa}} \sim MVN\left( {{{\bi{\mu}}_{\bi{\beta}}} = {\bf{Q}}{\bi{{\cal B}}}_{\bi{\theta}}^T{\bi{{{\cal C}}}_{\bi{\theta}}}{\Gamma _{{\bi{\nu}},{{\bi{\beta}}_{\rm{C}}},{\bi{\theta}}}}}{{\bf{W}}_{{\bi{\kappa}}}{{\bf{F}}^ \star },{\bf{Q}}} \right),$$

where ${\bf{Q}} = ({\bi{{\cal B}}}_{\bi{\theta}}^T\left[ {{{\bi{{\cal C}}}_{\bi{\theta}}}{\Gamma _{{\bi{\nu}},{{\bi{\beta}}_{\rm{C}}},{\bi{\theta}}}}{{\bf{W}}_{{\bi{\kappa}}}}{\Gamma _{{\bi{\nu}},{{\bi{\beta}}_{\rm{C}}},{\bi{\theta}}}}{{\bi{{\cal C}}}_{\bi{\theta}}}} \right]{{\bi{{\cal B}}}_{\bi{\theta}}} + \lambda {\bf{D}}{)^{ - 1}}$ and, due to our dropping of blocking, all the matrices in the square brackets are diagonal so can be multiplied rapidly. To create the smooth transition, we draw a realization from the, already estimated, tree-ring-based section of curve and condition on ${{\bi{\beta}}^R}$, the value of the three basis coefficients affecting the value of the spline at the joining knot. Then

$${{\bi{\beta}}^{ - R}}|{{\bi{\beta}}^R} = {\bi{{a}}}} \sim MVN({{\bar {\bi{\mu}}}_{\bi{\beta}}},{\bi{{\bar Q}}}}),$$

where ${{\bar {{\bi{\mu}}}}_{\bi{\beta}}} = {{\bi{\mu}}}_{\bi{\beta}}^{( - R)} - {{\bf{Q}}_{ - R,R}}{\bf{Q}}_{R,R}^{ - 1}({\bi{{a}}}} - {{\bi{\mu}}}_{\bi{\beta}}^R)$, and ${\bi{{\bar Q}}}}$ is the Schur complement of ${{\bf{Q}}_{R,R}}$ in Q. Here, for example ${{\bf{Q}}_{R,R}}$ corresponds to the elements of the covariance that relate to the realized tree-ring based spline coefficients we condition on.

MH updating the calendar ages ${\bi{\theta}}|{\bf{T}},{\bf{F}},{\bi{\beta}},{\bi{\nu}},{{\bi{\beta}}_{\bf{C}}},{\bi{\kappa}}$

This step is as for the update of the wiggle matched trees. Whether a dataset has a prior on their calendar ages or noisy observations does not affect the update. For each dataset ${\cal K}$ in turn we:

• Propose ${{\bi{\theta}}_{\cal K}}' \sim MVN({{\bi{\theta}}_{\cal K}},{\Sigma_{{\rm{prop}}}})$;

• Calculate ${{\bi{{\cal B}}}_{\bi{\theta}}}$ and ${\bi{{{\cal B}}}_{{\bi{\theta'}}}}$ (and similarly recalculate ${\Gamma _{{\bi{\nu}},{{\bi{\beta}}_{\bf{C}}},{\bi{\theta}}'}}$, ${C_{\bi{{\theta'}}}}$ and ${E_{{\bi{\theta'}}}}$ at ${{\bi{\theta}}_{\cal K}}'$). Accept according to Hastings Ratio:

$$HR = \min \left\{ {1,{{\pi ({{\bi{\theta}}_{\cal K}}'|{{\bf{T}}_{\cal K}},{\bf{F}},{\bi{\beta}},{\bi{\nu}},{\bi{\kappa}})} \over {\pi ({{\bi{\theta}}_{\cal K}}|{{\bf{T}}_{\cal K}},{\bf{F}},{\bi{\beta}},{\bi{\nu}},{\bi{\kappa}})}}} \right\},$$

where

$$\eqalign{\pi ({{\bi{\theta}}_{\cal K}}|{{\bf{T}}_{\cal K}},{\bf{F}},{\bi{\beta}},{\bi{\nu}},{\bi{\kappa}}) = \exp \left\{ { - {1 \over 2}\bigg[\left\{ {{{({{\bf{T}}_{\cal K}} - {{\bi{\theta}}_{\cal K}})}^T}{\bi{\Psi}}_{{\cal K},T}^{ - 1}({{\bf{T}}_{\cal K}} - {{\bi{\theta}}_{\cal K}})} \right\} + } \right. \cr \left\{ {{{({\bf{F}} - {\Gamma _{{\bi{\nu}},{{\bi{\beta}}_{C}},{\bi{\theta}}}}\left[ {{{\cal C}_{\bi{\theta}}}{{\cal B}_{\bi{\theta}}}\bi{\beta} + {{\bf{e}}_{\bi{\theta}}}} \right])}^T}{{\bf{W}}_{{\bi{\kappa}}}}({\bf{F}} - {\Gamma _{{\bi{\nu}},{{\bi{\beta}}_{\rm{C}}},{\bi{\theta}}}}\left[ {{{\cal C}_{\bi{\theta}}}{{\cal B}_{\bi{\theta}}}{\bi{\beta}} + {{\bf{e}}_{\bi{\theta}}}} \right])} \right\}\bigg]\bigg\} .}$$

Note that we only need to recalculate the elements in the likelihood of ${\bf{F}}$ (e.g. the elements of ${{\cal B}_{\bi{\theta}}}$) that relate to ${{\bi{\theta}}_{\cal K}}$. In the case of independent error, i.e. ${T_j} \sim N({\theta _j},\sigma _j^2)$, we use a proposal standard deviation ${\sigma _{{\rm{prop}}}}=2{\sigma _j}$ and we can accept/reject each proposal $j=1, \ldots ,{n_{\cal K}}$ independently. In the case of dependent error with a covariance matrix, i.e. ${{\bf{T}}_{\cal K}} \sim MVN({{\bi{\theta}}_{\cal K}},{{\bi{\Psi}}_{{\cal K},T}})$, or an equivalent prior, we update the whole of ${{\bi{\theta}}_{\cal K}}$ simultaneously by sampling a proposal ${\bi{\theta}}' \sim MVN({\bi{\theta}},{\Sigma}_{{\rm{prop}}})$, where ${\Sigma_{{\rm{prop}}}} \propto {{\bi{\Psi}}_{{\cal K},T}}$, to obtain higher acceptance rates.

MH updating the offsets ${{\bi \nu _{\boldcal K}}}|{\bf{F}},{\bi{\theta}},{\bi{\beta}},{\bi{\lambda}} ,{\bi{\xi}},{\bi{\tau}}$

Similarly, for each offset dataset:

• Propose a new constant dcf/MRA shift ${\nu _{\cal K}}' \sim N({\nu _{\cal K}},{\sigma _{{\rm{prop}}}})$;

• Calculate resultant ${F^{14}}{\rm{C}}$ multiplier ${\Gamma _{{\bi{\nu}}',{{\bi{\beta}}_{\rm{C}}},{\bi{\theta}}}}$. Accept according to Hastings Ratio:

$$HR = \min \left\{ {1,{{\pi ({\nu _{\cal K}}'|{\bf{F}},{\bi{\theta}},{\bi{\beta}},{{\bi{\nu}}^{ - {\cal K}}},{{\bi{\kappa}}})\varphi ({\nu _{\cal K}}';{\rho _{\cal K}},\omega _{\cal K}^2)} \over {\pi ({\nu _{\cal K}}|{\bf{F}},{\bi{\theta}},{\bi{\beta}},{{\bi{\nu}}^{ - {\cal K}}},{\bi{\kappa}})\varphi ({\nu _{\cal K}};{\rho _{\cal K}},\omega _{\cal K}^2)}}} \right\},$$

where

$$\pi ({\nu _{\cal K}}|{\bf{F}},{\bi{\theta}},{\bi{\beta}},{{\bi{\nu}}^{ - {\cal K}}},{\bi{\kappa}}) = \exp \left\{ { - {1 \over 2}[({\bf{F}} - {\Gamma _{{\bi{\nu}},{{\bi{\beta}}_{\rm{C}}},{\bi{\theta}}}}\left[ {{{\cal C}_{\bi{\theta}}}{{\cal B}_{\bi{\theta}}}{\bi{\beta}} + {{\bf{e}}_{\bi{\theta}}}} \right]{)^T}{{\bf{W}}_{{\bi{\kappa}}}}({\bf{F}} - {\Gamma _{{\bi{\nu}},{{\bi{\beta}}_{\rm{C}}},{\bi{\theta}}}}\left[ {{{\cal C}_{\bi{\theta}}}{{\cal B}_{\bi{\theta}}}{\bi{\beta}} + {{\bf{e}}_{\bi{\theta}}}} \right])]} \right\}$$

and $\varphi ( \cdot ;{\rho _{\cal K}},\omega _{\cal K}^2)$ is a normal density with mean ${\rho _{\cal K}}$ and variance $\omega _{\cal K}^2$.

As when updating the calendar ages for each dataset in turn, when calculating this Hastings Ratio we can restrict our attention to the dataset in question, as elements relating to other datasets will cancel, making this step considerably faster.

Gibbs updating the Cariaco Basin unvarved MRA ${{\bi{\beta}_C}}|{\bf{F}},{\bi{\theta}},{\bi{\beta}},{\bi{\lambda}},{\bi{\nu}},{\bi{\kappa}}$

Updating the spline MRA coefficients ${{\bi{\beta}}_C}$ for the Cariaco Basin is analogous to the update of the calibration curve coefficients ${\bi{\beta}}$. This can intuitively be seen by considering how these determinations are represented and swapping the offset multiplier ${\gamma _{\cal K}}({\theta _{\cal K}})$ with the calibration curve estimate $f(\theta )$, i.e. for determinations ${F_i}$ belonging to the Cariaco Basin unvarved record,

$${F_i} = {\gamma _{{\cal K}(i)}}({\theta _i})f({\theta _i}) + {\varepsilon _i} = f({\theta _i}){\gamma _{{\cal K}(i)}}({\theta _i}) + {\varepsilon _i},$$

and noting that, given ${\bi{\beta}}$ and ${\bi{\theta}}$, the curve estimates $f({\theta _i})$ are known and ${\gamma _{\cal K}}({\theta _i}) = \sum\nolimits_{k=1}^{30} {\beta _{{\rm{C}},k}}{B_{{\rm{C}},k}}({\theta _i})$, a linear combination of the MRA spline coefficients ${{\bi{\beta}}_C}$. More rigorously, this can be seen if we restrict our likelihood to solely the Cariaco determinations, and change the ordering of ${\Gamma _{{\bi{\nu}},{{\bi{\beta}}_{\rm{C}}},{\bi{\theta}}}}$ and $\left[ {{{\bi{{\cal C}}}_{\bi{\theta}}}{{\bi{{\cal B}}}_{\bi{\theta}}}{\bi{\beta}} + {{\bf{e}}_{\bi{\theta}}}} \right]$. Then, letting ${{\cal F}_{\cal K}} = {\rm{diag}}({{{\bi{\cal C}}}_{{{\bi{\theta}}_{\cal K}}}}{{{\bi{\cal B}}}_{{\bi{{\theta}}_{\cal K}}}}{\bi{\beta}} + {{\bf{e}}_{{{\bi{\theta}}_{\cal K}}}})$ be a diagonal matrix representing the current atmospheric calibration curve at the Cariaco calendar ages, our posterior for the MRA spline ${{\bi{\beta}}_C}$ becomes

$$[{{\bi{\beta}}_C}|{F_{\cal K}}][{{\bi{\beta}}_C}] \propto \exp \left\{ { - {1 \over 2}\big[({{\bf{F}}_{\cal K}} - {{\cal F}_{\cal K}}{{\cal B}_{C,{{\bi{\theta}}_{\cal K}}}}{{\bi{\beta}}_C}{)^T}{W_{{\bi{\kappa}}}}({{\bf{F}}_{\cal K}} - {{\cal F}_{\cal K}}{{\cal B}_{C,{{\bi{\theta}}_{\cal K}}}}{{\bi{\beta}}_C}) + {\lambda _{\rm{C}}}{\bi{\beta}}_{\rm{C}}^T{{\bf{D}}_{\rm{C}}}{{\bi{\beta}}_{\rm{C}}}\big]} \right\}.$$

Gibbs updating the smoothing parameter ${\bi{\lambda}} |{\bi{\beta}}$

As before, $\lambda |{\bi{\beta}} \sim Ga(A',B')$ where $A' = A + {\rm{rank}}({\bf{D}})/2$ and $B' = {\left( {{{{{\bi{\beta}}^T}{\bf{D}}{\bi{\beta}}} \over 2} + {1 \over B}} \right)^{ - 1}}$.

Gibbs updating the precision multipliers ${{\bi{\kappa}}}|{\bf{F}},{\bi{\beta}},{\bi{\theta}},{{\bi{\nu}}},{{\bi{\beta}}_{\rm{C}}},{\bi{\varrho}}$

For each datum,

$${\kappa _i}|{F_i},{\bi{\beta}},{\theta _i},{\gamma _{{\cal K}(i)}} \sim Ga({{{\varrho _{{\cal K}(i)}}} \over 2} + {1 \over 2},{{{\varrho _{{\cal K}(i)}}} \over 2} + {{R_i^2} \over {2\sigma _i^2}}),$$

where ${R_i} = {F_i} - {\gamma _{{\cal K}(i)}}({\theta _i})C({\theta _i})B{({\theta _i})^T}{\bi{\beta}} - E({\theta _i})$ is the current residual, and ${{{\varrho _{{\cal K}(i)}}} \over 2}$ is the current shape/rate prior for the dataset ${\cal K}$ to which datum $i$ belongs.

MH updating the expected tail behavior ${{\bi{\varrho}}_{\boldcal K}}|{\bi{\kappa}}$

For each dataset ${\cal K}$:

• Propose new parameter ${\varrho _{\cal K}}' \sim N({\varrho _{\cal K}},{\sigma _{{\rm{prop}}}})$;

• Accept according to Hastings Ratio:

$$HR = \min \left\{ {1,{{\pi ({\varrho _{\cal K}}'|{\bi{\kappa}})\pi ({\varrho _{\cal K}}')} \over {\pi ({\varrho _{\cal K}}|{\bi{\kappa}})\pi ({\varrho _{\cal K}})}}} \right\},$$

where

$$\pi ({\varrho _{\cal K}}|{\bi{\kappa}}) = \prod\limits_{j \in {\cal K}} {{{{{{{\varrho _{\cal K}}} \over 2}}^{{{{\varrho _{\cal K}}} \over 2}}}} \over {\Gamma ({{{\varrho _{\cal K}}} \over 2})}}\kappa _j^{{{{\varrho _{\cal K}}} \over 2} - 1}{e^{ - {{{\varrho _{\cal K}}} \over 2}{\kappa _j}}}\quad {\rm{and}}\ \pi ({\varrho _{\cal K}}) = {{B_\varrho ^{{A_\varrho }}} \over {\Gamma ({A_\varrho })}}\varrho _{\cal K}^{{A_\varrho } - 1}{e^{ - {B_\varrho }{\varrho _{\cal K}}}}.$$

Choice of Splines

The approach to creating the older part of the calibration curve is considerably more computationally intensive than for the predominantly dendrodated section of curve. In this older section, we have to update the calendar ages of the constituent determinations requiring recalculation of the value of the curve each time. This becomes less feasible as the number of knots increases. Consequently, for this older section, we are not able to choose such a large number of knots for our spline basis as used in the dendrodated section. To select our basis for the older section of curve, we therefore began by placing 400 knots at quantiles (with the addition of a small amount of random jitter) of the observed/prior calendar ages ${T_i}$. This meant that a knot was placed approximately every $5$th observed calendar age, similar to the dendrodated section of curve. However, in order to make the fitted placement of the three Bølling-Allerød floating tree-ring sequences more equitable, we overrode this knot choice in the period from 14.2–15.2 cal kBP. In this period, 50 knots were placed evenly (i.e. one every 20 calendar years) and all other knots were removed. Tree-ring sequences P305u and P317 in particular show considerable variability in $^{\rm14}{\rm{C}}$ over their course; see Adolphi et al. (Reference Adolphi, Muscheler, Friedrich, Güttler, Wacker, Talamo and Kromer2017) and Figure 7. If the knots used to model the calibration curve over the potential fitted range of these trees lie unevenly in calendar age, as would be the case if we maintained a knot placement based upon the observed quantiles of the other IntCal20 data, the highly variable $^{\rm14}{\rm{C}}$ nature of these tree-ring sequences may have meant that the method had a preference/bias to locate them where there were more knots as opposed to where they best fitted the other data. Conversely placing knots at even calendar age spacings aims to remove this potential bias. One knot every 20 calendar years is still sufficient to pick out required curve detail in this period. Finally, we added into the spline bases the 5 knots required for merging with the tree-ring section of curve and removed other knots in this period. In total, after these modifications, we were left with 444 knots meaning that the overall method still maintained computational feasibility. These knot locations then remained fixed throughout the sampler, i.e. did not change as the calendar ages ${\theta _i}$ were updated.

Updating with Parallel Tempering

Tempering our target likelihood for the calendar age terms $[{\bf{T}}|{\bi{\theta}}]$ (or alternatively the prior $[{\bi{\theta}}]$) concerning a particular dataset ${\cal K}$, we choose

$${[{{\bf{T}}_{\cal K}}|{{\bi{\theta}}_{\cal K}}]^{1/{\alpha _1}}} \propto \exp \left\{ { - {1 \over {2{\alpha _1}}}\Big[\left\{ {{{({{\bf{T}}_{\cal K}} - {{\bi{\theta}}_{\cal K}})}^T}{\bi{\Psi}}_{{\cal K},T}^{ - 1}({{\bf{T}}_{\cal K}} - {{\bi{\theta}}_{\cal K}})} \right\}\Big]} \right\}$$

This is equivalent to the untempered version but with an adjusted covariance ${{\bi{\Psi}}_{{\cal K},T}}({\alpha _1}) = {\alpha _1}{{\bi{\Psi}}_{{\cal K},T}}$. To run the chain at this temperature we simply adjust the calendar age covariance (for observed ${\bf{T}}$ or prior) accordingly. We can then perform all updates as described in Section 4.4.2 but with this modified calendar age covariance.

In tempering the likelihood for the observed ${{\bf{F}}_{\cal K}}$ for a particular dataset ${\cal K}$, we select

$$\eqalign{ {[{{\bf{F}}_{\cal K}}|{{\bi{\theta}}_{{\bi{\cal K}}}},{\bi{\beta}},{\bi{\nu}},{{\bi{\beta}}_{\rm{C}}},{\bi{\kappa}}]^{1/{\alpha _2}}} \propto \left( {\prod\limits_{j \in {\cal K}} {{\sqrt {{\kappa _j}} } \over {\sqrt {{\alpha _2}\ \sigma _j^2} }}} \right)\cr \quad \quad \times \exp \left\{ { - {1 \over {2{\alpha _2}}}{{({{\bf{F}}_{\bi{{\cal K}}}} - {\Gamma _{{\nu _{\bi{{\cal K}}}},{{\bi{\theta}}_{\bi{{\cal K}}}}}}\left[ {{{\bi{{\cal C}}}_{{{\bi{\theta}}_{\bi{{\cal K}}}}}}{{\cal B}_{{{\bi{\theta}}_{\bi{{\cal K}}}}}}{\bi{\beta}} + {{\bf{e}}_{{{\bi{\theta}}_{\bi{{\cal K}}}}}}} \right])}^T}{{\bf{W}}_{{{\bi{\kappa}}_{\bi{{\cal K}}}}}}({{\bf{F}}_{\bi{{\cal K}}}} - {\Gamma _{{\nu _{\bi{{\cal K}}}},{{\bi{\theta}}_{\bi{{\cal K}}}}}}\left[ {{{\bi{{\cal C}}}_{{{\bi{\theta}}_{\bi{{\cal K}}}}}}{{\cal B}_{{{\bi{\theta}}_{\bi{{\cal K}}}}}}{\bi{\beta}} + {{\bf{e}}_{{{\bi{\theta}}_{\bi{{\cal K}}}}}}} \right])} \right\},}$$

where for example, in a slight abuse of notation, ${\Gamma _{{\nu _{\cal K}},{{\bi{\theta}}_{\cal K}}}}$ denotes the diagonal offset multiplier matrix for the specific dataset ${\cal K}$, and will depend upon only ${\nu _{\cal K}}$. In the case of the Cariaco Basin unvarved record, we would instead have ${\Gamma _{{{\bi{\beta}}_{C}},{{\bi{\theta}}_{\cal K}}}}$. This is also equivalent to the untempered version but modifying all the ${F^{14}}{\rm{C}}$ variances within dataset ${\cal K}$ to be $\sigma _j^2({\alpha _2}) = {\alpha _2}\sigma _j^2$ for all $j \in {\cal K}$. To perform MCMC on the chain with this temperature, we make this modification and then all update steps remain as in Section 4.4.2 but with the temperature-adjusted F ¹⁴C variances.

The amount of tempering applied is dependent upon the datasets. We apply most tempering to the two floating Bølling-Allerød tree-ring sequences P305u and P317 since their internal $^{\rm14}{\rm{C}}$ determinations are highly variable and we wish to fully explore their potential fitting location; and also the Cariaco Basin unvarved record which has a large number of observations with high calendar age covariances where again we might otherwise be concerned about mixing. For these three datasets, we chose values of ${\alpha _2}$ that were equivalent, at the highest temperature of the four parallel chains, to increasing the uncertainty ${\sigma _i}$ on the ${F^{14}}{\rm{C}}$ determinations by $ \approx 36\%$ from the quoted laboratory values. For the other data, the highest ${\alpha _2}$ temperature was equivalent to increasing ${\sigma _i}$ by $ \approx 5\%$. We performed the same, small amount of tempering for the calendar ages as otherwise the acceptance rate for swapping between chains became very low. For all the data, the highest ${\alpha _1}$ temperature of the four chains was equivalent to increasing the calendar age covariance by $ \approx 10\%$.

Run Time and Convergence

We ran the four tempered chains in parallel, each for 250,000 iterations. Four independent swaps between the tempered chains were proposed every 5th iteration. Convergence was assessed by initializing the chains at several different starting points and then visually observing the robustness in the curve estimates; along with assessment of other outputs such as the posterior estimates for the calendar ages of the floating tree-ring sequences and the spline based MRA of the Cariaco Basin; and the overall model likelihood. For a more thorough discussion of convergence and the effect of tempering, see the Supplementary Information. The level of tempering chosen was based upon a trade-off between maintaining a reasonable acceptance rate of chain swaps, and a highest temperature that both allowed the Bølling-Allerød trees P305u and P317 to move and gave a “highest-temperature” curve estimate that had removed most fine detail; again see the Supplementary Information. In total, there were 4407 accepted chain swaps, with 559 between the two lowest temperatures, i.e. swaps that involved the untempered target of true interest. To create the final curve, we removed the first half of the target chain as burn-in and then thinned to every 50th realization leaving us with 2500 posterior realizations covering the older time-period from which to create the final complete curve.

Splitting the Hoffman Speleothem

The Hoffman speleothem (Hoffmann et al. Reference Hoffmann, Beck, Richards, Smart, Singarayer, Ketchmark and Hawkesworth2010) consists of two separated sections—a set of relatively young data for which there was some overlap with the trees, and then a set of much older measurements. These were split in curve creation and a separate dcf was applied to the older measurements. For the more recent section, the prior for the mean dcf ${\nu _{\cal K}}$ was based on the overlapping determinations, while an uninformative prior was placed on the mean dcf for the older measurements. Again this keeps the internal shape/structure but allows the dataset to jointly move up/down to fit best with the other material.