Counter-intuitive influence of Himalayan river morphodynamics on Indus Civilisation urban settlements

Urbanism in the Bronze-age Indus Civilisation (~4.6–3.9 thousand years before the present, ka) has been linked to water resources provided by large Himalayan river systems, although the largest concentrations of urban-scale Indus settlements are located far from extant Himalayan rivers. Here we analyse the sedimentary architecture, chronology and provenance of a major palaeochannel associated with many of these settlements. We show that the palaeochannel is a former course of the Sutlej River, the third largest of the present-day Himalayan rivers. Using optically stimulated luminescence dating of sand grains, we demonstrate that flow of the Sutlej in this course terminated considerably earlier than Indus occupation, with diversion to its present course complete shortly after ~8 ka. Indus urban settlements thus developed along an abandoned river valley rather than an active Himalayan river. Confinement of the Sutlej to its present incised course after ~8 ka likely reduced its propensity to re-route frequently thus enabling long-term stability for Indus settlements sited along the relict palaeochannel.

depending on the lithology, the whole drilling assembly was taken out and the drill bit was removed to pull the PVC pipe out of the drilling core tube barrel. The PVC pipe containing the sediment core was capped at both ends and labeled. The drill cuttings in the drill core bit were used for preparing field logs. Complete sediment cores are stored in the Indian National Core Archive facility at the Indian Institute of Technology -Kanpur, Kanpur (UP), India.
For sediment analysis, drill core pipes were split into two halves in a dark room preserving one half of the pipe for sampling for optically stimulated luminescence dating while the other half was used for detailed sedimentological logging and sample collection for analysis for dating and provenance. Colour, grain size, texture, mineral assemblages, inclusions such as nodules, rootlets, shells, carbonaceous material, and features like color mottling were documented. A qualitative estimate of grain size was made by following the grade scale proposed by Wentworth (ref. 5). The modal grain size was estimated by means of visual comparison to a reference card 6 . Texture is a combination of relative grain size, compactness and presence of any structures like lamination or lenses, or any inclusions like nodules, rhizoconcretions, shells, burrows. The core sediments were classified into sedimentary facies that are summarised in Supplementary Table 2. Three major sedimentary facies associations are identified: channel sand, floodplain fine-grained sediments and aeolian sand. Channel facies were identified mainly on the basis of sand grade particle size, poor to moderate sorting, thickness of sedimentary succession, and upward fining. Floodplain facies comprise finer-grained sediments, mostly silt and mud, in varying proportions along with iron and/or calcrete nodules, mottling and rhizoconcretions. Aeolian facies comprise yellowbrown, fine-grained, well-sorted sands (Supplementary Figure 1).

Optically Stimulated Luminescence Dating
Optically stimulated luminescence (OSL) dating determines the time elapsed since the last exposure of sediment grains to daylight (the burial age); this light exposure resets any prior latent OSL signal. During burial, mineral grains are exposed to environmental ionising radiation resulting in the creation of free charge. Some part of this free charge is trapped at defects and the net trapped charge is related to the total radiation exposure. In the case of quartz and feldspar, optical stimulation can release this trapped charge resulting in the emission of luminescence called OSL. The OSL signal can be calibrated in terms of absorbed dose (the equivalent dose, D e ; expressed in Gy). The burial age is then determined by dividing the D e by the independently measured environmental dose rate (expressed in Gy ka -1 ). We dated 52 samples extracted from seven individual sediment cores using the infra-red stimulated luminescence (IRSL) signal at 50˚C (henceforth IR 50 ) from multi-grain K-feldspar aliquots (2 mm) as well as the blue/green OSL signals from multi-grain (8 mm) and singlegrain quartz aliquots.
Fifty-three OSL samples (one modern sample from a Sutlej River sand bar and 52 samples from the seven cores) were prepared under subdued orange light conditions in the Nordic Laboratory for Luminescence Dating, DTU Risø Campus. The potentially light-exposed material from the ends of each selected core section was reserved for radionuclide concentration measurements, and the inner portion used for water content and OSL measurements. For the majority of samples (coarse to medium sand) the grain size fraction 180-250 µm was first separated by wet sieving and then treated with HCl (10%), H2O2 (10%), and HF (10%). This quartz-and feldspar-rich mixture was then density separated using aqueous solutions of 'Fastfloat' to give quartz-rich (>2.62 g cm -3 ) and K feldspar-rich (<2.58 g cm -3 ) extracts. The quartz-rich extract was again etched in HF (40%) to minimize any residual feldspar contamination, before re-sieving to >180 µm. In those cases where the samples were dominated by fine sand, the 90-180 µm fraction was processed in the same manner, but with final sieving to >90 µm (see Supplementary Table 7). Finally, the luminescence purity of the resulting quartz extract was examined by testing the IRSL sensitivity (a measure of feldspar contamination). Thirteen of the 38 quartz samples were found to have a significant IRSL signal which could not be removed by repeated HF treatment; nine of these samples were then measured using pulsed optical stimulation (see below).
Automated Risø TL/OSL readers (TL-DA-20, 7 ) were used for all OSL measurements. Optical stimulation of multi-grain aliquots used infra-red (IR, 870±40 nm; ~130 mW cm -2 ) and blue (470±30 nm; ~80 mW cm -2 ) light emitting diodes (LEDs) and pulsed optical stimulation (POSL) used an integrated pulsing option to control the LEDs 8 . All POSL curves were measured using an on-time (pulse width) and off-time (time elapsed between two consecutive pulses) of 50 µs; this setting has been shown to give the best separation of quartz and feldspar OSL signals 9 . Data were only acquired during the off-time, collection began after a delay of 5 µs after the end of the light pulse, and stopped 0.7 µs before the beginning of the next pulse. OSL signals were detected using EMI 9635QA photomultipliers in combination with various detection filters. Feldspar IRSL signals were detected using a filter pack consisting of 2 mm Schott BG-39 in combination with 4 mm Corning 7-59. Quartz OSL signals were detected through 7.5 mm Hoya U-340 filters. In situ irradiations were made using calibrated 90 Sr/ 90 Y beta sources providing dose rates ranging between ~0.1 and ~0.3 Gy s -1 to quartz mounted on stainless steel discs.
Automated Risø readers fitted with single grain laser attachments 10 were used for single-grain OSL measurements. The stimulation light source was a 10 mW Nd:YVO 4 solid-state diodepumped laser emitting at 532 nm, which was focused sequentially on a 10 by 10 array of grain holes spaced on a 600 µm square grid in special 10 mm diameter aluminium sample discs. The diameters of the single-grain holes were either 200 or 300 µm depending on the grain-size fraction measured. All single-grain discs were screened for contamination prior to measurement. Correcting the single grain data for beta source non-homogeneity using the approach described by Lapp et al. (ref: 11) did not result in significant changes to dose or dose distribution.

OSL measurements
The single-aliquot regenerative-dose (SAR) procedure 12 was used for equivalent dose determination. For K-rich feldspar measurements, optical stimulation for 200 s used IR LEDs with the sample held at 50˚C. The preheat and cutheat temperature were both 250˚C for 60 s 4 (ref. 13). A high-temperature IR bleach at 290˚C for 40 s was inserted in between the various SAR cycles to minimize recuperation. Initial pIRIR 290 measurements 14 were made using a preheat and cutheat of 320˚C for 60 s and a prior IR stimulation temperature of 50˚C. For quartz measurements, a double SAR procedure employing IR stimulation prior to blue or green stimulation (in both cases, sample at 125˚C, stimulation for 40 s) was used on all samples. A preheat temperature of 260˚C, 10 s and a cutheat temperature of 220˚C was used. A high-temperature blue bleach at 290˚C for 40 s was inserted in between each SAR cycle to minimize recuperation effects 15 .
The IRSL signals from 2 mm diameter aliquots of K-feldspar were summed over the first 2 s of stimulation and corrected for background using the last 20 s of stimulation (i.e., late background subtraction). The OSL signals from multi-grain quartz aliquots (8 mm) were summed over the first 0.64 s of stimulation and corrected for background using the signal detected between 1.12 and 1.76 s of stimulation (i.e., early background subtraction EBG 16 ). The effect of using EBG on quartz multi-grain aliquot data for these samples is shown in section 'Multi-grain quartz luminescence characteristics' below. The OSL signals from single grain measurements were collected during 1 s stimulations, and summed over the first 60 ms of stimulation and corrected for background using the last 0.25 s of stimulation.
The SAR dose response curves were based on a minimum of three regenerative dose points, a recuperation point and a recycling point. The doses reported were estimated by interpolation onto sensitivity-corrected dose-response curves fitted with a saturating exponential function where " # is the sensitivity-corrected OSL response to a given dose D (commonly denoted as the ratio L x /T x for laboratory dose and L n /T n for natural), & # is the saturation value and D 0 is a constant. Based on this interpolation model, the laboratory dose for which the response " # equals the corrected natural OSL signal , # is the so-called equivalent dose, D e . The uncertainties on D e estimated using individual aliquots or grains were assigned using "Analyst 4.00" (ref. 17); these consist of error propagation from counting statistics and curve fitting. For a significant number of quartz aliquots no bounded D e estimates could be derived.
Here we define a bounded dose estimate to be an estimate of D e for which . + " (2) was finite, where " is the standard error. Thus, an unbounded dose estimate is one where the natural sensitivity-corrected signal , # + 1 2 3 where 1 2 3 is the standard error of the natural sensitivity-corrected signal. Equivalent dose estimates were accepted if the uncertainty on the natural test dose response (i.e. the first test dose signal) is less than 20% and if the recycling ratio is within 15% of unity (multi-grain aliquots) or 2σ of unity (single grains). For single-grain analysis, dose estimates were also rejected if the IR depletion ratio 18 was not within 2 σ of unity, and if the recuperation dose was larger than 1 Gy.

Dosimetry
Radionuclide concentrations were estimated using high resolution gamma spectrometry, with ground and homogenized ~200 g samples mixed with wax to retain radon and cast in a reproducible geometry 19 . The resulting 238 U, 226 Ra, 232 Th and 40 K activity concentrations (Supplementary Table 7) were converted to dose rates following Guerin et al. (ref. 20). For K-feldspar extracts the internal beta dose rates from 40 K and 87 Rb were calculated assuming an effective potassium concentration of 12.5±0.5% and a rubidium concentration of 400±100 ppm (ref. 21).
Parent 238 U is not precisely measured in our gamma spectrometry facility; nevertheless, the arithmetic mean of the 226 Ra to 238 U activity ratio is 1.12±0.03 and the weighted mean is 0.94±0.02. We conclude that there is no evidence for significant secular disequilibrium in our samples in the first half of the uranium series. Nevertheless, for dosimetry calculations we conservatively assumed 234 U and 230 Th concentrations midway between the observed values of 238 U and 226 Ra with uncertainties equal to those of the 238 U analyses. We also assumed a 20±10% loss of 222 Rn compared to its parent 226 Ra. Supplementary Table 7 lists the resulting infinite matrix dry beta and gamma dose rates. Had we instead assumed uranium series equilibrium down to and including 226 Ra and based our uranium series activities only on 226 Ra analyses, then these dose rates would be, on average <1% different.
The total dose rates given in Supplementary Table 7 include a water content correction 22 based on the presumed lifetime average water content shown in column 6. These lifetime averages were based on laboratory-measured saturation water contents (with an associated uncertainty of ±4%) because most of our samples have lain below the water table for the entire burial period; even the shallowest samples will have been saturated for the majority of the burial period. Finally, a contribution from cosmic rays was added to the total dose rate 23 .

K-Feldspar analysis
We initially studied the behavior of both the IRSL measured at 50˚C (IR 50 ) and that of the subsequent IRSL measured at 290˚C (i.e. pIRIR 290 ), both following a preheat at 320˚C for 60 s. It rapidly became clear that the pIRIR 290 signal was very poorly bleached, and thus we focus here on the IR 50 signal. Figure 9a shows a typical natural IRSL stimulation curve measured at 50˚C (IR 50 ) for sample 101906, together with the regenerated stimulation curve resulting from a regeneration dose (in this case 143 Gy), both following a preheat at 250˚C for 60 s. The difference in decay shape between the two signals is presumably caused by anomalous fading in the natural signal 24 . Supplementary  Figure 9a shows a typical dose response curve for the same sample for the IR 50 signal. The dose response curve has been fitted using a single saturating exponential function and has a D 0 value of 245±7 Gy. The recuperation is small and the recycling value is 1.002±0.001. The natural dose is estimated to be 114±1 Gy. Supplementary Figure 9b is a histogram of the average apparent sensitivity change taking place in the SAR measurement for all natural Krich feldspar samples. The average total sensitivity change is 0.81±0.03 (n=37). The inset to We also measured this modern sample using the pIRIR 290 protocol 14 and measured a dose of 48±20 Gy (n=6); corresponding to an age of 13±5 ka. This age off-set is comparable with the ages of the younger samples and indicates that the pIRIR 290 signal was not well-bleached at deposition; we do not consider this signal further.

Multi-grain Quartz analysis
Multi-grain quartz luminescence characteristics Stimulation curve shape: Supplementary  Figure 9e shows a typical quartz dose response curve from sample 101906. The dose response curve has been fitted using a single saturating exponential function and has a D 0 value of 43±2 Gy. The recuperation is small and the recycling value is 0.97±0.02. The natural dose is estimated to be 82±11 Gy. The inset shows the natural OSL stimulation curve as well as the stimulation curve measured for the highest regeneration dose (in this case 197 Gy). These are very similar in shape, but in some aliquots the regenerated stimulation curve was observed to decay more slowly than the natural. The use of a background signal taken very close to the signal summation interval helps to make quartz OSL signals insensitive to such changes in shape by preferentially isolating the fast OSL component. To test whether or not this change in stimulation shape affects our large multi-grain aliquot quartz dose estimates the equivalent doses were repeatedly calculated using an initial signal interval consisting of the first 0.32 s and a background 0.48 s wide. This background was moved progressively down the curve from immediately adjacent to the signal interval to the final 39.52 to 40 s for four samples (091901, 091906, 111901 and 111915). The average ratio of the dose based on the most stringent summation interval (first 0.32 s minus the subsequent 0.48 s) to that based on the standard summation interval (given in Supplementary Methods: OSL Measurements) is 1.03±0.03 (n=4). It is clear that the average equivalent dose is independent of the background interval chosen for these samples and that all results are consistent with those from our standard intervals. It is concluded that our equivalent dose estimates are insensitive to differences in decay shapes between the natural and regenerated signals.

Sensitivity change
The SAR procedure is widely used in quartz OSL dating because of its ability to monitor and correct for sensitivity changes using the test dose correction. The sensitivity change as a function of SAR cycle number for the aliquot displayed in Supplementary Figure 9e is shown as an inset in Supplementary Figure 9f. This particular aliquot changes its sensitivity by a factor of ~22 and the average sensitivity change for all accepted aliquots for sample 101906 is 13.5±1.3 (n =15). In Supplementary Figure 9f the average total sensitivity change for all individual samples (n=37) is shown as a histogram. The average sensitivity change is about 5 times with the largest sensitivity change occurring for sample 111906 (27±2, n =21) and the smallest for 111905 (1.28±0.05, n =28). The remeasurement of a point on the dose response curve (i.e. a recycling point) gives information on the reproducibility of sensitivity-corrected luminescence signals; on average this recycling ratio was 1.008±0.009 (n=37) for all samples. Thus, although the sensitivity varies significantly during the measurement of the laboratory dose response curve, our SAR protocol is able to correct successfully for these changes.
Preheat plateau To investigate the dependence of the equivalent dose on our choice of preheat temperature we measured preheat plateaus for samples 101906, 111902 and 111903. The cutheat temperature (before the T x cycle) was kept at 40˚C below the preheat temperature. Supplementary Figure 9g illustrates these measurements for sample 101906 and 111902. Each point is the average of at least four aliquots. It appears that our equivalent doses are insensitive to preheat temperatures £ 280˚C. Here we choose a standard preheat of 260˚C for 10 s and a cutheat of 220˚C.

Multi-grain quartz dose recovery
To investigate whether a known dose can be recovered using the SAR protocol prior to any thermal treatment we undertook dose recovery experiments, where the aliquots were bleached twice for 100 s using blue light stimulation at room temperature, separated by a pause of 10 ks; this pause is intended to allow charge transferred to the 110˚C TL peak to decay to negligible levels. The aliquots were then given a beta dose in the range 15 to 115 Gy before subsequent measurement. Dose recovery measurements (see Supplementary Figure 12a) were undertaken for a total of 140 aliquots of 24 samples and we obtain an average dose recovery ratio of 1.03±0.02 (n=140), suggesting that a laboratory dose given before the first heating of these samples can be measured accurately using our protocol. This is of particular relevance to these samples given the very large sensitivity changes (up to 27 times) that are observed in both dose recovery and equivalent dose measurements. Murray 36)) from young or modern quartz samples from fluvial and colluvial environments and find that the average residual dose is ~2 Gy. These authors further point out that their data should be considered as providing upper limits to likely residual doses at deposition, because they were taken from samples still undergoing transport; such samples are expected to undergo further light exposure before final deposition and eventual preservation. Thus, in such environments, incomplete bleaching is likely to be of concern only for young samples (e.g., <10 ka).

Quartz modern analogue
The samples measured in this study are all (with the exception of samples 091905 and 111909) extracted from a major fluvial system, between ~50 km (SRH-5a) and ~300 km (GS cores) from the Himalayan mountain front onto a low gradient alluvial plain. Based on this geomorphological setting, we do not anticipate incomplete bleaching of quartz in this environment to be a significant problem. To test this hypothesis, we measured the quartz fraction of the modern analogue (sample 121915) taken from where the Sutlej river leaves the Siwalik foothills and enters the floodplain. For the multi-grain quartz measurements, a dose of 0.9±0.2 Gy (n=24) was measured. This corresponds to an age of ~ 400 years. We also measured 2,400 single grains and obtained a dose distribution with an average dose of 0.4±0.3 Gy (n=90) and a weighted dose of 0.06±0.01 Gy, corresponding to an age of ~25 years. Thus, if we assume that this modern analogue is appropriate to the buried samples being dated here, we can conclude that incomplete bleaching is not of significant concern for our samples.
Natural multi-grain quartz dose distributions Natural dose distributions were measured using large (8 mm) multi-grain quartz aliquots and the resulting dose distributions have relative standard deviations ranging between 16% (111911) and 59% (111913) with an average value of 33%. All dose response curves have been fitted using a single saturating exponential function and D 0 (exponential coefficient) values have been obtained for all aliquots. The average D 0 value for multi-grain quartz aliquots is 65±4 Gy (n=1197). A However, we were not able to derive bounded dose estimates for a significant fraction of the aliquots (~17% on average, using CW stimulation) which otherwise passed the rejection criteria outlined in section OSL Measurements. In the literature 37,38 , it has been pointed out that if the L' n value is close to the saturation value L' µ , the uncertainties on the associated dose estimate D will be large and asymmetric 37,38 and Wintle and Murray (ref. 39) stated that it seems prudent to ensure that accepted dose estimates are smaller than 2×D 0 . Although, the application of this 2×D 0 criterion may be appropriate when applied on a sample-by-sample basis, it will inevitably produce a bias to lower doses if applied on an aliquot-by-aliquot basis within a sample 40 . Because of this we have included all results which appear to intercept the laboratory dose response curve 41 .
IR contamination As stated in section Sample Processing, 13 of the 38 quartz samples (including the modern sample 121915) showed a significant reduction in OSL signal after stimulation with IR, which is normally assumed to indicate feldspar contamination. All samples were measured using a double SAR protocol 42,43 where an IR bleach is inserted prior to blue stimulation to reduce the likely contribution from feldspar to the quartz dose. However, a blue-stimulated feldspar signal remained after IR stimulation and it has been shown that the simple IR bleach (here at 125˚C, 100 s) in the double SAR protocol is not always sufficient to result in a pure quartz dose 30 . However, the luminescence lifetime of feldspar during blue light stimulation is significantly shorter than that of quartz [44][45][46] so it is possible to achieve an instrumental separation of the two signals using pulsed (POSL) stimulation. We applied a double SAR POSL protocol to 11 of the multi-grain quartz samples (9 samples with poor IR depletion ratios (111904, 05, 07, [19][20][21][22][23][24]).
In Supplementary Figure 13 the average CW D 0 values for individual samples are shown as a function of average IR depletion ratio (no IR depletion ratio has been measured for sample 111901, but the IR to blue ratio indicated significant IR contamination). If the IR signal arises from feldspar it is to be expected that samples with IR depletion ratios significantly smaller than unity would have substantially higher D 0 values, because the blue stimulated signal from feldspar has much higher D 0 values than those from the quartz fast component. However, Supplementary Figure 13 does not show any apparent correlation between D 0 and IR depletion ratio. We can also compare the D 0 values obtained for CW and POSL stimulation. If the stimulation curves from these IR contaminated samples are significantly influenced by the presence of contaminating feldspar, one would expect that the D 0 values would be larger in CW mode than in POSL: the D 0 ratio between CW and POSL modes is 1.23±0.17 (n=9) indistinguishable from unity. Only two of the samples (111904 and -05) have D 0 ratios significantly greater than unity, e.g., 2.3±0.3 and 1.7±0.2, respectively. For these samples alone, it would appear that at high doses, where the quartz signal is in saturation, feldspar contamination may govern the shape of the dose response curve. For those samples where both methods were used, the ratio of the CW to POSL doses is 1.01±0.09 (n=8; data not shown; no POSL dose estimate could be derived from sample 111907, see below) indicating that any blue-stimulated feldspar contamination remaining after the IR bleach at 125°C is insufficient to affect the CW dose estimate significantly.
Although the average doses and D 0 values seems more or less unaffected by any feldspar contamination, the proportion of aliquots for which we are unable to provide bounded dose estimates increases significantly when we apply POSL stimulation. The average relative number of unbounded dose estimates for the nine samples in CW mode is 12±5% (n=9) and in POSL mode 45±11% (n=9), e.g. for sample 111907 we are not able to derive bounded dose estimates in CW for 33% of the measured aliquots (n=16/48) due to saturation whereas in POSL mode the corresponding number is 100% (n=7/7). This suggests that many of these samples are indeed at or close to saturation on the quartz growth curve, but that under CW stimulation the feldspar signal (unimportant at low doses) does allow the dose response curve to continue to increase at high doses; this results in an incorrect and apparently finite quartz dose estimate. We use the POSL data for the 9 IR contaminated samples. For the four remaining samples with significant IR signals (091904, 101902, 111901 and 111908) we did not pursue the use of the POSL mode since the average doses were indistinguishable between the POSL and CW-OSL mode. For these samples, we only used CW-OSL, and the relative number of aliquots rejected due to saturation is 18, 47, 0 and 33%, respectively.
Supplementary Figure 9h shows a histogram of all measured D 0 values (CW and POSL); the average D 0 is 59.9±1.0 Gy (n=1161), although the distribution is significantly positively skewed with the highest D 0 value being 426±168 Gy. These data are used below when discussing bounded and unbounded estimates of equivalent dose.
Dose saturation One of the major drawbacks of using quartz as a natural dosimeter for old samples is that a fast-component dominated signal saturates at relatively low doses, e.g. ~100 Gy (ref. 39). Of the 37 core samples measured using multi-grain quartz aliquots (8 mm), only 6 samples (101901, -03, 111901, -02, -10, -12) had no aliquots with the natural sensitivity corrected signal (L n /T n ) appearing to be in saturation. These samples are all relatively young with OSL ages of ~ 6 ka, with the exception of sample 111902 which has an age of ~15 ka. The fact that many natural L n /T n values lie significantly above the saturation value for the majority of samples is of concern and implies that the measured natural signal did not grow up the same dose-response curve as is generated in the laboratory. This casts doubt on the accuracy of all equivalent dose estimates derived using multi-grain quartz aliquots; all aliquots presumably contain some grains for which the laboratory dose response curve does not describe the growth of the natural signal. Certainly, if a substantial number of aliquots appears to be in or above saturation it is likely that the derived average equivalent dose should be regarded as a minimum dose estimate. The relative number of aliquots rejected due to saturation is given in Supplementary Table 8 for each sample. Here, we have arbitrarily chosen to regard all multi-grain average equivalent dose estimates as minimum dose estimates if the proportion of aliquots in saturation is larger than 15%. Thus, we regard the multi-grain quartz D e estimates from 20 of the 37 samples as minimum dose estimates.

Single grain quartz analysis
To test if quartz and feldspar grains were well-bleached at deposition, quartz single-grain analysis was undertaken for 30 samples. Single grain OSL analysis is appropriate if significant incomplete bleaching is observed or anticipated 47,48 ; it may also be justified for samples which have experienced post-depositional mixing 49 , or as a means of rejecting grains with poor luminescence characteristics 50,51 -i.e., grains which do not accurately record the burial dose. We undertook single grain analysis using 30 (out of 37) core samples to investigate the unusually high rejection rate of multi-grain analyses due to saturation, and to further test whether incomplete bleaching may contribute to the scatter in inter-and intrasample dose and age estimates.
Single grain luminescence characteristics All single grain analyses were measured using a standard preheat of 260˚C for 10 s and a cutheat of 220˚C as these thermal treatments were determined to be appropriate using multi-grain aliquots. The luminescence characteristics of single grains of quartz have consistently been reported to be highly variable between different grains both in terms of intrinsic brightness, stimulation curve and dose response curve shape for the same measurement conditions 47,52,53 , and some selection criteria are required to identify grains suitable for further analysis.
Single grain rejection criteria We applied standard rejection criteria (see OSL Measurements) in a first attempt to obtain the most reproducible dose distributions. The effect of applying these rejection criteria to 30 single grain natural dose distributions on the normalized Central Age Model (CAM) 54 dose, the normalized over-dispersion and the normalized number of accepted grains are shown in Supplementary Figure 14a, b and c, respectively. In general, the application of these rejection criteria does not affect either the central dose or the observed scatter (quantified by the relative over-dispersion, OD) significantly; this conclusion applies both to natural dose distributions and to dose recovery dose distributions. Sample 111912 is the single exception to this statement. In this sample, although the central dose does not change significantly, the OD changes from 70±7% (s TN <20%) to 32±6% (all rejection criteria). However, the ratio of ODs determined with and without rejection criteria for all other samples are not statistically different from unity. A similar conclusion was reached by Thomsen  It would seem that the only significant effect of applying these rejection criteria is to significantly reduce the number of accepted dose estimates (by ~50%, on average). Nevertheless, for consistency with earlier published work, below we only consider the data which pass all the standard rejection criteria.
Intrinsic OSL brightness The natural OSL signals from 103,300 individual grains were measured from 30 samples. Of these grains, 3.6% gave both an acceptable test dose response (i.e. s TN <20%) and a natural sensitivity-corrected signal which intercepted the laboratory regenerated dose response curve (i.e. gave bounded dose estimates). After application of the remaining single grain rejection criteria, 1.9±0.4% (n=30 samples) of the measured grains were accepted into the respective natural dose distributions. Supplementary Figure 15 shows the cumulative light sum 57 for all grains from one sample from each of the six cores; 80% of the light sum is derived from 5 to 20% of the grains and ~40% of the grains do not contribute significantly to the light sum. In single grain analysis it is common for the sensitivity of individual grains to vary by over four orders of magnitude 58 but for these samples the average variation in sensitivity is only about two orders of magnitude with the majority of grains being only weakly luminescent. The inset to Supplementary Figure 15 shows the net OSL signal per Gy from the natural test dose signal for six representative samples for the grains accepted according to the selection criteria given in section OSL Measurements; the grains have been ranked from brightest to dimmest 58 . For comparison, the results from a Namibian fluvial sample (031304 (ref. 59)) and a Danish coastal deposit 55 are also given. The quartz samples from our study are in comparison significantly dimmer than published material. For the samples in this study the median of the natural test dose OSL response per grain is 1.78±0.09 Gy -1 (summed over the initial 60 ms of stimulation; 1966 grains, 30 samples). For comparison, the corresponding numbers for the fluvial sample (031304) from Namibia is 94 Gy -1 (525 grains) and for the recent Danish coastal deposit is 78 Gy -1 (358 grains).
Single grain OSL stimulation curve shape Multi-grain quartz measurements have shown that the quartz OSL signal consists of different components with differing photoionisation cross sections and thermal stabilities 51,60 . It is generally accepted that only samples dominated by the fast component are likely to give accurate estimates of burial dose 61,62 . As a consequence, many laboratories routinely use an early background subtraction (EBG) for multi-grain quartz measurements to minimize the contribution from slower decaying and potentially less bleached components. Supplementary Figure 16 shows five representative, normalized OSL stimulation curves for the samples measured in this study; there is a large variability in curve shape. Although such large variability is commonly observed, only a few single grain OSL studies have examined the impact of including dose estimates derived from grains with varying decay rates on the final dose distribution. Ballerini et al. (ref. 63) concluded that single grain dose estimates can be over-estimated if doses derived from slowly decaying OSL signals are accepted into the dose distributions and suggested that EBG ought to be applied. In a single grain quartz study by Reimann et al. (ref. 64) it was observed that using EBG instead of LBG (late background subtraction) apparently reduced the relative OD (from 34±2% to 30±3%) because the use of EBG preferentially leads to the rejection of low sensitivity quartz grains, but that this had no significant effect on the average burial dose estimate. This is in agreement with the observation by Thomsen et al. (ref. 55) who reported that doses from weakly luminescent single grains were significantly more over-dispersed than those from bright grains. The EBG approach was used by a number of authors [65][66][67] , but all concluded that for their samples no significant benefits could be observed. Madsen and Murray (ref. 68) and Durcan and Duller (ref. 69) suggested using the fast ratio F R to classify the decay shapes of individual stimulation curves. F R is given by (F-BG)/(M-BG), where F is the initial summation, M is a subsequent summation and BG is the late background. The fast ratio F R was applied by Duller (ref. 70) in his single grain dose recovery experiments and he observed a significant improvement in laboratory relative OD although the dose recovery ratios remained underestimated for given doses approximately > 100 Gy. In contrast, Feathers (ref. 71) observed both an improved dose recovery and a reduction in OD in laboratory and natural dose distributions. Thomsen et al. (ref. 55) also applied the F R approach to their single grain dose recovery data but only observed a small reduction in relative OD for some of their samples. Jacobs et al. (ref. 65) applied the F R approach to their natural samples and found no significant effects on either burial dose or relative OD.
We calculated the F R for six of our single grain dose recovery data sets (e.g., 111904, -10, values are all positively skewed with an average mode of 1.8±0.1. To assess if the rate of decay influences the estimated dose and the relative OD, we have calculated the average (arithmetic) dose, the CAM dose and the relative OD for grains with F R >3 and for all grains irrespective of F R value, respectively, and then derived the ratio of these values calculated from the two data sets. (On average only ~20% of the accepted grains have F R >3.) The weighted mean of this ratio is 1.04±0.04 (n=6) for the average dose and 1.07±0.04 (n=6) for the CAM dose. Thus, we do not see a significant change in estimated dose by applying the F R criterion. For the ratio of the relative OD, the weighted mean is 0.92±0.11 (n=5), and all ratios (except for 111915) are consistent with unity within 2 standard errors. Thus, rejecting grains based on the F R values does not seem to introduce a significant change in relative OD. However, for sample 111915, which has a relative OD of 26±4% (n=72) when all grains irrespective of F R value is included in the calculation, we observed that by rejecting grains with a F R value less than 3, the data set becomes under-dispersed. This does lend some support to the observation made by Duller (ref. 70), although it should be noted that the F R <3 data set only contains 9 dose points, and so the apparent reduction in the unexplained scatter may be simply a statistical fluctuation.
We also calculated the F R for 25 of our natural data sets. The weighted ratios for the arithmetic mean and the CAM are 1.10±0.04 (n=25, 96% consistent with unity within 2 standard errors) and 1.08±0.04 (n=25, 88% consistent with unity within 2 standard errors), respectively. Thus, if there is an effect on the estimated dose it is small. The weighted ratio for the relative OD is 0.79±0.07 (n=25, 96% consistent with unity within 2 standard errors), which would imply that we observe a significant reduction in OD. However, this ratio includes sample 111921, which has a relative OD of 55±11% (n=34), when all grains are included and a relative OD of 15±9% (n=11), if we only include those grains with a F R >3. If we were to eliminate this sample from the data set then the weighted ratio is 0.87±0.07 (n=24), which is again consistent with unity. Thus, it would appear that the decay rate of individual single grains does not seem to significantly affect either the estimated dose or the relative OD.
Single-grain OSL dose response curves In the literature it has been reported 53,70 that the shapes of the dose-response curves of individual grains are highly variable, e.g., some dose response curves saturate at high doses (high D 0 values) whereas others saturate at low doses (low D 0 values). In this study, the D 0 values of the dose response curves (all fitted with a saturating exponential function) also vary widely from grain to grain. Supplementary Figure  16b shows three measured dose response curves for three of the grains whose stimulation curves are shown in Supplementary Figure 16a. In contrast to Jacobs et al. (ref. 65), who stated that there was some support for the suggestion that grains with low D 0 values also had the slowest rates of OSL decay (e.g., low F R ratios), we see no correlation between D 0 and F R for these samples. Supplementary Figure 17b is a frequency histogram of the D 0 values known to better than 30% of the accepted grains from sample 111909. The average weighted D 0 value for this sample is 112±5 (n=229) Gy with a median value of 109 Gy. The average weighted D 0 for all samples is 94±6 Gy (n=30) which is significantly higher than that observed for multi-grain aliquots (60±1 Gy). For all samples, we were unable to derive bounded dose estimates for ~30% (880) of the grains which otherwise pass all standard rejection criteria. The number of grains rejected due to saturation varies between 4% (111912) and 63% (091905); the fraction of grains in saturation generally increases with depth (age), as might be expected. It is interesting to note that the D 0 values of grains rejected due to saturation are, on average, only 52±3% (n=30) of the D 0 values of the accepted grains. A similar observation was made by Thomsen et al. (ref. 55) in their dose recovery experiments. However, the average L n /T n ratio of the grains rejected due to saturation is 2.3±0.1 (n=30) times those of the accepted grains. Thus, grains rejected because of saturation have both smaller D 0 values and larger L n /T n ratios than accepted grains. Both these characteristics will individually lead to unbounded dose estimates. Again, we observe no correlation between D 0 value for these saturated grains and F R value (data not shown).
Single grain sensitivity change As found using multi-grain aliquots, the accepted single grains also showed significant sensitivity changes throughout the SAR measurement sequence. The average sensitivity change for all samples is 4.4±0.6 (n=30 samples). As in the multi-grain data, the largest average sensitivity change (11.5±1.3, n=59 grains) is observed for sample 111906 and the smallest (1.2±0.2, n=18 grains) for sample 111905. The CAM average of the sensitivity change observed in multi-grain measurements to that observed using single grain measurements is 1.33±0.15 (n=30 samples); in addition, there is no indication, in terms of sensitivity change, that rejected single grains behave differently from accepted grains. In contrast to Gliganic et al. (ref. 72), we see no correlation between the size of the sensitivity change and the rate of OSL decay (data not shown).

Single grain thermal transfer and recuperation
To assess the importance of thermal transfer and recuperation for the chosen thermal treatment (preheat of 260˚C, 10 s and cutheat of 220˚C) portions of samples 091903, -04 and 111904 were bleached for 1 hour at a distance of 80 cm in a daylight simulator (Hönle SOL2). A total of 4200 grains were measured and 39 grains were accepted into the dose distribution (data not shown), to give an average dose of 2.2±1.0 Gy (n=39), a CAM dose (assuming a normal distribution) of 0.6±0.2 Gy, and an over-dispersion of 0.6±0.3 Gy. Thus, thermal transfer and recuperation do not appear to be a significant problem for these samples for the chosen thermal treatment.  Figure 12b shows the dose recovery ratio as a function of given dose using the standard rejection criteria (filled symbols). We appear to be able to measure accurately a dose given in the laboratory prior to any thermal treatment. However, the variability observed in the individual dose distributions is significantly larger than expected from the assigned uncertainties (see section OSL Measurements), i.e. the average relative OD is 32±2% (n=10), significantly larger than the typical OD (~10%) reported in beta dose recovery experiments. In a series of dose recovery experiments performed on a recent coastal deposit from Denmark, Thomsen et al. (ref. 55) showed that the observed OD in their single grain dose distributions depends on the intrinsic brightness; i.e. the OD calculated from bright grains is smaller than that calculated from dim grains. Given the low inherent OSL sensitivity of the quartz samples investigated here, it is possible that at least a part of this additional variability arises from this phenomenon. Regardless of the origins of this OD, these dose recovery results indicate that although we can, on average, measure a given dose accurately, the precision is low; we cannot know an individual dose estimate to better than ~30%, regardless of the calculated uncertainty associated with the individual dose. Thus, we have added (in quadrature) a 30% uncertainty to that derived from counting statistics and fitting (unless otherwise stated) to all individual dose estimates before any subsequent statistical analyses.

Single grain dose recovery experiments
Single grain natural dose distributions Natural quartz single grain dose distributions were measured from 30 of the 37 samples discussed here. Individual D e estimates were obtained for between 1200 and 7200 grains from each of the 30 samples. The number of accepted individual dose estimate ranges between 18 and 229. Relevant statistical information is provided in Supplementary Table 9. Supplementary Figure 18 shows representative single grain dose distributions for four of the samples (111910, 091906, 111915 and 111918). On the left the natural test dose signal is shown as a function of the dose estimated from individual grains and the insets show the same data as a simple frequency histogram. Radial plots of the same data are shown on the right, with individual data points assigned a minimum uncertainty of 30% (in addition to the uncertainty arising from counting statistics and curve fitting errors). (Note that Supplementary Figure 17a shows an additional dose distribution for sample 111909; there the natural test dose signal is expressed as a function of dose.) All dose distributions include a wide range of individual dose estimates, e.g. the dose distribution for sample 111918 contains individual dose estimates ranging between -0.2±1.0 Gy and 208±65 Gy and for sample 111909 they range between -8±10 Gy and 560±100 Gy. These two samples were extracted from the cores at 5.85 m (GS-14) and 36.89 m (GS-10), respectively and we confidently reject the possibility that the zero dose estimates for these samples arise from young intruding grains. Others 73 have also reported the existence of such zero dose grains in old samples; they were also unable to explain these apparently very young grains in terms of their depositional context. However, as we have screened the individual grain holes for contamination (stuck grains) before use, we are confident that any zero dose grains that we record have not arisen from contaminated discs. It would therefore be incorrect to simply reject these grains, because of the inherent risk of incorrectly biasing the dose distributions towards higher doses. In the absence of a credible extrinsic explanation for these doses (e.g. post-depositional mixing) we must accept that such extreme outliers simply indicate the degree of scatter in our measurements -that is, they reflect the intrinsic OD. Nine of the 30 single grain dose distributions contain such negative dose estimates (111912, 101901, 111903, -09, 101903, -04, 091901, -04 and 111918).
All natural single grain quartz distributions have been obtained using a double SAR protocol with an IR diode bleach at 125˚C for 40 s prior to green laser stimulation. In addition to this, only grains with an acceptable IR depletion ratio have been accepted. Thus, we are confident that our dose distributions are not significantly affected by feldspar contamination. Nevertheless, the existence of grains with very low doses extracted from old samples might indicate signal instability; although the fast component OSL signal from multi-grain quartz aliquots has been shown to be stable (>10 8 years), little information concerning stability on a grain by grain basis is available. To investigate the short-term stability, we undertook simple fading measurement on 2,400 grains from sample 111909, for which we previously had determined the natural dose. The fading measurement sequence consisted of three prompt L x /T x measurements separated by two delays of 3 days duration. In this experiment, the ratio of the averaged delayed L x /T x to the average prompt L x /T x measurements was 0.97±0.02 (n=83) and we did not observe any correlation between this ratio and the measured dose. This confirms that at least this quartz sample is unlikely to suffer from anomalous fading, or other forms of short-term instability.
The relative ODs of the natural dose distributions (calculated with individual uncertainties assigned on the basis of counting statistics and curve fitting errors alone) range between 32±6% (n=73, sample 111912) and 105±19% (n=61, sample 101904), with an unweighted average of 64±3% (weighted 61±3% with an OD of 20±6%) for all samples. Only two samples (111910 and -12 from GS-7) have ODs consistent with the OD reported for the single grain beta dose recovery experiments. There appears to be no correlation between OD and dose.

U-Pb geochronologic analysis of detrital zircons
Zircon U-Pb analyses were performed by LA-ICPMS using a New Wave Nd:YAG 213 nm laser ablation system, coupled to an Agilent 7500 quadrupole mass spectrometer. Real-time data were processed using GLITTER v4.4 data reduction software (www.glitter-gemoc.com). Repeated measurements of the zircon Plesovice standard (TIMS reference age 337.13 ± 0.37 Ma (ref. 74) and NIST 612 silicate glass 75 were used to correct for instrumental mass bias and depth-dependent inter-element fractionation of Pb, Th and U. Reported ages are based on the 206/207 ratio for grains older than 1.1 Ga and the 206/238 ratio for younger zircons. Grains with discordance >10% were omitted. Examination of discordant grain ratios showed that in most cases the discordance was due to mixing of multiple growth stages, rather than simple lead loss. Details of locations of modern river and dune samples are given in Supplementary Table 3, and locations of cores together with sample depth and lithology are given in Supplementary Tables 1 and 4. U-Pb isotope data for detrital zircons are presented in Supplementary Database 1.

Ar/ 39 Ar geochronologic analysis of detrital mica
The location and details of sample depths in cores are provided in Supplementary Tables 1, 3, 5. Samples were subjected to magnetic separation. Mica was handpicked under a binocular microscope to ensure a pure separate. Approximately 150 crystals per sample were picked. After rinsing in de-ionised water and methanol, the grains were parceled into Cu packets and positioned within an Al holder for irradiation with International 40 Ar/ 39 Ar Age Standard Fish Canyon sanidine (FCs, 28.201 ± 0.023 Ma (ref. 76)). Samples were irradiated in two batches. The first batch of samples were irradiated for 10 hours in the Cd-lined facility at the McMaster Nuclear Reactor, Ontario, Canada. The second set were irradiated in the Cd-lined CLICIT facility TRIGA reactor at Oregon State University, USA. Monitors were analyzed by total fusion with a focused CO 2 laser.
Single grains (200-400 µm) of unknowns (i.e., muscovite) were loaded into a Cu planchette in an ultra-high vacuum laser cell with a doubly pumped ZnSe window. Using a CO 2 laser, the muscovite crystals were fused. All gas fractions were subjected to 180 s of purification with two SAES GP50 getters (one at room temperature the other at 450˚C) and a cold finger maintained at -95.5˚C using a mixture of dry ice (CO 2[s] ) and acetone. Argon isotope ratios (i.e., ion beam intensities) were measured using a MAP 215-50 mass spectrometer in peak jumping mode. The mass spectrometer has a measured sensitivity of 1.13 × 10 -13 moles volt -1 . The extraction and cleanup line was automated, as were the mass spectrometer peak jumping routines and data acquisition. Blanks (full extraction line and mass spectrometer) were measured after every two analyses of unknowns. The average blank ± standard deviation of n = 812 from the entire blank run sequence was used to correct raw isotope measurements of unknowns ( 40 Ar 1.02 × 10 -15 moles, 39 Ar 3.10 × 10 -17 moles, 38 Ar 1.90 × 10 -17 moles, 37 Ar 7.85 × 10 -17 moles, 36 Ar 1.38 × 10 -17 moles). Mass discrimination was monitored by analysis of air pipette aliquots after every ten sample analyses (n = 174, 7.21 × 10 -14 moles 40 Ar, 40 Ar/ 36 Ar = 289.61 ± 0.57). The Ar isotope data were corrected for backgrounds, mass discrimination, and reactor-produced nuclides and processed using standard data reduction protocols. The decay constants of Steiger and Jager (ref. 77) and atmospheric argon ratios of Nier (ref. 78) were employed.
The BGC software MassSpec was used for data regression. Raw Ar isotope data are presented in Supplementary Database 2. Table 8 shows all fading corrected ages obtained using the IR 50 signal; these are shown on stratigraphic panels in Figures 5, 7 and 8, and shown as a function of depth in Supplementary Figure 11a for cores GS-10 and GS-11. The relative uncertainty (1 standard error) on these ages ranges between 4 and 13% with an average of 5.0±0.3% (n=52). The feldspar ages are generally in stratigraphic order with the exception of two young samples (111901/101901 and 111912). In GS-10 we appear to have a small but significant age inversion at the very top of the core; sample 111901 at ~2.3 m gives an age of 6.3±0.3 ka whereas sample 101901 at ~4.3 m gives an age of 4.9±0.2 ka. In core GS-7 (where only 3 fine-grained samples were measured) sample 111912 at ~5.1 m has an age of 4.0±0.2 ka and may be younger than the two overlying samples 4.5±0.2 ka and 4.7±0.2 ka). In core SRH-5 we have measured nine samples taken between ~2.4 and 44 m. The ages are all in stratigraphic order down to the bottom of the core which is dated to 67±3 ka.

K-feldspar ages Supplementary
In core GS-10 a total of 11 samples have been measured. These samples were taken between ~2.3 and 37 m (see Supplementary Figure 11a). Apart from the age inversion observed in the very top of the sequence (see above) the feldspar ages increase monotonically until a depth of ~14 m, below which the ages are all consistent with an average of 66±2 ka (n=6) until beyond 32 m. The deepest sample (111909) at 37 m, of aeolian origin, gives an age of 150±6 ka. The age of the uppermost medium sand sample in GS-10 is 23±2 ka.
In core GS-11 we have measured a total of 12 samples taken from depths between ~ 0.5 m and 37 m (see Supplementary Figure 11b). The ages increase monotonically with depth and as in core GS-10 there is evidence for a rapid sedimentation period at ~65 ka (mean age of samples between ~17 and ~32 m is 64±2 ka; n=3). The deepest sample (091905) taken at ~37 m from an aeolian deposit is 152±8 ka, which is consistent with the age obtained from the aeolian horizon in core GS-10. The age of the uppermost medium sand sample is 23.7±1.0 ka.
We have measured three samples from core GS-13 and five samples from core GS-14. All the ages are in stratigraphic order and the ages of the uppermost medium sand samples are 25.4±1.0 ka and 23.0±1.1 ka, respectively. The average age for all available uppermost fine/medium sand samples from all cores on the GS-transect is 23.7±1.0 ka (n=4).
Nine samples were measured for core MNK-6. A remarkable hiatus of ~55 ka is detected between 14.5 and 17.5 m (Figure 3, Supplementary Table 8). Below this break, the ages increase with depth and range from ~65 ka at 17.5 m to ~95 ka at the bottom of the core (35.2 m). The top of the core is characterised by a period of rapid deposition at 8.7±0.3 ka (n=4). Table 8 summarizes the quartz multi-grain equivalent doses and resulting ages for all 38 samples. The uncertainty (1 se) on individual ages range between 5 and 25% with an average uncertainty of 9.0±0.6% (n=37). As discussed above, the majority of these samples give minimum ages; only 17 of the 37 core samples are regarded as giving bounded ages. In core SRH-5A, only sample 111923 (22.6±1.5 ka, n = 37) is regarded as a bounded age. However, all ages from a depth of ~3 m are consistent with an average age of 20.6±0.7 ka (n=5). The youngest age (not to be confused with modeled singlegrain minimum ages) of the uppermost fine-to medium-grained sand sample is >20 ka.

Quartz multi-grain ages Supplementary
In the very fine-grained section of core GS-7, all three samples have less than 3% of the aliquots in saturation and so all are bounded. As was seen in the feldspar data, the quartz multi-grain age for the deepest of these samples (111912, 4.9±0.3 ka) is significantly younger than the age of the overlying sample 111910 (6.0±0.4 ka); this suggests some mixing of well-bleached material, or an error in dose rate, rather than incomplete bleaching (because a small degree of incomplete bleaching in feldspar -as seen here -should imply complete bleaching of quartz 10 ).
In core GS-10, five samples give bounded age estimates (see Supplementary Figure 11c, closed symbols). At the top of the core, there is a small age inversion between samples 111901 (2.3 m; 8.0±0.6 ka) and 101901 (4.3 m; 6.3±0.3 ka). The same inversion is present in the feldspar data, and (as above) we deduce that this is unlikely to arise from incomplete bleaching. Sample 111901 showed a significant IR contamination, but unfortunately, we did not have sufficient quartz to re-measure this sample using POSL. Otherwise the ages are in stratigraphic order and the age of the uppermost coarse/medium sand sample is 20±2 ka.
In core GS-11 there are several stratigraphic inconsistencies for the ages obtained for samples taken from below 6 m (Supplementary Figure 11d). However, all of these ages are derived from samples with a significant number of unbounded dose estimates and we thus consider these to be minimum ages. The remaining five samples are in stratigraphic order. The age of the uppermost fine-medium-grained sand sample is 23±2 ka.
In core GS-13, two of the three multi-grain quartz ages are considered to be minimum ages, but they all appear to be in stratigraphic order. The minimum age of the uppermost fine-to medium-grained sand sample is >26 ka.
In core GS-14, the three samples are not in stratigraphic order, and the ages of the first two samples are both regarded as minimum ages. The age of the deepest sample is 24±2 ka and the minimum age of the uppermost fine-medium-grained sand sample is >19 ka. Bearing in mind that two of the four ages for the uppermost fine-medium-grained sand samples from the GS-transect are considered to be minimum ages, the average age of these samples is 22.3±2.3 ka (n=4), completely consistent with the corresponding result from multi-grain feldspar, of 23.7±1.0 ka (n=4).

Supplementary Note 2: Single grain Quartz ages
Single grain quartz ages have been estimated using a range of different approaches. In Supplementary Table 9 we summarize average ages, both weighted (CAM) and unweighted (arithmetic) averages, minimum ages derived using three different minimum age models, the most prominent finite mixture model ages (FMM prom ) 79 all derived from the data sets obtained using the standard rejection criteria (see section OSL Measurements).
In general, the rationale for carrying out single grain measurements is to enable the identification of minimum burial ages for samples expected to be incompletely bleached or to identify (and reject) grains likely to have undergone post-depositional transport. It has often also been stated that the single grain technique enables the OSL characteristics of each individual grain to be evaluated, so that only reliable dose estimates contribute to the age determination 50,73 . We suggested that a large fraction of the multi-grain quartz ages must be regarded as minimum age estimates because a significant fraction (>15%) of the individual aliquots were in saturation (see Supplementary Note 1). By undertaking single grain analysis, we might be able to come up with more accurate age estimates for all samples, because we are able to reject individual grains with L n /T n values in or above saturation.

Average single grain quartz ages
For well-bleached samples that have not undergone significant post-depositional mixing, the best estimate of the burial age should be based on an average (usually weighted) dose estimate, even in the presence of a spatially variable dose rate 80 . As stated above, nine of the single grain dose distributions contain negative dose 20 estimates. The original CAM (i.e. CAM log ) assumes that natural dose distributions are lognormal and this prevents the application of this model to distributions containing non-positive dose estimates. For these nine samples, we calculate CAM (i.e. CAM UL in the terminology of Arnold et al. (ref. 81)) ages assuming the dose distributions are normally distributed; this allows the inclusion of all values. However, we note that the assumption of normality as against log-normality has no significant effect on the estimated dose or OD. The ratio of CAM log dose (calculated after the rejection of negative dose estimates) and the CAM UL doses (including all accepted dose estimates) is 1.01±0.04 (n=9 samples) and the ratio between the corresponding ODs is 0.99±0.04 (n=9 samples). If we were to apply the L n >3 s BG (refs. 73,82) criterion, in an attempt to avoid the inclusion of aberrant "0 Gy" grains, then these ratios are 1.07±0.05 and 0.96±0.04, respectively. If we only include samples expected to be older than 10 ka (based on the KF ages for the same samples) then the ratios are 1.07±0.07 and 0.99±0.04 (n=6 samples). Thus, for the nine samples containing negative dose estimates, the CAM ages and ODs given in Supplementary Table 9 are, within uncertainties, insensitive to the underlying assumption of the shape of the dose distribution.
The CAM ages for the individual cores are broadly in stratigraphic agreement with each other. The three ages for SRH-5A are all consistent with an average age of 13.5±0.6 ka. The CAM age of the uppermost medium-coarse-grained sand sample (111922) is 14±2 ka (n=37).
Only two samples were measured for GS-7 (111910 and -12) and they give CAM ages of 6.1±0.4 ka (n=76) and 4.6±0.3 ka (n=73), respectively. This apparent age reversal is also seen in both the feldspar and the multi-grain quartz ages.
The CAM ages for GS-10 are all in stratigraphic order (with the exception of sample 111906), although the data becomes noisy below 15 m (see Supplementary Figure 11e). The CAM age from the uppermost medium-grained sand sample (111903) in this core is 16±2 ka (n=69).
In core GS-11, all ages are in stratigraphic order except for sample 091904, which appears to be significantly younger than the ages derived from the samples taken above and below (see Supplementary Figure 11f). The single grain dose distribution for this sample contains a significant number of grains consistent with zero, nine of which are negative. Thus, it is not surprising that the CAM age derived for this sample appears to be too young. Unfortunately, we have not been able to make single grain measurements for the uppermost medium sand sample in this core (101905) due to lack of sample material.
The average CAM age for the uppermost fine-to medium-grained sand samples from the GStransect is 17.3±1.4 ka (n=3).
In single grain measurements, there is often a large variability in the intrinsic luminescence sensitivity and consequently a large variability in the relative uncertainties assigned to each individual dose estimate. Thus, it seems prudent to weight each individual dose estimate with respect to its individual uncertainty as is done in the CAM. However, the large relative ODs observed for these samples, even in dose recovery experiments, raises doubts as to the validity of the weighting factors. Thus, we have also calculated average arithmetic single grain ages for the accepted grains. These ages follow the same stratigraphic order as the CAM ages but are on average 21±2% (n=30) older than the CAM ages and gives an average arithmetic age of the uppermost fine-to medium-grained sand samples on the GS-transect of 20.9±1.8 ka (n=3).

Minimum single grain ages
In our view, the main reason for undertaking single grain analysis is to establish whether incomplete bleaching is significantly affecting these samples. In the literature it has often been argued that well-bleached samples have a relative OD of 20% or less 58,83 ; from this it is deduced that samples with a larger OD are likely to be affected by incomplete bleaching and thus that the burial age should be determined using minimum age models 83 . If the minimum uncertainty of 30% determined from the single grain dose recovery experiments is taken into account then all dose distributions (derived using standard rejection criteria) have relative ODs significantly larger than 20% with the exception of the two young samples measured from core GS-7 (samples 111910 and -12, which have relative ODs of 20±3% and 12±2%). Bailey and Arnold (ref. 84) suggested a decision tree making use of the relative OD, the skewness and the kurtosis of individual dose distributions to decide which model to apply to extract an accurate burial dose. Although this decision tree was developed for fluvial systems, it has since been applied to many different depositional environments including fluvial, glacio-fluvial, colluvial, alluvial, coastal and aeolian (e.g., If we apply this decision tree to our fluvial samples then, a minimum age model (either MAM3, MAM4 or the lowest 5%) is chosen for all samples except 111912, 111913, and 111919 for which the CAM should be applied. Thus, from the literature it is to be expected that at least the majority of these samples are incompletely bleached and hence that minimum age models should be applied. Consequently, we have applied three different minimum age estimation methods: (i) the four parameter minimum age model (MAM4 (ref: 54 Table 9 are the lowest dose component derived by the FMM containing more than 10% of the grains in the dose population. That is, if the smallest dose component contained less than 10% of the dose population the next youngest component is given. The number of components that were used to fit the data was determined using the log likelihood (llik) and the Bayesian Information Criterion (BIC (ref. 93)). Because of the presence of negative dose estimates we could not directly apply the FMM to the nine samples given above. We can compare the FMM min doses with the IEU doses for 19 samples (i.e. dose distributions with no negative dose estimates). The FMM min doses tend to be larger than or consistent with the IEU 22 doses (average FMM min /IEU=1.7±0.2, n=19); only 74% of the dose estimates are consistent with each other at two standard errors. Because of the good agreement between the MAM and IEU doses, we consider the FMM min doses to be less reliable. In the following discussion, we choose to consider the IEU ages only, because we could determine these for the largest number of samples.
The minimum ages (determined using the IEU) result in significant age reversals in all cores. In core GS-11 an unrealistically young IEU age of 0.09±0.07 ka is obtained for sample 091901 taken from a depth of ~2 m. Of even more concern is the unrealistic young minimum age of 0.4±0.2 ka for sample 091904 (GS-11) taken at a depth of ~25 m. However, all the unrealistic young ages are all obtained from the dose distributions containing negative dose estimates. As stated above, we regard it as physically very unlikely that a significant fraction of modern grains were transported to a depth of 25 m in these cores.
A comparison of CAM and IEU ages gives useful insights towards understanding the meaning of minimum ages in these samples. In Supplementary Figure 19 we plot the (CAM -IEU) difference as a function of CAM age and it is clear that there is a correlation. All data can be well represented by a straight line with a slope of 0.62±0.04. This observation can be interpreted in two ways: (1) the absolute dispersion in the data gets bigger with age (e.g., for a constant relative dispersion), and the minimum age is just sampling an edge of such a distribution. This scenario is very likely for a situation such as ours in which all sources of dispersion have not been accounted for and therefore the uncertainty budget is not complete.
(2) Alternatively, the difference between the CAM and the IEU represents an off-set caused by presumed incomplete bleaching, but this premise leads to the absurd conclusion that the degree of incomplete bleaching depends on the subsequent time of burial of the sample. We therefore reject this latter premise and rather interpret the observed trend to be an artefact of the application of minimum age models to those wide distributions for which the leading edge does not arise from an extrinsic phenomenon such as the degree of bleaching, but rather is an inherent part of a distribution whose central tendency is well represented by CAM. The inset to Supplementary Figure 19 shows the normalized difference between CAM and IEU (normalized by the CAM age) as a function of the proportion of the presumed wellbleached grains identified by the IEU (i.e., the proportion of grains used to derive the IEU age). All points are consistent with a straight line with a slope of 0.009±0.001 supporting our conclusion that the IEU does not identify an independent well-bleached population, but rather just extracts a part of it.
Based on this analysis we do not regard minimum age modeling of the data from these samples as physically meaningful.

Most prominent finite mixture model (FMM prom ) ages
In the decision tree of Bailey and Arnold 84 it is implicitly assumed that natural dose distributions are normally distributed. However, Arnold et al. (ref. 81) argue that natural dose distributions are log-normal and that the skewness should be calculated using logged dose values. Arnold and Roberts (ref. 73) propose a dose estimation decision tree based on the relative OD as well as the log-weighted skewness. If we apply this decision tree to our samples we conclude that the FMM should be applied to all samples except 111910, and -12 for which the CAM should be applied (based on the llik and the BIC score the FMM predicts that the CAM is more likely to be applicable than the FMM). In the past decade the FMM has widely been applied to samples suspected to be affected by post-depositional mixing 94 . The FMM identifies discrete dose components within a dose distribution; the most prominent dose component (FMM prom ) is usually assumed to be that relevant to the burial age, and this is calculated by dividing FMM prom by the average dose rate 72,[95][96][97][98][99][100][101] . It has also been argued that such dose components may arise from discrete dose-rate components, and some have suggested dividing FMM prom by a modeled dose-rate component 102 The most prominent dose component contains on average only 65±2% (n=27, not including the three samples where the FMM predicts that CAM is a more appropriate model) of the total accepted grain population implying that we have considerable (and in our view, unrealistic) post-depositional mixing in the majority of these samples. The average ratio of the FMM prom to the CAM ages is 1.25±0.08 (n=27). In the following description of the apparent chronostratigraphy of the individual cores, we make use of the ages derived from FMM prom (or CAM for the samples from GS-7).
In core SRH-5A, the three FMM prom ages are all consistent with an average age of 20±3 ka. The age of the uppermost coarse/medium sand sample (111922) is 25±3 ka.
For both samples in core GS-7, the FMM predicts that the CAM is likely to provide the best age estimate (see above).
In core GS-10, the FMM prom ages increase smoothly with depth until about 20 m; below this the ages are ~50 ka (see Supplementary Figure 11e). In this core, we appear to have two age inversions. The first inversion occurs for sample 111903, where the FMM identifies two components and gives a FMM prom age of 10.1±1.1 ka (p=58±9%), which appears to be too young. However, if we choose the second component of 33.4±4.0 ka (p=42±9%), the inversion disappears. Thus, one could argue that the older component should be chosen based on stratigraphic consistency. In any case, the choice of the largest component is arbitrary. If one postulates gross post-depositional mixing between sedimentary units, there is no a priori reason that grains from the original unit must always be in the majority. Thus, it is possible to argue that the age of the uppermost medium sand sample (111903) is 33.4±4.0 ka. The second inversion occurs for sample 111906, where the FMM predicts a single dose component giving an age of 29±2 ka. As stated previously, we have added 30% additional uncertainty onto individual dose estimates to account for intrinsic sources of variability (see section Single-grain dose recovery experiments). However, if one were to add 15% instead of 30% the FMM predicts two components, each containing ~50% of the grains, giving ages of 20±3 ka (p=47±15%) and 41±4 ka (p=53±15%), respectively. The average age of the samples taken from depths >18 m is 48±3 ka, so choosing the larger component removes this inversion. However, none of the beta dose recovery experiments resulted in dose distributions with ODs of ~15%, so we do not think it is justified to change the additional uncertainty parameter to 15%. It has been suggested that the appropriate additional uncertainty parameter can be determined by optimizing the llik and BIC scores provided by the FMM (ref. 93). In the following we test this approach by varying the additional uncertainty (in addition to the uncertainty assigned to individual dose estimates based on Poisson statistics and curve fitting uncertainties) between 5 and 60%. For sample 111906 it is not straightforward to determine the optimum uncertainty parameter as the llik score fluctuates considerably and the BIC score decreases smoothly for an additional uncertainty larger than 5% (see Supplementary Table 11). If we choose the optimum scores to occur at 25% then the resulting FMM prom age is 36±5 ka (p=71±24%). This age is not significantly different from the result obtained using 30% and does not fall in stratigraphic order. In Supplementary Table  11 we show the FMM results for 16 of the samples obtained using a range of additional uncertainties. If we choose the additional uncertainty parameter based on the highest llik score and the lowest BIC score, the additional uncertainty parameters range between 15 and 50%. The ratio of the FMM prom age determined using the optimized uncertainty parameter to the FMM prom age using an additional uncertainty of 30% is 1.01±0.05 (n=16) implying that no significant benefit is obtained by optimization. Only for one sample (111909) is the individual ratio not consistent with unity at two standard errors. Here, the FMM predicts an additional uncertainty of 50% and two dose components with ages of 24.7±3.1 ka (p=51±10%) and 72.1±8.7 ka (p=49±10%; see Supplementary Table 11). Using an additional uncertainty of 30%, the FMM identifies three components (see Supplementary Table 10 and  Supplementary Table 11) and predicts an age of 42±7 ka. Thus, for this sample one could argue for an age ranging between ~25 and 70 ka despite the fact that this sample is of aeolian origin.
In core GS-11, the FMM prom ages (30% additional uncertainty) broadly appear to increase with depth although samples 101906 and 091906 appear to be too young giving ages of ~13 ka. However, for both samples the FMM identifies two components in almost equal proportions. If we choose the largest dose components we obtain ages of 34±11 ka (091906) and 38±6 ka (101906) and the FMM prom ages are all in stratigraphic order and increase smoothly with depth (see Supplementary Figure 11f). Sample 091905, which is aeolian and expected to be much older, gives a FMM prom age of 115±11 ka.
For cores GS-13 and GS-14, all FMM prom ages are again in stratigraphic order and the ages of the uppermost fine-to medium-grained sand samples from the GS-transect are 12±2 ka and 15±2 ka, respectively.
The average FMM prom age for the uppermost fine-to medium-grained sand samples from all GS cores is 12.4±1.3 ka (n=3), which is significantly younger than the average CAM age (17. 3±1.4 ka, n=3).
It appears that for more than half the samples the FMM prom gives ages consistent with the CAM ages. However, the results suggest that many of the samples (including those from the deepest aeolian sediments) have been mixed after deposition, so that the original sediment now only represents about half of the total; the model suggests that the other half usually come from sediment that is on average ~25 ka different in age without any significant contribution from intermediate horizons.

Supplementary Note 3: Comparison of quartz and feldspar ages following standard selection criteria
Given the results from the modern analogue sample and the comparatively large apparent age in most of our core samples, we consider it unlikely that incomplete bleaching is a significant issue in this study; it may affect some of the younger ages but as discussed above these are not relevant to the determination of the time when major river flow ceased.
Nevertheless, in the absence of independent age control, we do not know which of our luminescence ages are necessarily more accurate. However, we can use stratigraphic constraints and a comparison of independent chronometers (quartz and feldspar) to identify the more reliable ages. Supplementary Figure 20a we plot all the multi-grain quartz ages against the fading corrected feldspar ages obtained using the IR 50 signal. The closed symbols (n=17) represent quartz ages which we term 'reliable' (<15% of aliquots giving an unbounded dose) whereas the open symbols (n=20) represent minimum ages (>15% unbounded dose estimates). The dashed line represents the 1:1 line and the shaded band around it represents ±10%. The inset shows the same data for feldspar ages less than ~35 ka. Although there is considerable scatter around the 1:1 line, the average ratio of all quartz to feldspar ages is 1.04±0.04 (relative s=26%; n=37). The slope of the fitted (solid) line in Supplementary Figure 12a is 0.90±0.03. If we include only the 17 reliable quartz ages, the average ratio is 1.07±0.05 (relative s=19%; n=17). These two quartz/feldspar ratios are both indistinguishable from unity, although the standard deviation is slightly smaller when using the "reliable" quartz ages. The quartz fast component and the IR 50 feldspar signal are known to bleach at significantly different rates in daylight (~1 order of magnitude difference). This, taken together with the results from the modern analogue and the age agreement discussed above, strongly suggest that both minerals were bleached well enough at deposition; this conclusion is consistent with published observations 32 . This agreement gives additional support to our conclusion that minimum age models are not relevant for our data, despite the predictions of the Bailey and Arnold's (ref. 84) decision tree. Supplementary Figure 20b, we plot the single-grain CAM quartz ages as a function of the feldspar ages. In the absence of significant incomplete bleaching and post-depositional mixing the single grain CAM is normally expected to give the most accurate ages. However, the CAM ages significantly underestimate the feldspar ages: the average quartz CAM to feldspar age is 0.67±0.05 (n=30) and the slope of the line in Supplementary Figure 12b is 0.46±0.02. Only 7 of the 30 ratios are consistent with unity within 2 standard errors. Given that the feldspar ages on average are ~70% older than the CAM ages; it is very difficult to attribute this disagreement to incomplete bleaching of feldspar as any effect of bleaching should, on average, be additive, not proportional. Moreover, the two oldest samples (111909 and 091905) known to be of aeolian origin, also appear to show that single-grain CAM quartz age underestimates the feldspar age by >40%. One explanation could be over-correction for anomalous fading in feldspars. The average fading correction is ~40% (based on a g value of ~3.5 % per decade); this is a typical value for the IR 50 signal 30,106 . Even if we assume that these samples do not suffer from any anomalous fading (despite the laboratory measurements), the average slope for all the samples is only increased to 0.65±0.01 (n=30, data not shown). However, the evidence for anomalous fading is unambiguous and we do not consider it likely that such over-correction is the cause of the difference between the single grain CAM ages and the feldspar ages. From these data alone it is difficult to decide which of the two data sets is likely to be the more accurate, although the feldspar ages suffer less from stratigraphic inversions than the quartz ages. However, this analysis does not take into account the effect of bias resulting from the rejection of grains in dose saturation. Supplementary Figure 20c, we plot the single grain FMM prom quartz ages against the feldspar ages. The slope of the linear fit to these data (excluding the extreme outlier, sample 111909 open symbol) is 0.76±0.02 (n=29) and the average quartz to feldspar age ratio is 0.81±0.06 (n=29). There is significant scatter in these data and the agreement is not as good as for standard multi-grain aliquot quartz measurements. The overall tendency is that the FMM prom ages are smaller than the feldspar ages. This is particularly pronounced for the deepest sample in GS-10 (111909) which is of aeolian origin. The FMM identifies three components for this sample with the most prominent component based on 47±7% of the grains and giving an age of 42±7 ka. The component with the largest dose contains 26±7% of the grains and gives an age of 104±13 ka. The feldspar age for this sample is 150±6 ka, so even if one were to choose the larger dose component for age estimation (despite the fact that this would not be advocated in the literature), the FMM prom age would significantly underestimate the feldspar age. At this stage in the argument we cannot identify whether the feldspar age is in error or the single grain quartz FMM prom age underestimates the deposition age.

Single-grain (FMM) compared with feldspar and interpretation of FMM results In
The most widespread use of the FMM (refs. 73,94,99) is to identify discrete dose components arising as the result of post-depositional mixing. However, for >30% of the samples investigated here, the main dose component contains less than 60% of the grain population; this implies that more than 40% of the grains in the unit of interest have migrated from one or two discrete units often lying at a considerable distance (in some cases above, in some below) from the layer of interest.
Three examples serve to illustrate that this is stratigraphically extremely unlikely. Firstly, in core SRH-5A, the FMM gives an older age of 25±3 ka (54% of grains) for sample 111922 (429 cm depth). The younger component of the mixture is 7.2±1.4 ka (46% of grains). This sample is at the top of the grey fluvial sand succession and is immediately overlain by floodplain clay and silt ( Fig. 8 and Supplementary Fig. 5). Thus, the younger component can only be sourced in the floodplain sediment. To dilute the original fluvial sand such that today 46% of the deposit is made up of grains derived from floodplain sediment would require (a) a large-scale process of very differentiated grain transport preferentially selecting only the sand grains from the overlying clay-silt, and (b) destruction of the distinct grain size transition between the two sedimentary facies. It is implausible that there has been mixing between the fluvial sands and the fine-grained deposits, because the latter do not contain any admixed fluvial medium-to coarse grained sand. Secondly, in core GS-11, the FMM gives an older age of 56±7 ka (51% of grains) for sample 091904 (2495 cm). The younger component of the mixture is 12±2 ka (49% of grains). This sample was taken from a grey fluvial sand unit and also lies immediately below a floodplain clay and silt unit (Fig. 5). In this case, the younger component could in principle be sourced from any of the fluvial sand units overlying this floodplain unit (see minor components in GS-11 in Supplementary Table 10). However, sediment of consistent age (13±2 ka) is not a major component of any deposit until a depth of ~920 cm (sample 101906) (Fig. 5). This more likely source unit is separated from sample 091904 by >15 m and four stratigraphically distinct floodplain clay and silt units. Again, to dilute the original fluvial sand such that today 49% of the deposit is made up of the material sourced from at least ~9 m depth would require the destruction of a number of grain size transitions observed in the sedimentary succession between the two units. Mixing of grain components between these two depths without homogenization of the succession is highly implausible. Finally, in core GS-13, the FMM gives an older age of 20±3 ka (71% of grains) for sample 111914 (507 cm). The younger component of the mixture is 5±2 ka (29% of grains). This sample was taken from a grey very fine sand bed and also lies immediately below a floodplain clay and silt unit (Figs. 5 and 7; Supplementary Figs. 2 and 3). The youngest FMM age in the grey very fine sand (sample 111913; 245 cm depth) overlying this floodplain unit is 12±2 ka (68% of grains) which is too old to source the younger component of sample 111914. Thus, the potential source must lie closer to the surface, in the undated aeolian or red silty clay units. To dilute the original fluvial sand such that today 29% of the deposit is made up of the material sourced from <245 cm depth would require the destruction of at least three distinct grain size breaks separating the two units. In particular, a silty claystone unit that is interstratified between the sand units of samples 111913 and 111914 (and which represents about 1 m of section between 4 and 5 m depth bgl) does not contain any sand grains ( Supplementary Fig. 2), thus invalidating a post-depositional mixing explanation. Thus, the presence of distinct grain size breaks in the stratigraphic succession clearly argues against any large-scale sediment mixing in this core as there is no evidence of grain size homogenization.
Animal burrows are frequently mentioned in justifying such results, but such gross mixing is undifferentiated and cannot explain the result from SRH-5A and, in any case, would have completely destroyed the observed grain size transitions observed in all examples. In summary, we conclude that physical mixing of these stacked units can be safely dismissed as an explanation for the apparent dose mixtures observed in our data.
Another explanation for the presence of components in a FMM analysis is that they are a consequence of extreme dose rate heterogeneity. However, this explanation can be dismissed with confidence, as Guérin et al. (ref. 80) have pointed out, modeling has shown that discrete dose rate populations do not exist except in very unusual conditions, and these do not apply to well mixed and homogeneous sands such as found in our samples.
The FMM components could arise from intrinsic characteristics such as differing luminescence properties. Although at this point it is difficult to see why some quartz grains should result is such a completely different dose estimate from other grains, but the possibility that, for instance, sensitivity changes during the measurement of the natural signal may be systematically different between two groups of grains cannot be dismissed.
Finally, the various components may be a statistical artifact resulting from the application of the model to a very broad dose distribution with unrealistic estimates of the intrinsic overdispersion. In the Supplementary Table 11 we investigate the effects of using different estimates of the intrinsic over-dispersion on the FMM ages; from these data, it appears that the FMM predicts a single dose distribution if the OD is set to ~50%.
Given the unlikely results obtained using the Finite Mixture Model, we are unable at this stage in the argument to decide whether the most accurate ages are obtained from quartz single-grained, from quartz multi-grain aliquots or from multi-grain feldspar aliquots. In Supplementary Note 4 we present two additional single-grain selection criteria to address this problem.

Supplementary Note 4: Additional single grain rejection criteria
As discussed in section Single-grain OSL dose response curves, the D 0 values of the dose response curves for individual grains vary significantly between grains (see Supplementary  Figure 17b), and grains which are rejected due to saturation tend to have lower D 0 values than those which are accepted. The fact that there is a difference in D 0 values between grains rejected only due to saturation and accepted grains is potentially of serious concern. For instance, consider a group of grains all with a D 0 value of 75 Gy extracted from a sample with a dose of 100 Gy (see Supplementary Figure 17c). For ease of calculation we have assumed the saturation value A (see section B.4) to be unity. Suppose that the L n /T n values are normally distributed with a standard deviation of 10%. This distribution of L n /T n values will be centered on ~0.74. If we then extract 1,000 random L n /T n values from this normal distribution and interpolate them on to the dose response curve, we obtain a skewed dose distribution (grey histogram in Supplementary Figure 17c). The degree of skewness will vary depending on the relationship between the given dose and the D 0 value, but as Murray and Funder 37 pointed out, the median of this dose distribution will be closer to the given dose than the average. However, let us now suppose we have a group of grains all with dose response curves with a D 0 value of, say, 25 Gy. Thus, the L n /T n distribution will be centered on ~0.98 just 2% lower than the saturation value. Again, we assume that it is a normal distribution with a standard deviation of 10%. If we randomly select 1,000 values from this L n /T n distribution and interpolate them onto the dose response curve with a D 0 value of 25 Gy, we obtain dose estimates from only approximately 50% of the L n /T n values. These dose estimates are shown as the white histogram (and arise from the white part of the L n /T n distribution centered on 0.98). As can be clearly seen from Supplementary Figure 17c, this dose distribution is highly skewed with a mode of approximately 50 Gy. Thus, we deduce that it is entirely possible that at least some of the grains which have been accepted into our natural dose distribution actually come from a distribution of grains with low D 0 values, and so which have absorbed a saturation dose; inclusion of such grains will bias it to low doses. In an attempt to circumvent this problem, we suggest all grains with low D 0 values (compared to the natural dose) should be rejected. Gliganic et al. (ref. 72) suggested removing all dose estimates with a D 0 value less than 25 Gy, but this was suggested solely because of their inability to obtain acceptable dose recovery ratios for their samples. We suggest a more fundamental reason for using such a criterion. For the OSL signal to be a useful dosimeter it is important that it is both optically and thermally stable over geological time scale. The fast component in quartz OSL has been shown to derive from the 325 °C TL peak which has a life time of ~10 8 y (ref. 12). However, other less thermally stable traps may contribute to the OSL signal 51,60 but it is anticipated that most of these contributions are eliminated by the chosen thermal pre-treatment (preheat).
However, Moskva and Murray (ref. 107) reported that when they sensitized the fast component in their previously insensitive samples a large fraction of their resulting fast component OSL signal was unstable. An unstable rapidly decaying OSL component has also been reported in heated quartz 108 . In section B.1.3 we showed for large multi-grain aliquot measurements that a preheat of 260 °C, 10 s and a cutheat of 220°C was an appropriate thermal treatment for these samples. However, it is possible that this may not be true for every grain 109 . To investigate the thermal stability of individual grains we measured pulse anneal curves for 2,400 grains from sample 111909. Prior to these measurements, we measured the natural doses in these grains as described in section OSL Measurements. The pulse anneal curves were measured using the SAR procedure in which the entire signal is measured for each preheat temperature. The preheat temperatures ranged from 200 to 340 °C and were held for 10 s. After the final measurement at 340 °C four repeat measurements at 220, 260, 280 and 320 °C were made. Any sensitivity change was monitored using a test dose with a preceding cutheat of 220 °C. High temperature blue stimulation at 280 °C for 40 s was inserted after each L x and T x cycle to minimize recuperation effects. The regeneration and test doses were both 50 Gy in this experiment. The shapes of these pulse anneal curves varies considerably. In Supplementary Figure 17e we show representative pulse anneal curves for four selected single grains from sample 111909 to illustrate this variability. The pulse anneal curve for grain D24G43 (circles) appears to be similar in shape to published data from multigrain aliquots 51,110 , whereas the other three grains appear to be less thermally stable to varying degrees. The apparently least thermally stable grain (D2G100, squares) has the smallest estimate of equivalent dose of 46±18 Gy, whereas the most apparently thermally stable grain (D24G43, circles) has the highest equivalent dose of 254±13 Gy (as observed by Fan et al. (ref. 109)). Although it probably would be advantageous to measure pulse anneal curves for every grain to establish the thermal stability on a grain-by-grain basis, we recognize that this would be very labor and time consuming. Instead we propose to make a single measurement using the chosen thermal treatment (in this case a preheat temperature of 260 °C and a cutheat temperature of 220 °C) and the same doses for both L x and T x . Since the regeneration and test doses are the same, if the signal is stable one would expect a similar response during both the regeneration and test dose measurement, despite their different prior thermal treatment. Thus, we can subsequently reject grains with a low value of this L x /T x ratio (hence forth called the L ts /T ts ratio) due to the possibility of low thermal stability (although of course some of this difference could result from sensitivity change). In Supplementary Figure 17f we investigate the effect of rejecting grains based on the L ts /T ts ratio on the average (arithmetic) dose. If we apply this stability criterion to the dose distribution obtained using the standard rejection criteria (D 0 > 0, squares) then we see little effect of rejecting grains based on the L ts /T ts ratio (Supplementary Figure 17f, green squares). However, if we apply it to the dose distribution only containing grains with dose response curves with D 0 > 100 Gy, we see that the average dose increases steadily from ~140 Gy to ~200 Gy for L ts /T ts values greater than 0.8. If we tighten the additional D 0 criteria so to only include dose estimates with D 0 values greater than 200 Gy, we see a similar pattern but the average dose increases to ~330 Gy (CAM dose of ~310 Gy). The shaded band in Supplementary Figure 17f is the quartz dose predicted from measurement of feldspar (Supplementary Note 1). In Supplementary Figure 17a the effect on the dose distribution of applying the L ts /T ts > 0.8 criterion in combination with the D 0 >100 Gy (light and dark grey symbols) or D 0 >200 Gy (only dark grey symbols) criteria is shown. When these criteria are invoked we see that we primarily reject grains with low D e values. Of the original 229 dose estimates, we are left with just 40 and 9 dose estimates when applying the additional criteria for D 0 >100 and D 0 >200 Gy, respectively. We have carried out the same analysis for all 30 single-grain dose distributions using a L ts /T ts ratio > 0.8 and determined the D 0 cut-off value appropriate for each sample (see Supplementary Table 9). The ratio of the doses estimated using the two new additional rejection criteria to the doses estimated using the standard rejection criteria is 1.39±0.08 (n=30) for the arithmetic average and 1.46±0.09 (n=30) for the CAM average. In general, the younger the sample, the smaller the effect of applying these criteria, e.g. for samples with KF ages <10 ka (i.e.101901, -03, 111910, -12) the average ratio of the arithmetic mean is 1.03±0.01 (n=4) and 1.06±0.04 (n=4) for the CAM. In Supplementary Figure 11e and Supplementary Figure 11f we show the derived ages as a function of depth for cores GS-10 and GS-11, respectively.
Although the single-grain quartz dose distributions obtained using the additional rejection criteria generally contain only a few grains (between 4 and 139) it is still worthwhile considering the average effect on the relative OD. The average ratio of the OD calculated for the dose distributions obtained with the new and the standard rejection criteria is 0.70±0.05, i.e., the relative OD is on average decreased by ~30%.
We also applied these additional two rejection criteria to 6 of the dose recovery data sets (111910, -12, -15, 101903 and -06). For sample 111916 we are left with just 3 dose estimates and this data set is subsequently excluded from the following analysis. The average dose recovery ratio for the five samples using the additional rejection criteria is 1.09±0.04 and the average relative OD is 22±3%. These values are completely consistent with the results obtained using the standard rejection criteria for the same five samples, i.e. a dose recovery ratio of 1.01±0.06 and an average relative OD of 27±3%.
In Supplementary Figure 20d we show the average single grain quartz ages obtained using the two additional rejection criteria as a function of the feldspar ages. The slope of the linear fit to these data is 0.98±0.02 (n=30) and the average quartz to feldspar age ratio is 1.04±0.05 (n=30) with 85% of the ratios falling within 2 standard errors of unity. The corresponding numbers for the CAM ages using the additional rejection criteria are 0.89±0.02 and 0.90±0.04. Thus, it would appear that we obtain a good agreement with feldspar if we employ the two proposed rejection criteria, despite the fact that we are left with very few quartz grains in each dose distribution. This gives us sufficient confidence in our feldspar ages to use these as the best estimates of burial ages.

Supplementary Note 5: Comparison of OSL results
We used OSL to obtain burial ages for 52 samples extracted from seven individual sediment cores. We dated these samples using the IR 50 signal from multi-grain K-feldspar aliquots (2 mm) and the blue/green OSL signals from multi-grain (8 mm) and single-grain quartz aliquots.
The fading corrected IR 50 feldspar ages are in stratigraphic order except for the two younger fine-grained samples (111912 and 111901/101901) in cores GS-10 and GS-11, respectively. These apparent inversions are also observed in the quartz ages suggesting that these minor inversions of ~1.5 ka most likely result from dose rate measurement problems. The measured dose rates for samples 111912 and 101901 are ~25% higher than those measured for the samples immediately above (111911 and 111901, respectively). If one were to use the dose rates determined for the samples directly above the age inversions would disappear. However, given the small absolute off-set of these inversions (~1.5 ka) and the otherwise good stratigraphic relationship between the feldspar ages, we do not regard these inversions as important.
To assess the importance of any potential incomplete bleaching we dated a modern analogue taken from where the modern Sutlej river leaves the Siwalik foothills and enters the Indo-Gangetic plain and determined a burial age of ~0.7 ka. Thus, we would regard these samples as probably well-bleached although it must be recognized that the feldspar ages for the very young samples (<10 ka) may be slightly overestimated. Multi-grain quartz analysis of the modern analogue gave a residual dose of 0.9±0.2 Gy (n=24), corresponding to an apparent burial age of ~400 years, which also indicates that the sediments delivered by this fluvial system are likely to be relatively well-bleached. Given the significant difference in bleaching rates between quartz and feldspar (approximately one order of magnitude), we interpret this agreement as indicating that both minerals were well-bleached when deposited.
The ages derived from multi-grain quartz measurements are broadly in stratigraphic order, but with the exception of the two young samples mentioned above as well as samples 091906 and 111916, which both appear too old. However, a significant fraction (23%) of the measured multi-grain aliquots from many samples were in or above saturation on the laboratory dose response curves and thus no bounded dose estimates could be derived. We have arbitrarily chosen only to consider ages to be reliable if we can derive bounded dose estimates for more than 85% of the measured aliquots. The remainder are regarded as minimum ages; of the 37 samples measured only 17 samples are bounded (Supplementary  Table 8). For the latter, we obtain a broadly good agreement with the corresponding feldspar ages (Supplementary Figure 20a), although there are two out of 17 samples (i.e. samples 111906 and -23) for which the quartz to feldspar age ratio is more than 3 standard errors away from unity.
To investigate whether or not these discrepancies arose because of intrinsic problems such as poor luminescence characteristics and/or extrinsic processes such as incomplete bleaching or post-depositional mixing, single grain dose distributions were measured for 30 of the 37 samples (Supplementary Table 9). A comparison of all of the single grain CAM ages with feldspar ( Supplementary Figure 20b) gives us a systematic fractional underestimate (with no significant intercept on the feldspar axis); this is not consistent with the expected effects of incomplete bleaching. This conclusion is supported by the application of minimum age models to the single grain quartz data which gives unrealistically low ages with several very significant stratigraphic inversions (Supplementary Table 9). By comparison with the CAM ages we also observe an effect of apparent incomplete bleaching which is a constant fraction of age (~60%; see Supplementary Figure 19). This is clearly unrealistic and the minimum ages are dismissed as inappropriate and inaccurate.
It is difficult to explain the differences between the multi-grain aliquot quartz and K-feldspar ages in terms of post-depositional mixing, because such mixing should have affected both ages in a similar manner. Thus, post-depositional mixing may result in stratigraphic inconsistencies but these should not be visible in a comparison of the two methods such as the one shown in Supplementary Figure 20a. Nevertheless, it may be that presumably intrinsic characteristics have led to the discrepancies in Supplementary Figure 20a and it may be possible to identify these "good" and "bad" grains using the FMM. However, application of the FMM results in significantly larger dispersion in the data compared to feldspar than that observed in the multi-grain quartz data, although it must be recognized that in the multigrain analyses samples that had a large number of aliquots in saturation were rejected completely, whereas in the FMM only grains that did not give bounded dose estimates were rejected (this did not result in the rejection of any entire sample). Nevertheless, when using standard rejection criteria, it remains unclear whether the FMM results or the K-feldspar results are to be preferred, although the mixing implications of the FMM are physically unrealistic.
Finally, we tested two additional single grain quartz rejection criteria based on the shape of the dose response curve and apparent thermal stability. Application of these criteria identified the majority of grains as unable to record accurately doses in the range of interest. The ages based on the few remaining accepted grains compare well with feldspar ages (Supplementary Figure 20d) with no detectable intercept on the feldspar axis and an average ratio of 1.04±0.04 (CAM ratio of 0.90±0.04). This agreement strongly suggests that these single grain quartz and K-feldspar ages should be preferred as the best estimates of deposition age. While we reject the result of the FMM analysis as indicating unrealistic degrees of mixing, the agreement in Supplementary Figure 20d does not necessarily rule out entirely the possibility of some degree of post-depositional mixing which presumably could affect quartz and feldspar to the same degree. Unfortunately, there are not enough single quartz grains identified as reliable dosimeters in the range of interest to allow for the application of the FMM to the accepted dataset. Nevertheless, the stratigraphic consistency of the feldspar data suggests that post-depositional mixing is not likely to have affected the ages significantly.
Given the absence of any suggestion of a systematic difference between the accepted quartz ages (both single grain and multi-grain) and fading-corrected multi-grain feldspar ages, and the fact that residual doses in related modern sediment are negligible, we consider the feldspar ages as the most accurate and precise.
Supplementary Figure 1 Characteristics of sediments at base of GS cores. a, Photograph of core section at base of GS11 core at depth of ~35 m bgl showing abrupt transition from light yellow brown well sorted, fine-grained aeolian sand to dark grey, mica-rich, medium-coarse-grained fluvial sand. b, Detail of facies transition. Note the light grey, calcrete nodules at top of aeolian sand (arrowed). The abrupt transition to fluvial sands records an incursion of a fluvial system into the area.  All average values for CAM and OD are consistent with unity. "s TN <20%" includes all grains with a bounded dose estimate with uncertainties on the natural test dose signal of less than 20%. "Recyc." includes all grains from "s TN <20%" except the ones with recycling values not consistent with unity at 2s. "IR depl." includes all grains from "Recyc." except the ones with IR depletion ratios not consistent with unity at 2 s. "Recup" includes all grains from "IR depl." except the ones with recuperation doses statistically greater than 1 Gy.   Grey micaceous medium-fine sand. This highly mica rich horizon is ~5 m thick.

Modern river samples Yamuna River sand
Medium-fine sand collected from modern sand bar.
Sutlej River sand Medium-fine sand collected from modern sand bar. Summary of the grain size, estimated water content (w.c.), and radionuclide concentrations used to calculate the total quartz and feldspar dose rates given in the last two columns. "KF" is potassium feldspar, "MG Q" is multi-grain quartz and "SG Q" single-grain quartz. Uncertainties represent one standard error. † assumed value. A water content of 30% is used for the 1319 samples based on the mean saturation water content of all other samples in this study. The transition between fine-and medium-grained sand and very fine-grained facies lies between the two entries in bold in each core, except for core GS-13. Supplementary Table 9. Summary of single grain data. The transition between fine medium sand and very fine grained facies lies between the two entries in bold in each core, except for core GS-13. The single grain measurements were made in single grain discs with hole diameters Ø of either 200 or 300 µm. "# neg. De" is the number of negative dose estimates in each dose distribution. "N" is the total number of grains measured. "In Sat" is the number of grains rejected solely because bonded dose estimates could not be derived. "OD" is the relative over-dispersion calculated solely based on counting statistics and curve fitting errors. "n" is the number of grains accepted into each dose distributions. "Average age" is the age calculated using the arithmetic mean of all accepted dose estimates. "CAM age" is the age calculated using all dose estimates and the Central Age Model (52). The CAM ages calculated for samples containing non-positive dose estimates is based on an assumption of normality. "FMMprom age" is the age calculated based on the most prominent dose component identified by the Finite Mixture Model (FMM, (72)). All non-positive dose estimates were removed prior to application of the FMM. The FMM model was applied after an additional uncertainty of 30% was added in quadrature to individual dose estimates. "MAM age" is the age calculated based on the burial dose estimate identified by the Minimum Age Model (MAM, (52)). For the dose distributions containing dose estimated the MAM age has been derived by applying the original MAM scripts after doing a simple exponential transformation (66). The MAM model was applied after an additional uncertainty of 30% was added in quadrature to individual dose estimates. "MAM age (no neg)" are the MAM ages calculated after elimination of negative dose estimates. "IEU age" is the age calculated using the Internal/external consistency criterion (IEU, (58), (78)). "FMMmin age" is the age calculated based on the lowest dose component containing more than 10% of the grains identified by the FMM. All non-positive dose estimates were removed prior to application of the FMM. The FMM model was applied after an additional uncertainty of 30% was added in quadrature to individual dose estimates. "Average age" and "CAM age" using the two additional rejection criteria is based on the arithmetic mean and the CAM mean of the dose estimates left after rejecting grains with D0 < 100 Gy (or 200 Gy for samples 111909 and 091905) and Lts/Tts ratios < 0.8, respectively.  Detailed FMM results for 16 of the 30 natural single-grain dose distributions for the most prominent (FMM prom ) and second most prominent (FMM prom,2nd ) grain populations." n neg ." is the number of non-positive dose estimates removed prior to running the FMM. "llik" and "BIC" are the log likelihood and the Bayesian Information Criterion, respectively. "Add (%)" is the additional uncertainty added on to the uncertainties assigned to individual dose estimates based on Poisson statistics and curve fitting uncertainties. "k" is the total number of components. The "Age" of each component is given in units of ka at 1 standard error and "p" is the proportion of grains attributed to each dose component. "FMM prom " is the most prominent component (i.e. the component with the highest value of p), and "FMM prom,2nd " is the second largest component. The ages highlighted in boldface are the ages obtained by optimising llik and BIC.