Luminescence age calculation through Bayesian convolution of equivalent dose and doserate distributions: the D_{e}_D_{r} model
 ^{1}Archéosciences Bordeaux, UMR 6034 CNRS – Université Bordeaux Montaigne, Esplanade des Antilles, 33600 Pessac, France
 ^{2}Laboratoire de Mathématiques Jean Leray, UMR 6629 CNRS – Nantes Université, 44000 Nantes, France
 ^{3}Geography & Earth Sciences, Aberystwyth University, Llandinam Building, Penglais Campus, Aberystwyth, SY23 3DB, Wales, UK
 ^{1}Archéosciences Bordeaux, UMR 6034 CNRS – Université Bordeaux Montaigne, Esplanade des Antilles, 33600 Pessac, France
 ^{2}Laboratoire de Mathématiques Jean Leray, UMR 6629 CNRS – Nantes Université, 44000 Nantes, France
 ^{3}Geography & Earth Sciences, Aberystwyth University, Llandinam Building, Penglais Campus, Aberystwyth, SY23 3DB, Wales, UK
Correspondence: Norbert Mercier (norbert.mercier@ubordeauxmontaigne.fr)
Hide author detailsCorrespondence: Norbert Mercier (norbert.mercier@ubordeauxmontaigne.fr)
In nature, each mineral grain (quartz or feldspar) receives a dose rate (D_{r}) specific to its environment. The doserate distribution therefore reflects the microdosimetric context of grains of similar size. If all the grains were well bleached at deposition, this distribution is assumed to correspond, within uncertainties, with the distribution of equivalent doses (D_{e}). The combination of the D_{e} and D_{r} distributions in the D_{e}_D_{r} model proposed here would then allow calculation of the true depositional age. If grains whose D_{e} values are not representative of this age (hereafter called “outliers”) are present in the D_{e} distribution, this model allows them to be identified before the age is calculated, enabling their exclusion. As the D_{e}_D_{r} approach relies only on the D_{r} distribution to describe the D_{e} distribution, the model avoids any assumption about the shape of the D_{e} distribution, which can be difficult to justify. Herein, we outline the mathematical concepts of the D_{e}_D_{r} approach (more details are given in Galharret et al., 2021) and the exploitation of this Bayesian modelling based on an R code available in the R package “Luminescence”. We also present a series of tests using simulated D_{r} and D_{e} distributions with and without outliers and show that the D_{e}_D_{r} approach can be an alternative to available models for interpreting D_{e} distributions.
 Article
(1937 KB) 
Supplement
(3449 KB)  BibTeX
 EndNote
For luminescence dating of sediments, the development of equipment to perform optically stimulated luminescence (OSL) analyses at the singlegrain (SG) level (Duller et al., 1999a, b) has been a significant technological breakthrough, offering the possibility to produce a distribution of individual equivalent doses (D_{e}) for a given sample. This advance has also fostered the development of statistical approaches to analyse these D_{e} distributions (e.g. Galbraith et al., 1999; Roberts et al., 2000; Fuchs and Lang, 2001; Lepper and McKeever, 2002; Thomsen et al., 2007; Woda and Fuchs, 2008; Cunningham and Wallinga, 2012; Cunningham et al., 2015; Guibert et al., 2017; Guérin et al., 2017). Most of these statistical models target the component comprising the grains whose deposition is relevant for the event to be dated (i.e. the target population) and calculate a (believed) representative D_{e} value from this identified subpopulation. The latest proposed model (Li et al., 2021) follows the same strategy but allows identifying outliers not representative of the depositional event for several different reasons. Therefore, all these approaches focus only on the D_{e} distribution and require assumptions on how the individual D_{e} values are distributed. It is also worth recalling here that the mean environmental dose rate (D_{r}) representative for the grains constituting the selected subpopulation has to be determined with confidence.
In parallel to these developments, a series of investigations approached the dose rate as a cause of dispersion of the individual D_{e} values. These investigations were experimental (Kalchgruber et al., 2003; Cunningham et al., 2012) and/or numerical (Nathan et al., 2003; Mayya et al., 2006; Guérin et al., 2015). They all demonstrated that the spatial distribution of radionuclidebearing minerals such as Kfeldspars, but also micas or zircons, might become driving agents dominating the D_{e} distribution. In the literature, these microdosimetric effects are usually grouped and considered a significant source of unexplained variance (overdispersion, ext_OD). Another important source of external overdispersion is the presence of outlier grains (due, for instance, to sediment mixing or incomplete bleaching); this second source is in addition to the overdispersion caused by the D_{r} distribution inherent to the sample. To a lesser extent, the measurement process of the D_{e} values causes an additional dispersion. This component includes a purely experimental and a more theoretical part: the first refers mainly to the reproducibility of the measurement equipment, whereas the second relates to the fact that the protocol applied to determine individual D_{e} values is not best tailored to individual grains but represents a compromise of settings deemed optimal. The dispersion induced by these phenomena constitutes the internal overdispersion (int_OD,) which combines quadratically with the ext_OD.
Different experimental approaches (Rufer and Preusser, 2010; Romanyukha et al., 2017) have been proposed for quantifying the microdosimetric effects, whereas Martin et al. (2015a, b, 2018) and Fang et al. (2018) developed numerical sediment models to calculate the D_{r} distribution for a given granulometric fraction. Even though such experiments and applications remain rare to date, in this contribution, we want to put forward two questions: does the information characterizing the D_{r} distribution provide valuable data to calculate a luminescence age? Furthermore, if so, what would be the way to do it? Moreover, assuming that our contribution convincingly outlines an approach: how does such an approach help identify intrusive or poorly bleached grains potentially present in a D_{e} distribution?
2.1 Basics
Let us start with a thought experiment assuming the following setting: (1) one considers a series of grains of similar shape and size behaving similarly in terms of luminescence and dose response. (2) These grains are perfectly bleached and have no residual dose. (3) They are then mixed in a matrix rich in diverse radionuclidebearing mineral phases, generating a heterogeneous flux of alpha and beta particles. One also assumes (4) that the equipment used for their future analysis is perfectly reproducible. With these conditions, we expose each grain to a specific dose rate, D_{r}, which is the sum of a common gamma and cosmic dose contribution and heterogeneous alpha and beta doserate components. If we wait for 50 kyr and measure a massive number of D_{e} values from these grains, we would expect to obtain a D_{e} distribution with the same shape as the D_{r} distribution but offset by a factor of 50 000.
If the depositional setting was further complicated by supplementing the matrix of wellbleached grains of interest with a series of grains having nonzero residual doses, then superimposition of the D_{e} and D_{r} distributions could potentially identify these outliers. Consequently, our thought experiments show that thanks to the combination of the D_{e} and D_{r} distributions and without any assumption about the shape of these distributions, the depositional age can be determined even if outliers are present. The mathematical details are somewhat more cumbersome than thought experiments, and hence we will outline them in the following section.
2.2 Mathematical model
The main idea behind the D_{e}_D_{r} model is to combine the information from the D_{e} and D_{r} distributions in a Bayesian framework to detect outliers (i.e. grains not representative of the target population) automatically (if they are present) before discarding them and computing the depositional age.
2.2.1 General considerations
In real life, the number of D_{e} values measured for a sample is not extremely large. Even in cases in which thousands of grains are analysed, the low percentage of grains emitting light combined with applying a series of rejection criteria may lead to a final D_{e} distribution comprising at best a few hundred values. In contrast, when the D_{r} values are obtained by a numerical simulation of the sediment sample, for instance, their number is only limited by the lab resources in terms of computation power.
Another key difference between the D_{e} and D_{r} distributions concerns individual uncertainties: current numerical models do not report uncertainties for individual beta doserate values. This contrasts with the D_{e} values since each one has an error term related to the uncertainties associated with the luminescence signal and the process of its determination (fitting and interpolation). Nevertheless, the D_{r} distribution is not free of uncertainties: at least three terms (gamma, cosmic, and beta dose rates) must be considered, and at least two of them (gamma and cosmic dose rates) are characterized by a mean value and an associated error.
As the D_{e}_D_{r} model relies on the shape of the D_{r} distribution to describe the expected shape of the D_{e} distribution and identify outliers, the int_OD of the D_{e} distribution (such as measured with a dose recovery test – DRT) needs to be incorporated into the D_{r} distribution. To do this, individual D_{r} values (written ${\stackrel{\mathrm{\u0303}}{D}}_{\mathrm{r}}$ in Eq. 1 below) are transformed into internally overdispersed D_{r} values (D_{r}) using the following equation:
where D_{r} is a value comparable to any D_{e} value, int_OD the standard deviation characterizing the DRT distribution, and ε a Gaussian variable with uninformative mean and standard deviation (also denoted N(0,1)).
2.2.2 Mathematics underpinning the model
In this section, we reiterate the method used for detecting outliers in the frame of the hierarchical model introduced by Galharret et al. (2021). This Bayesian method can estimate an OSL age for a sample with both singlegrain equivalent dose values and simulated (or measured) doserate distributions.
We assume that the classical relation between the equivalent dose D_{e}, the corrected dose rate D_{r} (according to Eq. 1) and the OSL age A,
is satisfied but applies to the probability distributions. More precisely, we assume that the probability distribution of D_{e} is equal to the probability distribution of A×D_{r}.
To determine A, the first step of the process is to estimate the sample's D_{r} distribution when the internal overdispersion of the D_{e} distribution is incorporated, as described in Eq. (1). Because of the wide variety of possible distributions, we chose a Gaussian finite mixture with an unknown number of components. This is a very flexible class of distributions allowing us to catch symmetric, asymmetric, and multimodal distributions. Note that a Gaussian finite mixture model is a weighted sum of K Gaussian distributions $\sum _{k=\mathrm{1}}^{K}{\dot{p}}_{k}N\left({\dot{\mathit{\mu}}}_{k},{\dot{\mathit{\sigma}}}_{k}^{\mathrm{2}}\right)$. All the model parameters $\left(K{\dot{p}}_{\mathrm{1}}{\dot{p}}_{K}{\dot{\mathit{\mu}}}_{\mathrm{1}}{\dot{\mathit{\mu}}}_{K}{\dot{\mathit{\sigma}}}_{\mathrm{1}}{\dot{\mathit{\sigma}}}_{K}\right)$ can be easily estimated using an expectation–maximization (EM) algorithm (Dempster et al., 1977) and the optimal value of the number of components K selected according to the Bayesian information criterion (BIC). This method is implemented in the R package “mclust” (see Scrucca et al., 2016, for details on statistical and numerical aspects). After fitting the mixture parameters on the D_{r} distribution, we fix their values for the rest of the modelling. According to Eq. (2), the distribution of the D_{e} values is also approximated by a Gaussian finite mixture model with the following parameters.
The second step is to estimate A considering any outliers present and the measurement errors on the D_{e} values, which are assumed to be Gaussian with zero mean and known variance. Here, the main idea of the modelling is to associate each measured D_{e} with an individual age. We denote as ${a}_{\mathrm{1}},\mathrm{\dots},{a}_{n}$ these individual ages which are related to age A as follows:
where ${\mathit{\epsilon}}_{i},\mathrm{\dots},{\mathit{\epsilon}}_{n}$ represent independent Gaussian distributions with a zero mean. In the absence of outliers, we can assume that these errors have a common variance. The density of the prior variance is then
(see Galharret et al., 2021). This probability distribution is named a shrinkage distribution with parameter ${s}_{\mathrm{0}}^{\mathrm{2}}$. This is a usual choice of prior for the variance parameter for metaanalysis models (see Spiegelhalter et al., 2004). The parameter ${s}_{\mathrm{0}}^{\mathrm{2}}$ allows controlling the dispersion of the individual ages ${a}_{\mathrm{1}},\mathrm{\dots},{a}_{n}$ around A. Note that a preliminary estimate of individual ages is necessary to get an order of magnitude of the age errors. To do that, we consider the shrinkage parameter to be the harmonic mean of the variance of the individual ages. This choice ensures that errors on individual ages are not favoured over the dispersion of the individual ages ${a}_{\mathrm{1}},\mathrm{\dots},{a}_{n}$ and vice versa. In other words, neither is assumed to be negligible relative to the other, both having the same weight under the prior information.
At this step, we may refer to this model as a Bayesian central age model (BCAM) because it can be viewed as a Bayesian version of the seminal central age model (Galbraith et al., 1999) even though differences exist, the most important being the absence of any predefined function representing the D_{e} distribution. However, this model is not robust to the presence of outliers. Hence, before estimating A, we add an additional step to detect and remove the outliers if they are present in the D_{e} distribution.
In this additional step, we adapt the BCAM in including individual random effects. This is the same principle as applied in the event model introduced by Lanos and Philippe (2017, 2018). It amounts to the assumption that the errors ${\mathit{\epsilon}}_{i},\mathrm{\dots},{\mathit{\epsilon}}_{n}$ have individual variances ${\mathit{\sigma}}_{\mathrm{1}}^{\mathrm{2}},\mathrm{\dots},{\mathit{\sigma}}_{n}^{\mathrm{2}}$ independently and identically distributed from the same shrinkage distribution as previously chosen for the BCAM. While the event model can be used to estimate A, it suffers from a lack of precision due to the summation of individual variances. Thus, in our approach, we use the posterior distribution of individual variances ${\mathit{\sigma}}_{\mathrm{1}}^{\mathrm{2}},\mathrm{\dots},{\mathit{\sigma}}_{n}^{\mathrm{2}}$ for constructing a decision rule to detect outliers. Indeed, these parameters measure the dispersion of individual ages around the central age. Therefore, if an equivalent dose is detected as an outlier, its corresponding individual age will take large values with respect to the prior information on ${\mathit{\sigma}}_{\mathrm{1}}^{\mathrm{2}},\mathrm{\dots},{\mathit{\sigma}}_{n}^{\mathrm{2}}$. Thus, a D_{e} value is identified as an outlier if the posterior distribution of its individual variance is stochastically greater than its prior distribution. To do that, we use quantiles and compare the prior and posterior distributions. More precisely, we fix a probability 1−α close to 1 (for instance, $\mathrm{1}\mathit{\alpha}=\mathrm{0.95}$): if the posterior (1−α) quantile is greater than the prior (1−α) quantile, the associated D_{e} is tagged as an outlier and removed from the D_{e} distribution (Fig. 1).
When this selection is completed, the age A is estimated with BCAM from the D_{e} distribution with the outliers removed, while the posterior distributions are approximated from Markov chain Monte Carlo (MCMC) samples. In practice, we use the Gibbs sampler JAGS (Plummer, 2003) through the associated R (R Core Team, 2021) package “rjags” (Plummer, 2019).
2.2.3 Original data and structure of the model
Input data for the model are values from the ${\stackrel{\mathrm{\u0303}}{D}}_{\mathrm{r}}$ and D_{e} distributions. The D_{e} distribution is a series of central values with associated errors, whereas the ${\stackrel{\mathrm{\u0303}}{D}}_{\mathrm{r}}$ distribution represents the probability of each doserate value. Additionally, the internal overdispersion (int_OD) obtained from the DRT experiment is required. This parameter is used to modify the D_{r} distribution to be the same shape as the expected D_{e} distribution.
In summary (Fig. 2), the mathematical model to combine the D_{e} and D_{r} distributions consists of four steps.

Each ${\stackrel{\mathrm{\u0303}}{D}}_{\mathrm{r}}$ value is transformed according to Eq. (1), considering the int_OD value.

The D_{r} distribution is fitted with a weighted sum of normal (Gaussian) densities. The number of functions and their height and width are automatically adjusted to maximize the likelihood function (Fig. 3).

After a rough estimation of the individual ages (corresponding to the D_{e} values divided by the mean dose rate), the “distance” of each D_{e} value and its uncertainty with the model are computed using an MCMC process and compared to a fixed threshold set to 5 %. D_{e} values scoring lower than 95 % are considered outliers (Fig. 4).

Finally, D_{e} values corresponding to the identified outliers are removed from the D_{e} distribution, and the age is computed by the Bayesian central age model from this new D_{e} distribution. The cumulative probability distribution of the resulting model is then compared with this new D_{e} distribution and the original data (Fig. 5).
2.3 The implementation of BCAM in R
The mathematical model was implemented in R and is available in the package
“Luminescence” (Kreutzer et al., 2012) version ≥ 0.9.16
(Kreutzer et al., 2022) under the function name combine_De_Dr()
. The D_{e} and D_{r} distributions can be imported
directly from an Excel™ spreadsheet or CSV file or simply passed as a data.frame
(a data object in R, comparable to a spreadsheet) imported through other formats. The other values are directly passed to the function as parameters.
The function combine_De_Dr()
returns four plots (Figs. S2–S3 in the Supplement): the first two figures are related to detecting outliers and illustrate the variation of the individual standard deviation of the posterior age distributions. The last two figures show a kernel density plot of the posterior ages and the empirical cumulative distribution function plot. This last figure compares the cumulative D_{e} distributions (with or without the identified outliers) with the modelled D_{e} distribution (A×D_{r}). We provide a simple example with R code in the Supplement.
Our tests rely on simulated numerical data (Supplement S2). Complex D_{r} distributions were built with a series of values (at least 1000 per series) randomly sampled from normal and/or lognormal distributions. From each obtained D_{r} distribution, 100 values were randomly drawn and multiplied by 50 to represent individual D_{e} values (the D_{r} values vary around 1 Gy ka^{−1}, and the D_{e} values are then around 50 and expressed in Gy). Each D_{e} value was then associated with an uncertainty randomly sampled from a normal distribution of relative uncertainties N(0.1,0.05).
In cases in which outliers were added to the initial D_{e} distribution, their values were randomly determined from a normal or lognormal distribution, and uncertainties were defined as mentioned in the previous paragraph.
3.1 Tests without outliers
Table 1 lists the results of tests performed using four different D_{r} distributions: (1) a single normal distribution, (2) a sum of two normal distributions, (3) a single lognormal distribution, and (4) a sum of two lognormal distributions. For each D_{r} distribution, five runs were computed, and D_{e}_D_{r} ages were calculated using the combine_De_Dr()
function. Figure 6 shows an example abanico plot (Dietze et al., 2016) of a D_{e} distribution for each type of simulated D_{r} distribution.
For the four distribution types considered in these tests, the D_{e}_D_{r} age is very close to the given age, i.e. 50 ka, demonstrating the efficiency of the D_{e}_D_{r} model. It is also worth mentioning that although we did not add outliers to the initial D_{e} distribution, a few values were identified by the model as outliers and were then discarded before the final age was calculated. However, this is not surprising and can be explained by the stochastic nature of the sampling process of the D_{e} values, which each had an associated random uncertainty.
3.2 Tests with outliers
Table 2 reports results for D_{e} distributions, including 20 outliers in addition to the original 100 values. As indicated in this table, if the initial D_{e} values were sampled from a normal distribution, the outlier values were also sampled from a normal distribution (for instance, ${X}_{i,k}\sim N(\mathrm{1.3},\mathrm{0.05})$ for $j:=\mathit{\{}\mathrm{1},\mathrm{\dots},\mathrm{50}\mathit{\}}$, $k:=\mathit{\{}\mathrm{1},\mathrm{\dots},\mathrm{20}\mathit{\}}$). Furthermore, when the 100 D_{e} values were sampled from a lognormal distribution, the 20 outlier values were also sampled from a lognormal distribution (for instance, ${X}_{i,k}\sim \mathrm{log}N(\mathrm{1.3},\mathrm{0.05})$ for $j:=\mathit{\{}\mathrm{1},\mathrm{\dots},\mathrm{50}\mathit{\}}$, $k:=\mathit{\{}\mathrm{1},\mathrm{\dots},\mathrm{20}\mathit{\}}$).
The D_{e}_D_{r} age is slightly higher than 50 ka in all cases because a few outlier values overlap randomly with the initial D_{e} distribution and were therefore not identified as outliers by the D_{e}_D_{r} approach. However, this overestimation remains low (< 5 % of the true age), whereas outliers represent almost 17 % ($\mathrm{20}/\mathrm{120}$) of the D_{e} values. Examples are illustrated in Fig. 7 as abanico plots (Dietze et al., 2016).
To test the model's performance to identify outliers when their values are close to the initial D_{e} values, we simulated normal and lognormal distributions with outliers that followed our setting from above: ${X}_{i,k}\sim N(\mathrm{1.3},\mathrm{0.05})$ for $i:=\mathit{\{}\mathrm{1},\mathrm{\dots},n$, $k:=\mathit{\{}\mathrm{1},\mathrm{\dots},\mathrm{20}\mathit{\}}$ and ${X}_{i,k}\sim \mathrm{log}N(\mathrm{1.3},\mathrm{0.05})$ for $i:=\mathit{\{}\mathrm{1},\mathrm{\dots},n$, $k:=\mathit{\{}\mathrm{1},\mathrm{\dots},\mathrm{20}\mathit{\}}$, where this time n varied from 0 to 50 (then representing between 0 % and 33 % of the initial D_{e} distribution). The results are given in Table 3 and displayed in Fig. 8. The D_{e}_D_{r} ages increase with the percentage of outliers, but the overestimation remains below 10 % of the true age in all cases. This result is particularly interesting because these simulations represent cases in which a series of poorly bleached grains (i.e. the outliers) whose D_{e} values are not significantly different from the mean D_{e} have been measured in addition to wellbleached grains (initial D_{e} values).
Results show that the D_{e}_D_{r} model works well for D_{e} distributions without outliers. It also gives satisfactory results when the D_{e} values of the outliers are significantly different from the individual D_{e} values composing the target population. On the other hand, the existence of values defined as outliers but very close to the target population may be assimilated by the model to the target population and thus not be identified as outliers. This is related to the fact that the D_{e}_D_{r} model is a “majority rule” model.
This notion of majority is vital because it sets the limits of the model's applicability. If the number of outliers is significantly larger than the number of D_{e} values representing the target population, the D_{e}_D_{r} model will combine the D_{e} and D_{r} distributions as best as possible so that a maximum of D_{e} values corresponds to (A×D_{r}) values. A visual examination of the distributions calculated by the model (e.g. Fig. 5) is therefore indispensable, as is a visual examination of the outliers identified within the distribution of individual ages (e.g. Fig. 4).
On the plus side, it is also important to recall that the D_{e}_D_{r} model does not require a predefined function to represent the D_{e} distribution. For the tests carried out previously, the type of distribution (normal, lognormal, or a mixture of those distributions) was fixed to randomly draw the simulated values only. However, the type of chosen distribution and the parameters characterizing it (mean and variance) were not supplied to the model. In other words, the D_{e}_D_{r} model did not know about those parameters.
Nevertheless, the D_{e}_D_{r} model does require a precise determination of the D_{r} distribution. To date, this distribution can be obtained either from numerical sediment models considering bulk density, grain size composition, mineralogy and the spatial distribution of radioelements (for a possible 2D approach, see Dietze et al., 2021), or it can be obtained experimentally using nuclear detectors (e.g. Romanyukha et al., 2017; Fu et al., 2022). Unfortunately, at present, such experiments are scarce and remain relatively difficult to implement. Supposing they become more common, systematic comparisons between the D_{e}_D_{r} model, which provides the most probable age, and other models leading to the D_{e} value most representative of the event to be dated will become possible in a future contribution. Moreover, perhaps cases will be observed in which the D_{r} distributions do not follow a simple distribution (typically lognormal) as already suggested by Martin et al. (2015b).
One output of the model is the posterior distribution of the A defined
through a simulated Markov chain. The highest posterior density interval
(HPDI), a region of the density curve encompassing a particular credible
interval (e.g. 68 % or 95 %), can be calculated from this distribution. The HPD, the HPDI, the mean, $\stackrel{\mathrm{\u203e}}{A}$, and the standard deviation, SD_{A}, of the posterior distribution can be calculated with the function Luminescence::plot_OSLAgeSummary()
(see Supplement S1). If the posterior distribution of the age is a Gaussian distribution, the HPDI coincides with the interval $\left[\stackrel{\mathrm{\u203e}}{A}\pm {\mathrm{SD}}_{A}\right]$ at a 68 % credible level ($\left[\stackrel{\mathrm{\u203e}}{A}\pm \mathrm{2}{\mathrm{SD}}_{A}\right]$ at 95 %). The lower and upper end of the
value of the HPDI can be supplemented to facilitate systematic errors
associated with the average total dose rate and the source dose rate of the
equipment used for the D_{e} measurements. Two approaches are feasible. (1) The typical approach consists of modifying the precision, with the standard deviation SD_{A} being replaced by $\sqrt{{\mathrm{SD}}_{A}^{\mathrm{2}}+{\stackrel{\mathrm{\u203e}}{A}}^{\mathrm{2}}{p}^{\mathrm{2}}}$, where p denotes the relative error (in %) associated with the systematic error. (2) To preserve the HPD region that
considers the possible asymmetry of the posterior distribution, the
systematic error can be modelled by $\stackrel{\mathrm{\u0303}}{A}=A(\mathrm{1}+p\mathit{\u03f5})$, where
ϵ is a standard Gaussian variable independent of the age. One
can then easily sample from the corrected age and update the HPD region. If
the probability distribution of $\stackrel{\mathrm{\u0303}}{A}$ is a Gaussian distribution, the two
approaches are equivalent.
To date, the D_{e}_D_{r} model is thus the first model that allows considering the information from the equivalent doses and dose rates simultaneously, thus offering a substantial paradigm change compared to existing approaches.
The D_{e}_D_{r} model is an alternative to statistical models to determine the target population from a D_{e} distribution. Combining the information associated with the equivalent doses and dose rates experienced by the grains during burial, the model offers the possibility to determine the age of the target population without any predefined function representing the D_{e} distribution.
Future work should focus on tests carried out on welldated samples (typically crosschecked with ^{14}C dating) to validate the D_{e}_D_{r} model experimentally. This would, however, first necessitate access to accurately and precisely determined D_{r} distributions.
The source code of the model is part of the R package “Luminescence” (≥ v0.9.16) and available openaccess under GPL3 licence conditions (https://CRAN.Rproject.org/package=Luminescence, last access: 8 September 2021; https://doi.org/10.5281/zenodo.6345291, Kreutzer et al., 2022).
The supplement related to this article is available online at: https://doi.org/10.5194/gchron42972022supplement.
NM and CT initiated the work by writing the first paper draft and an initial R script. JMG and AP developed the mathematical basis for the model. SK implemented the model in the R package “Luminescence” and ran additional tests. All authors equally contributed to the discussion and the final paper writeup.
The contact author has declared that neither they nor their coauthors have any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Sebastian Kreutzer is grateful to the European Union for providing a Marie SkłodowskaCurie grant.
This research has been supported by the Agence Nationale de la Recherche (grant no. ANR10LABX52) and the European Union's Horizon 2020 research and innovation programme (CREDit, grant no. 844457).
This paper was edited by Julie Durcan and reviewed by Alastair Cunningham and one anonymous referee.
Cunningham, A. C. and Wallinga, J.: Realizing the potential of fluvial archives using robust OSL chronologies, Quat. Geochronol., 12, 98–106, https://doi.org/10.1016/j.quageo.2012.05.007, 2012.
Cunningham, A. C., DeVries, D. J., and Schaart, D. R.: Experimental and computational simulation of betadose heterogeneity in sediment, Radiat. Meas., 47, 1060–1067, https://doi.org/10.1016/j.radmeas.2012.08.009, 2012.
Cunningham, A. C., Wallinga, J., Hobo, N., Versendaal, A. J., Makaske, B., and Middelkoop, H.: Reevaluating luminescence burial doses and bleaching of fluvial deposits using Bayesian computational statistics, Earth Surf. Dynam., 3, 55–65, https://doi.org/10.5194/esurf3552015, 2015.
Dempster, A. P., Laird, N. M., and Rubin, D. B.: Maximum Likelihood from Incomplete Data Via the EM Algorithm, J. R. Stat. Soc. Ser. B Methodol., 39, 1–22, https://doi.org/10.1111/j.25176161.1977.tb01600.x, 1977.
Dietze, M., Kreutzer, S., Burow, C., Fuchs, M. C., Fischer, M., and Schmidt, C.: The abanico plot: visualising chronometric data with individual standard errors, Quat. Geochronol., 31, 12–18, https://doi.org/10.1016/j.quageo.2015.09.003, 2016.
Dietze, M., Kreutzer, S., Fuchs, M. C., and Meszner, S.: sandbox – Creating and Analysing Synthetic Sediment Sections with R, Geochronology Discuss. [preprint], https://doi.org/10.5194/gchron202139, in review, 2021.
Duller, G. A. T., BøtterJensen, L., Murray, A. S., and Truscott, A. J.: Single grain laser luminescence (SGLL) measurements using a novel automated reader, Nucl. Instrum. Meth. B, 155, 506–514, https://doi.org/10.1016/S0168583X(99)004887, 1999a.
Duller, G. A. T., BøtterJensen, L., Kohsiek, P., and Murray, A. S.: A HighSensitivity Optically Stimulated Luminescence Scanning System for Measurement of Single SandSized Grains, Radiat. Prot. Dosimetry, 84, 325–330, https://doi.org/10.1093/oxfordjournals.rpd.a032748, 1999b.
Fang, F., Martin, L., Williams, I. S., Brink, F., Mercier, N., and Grün, R.: 2D modelling: A Monte Carlo approach for assessing heterogeneous beta dose rates in luminescence and ESR dating: Paper II, application to igneous rocks, Quat. Geochronol., 48, 195–206, https://doi.org/10.1016/j.quageo.2018.07.005, 2018.
Fu, X., Romanyukha, A. A., Li, B., Jankowski, N. R., Lachlan, T. J., Jacobs, Z., George, S. P., Rosenfeld, A. B., and Roberts, R. G.: Beta dose heterogeneity in sediment samples measured using a Timepix pixelated detector and its implications for optical dating of individual mineral grains, Quat. Geochronol., 68, 101254, https://doi.org/10.1016/j.quageo.2022.101254, 2022.
Fuchs, M. and Lang, A.: OSL dating of coarsegrain fluvial quartz using singlealiquot protocols on sediments from NE Peloponnese, Greece, Quaternary Sci. Rev., 20, 783–787, https://doi.org/10.1016/S02773791(00)000408, 2001.
Galbraith, R. F., Roberts, R. G., Laslett, G. M., Yoshida, H., and Olley, J. M.: Optical dating of single and multiple grains of Quartz from Jinmium Rock Shelter, Northern Australia: Part I, Experimental design and statistical models, Archaeometry, 41, 339–364, https://doi.org/10.1111/j.14754754.1999.tb00987.x, 1999.
Galharret, J.M., Philippe, A., and Mercier, N.: Detection of outliers with a Bayesian hierarchical model: application to the singlegrain luminescence dating method, Electron. J. Appl. Stat. Anal., North America, 14, 318–338, http://sibaese.unisalento.it/index.php/ejasa/article/view/22661, last access: 26 October 2021.
Guérin, G., Jain, M., Thomsen, K. J., Murray, A. S., and Mercier, N.: Modelling dose rate to single grains of quartz in wellsorted sand samples: The dispersion arising from the presence of potassium feldspars and implications for single grain OSL dating, Quat. Geochronol., 27, 52–65, https://doi.org/10.1016/j.quageo.2014.12.006, 2015.
Guérin, G., Christophe, C., Philippe, A., Murray, A. S., Thomsen, K. J., Tribolo, C., Urbanova, P., Jain, M., Guibert, P., Mercier, N., Kreutzer, S., and Lahaye, C.: Absorbed dose, equivalent dose, measured dose rates, and implications for OSL age estimates: Introducing the Average Dose Model, Quat. Geochronol., 41, 163–173, https://doi.org/10.1016/j.quageo.2017.04.002, 2017.
Guibert, P., Christophe, C., Urbanova, P., Guérin, G., and Blain, S.: Modeling incomplete and heterogeneous bleaching of mobile grains partially exposed to the light_Towards a new tool for single grain OSL dating of poorly bleached mortars, Radiat. Meas., 107, 48–57, https://doi.org/10.1016/j.radmeas.2017.10.003, 2017.
Kalchgruber, R., Fuchs, M., Murray, A. S., and Wagner, G. A.: Evaluating doserate distributions in natural sediments using αAl_{2}O_{3}:C grains, Radiat. Meas., 37, 293–297, https://doi.org/10.1016/S13504487(03)00012X, 2003.
Kreutzer, S., Schmidt, C., Fuchs, M. C., Dietze, M., Fischer, M., and Fuchs, M.: Introducing an R package for luminescence dating analysis, Anc. TL, 30, 1–8, http://ancienttl.org/ATL_301_2012/ATL_301_Kreutzer_p18.pdf (last access: 8 September 2021), 2012.
Kreutzer, S., Burow, C., Dietze, M., Fuchs, M. C., Schmidt, C., Fischer, M., Friedrich, J., Mercier, N., Smedley, R. K., Christophe, C., Zink, A., Durcan, J., King, G. E., Philippe, A., Guérin, G., Riedesel, S., Autzen, M., Guibert, P., Mittelstrass, D., Gray, H. J., and Galharret, J.M.: Luminescence: Comprehensive Luminescence Dating Data Analysis (0.9.19), Zenodo [code], https://doi.org/10.5281/zenodo.6345291, 2022.
Lanos, P. and Philippe, A.: Hierarchical Bayesian modeling for combining dates in archeological context, J. Société Fr. Stat., 158, 72–88, 2017.
Lanos, P. and Philippe, A.: Event date model: a robust Bayesian tool for chronology building, Commun. Stat. Appl. Methods, 25, 131–157, https://doi.org/10.29220/csam.2018.25.2.131, 2018.
Lepper, K. and McKeever, S. W. S.: An objective methodology for dose distribution analysis, Radiat. Prot. Dosimetry, 101, 349–252, https://doi.org/10.1093/oxfordjournals.rpd.a005999, 2002.
Li, B., Jacobs, Z., and Roberts, R. G.: Bayesian analysis of D_{e} distributions in optical dating: Towards a robust method for dealing with outliers, Quat. Geochronol., 67, 101230, https://doi.org/10.1016/j.quageo.2021.101230, 2021.
Martin, L., Incerti, S., and Mercier, N.: DosiVox: Implementing Geant 4based software for dosimetry simulations relevant to luminescence and ESR dating techniques, Anc. TL, 33, 1–10, http://ancienttl.org/ATL_331_2015/ATL_331_Martin_p110.pdf (last access: 8 September 2021), 2015a.
Martin, L., Mercier, N., Incerti, S., Lefrais, Y., Pecheyran, C., Guérin, G., Jarry, M., Bruxelles, L., Bon, F., and Pallier, C.: Dosimetric study of sediments at the Beta dose rate scale: characterization and modelization with the DosiVox software, Radiat. Meas., 81, 134–141, https://doi.org/10.1016/j.radmeas.2015.02.008, 2015b.
Martin, L., Fang, F., Mercier, N., Incerti, S., Grün, R., and Lefrais, Y.: 2D modelling: A Monte Carlo approach for assessing heterogeneous beta dose rate in luminescence and ESR dating: Paper I, theory and verification, Quat. Geochronol., 48, 25–37, https://doi.org/10.1016/j.quageo.2018.07.004, 2018.
Mayya, Y. S., Morthekai, P., Murari, M. K., and Singhvi, A. K.: Towards quantifying beta microdosimetric effects in singlegrain quartz dose distribution, Radiat. Meas., 41, 1032–1039, https://doi.org/10.1016/j.radmeas.2006.08.004, 2006.
Nathan, R. P., Thomas, P. J., Jain, M., Murray, A. S., and Rhodes, E. J.: Environmental dose rate heterogeneity of beta radiation and its implications for luminescence dating: Monte Carlo modelling and experimental validation, Radiat. Meas., 37, 305–313, https://doi.org/10.1016/s13504487(03)000088, 2003.
Plummer, M.: JAGS: A program for analysis of Bayesian graphical models using Gibbs sampling, 3rd International Workshop on Distributed Statistical Computing (DSC 2003), 20–22 March 2003, Vienna, Austria, DSC Work. Pap., 124, 1–10, 2003.
Plummer, M.: rjags: Bayesian Graphical Models using MCMC, R package version 410, https://CRAN.Rproject.org/package=rjags (last access: 8 September 2021), 2019.
R Core Team: R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria, https://www.rproject.org, last access: 8 September 2021.
Roberts, R. G., Galbraith, R. F., Yoshida, H., Laslett, G. M., and Olley, J. M.: Distinguishing dose populations in sediment mixtures: a test of singlegrain optical dating procedures using mixtures of laboratorydosed quartz, Radiat. Meas., 32, 459–465, https://doi.org/10.1016/S13504487(00)001049, 2000.
Romanyukha, A. A., Cunningham, A. C., George, S. P., Guatelli, S., Petasecca, M., Rosenfeld, A. B., and Roberts, R. G.: Deriving spatially resolved beta dose rates in sediment using the Timepix pixelated detector, Radiat. Meas., 106, 483–490, https://doi.org/10.1016/j.radmeas.2017.04.007, 2017.
Rufer, D. and Preusser, F.: Potential of Autoradiography to Detect Spatially Resolved Radiation Patterns in the Context of Trapped Charge Dating, Geochronometria, 34, 1–13, https://doi.org/10.2478/v1000300900144, 2010.
Scrucca, L., Fop, M., Murphy, T. B., and Raftery, A. E.: mclust 5: Clustering, Classification and Density Estimation Using Gaussian Finite Mixture Models, R J., 8, 289–317, https://doi.org/10.32614/RJ2016021, 2016.
Spiegelhalter, D. J., Abrams, K. R., and Myles, J. P.: Bayesian approaches to clinical trials and healthcare evaluation, Wiley, Chichester, ISBN 9780470092590, 2004.
Thomsen, K. J., Murray, A. S., BøtterJensen, L., and Kinahan, J.: Determination of burial dose in incompletely bleached fluvial samples using single grains of quartz, Radiat. Meas., 42, 370–379, https://doi.org/10.1016/j.radmeas.2007.01.041, 2007.
Woda, C. and Fuchs, M.: On the applicability of the leading edge method to obtain equivalent doses in OSL dating and dosimetry, Radiat. Meas., 43, 26–37, https://doi.org/10.1016/j.radmeas.2007.12.006, 2008.