USING SEDIMENTOLOGICAL PRIORS TO IMPROVE ¹⁴C CALIBRATION OF BIOTURBATED SEDIMENT ARCHIVES

Abstract

Radiocarbon (¹⁴C) dating is often carried out upon multi-specimen samples sourced from bioturbated sediment archives, such as deep-sea sediment. These samples are inherently heterogeneous in age, but existing ¹⁴C calibration techniques were originally developed for age homogeneous material, such as archaeological artifacts or individual tree rings. A lack of information about age heterogeneity leads to a systematic underestimation of a sample’s true age range, as well as the possible generation of significant age-depth artifacts during periods of the Earth’s history coinciding with highly dynamic atmospheric Δ¹⁴C. Here, a new calibration protocol is described that allows for the application of sedimentological priors describing sediment accumulation rate, bioturbation depth and temporally dynamic species abundance. This Bayesian approach produces a credible calibrated age distribution associated with a particular laboratory ¹⁴C determination and its associated sedimentological priors, resulting in an improved calibration, especially in the case of low sediment accumulation rates typical of deep-sea sediment. A time-optimized computer script (biocal) for the new calibration protocol is also presented, thus allowing for rapid and automated application of the new calibration protocol. This new calibration protocol could be applied within existing age-depth modeling software packages to produce more accurate geochronologies for bioturbated sediment archives.

INTRODUCTION

Radiocarbon (¹⁴C) analysis is routinely used to determine the age of marine sediment archives up to ∼50 ka in age, and has been fundamental in increasing our understanding of the spatio-temporal development of palaeoclimate during the last glacial and the Holocene. However, due to ¹⁴C being a very rare radioisotope in the environment (approximately one in 10¹² carbon dioxide molecules in the atmosphere is ¹⁴CO₂), it is more difficult to measure than more common, stable carbon isotopes. From a practical standpoint, this rarity results in a requirement of relatively large sample sizes to attain a sufficient measurement signal using, e.g., accelerated mass spectrometry (AMS). In the case of, e.g., deep-sea sediment archives, many tens of single microfossil specimens from a discrete core depth are often pooled into a single sample for measurement. However, AMS will only report a mean ¹⁴C activity (along with a measurement uncertainty), so any information about the ¹⁴C heterogeneity of the sample will be lost.

Systematic bioturbation of deep-sea sediment causes discrete downcore intervals of deep-sea sediment to have an age distribution that is characterized by an exponential probability density function with a long tail towards older ages (Berger and Heath 1968). This age distribution is mainly governed by the sediment accumulation rate (SAR) and bioturbation depth (BD), the latter of which is typically around 10 cm (Trauth et al. 1997; Boudreau 1998). The presence of the aforementioned age distribution is supported by studies of, e.g., particle mixing, stable isotopes, ¹⁴C, species abundance and tephras (Bramlette and Bradley 1942; Nayudu 1964; Ruddiman and Glover 1972; Peng et al. 1979; Hutson 1980; Pisias 1983; Schiffelbein 1984; Andree 1987; Bard et al. 1987; Wheatcroft 1992; Trauth et al. 1997; Henderiks et al. 2002; Löwemark and Grootes 2004; Sepulcre et al. 2017; Lougheed et al. 2018; Abbott et al. 2018; Missiaen et al. 2020; Dolman et al. 2021).

In the case of a wholly non-bioturbated sediment archive (such as laminated sediment retrieved from an anoxic environment), a SAR of 5 cm ka^–1 can be assumed to have a temporal resolution of 1000/5 = 200 yr cm^–1. However, in the case of bioturbated sediment typical of the oxygenated deep-sea, the 1σ age value of a 1 cm slice of sediment with a typical deep-sea SAR of 5 cm ka^–1 and BD of 10 cm can be approximated as 10/5×1000 = 2000 yr (Berger and Heath 1968). Somewhat counter-intuitively, that same bioturbated sediment archive will also exhibit a downcore increasing mean age of 200 yr cm^–1, which can deceptively mask the fact that the sediment is bioturbated. In essence, downcore increase in mean age is not the same concept as the discrete-depth age variance. This concept is visualized in Figure 1 by a 5 cm ka^–1 single particle sediment simulation, where it can be seen that the discrete-depth median age of the single particle population increases by ∼200 yr per cm, whereas the actual age range contained in each discrete 1 cm depth is actually much greater, and characterized by an exponential distribution with a long tail towards older ages (which can be seen in Figure 1 as the decreasing density of single particles towards older ages). Hence, the true temporal resolution of the 5 cm ka^–1 discrete-depth archive is much greater than 200 yr, and failing to consider this point risks leading to false high-precision in age-depth chronologies.

Figure 1 5 cm ka^–1 sediment simulation of single particles using a global average BD of 10 cm (Trauth et al. 1997; Boudreau 1998) and a best-case 10⁴ simulated particles per cm. Shown also is the discrete 1 cm depth median age, as well as the associated 68.2% and 95.4% age range. Also shown are the calibrated age distributions that would result if one were to use the existing state of the art to calibrate the mean ¹⁴C activity resulting from all the particles contained in each 1 cm discrete depth. The single particle simulation is carried out using SEAMUS (Lougheed 2020), using the IntCal20 calibration curve (Reimer et al. 2020) and assuming no reservoir affect. Calibration is carried out using MatCal (Lougheed and Obrochta 2016).

The current state of the art in palaeoclimate includes no information about bioturbation when ¹⁴C calibrating multi-specimen samples retrieved from deep-sea (or lacustrine) sediment archives. In essence, the current state of the art considers only ¹⁴C-centric priors (calibration curve and reservoir effect) when estimating (calibrating) true age. Such an approach incorrectly treats deep-sea sediment as having discrete age increments, similar to non-bioturbated archives such as tree rings, speleothems and/or varves. By ignoring bioturbation when calibrating ¹⁴C measurements from sediment archives, one essentially assumes that a BD of 0 cm and SAR of 1000 cm ka^–1. Hence, the current lack of method for including correct sedimentological priors when applying the ¹⁴C method to sediment archives can lead to an underestimation of the full age uncertainty. This underestimation is also illustrated in Figure 1, where the existing state of the art in ¹⁴C dating and calibration is virtually applied to a simulated bioturbated sediment core created using the SEAMUS single foraminifera simulator (Lougheed 2020). In can further be seen in Figure 1 that the current calibration method can produce significant age-depth artifacts when applied to bioturbation, which is due to the mixing of single elements (e.g., foraminifera) from periods of past dynamic Δ¹⁴C into the same discrete depth interval (Lougheed et al. 2020).

METHOD

The new calibration protocol presented here involves complementing the traditional ¹⁴C priors (past Δ¹⁴C from a calibration curve, reservoir effect) with sedimentological priors (SAR, BD and temporal changes in species abundance). This improved calibration protocol for sediment archives allows us to estimate an improved age distribution from the ¹⁴C activity measurement carried out on a given bioturbated sample (Figure 2).

Figure 2 A flow chart demonstrating the principle of complementing ¹⁴C priors with sedimentological priors to produce a more accurate calibrated age distribution for bioturbated sediment archives.

Establishing a Prior Distribution for Calendar Age

In order to calibrate ¹⁴C activity measurements carried out upon heterogeneous samples retrieved from bioturbated sediment, the following sedimentological priors are defined:

• s = estimated sediment accumulation rate (SAR), in cm yr^–1

• m = bioturbation (mixing) depth (BD), in cm

• k = the fraction of the analyzed microfossils that are fragmented (a value between 0 and 1)

• a = time series of abundance of the analyzed species relative to itself (values between 0 and 1)

Both SAR and BD are considered here as a constant value, i.e., not as a time series of temporally variable values. These inputs are kept constant foremost to reduce computation time, and also because temporal changes in, e.g., SAR (the relationship between mean age and depth) are not known when an age-depth chronology has yet to be developed. In short, applying detailed information about temporal changes in SAR when the age-depth relationship of the sediment is not yet known would constitute circular thinking.

Prior information is often applied within Bayesian analysis to construct an expected prior probability distribution based on established understanding of physical processes. In this case, we use SAR and BD priors to construct a prior distribution of relative age for the sample being calibrated, based on theoretical understanding of the influence of bioturbation upon the age distribution of sediment. Following Berger and Heath (1968), the age distribution for a given depth of fully bioturbated sediment core can be represented by an exponential probability distribution, which can be considered the basis of the prior probability distribution for a sample’s calibrated age:

(1) $${p_{prior}}\left( {{r_1},{r_2}, \ldots ,{r_n}} \right) = exp\left( {{{ - \left( {{r_1},{r_2}, \ldots ,{r_n}} \right)s} \over m}} \right)$$

where r is the relative age (starting at 1 yr) within P _prior . The low-probability long tail of an exponential probability function continues to infinity, which obviously cannot be stored in computer memory. The prior distribution is therefore limited to the age equivalent value of five bioturbation depths, i.e., a relative age of r _limit = 5m/s, which is rounded to the nearest whole year.

When picking microfossils for ¹⁴C analysis, palaeoceanographers generally prefer to pick whole and/or pristine specimens. The fragmented and/or dissolved microfossils that are not picked have been resident in the bioturbation depth for a longer time and have been exposed to more bioturbation cycles, and as such represent the oldest fraction of the sample (Rubin and Suess 1955; Ericson et al. 1956; Emiliani and Milliman 1966; Barker et al. 2007). There is therefore a benefit in not picking the older, broken foraminifera, as it results in a more constrained age distribution (the long tail of the age distribution is shortened). Information regarding the fact that the oldest/broken foraminifera are not picked can be incorporated into the prior distribution. The estimated fraction of fragmented microfossils (k) can be related to the cumulative expression of Eq. (1):

(2) $$1 - k = 1 - {\rm{exp}}\left( {{{ - rs} \over m}} \right)$$

Eq. (2) can be solved to attain r(k), the threshold age for fragmented foraminifera:

(3) $$r\left( k \right) = {{ - m.{\rm{ln}}\left( k \right)} \over s}$$

Regions of the prior probability distribution (p _prior ) older than r(k) can, therefore, be considered to consist of fragmented microfossils that are not picked by palaeoceanographers. When r(k) < r _limit , p _prior is truncated at the discrete relative age r(k) to incorporate prior information from the picking process. When r(k) ≥ r _limit , r(k) is approximated to r _limit . All discrete probability values in p_prior are subsequently normalized such that they sum to 1.

Establishing a Distribution for ¹⁴C Activity

Please note that, to avoid ambiguity, throughout this text the use of the term “age” refers exclusively to true/calibrated age, while ¹⁴C activity is always referred to as ¹⁴C activity, i.e., not as “¹⁴C age”.

The new calibration protocol must incorporate the full uncertainty regarding ¹⁴C activity, which includes uncertainties regarding the laboratory ¹⁴C activity determination, the calibration curve ¹⁴C activity, and the ¹⁴C activity depletion as a result of the reservoir effect. These are expressed here as follows:

• A _det = The laboratory ¹⁴C activity determination of the sample (in ¹⁴C yr BP).

• σ _det = The measurement uncertainty associated with A _det (in ¹⁴C yr).

• A _cc (t) = The ¹⁴C activity (in ¹⁴C yr BP) predicted by the calibration curve for a discrete age t.

• σ _cc (t) = The uncertainty (in ¹⁴C yr) associated with A _cc (t).

• R(t) = The predicted ¹⁴C activity depletion (in ¹⁴C yr) of A_det relative to the calibration curve at discrete age t, due to a local reservoir effect (Stuiver et al. 1986). R(t) can be substituted with ΔR(t) in the case of a marine calibration curve.

• σ _R (t) = The uncertainty (in ¹⁴C yr) associated with R(t) (or ΔR(t)).

Activity depletion due to R(t) is considered here by incorporating it into the calibration curve ¹⁴C activity. This approach to handling R(t) allows, if desired, for temporally dynamic R(t) to be correctly incorporated (Waelbroeck et al. 2019). The calibration curve is adjusted as follows, for each discrete calendar age t:

(4) $${A_{ccR}}\left( t \right) = {A_{cc}}\left( t \right) + R\left( t \right)$$

Uncertainties pertaining to calibration curve ¹⁴C activity and the ¹⁴C reservoir effect (σ _cc (t) and σ _R (t)) are both Gaussian, so they can be easily propagated into one term, for each discrete calendar age t:

(5) $${\sigma _{ccR}}\left( t \right) = \sqrt {\left( {\sigma _{cc}^2\left( t \right) + \sigma _R^2\left( t \right)} \right)} $$

Before proceeding, all of the aforementioned ¹⁴C-related values are first converted into F¹⁴C space to facilitate more accurate calculations that take isotope mass balance into account, which is especially relevant in the case of wide range of ¹⁴C activity (Erlenkeuser 1980; Bronk Ramsey 2008; Keigwin and Guilderson 2009), such as is the case with bioturbated sediment archives.

A sequence of probabilities can describe the closeness of a sequence of ¹⁴C activities predicted for all discrete ages t (represented as T) available within the calibration curve (i.e., A _ccR (T)), to a single ¹⁴C activity predicted by the calibration curve for a discrete age t (i.e., A _ccR (t)). This closeness, which includes a quantification of calibration curve and reservoir effect uncertainties, can be evaluated using a normal distribution for each instance of t, summing through all n values available in T to give the total relative ¹⁴C probability for each t:

(6) $${p_{14C}}\left( {T{\rm{|}}t} \right) = \sum\nolimits_{{T_1}}^{{T_n}} {\left( {{1 \over {{\sigma _{ccR}}\left( t \right)\sqrt {\left( {2\pi } \right)} }}{\rm{exp}}\left( {{{ - {{\left( {{A_{ccR}}\left( T \right) - {A_{ccR}}\left( t \right)} \right)}^2}} \over {2\sigma _{ccR}^2\left( {t} \right)}}} \right)} \right)} $$

The Prior Calibration Process

The prior calibration process involves moving the p _prior distribution along a sliding window of calendar ages and each time computing the the hypothetical laboratory mean ¹⁴C activity determination (h _det ) that would result from each p _prior placed at a sliding window starting at each t:

(7) $${h_{det}}\left( t \right) = \mathop \sum \nolimits_{r = 1}^{r\left( k \right)} \left( {{A_{ccR}}\left( {t + r - 1} \right) \cdot {p_{14C}}\left( {T{\rm{|}}t + r - 1} \right) \cdot {p_{prior}}\left( r \right) \cdot a\left( {t + r - 1} \right)} \right)$$

Subsequently, it is possible to evaluate the single probability value of each h _det (t) as a function of its closeness to the normal distribution of the sample’s observed laboratory determination A _det ± σ _det :

(8) $${p_{{h_{det}}}}\left( t \right) = {1 \over {{\sigma _{det}}\left( t \right)\sqrt {\left( {2\pi } \right)} }}{\rm{exp}}\left( {{{ - {{\left( {{h_{det}}\left( t \right) - {A_{det}}} \right)}^2}} \over {2\sigma _{det}^2}}} \right)$$

For each sliding window placed at each t, a vector of calibrated age probabilities is calculated, corresponding to each discrete age in the sliding window:

(9) $${p_{cal}}\left( t \right) = {p_{{h_{det}}}}\left( t \right) \cdot \left( {{p_{prior}}\left( {r,r + 1, \ldots ,r\left( k \right)} \right) \odot a\left( {t,t + 1, \ldots ,t + r\left( k \right) - 1} \right)} \right)$$

Subsequently, each p _cal (t) is sorted into a large matrix, referred to here as M _cal (T):

(10) $${M_{cal}}(T) = \left[ {\matrix{ {{p_{cal}}{{\left( {{t_1}} \right)}_{\rm{1}}}} & {{p_{cal}}{{\left( {{t_1}} \right)}_{\rm{2}}}} & \ldots & {{p_{cal}}{{\left( {{t_1}} \right)}_{\rm{n}}}} &0 & 0 \cr 0 & {{p_{cal}}{{\left( {{t_2}} \right)}_{\rm{1}}}} & {{p_{cal}}{{\left( {{t_2}} \right)}_{\rm{2}}}}&\ldots & {{p_{cal}}{{\left( {{t_2}} \right)}_{\rm{n}}}} &0 \cr 0& {\ddots} & {\ddots} & {\ddots} & {\ddots} & 0 \cr 0& 0& {{p_{cal}}{{\left( {{t_n}} \right)}_{\rm{1}}}} & {{p_{cal}}{{\left( {{t_n}} \right)}_{\rm{2}}}} & \ldots& {{p_{cal}}{{\left( {{t_n}} \right)}_{\rm{n}}}} \cr } } \right]$$

The final credible calibrated probability distribution corresponding to all ages T can be calculated simply by summing all rows in M _cal (T):

(11) $${p_{cal}}\left( T \right) = \mathop \sum \nolimits_{i = 1}^n {M_{cal}}{\left( T \right)_{ij}}$$

All elements in the resulting vector p _cal (T) are subsequently normalized such that they sum to 1.

Script for Automated Calibration (biocal)

Here, a fully documented Matlab function (biocal.m) is provided for automated calculation of the calibration protocol outlined in this study, with full compatibility in Octave. Other programming language versions of the script (e.g., Python, Julia, R) are forthcoming and will be uploaded to the same software repository upon completion. The biocal script takes full advantage of computer memory to carry out calculations using vectorized programming, thus resulting in a time-optimized routine. In the calibration protocol described in the previous section, it is assumed that it is possible to calculate P _prior sliding windows along the the entire history covered by the calibration curve. However, as it would be computationally prohibitive to calibrate for the entire history of the calibration curve, biocal restricts its P _prior sliding window calculations to an interval of the calibration curve covering a 3σ distance in each direction from the laboratory ¹⁴C determination, with added padding to accommodate a long tail of P _prior sitting at +3σ distance. In future, when computer memory and processor power increases by another order of magnitude, it will be possible to compute sliding windows across the entire calibration curve, assuming that would ever be deemed necessary. For now, as long as the tails of the final calibrated probability distribution gradually fall to very small values near to zero, we can know that a sufficient interval of the calibration curve has been considered.

The calculation time and memory usage for biocal increases with decreasing SAR, increasing BD, increasing ¹⁴C measurement uncertainty, increasing calibration curve uncertainty and increasing reservoir effect uncertainty. Testing using Matlab 2020a on a Linux system with an Intel i7-9700 CPU resulted in the following times and memory usage: a Younger Dryas aged sample with SAR of 4 cm ka^–1 and BD of 10 cm required 1.7 s calculation time and 2GB memory; the same sample, but with a SAR of 20 cm/ka^–1, required 0.2 s to calculate and used 100 MB of memory.

GROUND-TRUTH EVALUATION

Evaluating Calibration Using Sedimentological Priors

Here, a test is carried out to determine if the calibration protocol incorporating sedimentological priors results in an improved calibration process (i.e., a better estimation of the true age distribution of the measured sample) for a number of SAR scenarios and using a globally representative BD of 10 cm. First, the established understanding of bioturbation’s effect upon age-depth models on geological timescales (Berger and Heath 1968; Guinasso and Schink 1975; Peng et al. 1979; Trauth et al. 1997; Trauth 1998; Dolman and Laepple 2018; Lougheed 2020) is used to calculate the associated annualized age distribution that would be expected for a discrete-depth, 1 cm sediment sample. For all scenarios, the mean value of the age distribution is set at 12 ka, and it is assumed that the oldest 10% of the foraminifera are broken foraminifera that are not picked and, therefore, not included in the distribution. These age distributions represent the ground-truth age distribution of our virtual sample (represented as solid blue lines in Figure 3), the target age distribution that a calibrated age distribution can be judged against.

Figure 3 Comparing the new ¹⁴C calibration protocol to the existing ¹⁴C calibration method, in the case of samples with a mean age of 12 ka, constant species abundance and various sedimentological prior scenarios. Shown in all panels: the ground-truth age distribution (solid blue line); the age distribution estimated using the new ¹⁴C calibration protocol with sedimentological priors (dashed orange line); the age distribution estimated using the traditional ¹⁴C calibration method (filled yellow area). A_det is the expected mean ¹⁴C activity determination resulting from the ground-truth age distribution according to IntCal20. The following scenarios are considered as sedimentological priors: Panel A: SAR 4 cm ka^–1, BD 10 cm; Panel B: SAR 6 cm ka^–1, BD 10 cm; Panel C: SAR 8 cm ka^–1, BD 10 cm; Panel D: SAR 10 cm ka^–1, BD 10 cm; Panel E: SAR 12 cm ka^–1, BD 10 cm; Panel F: SAR 14 cm ka^–1, BD 10 cm; Panel G: SAR 16 cm ka^–1, BD 10 cm; Panel H: SAR 18 cm ka^–1, BD 10 cm; Panel I: SAR 20 cm ka^–1, BD 10 cm. (Please see electronic version for color figures.)

Subsequently, we can carry out a “virtual AMS analysis” upon the ground-truth distribution by using the IntCal20 (Reimer et al. 2020) calibration curve to determine the mean ¹⁴C activity that could be expected, in a best-case scenario, to result from the aforementioned age distribution. For simplicity’s sake, no reservoir effect is included in this demonstration, and it is assumed that the mean ¹⁴C activity reported by IntCal20 perfectly represents the ¹⁴C activity recorded by the sediment archive, with linear interpolation applied to IntCal20 where necessary to achieve annual resolution.

Assuming an appropriate ¹⁴C measurement uncertainty of ±80 ¹⁴C yr, the mean ¹⁴C activity can then be calibrated in two ways, which can subsequently be compared to each other: (1) using IntCal20 and Matcal 3.1 (Lougheed and Obrochta 2016) to carry out the existing, standard ¹⁴C calibration procedure following, e.g., Bronk Ramsey (2008), shown in Figure 3 as filled yellow areas; (2) using the aforementioned biocal in combination with IntCal20, supplemented by the SAR and BD priors associated with each scenario, to carry out the new calibration protocol outlined in this study, which is represented in Figure 3 as broken orange lines.

As could be expected, the calibration protocol using sedimentological priors outperforms the standard calibration procedure in estimating the ground-truth age distribution, as shown in Figure 3 for a number of SAR scenarios ranging between 4 and 20 cm ka^–1, with a BD of 10 cm and constant temporal species abundance. In such use case scenarios, using the calibration protocol with sedimentological priors demonstrably leads to a more accurate calibrated age distribution, which would be ideal for improving age-depth modeling of low SAR sediment archives.

In Figure 4, we repeat the same SAR scenarios as previously, but in the case of a much older ground-truth scenario (mean age of 32 ka), whereby Gaussian uncertainties associated with both the sample ¹⁴C activity (±300 ¹⁴C yr assumed here) and the ¹⁴C calibration curve are both markedly increased. In Figure 4(e–i), it can be seen that these larger uncertainties, when combined with increasing SAR, lead to the sedimentological priors becoming overwhelmed by the Gaussian ¹⁴C uncertainties and, consequently, the calibrated age distribution determined by the procedure starts to approach a normal distribution. In these use case scenarios, the new calibration protocol using sedimentological priors does not necessarily offer any advantage over the traditional calibration method.

Figure 4 Comparing the new ¹⁴C calibration protocol to the existing ¹⁴C calibration method, in the case of samples with a mean age of 32 ka, constant species abundance and various sedimentological prior scenarios. Shown in all panels: the ground-truth age distribution (solid blue line); the age distribution estimated using the new ¹⁴C calibration protocol with sedimentological priors (broken orange line); the age distribution estimated using the traditional ¹⁴C calibration method (filled yellow area). A_det is the expected mean ¹⁴C activity determination resulting from the ground-truth age distribution according to IntCal20. The following scenarios are considered as sedimentological priors: Panel A: SAR 4 cm ka^–1, BD 10 cm; Panel B: SAR 6 cm ka^–1, BD 10 cm; Panel C: SAR 8 cm ka^–1, BD 10 cm; Panel D: SAR 10 cm ka^–1, BD 10 cm; Panel E: SAR 12 cm ka^–1, BD 10 cm; Panel F: SAR 14 cm ka^–1, BD 10 cm; Panel G: SAR 16 cm ka^–1, BD 10 cm; Panel H: SAR 18 cm ka^–1, BD 10 cm; Panel I: SAR 20 cm ka^–1, BD 10 cm.

Additionally, it is also possible to revisit the 5 cm ka^–1 scenario from Figure 1, where it was shown that the traditional calibration method would misrepresent the age distribution of bioturbated (deep-sea) sediment. The new calibration protocol using sedimentological priors is applied to the same simulated sediment core (Figure 5), resulting in a much-improved calibration, whereby the 95.4% age interval predicted by the new calibration protocol provides an almost complete overlap with the actual 95.4% age interval of the single particle population. There remain some minor age-depth artifacts which result from single particles during periods of highly dynamic Δ¹⁴C (e.g., the last deglaciation) being mixed into the same discrete depths. This is an unavoidable fact of ¹⁴C dating of bioturbated sediment records, so researchers should remain vigilant when interpreting apparent SAR changes during periods of highly dynamic Δ¹⁴C. However, when one uses the new calibration protocol detailed here, the relative effect of these age-depth artifacts is reduced due to the much more realistic and wider calibrated age confidence intervals.

Figure 5 5 cm ka^–1 sediment simulation of single particles using a global average BD of 10 cm (Trauth et al. 1997; Boudreau 1998) and best-case 10⁴ particles per cm. Shown also is the discrete 1 cm depth median age, as well as the associated 68.2% and 95.4% age range. Also shown are the calibrated age distributions that would result if one were to use the new calibration protocol outlined in this manuscript to calibrate the mean ¹⁴C activity resulting from all the particles contained in each 1 cm discrete depth. Specifically, the biocal routine is applied to the mean ¹⁴C age of each discrete depth, with a SAR prior of 5 cm ka^–1 and a BD prior of 10 cm. The single particle simulation is carried out as in Figure 1.

Evaluating Calibration Using Sedimentological and Abundance Priors

Temporal changes in species abundance (e.g., of foraminifera) will affect the shape of the species’ age distribution for a given discrete depth. Here, a sine wave with a wavelength of 2000 yr is used, purely for demonstrational purposes, as a theoretical temporal abundance function (Figure 6). In Figure 7, the same SAR scenarios as in Figure 3 are analyzed, but this time with the application of the abundance aspect. Firstly, the aforementioned sinusoidal temporal abundance function is applied to the ground truth distribution. Subsequently, the same abundance function is used as an additional prior input when running biocal, to complement the sedimentological priors. The results in Figure 4 demonstrate how known information about temporal changes in species abundance can be used to produce better informed calibrated age estimations for bioturbated sediment archives. In Figure 8, the 2000-year wavelength abundance function is also applied to in the case of an older ground-truth distribution, demonstrating that abundance priors can also be used as a tool to better constrain ¹⁴C analysis of older samples that have greater uncertainty.

Figure 6 Visualization of the theoretical species abundance function used in this study to demonstrate the incorporation of prior information about species abundance in the ¹⁴C calibration protocol developed in this study. The abundance function is implemented as a sine wave with a wavelength of 2000 yr.

Figure 7 Comparing the new ¹⁴C calibration protocol to the existing ¹⁴C calibration method, in the case of samples with a mean age of 12 ka, temporally dynamic species abundance and various sedimentological prior scenarios. Shown in all panels: the ground-truth age distribution (solid blue line); the age distribution estimated using the new ¹⁴C calibration protocol with sedimentological priors (broken orange line); the age distribution estimated using the traditional ¹⁴C calibration method (filled yellow area). A_det is the expected mean ¹⁴C activity determination resulting from the ground-truth age distribution according to IntCal20. The following scenarios are considered as sedimentological priors: Panel A: SAR 4 cm ka^–1, BD 10 cm; Panel B: SAR 6 cm ka^–1, BD 10 cm; Panel C: SAR 8 cm ka^–1, BD 10 cm; Panel D: SAR 10 cm ka^–1, BD 10 cm; Panel E: SAR 12 cm ka^–1, BD 10 cm; Panel F: SAR 14 cm ka^–1, BD 10 cm; Panel G: SAR 16 cm ka^–1, BD 10 cm; Panel H: SAR 18 cm ka^–1, BD 10 cm; Panel I: SAR 20 cm ka^–1, BD 10 cm.

Figure 8 Comparing the new ¹⁴C calibration protocol to the existing ¹⁴C calibration method, in the case of samples with a mean age of 32 ka, temporally dynamic species abundance and various sedimentological prior scenarios. Shown in all panels: the ground-truth age distribution (solid blue line); the age distribution estimated using the new ¹⁴C calibration protocol with sedimentological priors (broken orange line); the age distribution estimated using the traditional ¹⁴C calibration method (filled yellow area). A_det is the expected mean ¹⁴C activity determination resulting from the ground-truth age distribution according to IntCal20. The following scenarios are considered as sedimentological priors: Panel A: SAR 4 cm ka^–1, BD 10 cm; Panel B: SAR 6 cm ka^–1, BD 10 cm; Panel C: SAR 8 cm ka^–1, BD 10 cm; Panel D: SAR 10 cm ka^–1, BD 10 cm; Panel E: SAR 12 cm ka^–1, BD 10 cm; Panel F: SAR 14 cm ka^–1, BD 10 cm; Panel G: SAR 16 cm ka^–1, BD 10 cm; Panel H: SAR 18 cm ka^–1, BD 10 cm; Panel I: SAR 20 cm ka^–1, BD 10 cm.

ADVICE FOR DETERMINING PRIOR VALUES

In order to carry out the calibration protocol detailed here, prior values for SAR, BD, fraction broken foraminifera, temporal species abundance and temporal reservoir effect are required. A first order estimate for the sediment accumulation rate can be ascertained by examining the general relationship between age-depth determinations (including ¹⁴C-derived age estimates based on existing calibration methods without sedimentological priors). This approach does represent a Catch-22 situation, however: we need an approximate indication of the age-depth relationship to determine the SAR prior, but the combination of SAR and Δ¹⁴C history can influence the ¹⁴C age distribution shape for a particular sediment interval, and hence apparent age, of the sediment archive. It would be prudent, therefore, to test a number of realistic SAR priors and examine the consequences for geochronological interpretation.

It is possible to use an approximate prior for BD using an estimate based on globally representative values (generally between 8 and 12 cm) (Trauth et al. 1997; Boudreau 1998). One could also directly estimate for the sediment archive itself based on ¹⁴C investigations of the core top (Peng et al. 1979; Trauth et al. 1997; Henderiks et al. 2002), or by using ¹⁴C measurements on single foraminifera (Lougheed et al. 2018) or, more accessibly, by measuring ¹⁴C on a number of samples with low numbers of foraminifera and using a statistical analysis of the sample variation to infer downcore bioturbation depth (Dolman et al. 2021).

The fraction of unpicked, fragmented microfossils can be estimated by simply investigating the sample material (Le and Shackleton 1992). There is a risk, however, that the very oldest microfossils of the original population are completely dissolved and are therefore no longer present in the sample material as broken material (Ruddiman and Heezen 1967), which could affect assumptions regarding the p _prior age distribution. In any case, one can take into account the susceptibility of a particular species to breakage (Boltovskoy 1991; Boltovskoy and Totah 1992) in combination with knowledge of bottom water chemistry (Ruddiman and Heezen 1967; Parker and Berger 1971), as well as the average residence time in the bioturbation zone, itself a function of SAR and BD (Lougheed et al. 2020).

Additional challenges are associated with determining temporal changes in species abundance, seeing as the abundance record estimated from the depth domain (i.e., the downcore, discrete-depth record) is itself modified by bioturbation (Lougheed 2020), and therefore does not reflect the original species abundance signal in the time domain. Species abundance in the time domain, which is called for in the calibration protocol outlined here, could be based on an estimate from, e.g., a transient palaeoclimate model run linked to an ecological model (Lombard et al. 2011; Morard et al. 2013; Roche et al. 2018; Metcalfe et al. 2020), although estimating relative temporal abundance of a species using such an approach remains a challenging task. Temporal reconstructions of abundance represent an inherent difficulty for the interpretation not just of ¹⁴C chronological data, but downcore, multi-specimen microfossil records in general (Hutson 1980; Boyle 1984; Bard 2001; Löwemark and Grootes 2004; Löwemark et al. 2008; Lougheed 2020). If one is simply not aware of the temporal abundance history at a site, a suitable approach could involve applying multiple plausible abundance scenarios when calibrating ¹⁴C dates using the calibration protocol outlined here and examining if the spread of calibrated age outcomes significantly affects the geochronological interpretation. Such an approach is similar to the current state of the art, when one might reasonably experiment with multiple reservoir effect scenarios or calibration curve versions.

CONCLUSION

Current ¹⁴C calibration workflows for sediment archives do not incorporate information about sedimentological processes such as SAR and BD, meaning that current ¹⁴C-based geochronologies systematically underestimate the total age range of a multi-specimen sample, and potentially also contain age-depth artifacts. By taking into account sedimentological processes in addition to ¹⁴C uncertainties, a more credible calibrated age distribution can be ascertained using the protocol outlined here. This new calibration protocol offers most improvement in the case of lower SAR typical of deep-sea sediment archives. It should be noted, however, that SAR itself can influence the age distribution (and hence ¹⁴C activity distribution) of a sample, but in order to determine the SAR prior accurately one needs to know the approximate age-depth relationship of the sediment. This Catch-22 type situation inherently limits high-temporal resolution geochronological analysis of deep-sea sediment, so an exploratory approach involving a range of plausible scenarios could help understand consequences for geochronological interpretation and allow researchers to test the effect of their assumptions. Such an approach can be facilitated by the computerized implementation (biocal) of the calibration protocol presented here, allowing for many scenarios to be rapidly explored. This time-efficient, vectorized computer script could be ported to and included in existing geochronological software packages typically applied to sediment archives (Bronk Ramsey 1995; Haslett and Parnell 2008; Parnell et al. 2008; Blaauw 2010; Blaauw and Christen 2011; Lougheed and Obrochta 2019), thus leading to improved age-depth chronologies, and ultimately improving the accuracy of geochronological interpretation of sediment archives.