Luminescence age calculation through Bayesian convolution of equivalent dose and dose-rate distributions: the D_e_D_r model

In nature, each mineral grain (quartz or feldspar) receives a dose rate (D_r) specific to its environment. The dose-rate distribution therefore reflects the micro-dosimetric 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.

1 Introduction

For luminescence dating of sediments, the development of equipment to perform optically stimulated luminescence (OSL) analyses at the single-grain (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 sub-population. 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 sub-population 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 radionuclide-bearing minerals such as K-feldspars, but also micas or zircons, might become driving agents dominating the D_e distribution. In the literature, these micro-dosimetric 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 micro-dosimetric 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 Convolution of D_e and D_r distributions

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 radionuclide-bearing 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 dose-rate 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 well-bleached grains of interest with a series of grains having non-zero 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 dose-rate 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{̃}}{D}}_{\mathrm{r}}$ in Eq. 1 below) are transformed into internally overdispersed D_r values (D_r) using the following equation:

$$\begin{array}{}\text{(1)}& D_{\mathrm{r}}={\stackrel{\mathrm{̃}}{D}}_{\mathrm{r}}(\mathrm{1}+\mathit{int}\mathit{\_}\mathit{OD}\mathit{ϵ}),\end{array}$$

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 single-grain equivalent dose values and simulated (or measured) dose-rate 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,

$$\begin{array}{}\text{(2)}& D_{\mathrm{e}}=A\times D_{\mathrm{r}},\end{array}$$

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.

$$\sum _{k=\mathrm{1}}^K{\dot{p}}_{k}N\left(A{\dot{\mathit{\mu }}}_k,A^{\mathrm{2}}{\dot{\mathit{\sigma }}}_k^{\mathrm{2}}\right)$$

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:

$$\begin{array}{}\text{(3)}& a_i=A+{\mathit{ϵ}}_i,\end{array}$$

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

$$\begin{array}{}\text{(4)}& p\left(x\right)={\displaystyle \frac{s_{\mathrm{0}}^{\mathrm{2}}}{{\left(s_{\mathrm{0}}^{\mathrm{2}}+x\right)}^{\mathrm{2}}}}\end{array}$$

(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 meta-analysis 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 pre-defined 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.

Figure 1 Comparison of prior and posterior cumulative distribution functions of individual variance and their 95 % credible interval (bottom horizontal lines): the corresponding equivalent dose is detected as an outlier (a) or not (b).

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{̃}}{D}}_{\mathrm{r}}$ and D_e distributions. The D_e distribution is a series of central values with associated errors, whereas the ${\stackrel{\mathrm{̃}}{D}}_{\mathrm{r}}$ distribution represents the probability of each dose-rate 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.

Figure 2 Diagram representing the different steps of the estimation method.

In summary (Fig. 2), the mathematical model to combine the D_e and D_r distributions consists of four steps.

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

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

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

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

Figure 3 Approximation of the D_r distribution with a mixture of normal (Gaussian) functions.

Figure 4 Characterization of the D_e values: the values in red are identified as outliers.

Figure 5 Comparison of the cumulative distribution functions: A×D_r (red line), D_e (dotted blue line), and reduced – after removal of the detected outlier values – D_e (dashed green line).

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.

3 Model tests

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 log-normal 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⁻¹, 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 log-normal 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 log-normal distribution, and (4) a sum of two log-normal 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.

Table 1 Results of tests without added outliers. Tests were performed with four different shapes of D_r distributions: Norm(n, m, SD) and log-Norm(n, m, SD) indicate normal and log-normal distributions, respectively, where (n) is the number of random values, (m) the mean of the distribution, and (SD) the standard deviation. The number of Gaussian components identified by the model when fitting the D_r distribution is given, as is the number of points identified as outliers. Numbers in bold represent average values.

Figure 6 Abanico plots of the D_e distributions (100 values) without outliers. (a) Single normal distribution, (b) sum of two normal distributions, (c) single log-normal distribution, and (d) sum of two log-normal distributions.​​​​​​​

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 log-normal distribution, the 20 outlier values were also sampled from a log-normal distribution (for instance, $X_{i,k}\sim \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{\}}$).

Table 2 Results of tests with 20 outliers added to the original D_e distribution. Their values were determined following the function indicated in the second column. Notice that the (m) parameter of these functions varied from 1.3 to 1.6, leading to outlier values which, on average, increased, as is observable on the abanico plots (Fig. 7). Numbers in bold represent average values.​​​​​​​

Figure 7 Abanico plots of the D_e distributions (100 initial values sampled from the D_r distribution) to which 20 outlier values have been added: red dots indicate the values identified as outliers (not the values added as outliers). (a–d) Initial D_e distribution = normal distribution; (e–h) initial D_e distribution = log-normal distribution. Outlier values (20) added as follows: (a) 50 × Norm(20, 1.3, 0.05); (b) 50 × Norm(20, 1.4, 0.05); (c) 50 × Norm(20, 1.5, 0.05); (d) 50 × Norm(20, 1.6, 0.05) and (e) 50 × log-Norm(20, 1.3, 0.05); (f) 50 × log-Norm(20, 1.4, 0.05); (g) 50 × log-Norm(20, 1.5, 0.05); (h) 50 × log-Norm(20, 1.6, 0.05).

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

Table 3 For each type of D_r distribution (normal or log-normal), outlier values were added following either the function 50 × Norm(n, 1.3, 0.05) or the function 50 × log-Norm(n, 1.3, 0.05). The number of outliers (n) varied from 0 to 50 (then representing between 0 % and 33 % of the initial D_e distribution). The age error is the 95 % credible interval.

Figure 8 D_e_D_r age as a function of the number of outliers added to the initial D_e distribution (the expected age is 50 ka, indicated by the red line). Norm and log-Norm represent the functions from which the initial D_e distributions (comprising 100 values) were built. Error bars represent 95 % credible intervals. Dotted lines are ±10 %.

To test the model's performance to identify outliers when their values are close to the initial D_e values, we simulated normal and log-normal 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 \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 well-bleached grains (initial D_e values).

4 Discussion

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 pre-defined function to represent the D_e distribution. For the tests carried out previously, the type of distribution (normal, log-normal, 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 radio-elements (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 log-normal) 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{‾}}{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{‾}}{A}\pm {\mathrm{SD}}_A\right]$ at a 68 % credible level ($\left[\stackrel{\mathrm{‾}}{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{‾}}{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{̃}}{A}=A(\mathrm{1}+p\mathit{ϵ})$, 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{̃}}{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.

5 Conclusion

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 pre-defined function representing the D_e distribution.

Future work should focus on tests carried out on well-dated samples (typically cross-checked with ¹⁴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.