Challenges of dating quartz OSL samples with saturated grains

A suite of samples from an extensive aeolian sandscarp near Victoria Falls, Zambia was used to explore several different methods of calculating optically stimulated luminescence (OSL) ages that account for the effects of saturated quartz grains. Beta dose rate heterogeneity and early OSL signal saturation of the samples exacerbate the impact that saturated grains have on the equivalent dose (D e ) values calculated. Saturated grains that cannot calculate D e values are often rejected but the minimum burial dose information they contain can have a significant impact on a sample ’ s average D e value. This study compares multiple techniques for combining luminescence measurements that enables inclusion of this data and their sensitivity to a criterion that rejects grains with early OSL signal saturation. The methods tested are found to have different advantages and disadvantages, but reasonable agreement between the D e values they calculate suggests that including data from saturated grains makes a more significant difference to D e values calculated than the specific method used to combine the data.


Introduction
Optically stimulated luminescence (OSL) dating measurements of quartz grains, especially from older samples, often include data from grains with saturated palaeodoses (e.g.Thomsen et al., 2016;Colarossi et al., 2020).The proportion of saturated grains can vary substantially between samples and is influenced not only by the magnitude of each grain's equivalent dose (D e ), but also by the dose at which each grain's OSL signal saturates.Dose recovery experiments suggest that measurements with more than 50% of saturated grains may not provide accurate dose estimates (Duller 2012), and some studies of natural samples have attributed underestimations of single grain ages relative to multi-grain ages to the effects of saturated grains (e.g.Thomsen et al., 2016;Colarossi et al., 2020).If these saturated grains contain useful chronological data, more accurate D e values, and hence ages, could potentially be obtained if these grains could be incorporated into the analysis.
The impact of saturated grains on D e value calculations is likely to be more significant for samples with wide palaeodose distributions, where truncation of a small subset of grains with large palaeodoses can make a substantial difference to the average D e value calculated for the sample.
It is also more likely to be significant for samples with a larger proportion of saturated grains, which can be caused by either a substantial number of grains with large palaeodoses or early OSL signal saturation.
This study aims to compare several different methods for calculating D e values that attempt to account for saturated grains on single grain distributions for a suite of samples collected from an extensive aeolian sandscarp near Victoria Falls, Zambia.This comparison will consider the strengths and weaknesses of the different approaches, assess their ability to recover doses with saturated grains, and compare D e values calculated on the same datasets.

Study area and sample preparation
Two suites of samples from a ~34 m thick sandscarp in Mosi oa Tunya National Park (17 • 49.4 ′ S, 25 • 47.9 ′ E), near Victoria Falls and the associated town of Livingstone, Zambia were collected as part of the 'Deep Roots of Humanity' archaeological project investigating the transition from the Early to Middle Stone Age in south-central Africa.The ZS section (samples 232/ZS17-XX), sampled in 2017, is further North and adjacent to the airport runway, whilst the NP section (samples 232/NP18-XX), sampled in 2018, lies ~1 km to the south-east, closer to the town of Livingstone.This study focuses on three samples from the ZS section and six samples from the NP section although all nine samples are discussed together and plotted by depth beneath the modern plateau surface.
The extensive sandscarp has a distinctive deep red colour (Fig. 1) that is in sharp contrast to the white modern fluvial deposits of the nearby Zambezi River.The deposit is massive, lacks defining sedimentary structures, and has a bimodal grain size distribution with peaks at ~150 and ~500 μm.The sandscarp is interpreted to be an aeolian sandsheet that accrued gradually through time burying artefacts and surface grains under deeper and deeper layers of sand.Middle Stone Age archaeological artefacts are prolific on the modern surfaces of the sandscarp and have also been found in-situ buried within the sands.These dating samples will provide chronological context for the archaeological findings and contribute to filling a large regional gap in our understanding of the timing of the transition from the Early to Middle Stone Age in this part of Africa (Barham et al., 2015).This period encompasses the evolution of Homo sapiens and transformative technological innovations (Coe et al., 2022), with the work at Victoria Falls providing a chronological framework to compare with the eastern and southern African archaeological records.
Samples were collected from exposed sections, associated with archaeological finds, by hammering opaque plastic tubes (~25 cm in length and ~8 cm in diameter) into the sediment.Material from the inner part of the tubes was selected to extract quartz for OSL measurements to minimise the potential for exposure during sampling.The samples were treated with H 2 O 2 and HCl to remove organics and carbonates respectively, dried, and then sieved to obtain grains 180-212 μm in diameter.This grain size was then separated using sodium polytungstate solutions at densities of 2.62 and 2.70 g/cm 3 to remove any feldspars and heavy minerals.The remaining quartz-rich material was then etched in 40% hydrofluoric acid for 40 min to remove the alpha irradiated outer layer of the quartz.Samples were rinsed thoroughly in concentrated HCl following HF etching to prevent the formation of fluorides.

Equipment and protocols
Luminescence measurements were made on a Risø TL/OSL DA15 reader equipped with a ~0.075Gy/s Sr/Y beta source.The prepared quartz grains were mounted in single grain discs containing 100 holes and stimulated using a single grain attachment based around a 10 mW Nd:YVO 4 laser emitting at 532 nm (Bøtter-Jensen et al., 2003).Luminescence signals integrated the initial 0.2 s with a late background of the final 0.2 s.An instrumental uncertainty of 2.5% was used when calculating measurement results.
A single aliquot regenerative dose (SAR) protocol incorporating an OSL IR depletion ratio test (Duller 2003) was used for all measurements (Table 1).For dose recovery tests, grains were bleached using the green laser for 1 s on each position whilst holding the sample at room temperature, left for 10,000 s and then bleached for a further 1 s on each grain.Dose recovery preheat plateau tests were performed on samples 232/NP18-06 and 232/ZS17-04 with regenerative dose preheat temperatures from 180 to 280 • C at 20 • C intervals for 10 s, and a constant 10 s heating of 160 • C used for the test dose.The given dose of 25 Gy could be recovered within 10% for preheats from 200 to 260 • C, and a preheat of 220 • C for 10 s was used for dating all samples.
For each sample, at least 2 single grain discs (200 grains) were used to undertake a dose recovery test using the 220 • C preheat.The given laboratory dose for this test was adjusted for each sample to be similar to its measured equivalent dose (D e ).Equivalent doses were measured on 12-18 single grain discs (1200-1800 grains) for each sample.The same 10 regenerative doses (ranging from 0 to 506 Gy) and test dose (2 Gy) were measured for all the grains so that analyses of dose response curves and the impact of luminescence saturation could be studied consistently.Results for each grain were only accepted if they passed basic rejection criteria including: a natural test dose signal that is more than three standard deviations above its background, a relative error of less than 20% on the test dose signal following the natural, recycling ratios within two uncertainties of 10% of unity, and an IR depletion ratio within one uncertainty of 10% of unity.The rejection threshold for the IR depletion ratio was stricter because the measurements being compared had been measured consecutively whilst the recycling ratio compared measurements from the beginning and end of the protocol sequence.Additionally, a very small number of grains (a total of 5 grains found in three of the samples) were rejected for having negative D e values.All five of these grains also recorded a significant drop in luminescence sensitivity between the natural and first regenerative dose.

Dosimetry
The environmental dose rate of each of the samples can be broken down into different components: cosmic dose rate, external gamma dose rate, external beta dose rate, and internal alpha dose rate.The contribution from any external alpha dose rate is assumed to have been  Green laser for 1 s at 125 C 7 Repeat steps for next regenerative dose a IR diodes only used for last regenerative dose as part of IR depletion ratio test (Duller 2003).
M.S. Chapot et al. removed by the HF etching of the quartz grains.The internal alpha dose rate is estimated to be 0.026 ± 0.010 Gy/ka based on thick source alpha counting of quartz separates from similar quartz rich sediments in Botswana (Tooth et al., 2022).Cosmic dose rates were calculated using the equations in Prescott and Hutton (1994) based on longitude, latitude, altitude (980 m), overburden density (1.8 g/cm 2 ) and burial depth (sample specific).As the sample specific burial depths have likely increased over time, the cosmic dose rates were calculated assuming an average burial depth of half the modern depth.External beta dose rates were calculated from thick source beta counting data (using a GM-25 beta counter, Bøtter-Jensen and Mejdahl, 1988) on subsamples of the material collected from the end of the sample tubes which had previously been dried and ground in a TEMA mill to obtain a homogenised powder.External gamma dose rates for most of the samples were measured in situ using an Ametek DigiDart portable gamma spectrometer with a 2-inch diameter NaI detector.The threshold method was used for analysis of the resulting spectrum (e.g.Mercier and Falguères, 2007) and the response of the detector was calibrated using the concrete blocks in Oxford (Rhodes and Schwenninger 2007) and checked for consistency using an internal standard housed in the Aberystwyth Luminescence Research Laboratory.These field measurements were also adjusted for average burial moisture content by un-attenuating the measured values based on the measured moisture content of the sample (1-3%) and re-attenuating the dry value by the assumed average burial moisture content (5 ± 2%).Two of the samples, ZS17-05 and ZS17-06, did not have field gamma measurements and the external gamma dose rates for these samples were calculated using high resolution gamma spectrometry data measured at the Hannover Luminescence Laboratory using a Sakura detector on 50 g subsamples.

Beta dose rate heterogeneity
Previous studies have suggested that samples with lower dose rates are more susceptible to beta dose rate heterogeneity (Mayya et al., 2006;Smedley et al., 2020) as significant beta emitters may be localised in rare heavy mineral hotspots (Burrough et al., 2019;Jankowski and Jacobs, 2018).The average dose rate for the samples in this study is only 0.57 ± 0.02 Gy/ka and external beta dose rates (0.18-0.14 Gy/ka) comprise 25-30% of the total dose rates.As less than 0.05% potassium was measured for these samples, no more than 0.04 Gy/ka beta dose rates arise from K. The remaining 75-80% of the external beta dose rate originates from U and Th presumed to be isolated in heavy mineral grains and this dose rate component may vary considerably on a grain-to-grain scale.While most of the quartz grains are likely reasonably far from a heavy mineral, the small proportion of grains adjacent to a zircon or other significant beta emitter will have a much higher dose rate.

Defining the limits for obtaining equivalent doses
Each of the nine samples have accepted D e measurements for between 190 and 442 grains (Table 3).However, all the samples have a significant proportion of grains at or beyond OSL dating range, with the exact proportion depending on how that limit is defined (Table 3).The most obviously saturated grains are those for which the L n /T n ratio lies above the maximum asymptote (I max ) of the dose response curve (DRC); such grains will be referred to as 'saturated' in this paper.There are also grains for which the upper uncertainty of the L n /T n ratio is above I max , hereafter termed 'uncertainty saturated' grains.However, even grains with natural OSL signals below I max may not yield accurate D e values, as shown by Chapot et al. (2012) and other studies (e.g.Colarossi et al., 2020).Wintle and Murray (2006) suggested a prudent upper limit of 2D 0 , where D 0 is the characteristic dose of a single saturating exponential DRC.Mathematically, 2D 0 corresponds to the dose with a L x /T x value at 86.5% of I max .Many of the DRCs in this study were fit with the sum of two saturating exponentials and we consider an upper limit dose (UL dose ) as the dose where I is equal to 86.5% of I max .
Upper limit doses vary greatly between grains for each sample and for many grains this value is 80 Gy or less (grey histograms in Fig. 2).In dose recovery data sets where all grains were given the same known laboratory dose, all the uncertainty saturated grains (blue histograms in Fig. 2) have low UL doses (all are <80 Gy, and most <40 Gy) -these grains are incapable of recording the 41.4 Gy radiation dose administered (Fig. 2a).However, the uncertainty saturated grains from D e datasets are much more spread out and similar to the distribution of UL dose values for the sample (Fig. 2b).This difference likely reflects the range of  palaeodoses recorded by the grains and may influence the effects of signal saturation on average dose values calculated for the datasets.Initial analyses on the single grain datasets in this study focused on accepted grains for which Analyst (Duller 2015) could calculate D e values, meaning that grains were not automatically excluded if they had L n /T n values that were more than 86.5% of I Max , but 'uncertainty saturated' grains were excluded.The resulting D e distributions were broad with slight positive skew (see Fig. S1 for Abanico plots), which were interpreted as resulting from beta dose rate heterogeneity.
Overdispersion values calculated using the central age model (CAM, Galbraith et al., 1999) are primarily 60-71% at the ZS site and 38-57% at the NP site but increase to ~100% for the two uppermost samples, NP18-08 and NP18-09 (Table S1).The higher percentage overdispersion for these upper samples originate from grains with high D e values and is in part because the average D e values of the samples are lower, making the percentage overdispersions higher.
Central D e values for each sample (Table S1) were calculated using the Average Dose Model (ADM, Guérin et al., 2017).This model calculates an arithmetic mean for the SG D e distribution and is appropriate to divide by the dose rate for each sample, which is an arithmetic mean of a suspected similar distribution of grain-specific dose rates.A σ m parameter of 10% (representing intrinsic variability) was used for all ADM D e calculations based on the overdispersion values of laboratory dose recovery tests (Table 4).

Dose recovery results
Dose recovery tests were undertaken for each sample with the given dose ranging from 5.1 to 95.0 Gy.Most of the samples had 200 grains measured for the dose recovery test but two of the samples (NP18-06 and NP18-05) had 1000 grains measured.As might be expected, for the three lowest given doses (5.1-15.8Gy) no uncertainty saturated grains were observed (Table 4).For given doses from 26.1 to 95.0 Gy, the proportion of uncertainty saturated grains increases up to 43% of the accepted grains.Most of the samples had excellent dose recovery, even when 43% of the accepted grains were uncertainty saturated (Table 4).However, the dose recovery for sample NP18-05 (67.8 Gy given dose, 19% uncertainty saturated), was only just within one uncertainty of 10% of unity (0.88 ± 0.02) and may have been negatively affected by the exclusion of the uncertainty saturated grains.

Accounting for saturated grains
Unlike the dose recovery measurements, in the equivalent dose datasets all the samples have a significant proportion of grains that are saturated, uncertainty saturated, or beyond their grain-specific UL dose (Table 3).This is due to the combination of the grains recording a wide range of palaeodoses and many of the grains having relatively low UL doses (e.g.Fig. 2b).This section compares five different methods of calculating D e values that account for the presence of these grains.Each of the methods are described below and results for the D e datasets and the dose recovery dataset of NP18-05 are discussed and shown in Fig. 3 and S2.
Fig. 3 is comprised of 4 subfigures showing values calculated for three samples' equivalent dose datasets and the dose recovery dataset of NP18-05 using the five different methods.Dose estimates calculated with four of these methods are plotted against minimum UL dose on the bottom X-axis.Leftmost datapoints are plotted at the lowest UL dose of any accepted grain for the sample and are the dose values calculated when no minimum UL dose threshold is applied.As increasingly strict minimum UL dose thresholds are applied, the number of grains included is reduced.Solid black reference lines in each plot depict thresholds that might be chosen based on the dose estimate calculated, either where the sample average dose is equal to the minimum UL dose or where it is equal to half the minimum UL dose .The fifth method for calculating sample average doses, the L n /T n Method, is plotted against the UL dose of a standardised DRC for a grouping of grains on the top X-axis.For these  a Measured doses were calculated using ADM.
datapoints, the number of grains included in the analysis does not necessarily decrease with increasing values.Figure S2 is comprised of similar subfigures for every sample discussed in this study and also includes plots of the number of grains included in each calculation.

Synthetic super-aliquots
In multiple-grain measurements, the data from all the grains on an aliquot are combined because they all luminesce at the same time, summing their individual light intensities for each measurement, leading to a single dose response curve being constructed for the aliquot.The data from single grain measurements can be similarly combined by summing the luminescence signals of different grains to construct synthetic aliquots (e.g.Henshilwood et al., 2002).Synthetic aliquots are most often constructed by summing the luminescence from all the grains on an individual single grain disc (e.g.Colarossi et al., 2020;Henshilwood et al., 2002).However, they can also be constructed by summing the luminescence of grains from multiple single grain discs that were measured with the same protocol, and this summation produces a synthetic super-aliquot.
In this study, we investigated three ways of constructing synthetic super-aliquots, comparing the results when summing (1) all the grains measured, (2) only the grains that passed rejection criteria, and (3) only the grains that both passed rejection criteria and could calculate a D e value in Analyst.In most cases, data from synthetic super-aliquots of all grains measured and those of all grains that passed rejection criteria are very similar (e.g.Fig. 4).This is not surprising given that most of the rejected grains are rejected for having low signal intensity and they therefore contribute very little luminescence signal to the synthetic aliquots.Although individual grains exhibit a wide spectrum of UL doses (Fig. 2), averaging effects result in extremely similar synthetic super-aliquot DRCs for each of the samples with UL doses of ~120 Gy (Fig. S3); however, the synthetic super-aliquot DRC for sample 232/ NP18-04 has a slightly higher I max value relative to the other samples (Fig. S3).
Equivalent dose values of synthetic super-aliquots that excluded 'uncertainty saturated' grains (shown as pink X's in Fig. 4) are significantly lower than those calculated for synthetic super-aliquots including  all grains or all accepted grains (shown as black and blue circles in Fig. 4).This supports the hypothesis of Colarossi et al. (2020) that differences between single grain D e values and those from multiple-grain aliquots and synthetic aliquots of the same sample are likely the result of excluding uncertainty saturated grains.Synthetic super-aliquot D e values calculated using all grains that passed rejection criteria are shown as white circles in Fig. 3 and are generally insensitive to a UL dose criterion.For the deepest sample (NP18-04, Fig. 3b), discrete D e and uncertainty values were able to be calculated, but the L n /T n ratio of the synthetic super-aliquot is 91.3% of I max .The dose recovery ratio for sample NP18-05 calculated using the synthetic super-aliquot approach (0.90 ± 0.05) is marginally closer to unity than that calculated in Section 6.
One potential weakness of calculating D e values with a synthetic super-aliquot approach is the significant weighting towards the brightest grains.For example, in sample ZS17-05, out of 319 grains accepted for the analysis, fifteen grains provide more than 50% of the summed T n signal and removing them one by one from inclusion in a synthetic aliquot analysis can make a considerable difference to the D e value calculated (Fig. 5).The influence of these grains is present but significantly muted in calculations where a D e is calculated for each grain, such as the ADM, in which they are weighted more similarly to the other grains in the dataset.
It is interesting to consider that this significant signal intensity weighting underlies all multi-grain analyses and may account for a significant proportion of inter-aliquot variability (shown in grey in Fig. 4).Such significant weighting toward a few specific grains is likely to be of less importance in dose recovery tests or for natural samples where all the grains have palaeodose values within uncertainty of one another.However, for samples with significant dose rate heterogeneity, it is critical that the sample average D e value is not dominated by an exceptionally bright quartz grain that happened to be situated next to a zircon.

Minimum upper limit dose threshold
Thomsen et al. ( 2016) proposed a method to help resolve potential biases caused by the exclusion of saturated grains.In their method, the D e value was calculated for each sample after removing grains whose D 0 value was below a given threshold and repeating this process at different values for the D 0 threshold.As the D 0 threshold is increased, this approach reduces the proportion of grains excluded due to saturation by rejecting all grains with low D 0 values regardless of whether they were saturated.In this study, we modified the approach of Thomsen et al. (2016) to make it suitable for dose response curves fitted with the sum of two exponentials by having the threshold not defined by D 0 , but by UL dose (see Fig. S6 for relationships between grains' D e values and UL dose values for each sample).
For young samples (e.g.NP18-07, red circles Fig. 3a and S2) applying a minimum UL dose rejection criterion seems to make little difference to the ADM D e value calculated.However, other samples in this study are impacted by the criterion and the ADM D e values increase as increasingly high thresholds are applied (Fig. 3b-d and S2).Many of the samples with increasing ADM D e values have clear plateau regions where the calculated ADM D e remains consistent as the threshold value is increased (Fig. 3b and d and S2) but the deepest sample analysed in this study (Fig. 3c) rejects all measured grains before a plateau can be reached.Thomsen et al. (2016) recommended using a D 0 threshold equal to the sample's average D e value, which corresponds to a UL dose threshold of two times the sample's average D e value.Lines representing this threshold, as well as a threshold of UL dose equal to the sample's average D e value (necessary for all grains to be capable of recording the sample average burial dose) have been included in Fig. 3 and Fig. S2 for reference.The minimum UL dose criterion improves the dose recovery ratio for sample NP18-05 to within commonly accepted values (0.93 ± 0.02 for a plateau at a threshold of ~80 Gy, Fig. 3d).
This method attempts to mitigate the influence of rejecting saturated and uncertainty saturated grains by progressively restricting analysis to just those grains with high UL doses which are capable of recording higher doses.This involves also rejecting grains that calculate D e values but have early signal saturation.This means that instead of increasing the number of included grains by adding data from saturated grains, we are decreasing the number of grains by rejecting otherwise accepted grains, thereby increasing uncertainties.Additionally, this method has no mechanism of accounting for the effect of saturated grains with high UL doses , which may explain why it was ineffective at resolving differences between single and multiple grain analyses in the study of Colarossi et al. (2020).

Including saturated grains
Although D e values cannot be calculated for saturated grains, they still contain information about the minimum dose they were likely to have received during burial.Uncertainty saturated grains can calculate D e values but have infinite upper uncertainties suggesting they also record minimum burial dose estimates.Grains with natural OSL signals exceeding a prudent upper limit may be able to calculate discrete D e values and uncertainties but such values are likely to have large and asymmetric uncertainties and may be unreliable (Chapot et al., 2012).Even still, when the doses recorded by these three types of grains are greater than the average dose recorded by the other grains, the data they contain may be critical to calculating a more accurate age for a sample.This is especially true for samples with significantly heterogeneous micro-dosimetry or when saturated grains are a significant proportion of a sample's single grain dataset.
In this study we investigated the impact of including the UL doses of these three types of grains in ADM calculations as if they were D e values for the grains.Uncertainties for the UL doses were approximated by using the relative uncertainty of the grain's L n /T n value.This relative uncertainty was applied to 0.865*I max and interpolated onto the grain-specific DRC to calculate the grain's lower and upper bound UL dose uncertainties.Unfortunately, calculating UL dose uncertainties in this manner will result in infinite upper bound uncertainties for some grains and significantly asymmetric uncertainties for several others.
In order to have discrete symmetrical uncertainties to include in our age model calculations, we used the lower bound upper limit dose uncertainties as if they were symmetrical.Due to the shape of DRCs, UL dose uncertainties calculated in this manner are over precise to some degree, thereby giving these grains some degree of over-weighting in models M.S. Chapot et al. calculating sample D e values.However, it is worth considering that for a saturated grain, the UL dose is a minimum estimate and if the UL dose is greater than the sample average dose, the overweighting will shift the sample's age slightly older, potentially increasing the accuracy of the D e value calculated for the sample.
Including data from these three types of grains in the manner discussed here requires decisions such as whether to use Analyst D e values or the UL dose for grains with natural OSL signals greater than 0.865*I max and whether symmetrical lower bound uncertainties should be used for all grains of these types or only in the few instances where the upper bound uncertainty is infinite.However, initial investigations suggest these different approaches have little impact, well within uncertainties, to the ADM D e value calculated (section S1).For the values reported in this paper, ADM calculations on single grain datasets including saturated grains were comprised of Analyst calculated D e and uncertainty values for grains with L n /T n values less than 0.865*I max and the UL dose with symmetrical lower bound uncertainties for grains with L n /T n values greater than 0.865*I max .
Including the UL dose as the D e value for grains that would have otherwise been excluded results in a significant difference to the ADM D e values calculated (black plus symbols, Fig. 3).In almost all cases, ADM D e values calculated including these grains are higher than those excluding them.The main exception to this is for older samples (e.g.NP18-04, Fig. 3c and NP18-05, Fig. S2g) at low minimum UL dose thresholds.In these cases, the large number of grains with low UL doses brings the sample D e value below that calculated if saturated grains are excluded, but this quickly reverses as the minimum UL dose threshold increases.
For most of the samples, applying a minimum UL dose criterion to ADM D e calculations which include UL doses for grains with L n /T n values greater than 0.865*I max results in either an initial D e value plateau (Fig. 3a) or an increase in D e value leading to one or more plateaus (Fig. 3b and d), the weighted average and uncertainty of D e values comprising these plateaus are referred to as 'inclusive plateau' D e values (shown as red horizontal lines on Fig. 3 and S2).However, even with this approach, ADM D e values for the deepest sample do not reach a plateau (Fig. 3c).The 'inclusive plateau' dose recovery ratio of 0.93 ± 0.02 for sample NP18-05 is the same result calculated using only grains that provide Analyst D e values, suggesting that the inclusion of UL doses may have a more significant effect on natural samples than laboratory irradiations.

Bayesian statistics using the BayLum R package
The ADM uses frequentist statistics and cannot cope with asymmetrical or infinite uncertainties.Bayesian statistics provides a way around these difficulties.Combès and Philippe (2017) presented Bayesian models for calculating central equivalent doses which were later developed into the BayLum R package (Christophe et al., 2017;Philippe et al., 2019).Within the BayLum package, the Palaeodose_Computation function enables D e calculations using Cauchy, Gaussian, or lognormal distributions.Heydari and Guérin (2018) conducted dose recovery tests to compare these different distribution types for Bayesian models against the ADM and CAM D e calculations for quartz grains near saturation or with variable absorbed doses.The results of their study suggested that the Gaussian and lognormal-average Bayesian models were the most accurate and did not require a minimum D 0 rejection criterion to be applied.When a minimum D 0 rejection criterion was used, the ADM also provided accurate results in their experiment.
Calculating D e values using the Palaeodose_Computation function involves uploading the BIN files containing all the OSL measurements for each grain and informing the program which grains to include in the analysis.Each grain's DRC is fitted by the program as part of the function but neither the grain-specific D e value, nor the DRC equations fitted, are provided as an output numerically or visually.We included all grains that passed rejection criteria regardless of saturation level and fit Gaussian distributions to the data.Whilst the majority of the grains' DRCs were best fit with the sum of two saturating exponential functions in Analyst, this type of function is not an available option for the Palaeodose_Computation function, which also needs to fit the same type of function to all the included grains.Therefore, in this analysis, every grain's DRC was fit with a single saturating exponential plus linear equation that was not forced through the origin.
Equivalent dose values calculated with this method tend to be greater but within uncertainties of other D e values that were calculated when including saturated grains (blue triangles, Fig. 3 and S2) and most of the samples have stable BayLum D e values with increasingly strict minimum UL dose thresholds in agreement with the findings of Heydari and Guérin (2018).Unfortunately, it is not possible to determine whether differences in D e values between this method and others including saturated grains originate from the Bayesian analyses or from the individual grain DRCs being fit with a different type of equation, which may have a significant impact for grains with natural OSL signals near saturation.However, the dose recovery ratio of 0.92 ± 0.02 for sample NP18-05 is similar to the value obtained using other methods which accounted for saturated grains.The BayLum Palae-odose_Computation also produces the oldest age estimate for the deepest sample (Fig. 6).

L n /T n Method
Another method for including data from saturated grains was proposed by Li et al. (2017) and involves applying a frequentist age model (e.g.weighted-mean, ADM, CAM) to the L n /T n distribution of accepted grains and interpolating the results onto a re-normalised standardised DRC for the sample.One of the greatest challenges of applying this L n /T n method to quartz grains is the enormous variety in DRC shape between grains, even within the same sample (Fig. 2).In their publication testing this L n /T n method on simulated datasets, Li et al. (2020) recommend grouping the accepted quartz grains by applying a finite mixture model (FMM) to ratios of L x /T x values from each grain's fitted DRC.For example, the ratio of the L x /T x value at 400 Gy divided by the L x /T x value at 50 Gy.A central L n /T n value, a standardised DRC, and a D e value can then be determined for each grouping.
We applied the L n /T n method by grouping the accepted grains of each sample based on FMM results for the ratio of L x /T x values at 506 Gy and 32 Gy.This ratio will be close to unity for an early saturating grain and significantly larger for grains with high UL doses (see Fig. S4 for a plot of UL dose vs this L x /T x ratio for sample ZS17-05).Starting with two components, we repeatedly ran the FMM with an increasing number of components until the Bayes Information Criterion was no longer decreasing or one of the components was repeated.This process divided the accepted grains for each of the samples into 5-7 groupings per sample containing 2-143 grains (Table S1).L n /T n values and DRCs of the grains in each grouping were scaled and re-normalised using the LSnormalisation procedure (Li et al., 2017).Equivalent dose values were calculated for all groupings with at least 25 grains (Table S2) using the group average re-normalised L n /T n value, as calculated using the ADM with a σ m value of 10%, interpolated onto the group-specific re-normalised DRC.
The large number of components calculated with the FMM for each dataset reflects the large spectrum of UL doses of the grains and meant that each grouping had a relatively small number of included grains.Multiple D e values were calculated for each sample (one for each grouping) and are plotted in green on Fig. 3 and S2 against the UL dose of the standardised DRC of each component (top axis label).As the groupings were constructed using L x /T x ratios rather than UL doses and were sorted using the FMM, the UL doses of grains comprising each grouping are more complex than a simple roving window, but unlike the other D e values plotted on these figures, the L n /T n method component D e values include grains with UL dose values lower than where they are plotted and do not include all grains with higher UL dose values.
Equivalent doses calculated with the L n /T n method agree with other D e values calculated for the samples in this study but tend to be lower than other methods accounting for saturated grains (Fig. 3).In some cases, the scatter between D e values calculated for different groupings of the same sample can be significant and could reflect that the groupings do not each contain a sufficient number of grains to fully describe the complex equivalent dose distributions of these samples.However, for some of the datasets, including the NP18-05 dose recovery dataset (0.95 ± 0.04), D e values for L n /T n method groupings with at least 25 grains agree with one another.
Two of the samples with the greatest scatter in D e values for different L n /T n method groupings are ZS17-04 and ZS17-05 (green squares Fig. S2d and S2f).For sample ZS17-04, L n /T n method D e values increase with UL dose values of the standardised DRCs of each component (see Fig. S6 for a plot of D e value vs UL dose value for grains in each sample).Whereas for sample ZS17-05, a significant proportion of the scatter in D e values originates from the larger D e value calculated for its L n /T n method component 3. Abanico plots of the L n /T n distributions underlying the different groupings (Fig. S7) reveal that this component has a significant positive tail that shifts the ADM L n /T n average higher than would be calculated with a different average model such as CAM; however, component 3 is only the third most populous grouping for this sample and it's possible that with additional measurements the positive tail may become less prominent similar to components 1 and 2, which both have ~100 grains.

Ages and stratigraphy
The stratigraphic relationships between the samples collected from the aeolian sandscarp provides an initial test of the integrity of the luminescence ages.When plotted by depth below the modern surface, each of the different methods investigated produced age estimates in stratigraphic order (Fig. 6, Table S3).The ADM ages shown in pink on Fig. 6 were calculated using all accepted grains for which Analyst calculated D e values and is an example of D e values that ignore uncertainty saturated grains.L n /T n method results are shown as the weighted average of D e values calculated by components with at least 25 grains, and these values tend to be only slightly higher than the pink ADM ages.Ages calculated using inclusive plateau, synthetic super-aliquot, and Gaussian BayLum approaches all use data from the same grains and are generally in agreement suggesting that the inclusion of saturated grains in a D e calculation is more significant than the method used to include them.
Ages for the two uppermost samples agree with each other as well as the next sample ~3.5 m deeper (Fig. 6).This apparent age stagnation is likely related to the significantly higher overdispersion values of the two uppermost samples and provides an interesting opportunity to consider the ability of these different approaches to calculate minimum estimates.The synthetic super-aliquot and BayLum calculations have no clear mechanism for a minimum age approach.However, frequentist minimum age models could be applied to either the 'inclusive plateau' method or the L n /T n values of a grouping of the L n /T n method.The minimum age model (MAM) presented by Galbraith et al. (1999) calculates a central value from a truncated normal distribution based on a user-defined σ b value of the expected overdispersion for the sample.
Inclusive plateau MAM ages for the two uppermost samples were calculated using a value for σ b of 53% (the average overdispersion of the other samples at the NP site) and are shown on Fig. 6 as inverted light grey triangles.However, this model uses a geometric rather than an arithmetic mean for calculating the central value, which is used by the ADM.As the geometric mean itself will be lower than an arithmetic mean, using CAM ages (geometric mean) for the uppermost two samples also shifts their ages younger whilst not assuming a truncated distribution (shown on Fig. 6 as inverted dark grey triangles).
The significant scatter in ages for the deepest two samples, ZS17-06 and NP18-04, is what one might expect for samples at the limit of quartz OSL dating.Sample ZS17-06 has three inclusive plateau ages calculated due to a series of plateaus on an otherwise increasing trend, whilst no inclusive plateau age is calculated for NP18-04 as its ADM D e values only continued to increase with increasingly strict minimum UL dose thresholds without ever reaching a plateau.Synthetic super-aliquots for these two samples had natural signals of 67.6% (ZS17-06) and 91.3% (NP18-04) of I max .If we consider the synthetic super-aliquot ages to be accurate or potential minimum estimates, the ages calculated for the samples in this study suggest there may be a depositional hiatus or significant change in accumulation rate between samples NP18-05 and ZS17-06/ NP18-04 at ~11.6-12.9m below the modern surface.
Ages including data from saturated grains are suggested to be the more accurate ages for these samples based on the interpretation that the overdispersion and skew of the SG distributions primarily originates from beta dose rate heterogeneity.

Conclusions
Saturated grains, uncertainty saturated grains, and grains beyond their grain-specific UL dose can contain critical information for calculating more accurate D e values.There are many ways of accounting for these types of grains, some of which are discussed here.Applying a minimum UL dose criterion or calculating D e values with the L n /T n method can significantly reduce the number of grains included in a D e calculation and may require measurement of many grains to sufficiently represent complex single-grain distributions.Synthetic aliquots, inclusive plateau D e values, and BayLum palaeodose computations can incorporate data from saturated grains, but each method has weaknesses: Synthetic aliquots are sensitive to extreme D e values of very bright grains, multiple plateaus may be identified with an inclusive minimum UL dose analysis, and BayLum palaeodose computations don't provide output of individual grain DRC fits.Future research on sites with excellent stratigraphy and independent age control may be able to further test the accuracy of these different methods, but this is beyond the scope of this paper.Even still, good agreement between ages calculated using these different methods for this suite of samples with significant overdispersion and early OSL signal saturation suggests that including data from saturated grains has a more significant impact on D e calculation than which method is used to combine the data.
When dealing with equivalent dose datasets with saturated grains, it is important that such grains are not systematically rejected without consideration of the maximum dose that the grain can record.If

Fig. 1 .
Fig. 1.Photograph of the red sandy deposit where sample 232/NP18-07 was collected, which is representative of the remaining sample locations.

Fig. 2 .
Fig. 2. Histograms of UL doses for single grain OSL measurements of ZS7-04 for a) dose recovery measurements (given dose 41.4 Gy) and b) equivalent dose measurements.

Fig. 3 .
Fig.3.Equivalent doses calculated for some of the samples in this study using various methods and showing the influence of a minimum UL dose criterion.Note that while most of the points include all grains with D e values above the minimum threshold, the equivalent dose values for the L n /T n method (green) were calculated for groupings and refer to the top axis.The solid black reference lines in each plot depict potential minimum UL dose thresholds based on the average D e for the sample.Dashed reference lines in subfigure d) show unity and ±10% for the dose recovery test.

Fig. 4 .
Fig. 4. Synthetic DRCs for sample ZS17-06.Fine grey lines show synthetic aliquots summing all the grains on each of the 18 individual single grain discs measured.Datapoints and bold lines show synthetic super-aliquot DRCs summing grains from all 18 single grain discs measured after applying the criteria described in the legend.

Fig. 5 .
Fig. 5. Synthetic super aliquot D e values and ADM D e values for the same grains of sample ZS17-05 removing the grains with the brightest T n signals one at a time.

Fig. 6 .
Fig. 6.Ages calculated using the various approaches applied in this study plotted against depth.Filled points are analyses including saturated grains whilst open points were calculated only using grains with Analyst D e values.

Table 1
SAR Protocol used for Equivalent Dose and Dose Recovery Measurements.

Table 2
Environmental dose rates.
Water content of 5 ± 2% used for all samples.Internal alpha dose rate of 0.026 ± 0.010 Gy/ka used for all samples.aDepthsare below modern surface, half this value was used for cosmic dose rate calculation.b180-212 μm grain size used for all samples.

Table 3
Proportions of grains with burial doses at or beyond their respective OSL dating ranges.

Table 4
Dose Recovery test results calculated using all accepted grains (N) that gave discrete values in Analyst.