Maryam Heydari

Guillaume Guérin

Nikolay Sirakov

Philippe Fernandez

Catherine Ferrier

Aleta Guadelli

Jean-Claude Leblanc

Stanimira Taneva

Svoboda Sirakova

Jean-Luc Guadelli

The last 30,000 to 700,000 years ago: unravelling the timing of human settlement for the Palaeolithic site of Kozarnika

Kozarnika cave is a renowned prehistoric site in the Balkans, which contributes significantly to our understanding of the human past due to its rich assemblages associated with the Lower to Upper palaeolithic periods. To unravel the timing of human occupations in Kozarnika, various dating methods have been employed in the past. Radiocarbon dating was used to unfold the timeframe for the Kozarnikian lithic tradition uncovered in the Upper Palaeolithic sequence of the cave, and paleomagnetic dating assigned the Brunhes-Matuyama reversal to the layer beneath the Lower Palaeolithic assemblages.

In this study, we employed luminescence-dating methods including a body of different signals to date sediment layers from the top to the bottom of the sequence covering the period of ca 30 to 700 ka. Our investigations revealed that the Kozarnikian tradition in layers 5a-c falls between 30-35 ka. Following that, we suggest that the Middle Palaeolithic period initiated between 250-309 ka and lasted until 40-53 ka. More importantly, we have updated the age of the Neanderthal radius discovered in the Mousterian assemblages to 201±17 ka.

Introduction

The Balkans are renowned for witnessing different hominin and human species of inhabitants due to the location between the Near East and the Western Europe during the long period of the Palaeolithic. Human remains associated with the Neanderthals or anatomically modern humans (AMH) have been discovered in the Middle and Upper palaeolithic assemblages in various sites (e.g., [START_REF] Galanidou | Early Hominids in the Balkans[END_REF][START_REF] Glen | Human remains. In: Excavation in the Bacho Kiro Cave (Bulgaria)[END_REF][START_REF] Hoffecker | The spread of modern humans in Europe[END_REF]Hublin et al., 2020;[START_REF] Strait | The Human Fossil Record of Bulgaria and the Formulation of Biogeographic Hypotheses[END_REF][START_REF] Tillier | Evidence of Neanderthals in the Balkans: The infant radius from Kozarnika Cave (Bulgaria)[END_REF]. Additionally, the high possibility of interaction between two human species [START_REF] Fu | An early modern human from Romania with a recent Neanderthal ancestor[END_REF][START_REF] Hajdinjak | Initial Upper Palaeolithic humans in Europe had recent Neanderthal ancestry[END_REF] renders this region unique for further investigations.

Establishing a reliable temporal framework disclosing the timing of hominin species' settlements, the demise of one (e.g., Neanderthal), and domination of another (e.g., AMH) is a fundamental study, and it can be revealed through applying numerical dating methods. Among the whole body of techniques, radiocarbon dating plays a major role to track dynamics of the human peopling during the Late Middle to Upper Palaeolithic periods generally for Europe (e.g., [START_REF] Higham | The timing and spatiotemporal patterning of Neanderthal disappearance[END_REF] and in particular for the Balkans (e.g., [START_REF] References Alex | Radiocarbon chronology of Middle and Upper Paleolithic sites in Serbia, Central Balkans[END_REF][START_REF] Devièse | Direct dating of Neanderthal remains from the site of Vindija Cave and implications for the Middle to Upper Paleolithic transition[END_REF][START_REF] Fewlass | A 14C chronology for the Middle to Upper Palaeolithic transition at Bacho Kiro Cave, Bulgaria[END_REF][START_REF] Higham | Revised direct radiocarbon dating of the Vindija G1 Upper Paleolithic Neandertals[END_REF][START_REF] Trinkaus | An early modern human from the Pestera cu Oase, Romania[END_REF].

Whilst numerous temporal frameworks have been constructed to date the humans' fossils or the artefacts attributed to the Middle/Upper Palaeolithic, the number of available datasets is relatively scarce for the Lower to Middle Palaeolithic. This could be explained due to the shortage of the well-preserved Lower Palaeolithic sequences in eastern Europe [START_REF] Doronichev | The Lower Paleolithic in Eastern Europe and the Caucasus: A Reappraisal of the Data and New Approaches[END_REF]. Having said that, the Palaeolithic site of Kozarnika has been a special place all along not only because of its rich archaeological remains related to the Lower Palaeolithic but also for offering the first Bulgarian Neanderthal's remain, a radius associated with a rich Mousterian assemblage [START_REF] Tillier | Evidence of Neanderthals in the Balkans: The infant radius from Kozarnika Cave (Bulgaria)[END_REF]. Furthermore, the Upper palaeolithic industry found in the cave stands out from the other lithic productions and is known as Kozarnikian [START_REF] Guadelli | Une séquence du Paléolithique inferieur au Paléolithique récent dans les balkans : La grotte kozarnika à oreshets (nord-ouest de la bulgarie)[END_REF][START_REF] Tsanova | Les débuts du Paléolithique supérieur dans l'Est des Balkans. Réflexion à partir de l'étude taphonomique et techno-économique des ensembles lithiques des sites de Bacho Kiro (couche 11)[END_REF][START_REF] Sirakov | Un nouveau faciès lamellaire du début du Paléolithique supérieur dans les Balkans[END_REF].

During the last decades, different dating methods have been employed to establish reliable chronologies for the sequences in the cave. The Upper Palaeolithic sequence was dated using the radiocarbon method (Guadelli et al., 2005, in press). Moreover, preliminary luminescence ages attributed to the neanderthal radius were previously published [START_REF] Tillier | Evidence of Neanderthals in the Balkans: The infant radius from Kozarnika Cave (Bulgaria)[END_REF]. In another study, a palaeomagnetism sequence was established for stratified layers in Kozarnika [START_REF] Muttoni | An early Brunhes (<0.78 Ma) age for the Lower Paleolithic toolbearing Kozarnika cave sediments, Bulgaria[END_REF]. Finally, a timeframe for the Lower and Middle Palaeolithic assemblages was demonstrated on mammal biochronology (Guadelli et al., 2005, in press;Fernandez and Cregut-Bonnoure, 2007;[START_REF] Popov | An almost one million year long (Early to Late Pleistocene) small mammal succession from the archaeological layers of Kozarnika Cave in Northern Bulgaria[END_REF][START_REF] Sirakov | An ancient continuous human presence in the Balkans and the beginnings of human settlement in western Eurasia: A Lower Pleistocene example of the Lower Palaeolithic levels in Kozarnika cave (North-western Bulgaria)[END_REF].

In this study, we applied multiple luminescence-dating methods not only to elucidate the timeframe of humans resided during the Middle and Upper Palaeolithic but to shed light on the timing of the inhabitants at the Lower Palaeolithic and provide a full chronology covering the complete sequence through sediment samples. This will allow us to improve our understanding about the timing of humans' settlements, particularly for Kozarnika and more widely for the Balkans.

Site location

Kozarnika cave is located in the northwest of Bulgaria, 30 km away from the Serbian border (Fig. 1). The site is situated in the northern slopes of the Balkan Mountains, merging to the Danubian Plain, which is extended to the Danube River to the north. The cave is positioned at an altitude of 481 m on the mountainside of Bashovishka Bara valley. The streams in the valley drain into the Skomlia River irrigating the region, which feeds the Danube River ca 30 km away.

The first report about Kozarnika dates back to 1933 [START_REF] Popov (попов ; Попов | Пещерата Миризливка. Приносъ къмъ дилувиалната фауна и културата на дилувиалния човекъ въ България. National Archaeological Museum[END_REF] and years after in 1996, the site was excavated systematically by Bulgarian-French researchers [START_REF] Guadelli | Une séquence du Paléolithique inferieur au Paléolithique récent dans les balkans : La grotte kozarnika à oreshets (nord-ouest de la bulgarie)[END_REF]. Following years of collaborative work, 21 geological layers were discovered in a sequence of about 9 m containing assemblages attributed to the Lower, Middle, Upper Palaeolithic as well as Holocene.

Stratified geological layers including archaeological discoveries

The Bashovishka Bara valley has cut 185 m deep into the northern foothills of the Balkan through three main geological formations. The top formation belongs to the Early Jurassic which at the middle is replaced by the Middle Jurassic and finally at the bottom is superimposed on the Late Jurassic. The Kozarnika cave consists mainly of limestone, which is formed in the Upper Jurassic geological formation [START_REF] Guadelli | Une séquence du Paléolithique inferieur au Paléolithique récent dans les balkans : La grotte kozarnika à oreshets (nord-ouest de la bulgarie)[END_REF][START_REF] Sirakov | An ancient continuous human presence in the Balkans and the beginnings of human settlement in western Eurasia: A Lower Pleistocene example of the Lower Palaeolithic levels in Kozarnika cave (North-western Bulgaria)[END_REF]. The cave is part of a larger karst landscape with other caves discovered in the north-western region of Bulgaria.

The Pleistocene layers in the cave are generally well-preserved. However, we observed truncations in the cave for some layers containing Middle and Upper Palaeolithic assemblages [START_REF] Guadelli | Une séquence du Paléolithique inferieur au Paléolithique récent dans les balkans : La grotte kozarnika à oreshets (nord-ouest de la bulgarie)[END_REF][START_REF] Sirakov | An ancient continuous human presence in the Balkans and the beginnings of human settlement in western Eurasia: A Lower Pleistocene example of the Lower Palaeolithic levels in Kozarnika cave (North-western Bulgaria)[END_REF]. These truncations are likely stemming from the cave's subsequent use in modern times, when it was used as a mushroom farm and goat farm, and finally as a bunker in the 1960s.

The geological layers in the cave were categorised into four main groups covering the top to the bottom of the excavation profile (Fig. 2). Their general characteristics are explained below.

Fig. 2. 1.3 Geological layers containing the archaeological finds

Layers 3a-4. These layers with a thickness of 1 m to 1.4 m form the top of the sequence. Most of the sediment is of aeolian origin. Nevertheless, a considerable amount of calcareous clasts gravels in colours from light brown to whitish silty medium can be found in this part in addition to debris from the ceiling and the walls of the cave. The archaeological assemblages belong to the middle and late Kozarnikian tradition (the indigenous eponymous bladelets industry mentioned above) associated with the Upper Palaeolithic [START_REF] Guadelli | Une séquence du Paléolithique inferieur au Paléolithique récent dans les balkans : La grotte kozarnika à oreshets (nord-ouest de la bulgarie)[END_REF].

Layer 5a-10a. These layers show thicknesses of 1.2-1.5 m and were deposited through a combination of aeolian transport (recognised by beige silty-fine grains), and sediment fragments from the ceiling and walls (identified with coarse grains and calcareous gravel and pebble). Both groups of sediments are distributed irregularly over the section. Layer 5b contains some volcanic ashes superimposed on layers 6 and 7, which both were composed of silty sediment varying from light grey (6) to light brown (7). Due to the similarity between the origin of these two layers, they both were considered as one layer named 6/7. The calcareous gravels and charcoal also were sparsely observed in this layer [START_REF] Ferrier | Kozarnika. Etude geologique et sedimentologique[END_REF]. Generally speaking, the formation of layers 5a, 6/7, and 10a is dominated by aeolian processes. Layers 10a-9a contain archaeological elements related to the Middle Palaeolithic period. The archaeological finds in layer 6/7 are consistent with the Middle and Initial Upper Palaeolithic assemblages. The upper layer (5a) contains the blade industry known as the Kozarnikian tradition.

Layer 10b-10c. These layers originated from roughly weathered calcareous gravel and pebbles, with a combination of loamy sand with dark brown colour that could be attributed to human activity in the past. The sediment thicknesses of these layers vary from 15 cm to 50 cm To be submitted to the journal of Quaternary Science Reviews and contain archaeological artefacts associated with the very early Middle Palaeolithic with the Mousterian characteristics. Layer 11a-13. These layers were located at the bottom of the sequence having a thickness of ca 2.5 m. The sediments in this section consist of orange loamy fine and coarse fractions with flint, limestone, and debris from the ceiling and the walls of the cave. In particular calcareous pebbles and blocks with various concentrations were observed in layer11a as a common feature of this layer in a different part of the cave. The same pattern with a combination of a dark fraction of silty material with quartz origin applies for layers 11b, 11c, 12, and 13. It is likely that this sequence of layers is a result of mudslides according to [START_REF] Ferrier | Kozarnika. Etude geologique et sedimentologique[END_REF]. All of these layers comprise findings from an industry affiliated with the Lower Palaeolithic period.

Material and Methods

Luminescence dating estimates energy (radiation dose) deposited in natural minerals after the last sunlight exposure. This quantity is known as the equivalent dose (estimated in the laboratory), which is built up by energy emitted by natural radioelements in the sediment over the burial time (dose rate). Dividing the equivalent dose by the corresponding dose rate, results in the burial time (age). We used various techniques to determine the equivalent doses of the samples. First, we began with the blue optically stimulated luminescence (OSL) for quartz [START_REF] Huntley | Optical dating of sediments[END_REF], followed by post infrared stimulated luminescence at the elevated temperature ( Post-IR IRSL, [START_REF] Thomsen | Laboratory fading rates of various luminescence signals from feldspar-rich sediment extracts[END_REF] to completely cover the Middle to the Upper palaeolithic sequences. In addition, the infrared stimulated luminescence (IRSL, [START_REF] Hütt | Optical dating: K-feldspars optical response stimulation spectra[END_REF] signals for some samples were measured at 50 ˚C to crosscheck with the two other age series.

It is well-known that signal saturation is a key parameter hampering the application of luminescence dating for the Lower to Middle Palaeolithic. Given that, we attempted to tackle this problem by combining various luminescence signals with different saturation degrees to extend the age limit. We switched to violet stimulated luminescence (VSL) [START_REF] Ankjaergaard | Violet stimulated luminescence dating of quartz from Luochuan (Chinese loess plateau): Agreement with independent chronology up to ∼600 ka[END_REF][START_REF] Ankjaergaard | Towards dating Quaternary sediments using the quartz Violet Stimulated Luminescence (VSL) signal[END_REF][START_REF] Jain | Extending the dose range: Probing deep traps in quartz with 3.06eV photons[END_REF] on quartz accompanied by infrared radiofluorescence (IR-RF) (Trautmann et al., 1999b(Trautmann et al., , 1999a[START_REF] Trautmann | Investigations of feldspar radioluminescence: potential for a new dating technique[END_REF], for review see [START_REF] Murari | Infrared radiofluorescence (IR-RF) dating: A review[END_REF] on K-feldspar to estimate ages for the older samples.

The dose rates were estimated by means of both in situ and laboratory measurements for all samples (further details below).

Sample preparation

We present ages for 18 sediment samples extracted from the profile, using opaque tubes (see Table 1). Our sample preparation followed routine procedures applied in luminescence dating [START_REF] Preusser | Luminescence dating: basics, methods and applications[END_REF]. After conducting the grain size analysis (see Fig. S1 in the supplementary material) coarse grains with 180 µm to 250 m diameters were extracted using wet sieving. To eliminate carbonates, the obtained fraction of sediments was rinsed in HCl (30%) over a couple of hours until no further reaction was observed. Following that, they were soaked in H2O2 (30%) to remove organic compounds. Later, density separation was carried out through lithium hetero-polytungstates (LST) heavy-liquid. First, quartz grains were separated from the bulk sample using a density of 2.62 g cm -3 . Then, K-feldspar grains were extracted from the remaining material by applying the LST with a density of 2.58 g cm -3 . We removed the outer layers of quartz and K-feldspar grains affected by alpha particles by etching them with HF; (40% for 40 min on quartz and 10% for 10 min on K-feldspar). Etching was followed by resting the grains in HCl (15%) overnight to eliminate unwanted fluorides. Finally, all grains were sieved for the second time with the smaller mesh size (<180 µm) to remove any remaining fragmented grains.

Table 1

Instrumentation

Various instruments were employed for measuring luminescence signals using different minerals and techniques. Blue OSL was measured on Freiberg Instrument lexsyg SMART TL/OSL systems [START_REF] Richter | Lexsyg smarta luminescence detection system for dosimetry, material research and dating application[END_REF] equipped with five blue LEDs (458 Δ 3 nm, max. 100 mW cm -2 ) and five infrared LEDs (850 Δ 3 nm, max. 300 mW cm -2 ) as well as on a Freiberg Instrument lexsyg research reader (hereafter L1, [START_REF] Richter | lexsyga new system for luminescence research[END_REF] equipped with five blue LEDs (458 Δ 10 nm, max. 100 mW cm -2 ) and five infrared LEDs (850 Δ 3 nm, maximum 300 mW cm -2 ). The stimulation powers for lexsyg research and lexsyg SMART systems were set to 100 mWcm -2 and 40 mWcm -2 , respectively. Signals were detected using a UV filter combination (Hoya U340, 2.5 mm and Delta BP 365/50 EX) in both systems.

Post-IR IRSL (pIR50IRSL290) signal was measured using a blue-violet filter combination Schott BG 39 (3 mm) and Chroma D410/30x in Freiberg Instruments lexsyg research reader with a stimulation power of 100 mWcm -2 . IRSL signals were measured using Risø TL/OSL DA-15 and DA-20 readers equipped with five blue LEDs (470 Δ 20 nm, max. 105 mW cm -2 ) and five infrared LEDs (870 Δ 40 nm, max. 145 mW cm -2 ). Signals were detected through a combination of Schott BG39 and Corning 7-59 [START_REF] Bøtter-Jensen | Developments in radiation, stimulation and observation facilities in luminescence measurements[END_REF].

To be submitted to the journal of Quaternary Science Reviews VSL signals were measured using Freiberg Instruments lexsyg research reader containing violet-emitting LEDs (405 Δ 3nm, max. 100 mW cm -2 ). The stimulation power was set to 100 mWcm -2 and the signal was detected using a combination of a Hoya U340 (2.5 mm) and an AHF-BL HC 340/26 interference filter.

All systems were operated with a 90 Sr/ 90 Y-source which was calibrated using Risø calibration quartz [START_REF] Hansen | A new irradiated quartz for beta source calibration[END_REF]. The dose rates available for each sample in lexsyg research and SMART as well as Risø TL/OSL systems were ca 8 Gy min -1 , 11 Gymin -1 , and 7 Gy min -1 at the measurement dates, respectively. IR-RF was measured through another Freiberg Instrument lexsyg research system, which was equipped with a Hamamatsu H7421-50 PMT in conjunction with a Chroma D850/40 interference filter. A calibrated 90 Sr/ 90 Y-source delivered a dose rate of ca 3 Gy min -1 . For the bleaching of the natural IR-RF, we used the internal solar simulator of the machine. All samples were prepared using stainless steel cups with a diameter of 8 mm.

Approximately, a hundred grains were mounted in diameter of 3 mm using silicon oil for measuring blue OSL, pIR50IRSL290, and IRSL. The number of grains was increased for VSL and IR-RF; in these cases, the cups were fully covered by grains, representing several thousand grains.

Dose rate

The dose rate available to quartz and feldspar consists of external and internal parts due to the activity of the radioelements outside and inside of the grains. For the external part, we have only estimated beta and gamma dose rates because the etching removed the outer layer of the grains affected by alpha-radiation.

We estimated the radionuclide concentrations of U, Th, and K using high-resolution low background gamma-ray spectrometry [START_REF] Guibert | TL dating: Low background gamma spectrometry as a tool for the determination of the annual dose[END_REF]. The obtained concentrations were then turned into beta-dose rates using conversion factors after [START_REF] Guérin | Dose-rate conversion factors: update[END_REF] and corrected using grain size attenuation factors following Guérin et al. (2012). The effect of water on dose rate was corrected after [START_REF] Nathan | On the dose-rate estimate of carbonate-rich sediments for trapped charge dating[END_REF]. We estimated the saturation water content based on the grain size distribution as suggested in [START_REF] Nelson | Using grain-size characteristics to model soil water content: Application to dose-rate calculation for luminescence dating[END_REF]. The sample-specific percentages of the sand, silt, and clay fraction are displayed in Table S1 in the supplementary material. We have used the fraction of 0.5 0.2 (of the saturation value) as a representative value to approximate the past water content. We increased the obtained value to 20 8% for samples (Koz16 to Koz21) due to the mudflow origin of layer 11b downwards.

To be submitted to the journal of Quaternary Science Reviews Gamma-dose rates were estimated using artificial dosimeters (Al2O3:C chips; see Kreutzer et al., 2018 for further details) for all samples. The chips were first heated up for 10 min at 350° C. Three chips were inserted in each aluminium tube and remained in the sediment for 271 days. After retrieving, the stored energy was measured using Daybreak 2200 OSL reader system [START_REF] Bortolot | A new modular high capacity OSL reader system[END_REF]. The chips were stimulated through green light (Nichia NSPG10) and the emissions were recorded through a 7.5 mm Hoya U-340 filter. For comparison, the gammadose rate was also estimated using radioelement concentrations derived from the gamma-ray spectrometry. The concentrations were converted to the dose rate after convention factors tabulated in [START_REF] Guérin | Dose-rate conversion factors: update[END_REF] and corrected for water content following (Guérin and Mercier, 2012).

The cosmic-dose rate was calculated using the coordinate of the site and the overburden depth after the procedure in [START_REF] Prescott | Cosmic ray contributions to dose rates for luminescence and ESR dating: Large depths and long-term time variations[END_REF][START_REF] Prescott | The contribution of cosmic radiation to the environmental dose for thermoluminescence dating. Latitude, altitude and depth dependences[END_REF]. The uncertainty of the cosmic dose rate was estimated at 10% after [START_REF] Prescott | Cosmic ray contributions to dose rates for luminescence and ESR dating: Large depths and long-term time variations[END_REF].

The internal dose rate was estimated using K and Rb concentrations of 12.5 0.5% [START_REF] Huntley | The K content of the K-feldspars being measured in optical dating or in thermoluminescence dating[END_REF] and 460 20 ppm after [START_REF] Mejdahl | Internal radioactivity in quartz and feldspar grains[END_REF] respectively. The internal dose rate for quartz coarse grains was taken from [START_REF] Mejdahl | Internal radioactivity in quartz and feldspar grains[END_REF].

Luminescence signal measurements

We applied several protocols to both quartz and feldspar grains and accordingly the corresponding signals were used for dating. We have measured blue-OSL and VSL for quartz grains. In addition, IRSL at 50°C, and pIR50IRSL290 at an elevated temperature (290 °C) were measured to establish ages based on feldspar grains. Nevertheless, this approach potentially could disclose the history of grains' bleaching, which is a worthy matter of consideration since the accuracy of the obtained ages heavily depends on it. Below, the setting for each protocol is explained.

Blue optically stimulated luminescence (blue-OSL)

We used the single aliquot regenerative (SAR) dose protocol [START_REF] Murray | Luminescence dating of quartz using an improved singlealiquot regenerative-dose protocol[END_REF] to estimate the equivalent dose (De) for each sample.

First, the dominance of the fast-decaying signal favourable for relatively young sediment samples was tested. We simply compared the blue-OSL signal from one of Kozarnika's samples with calibrated quartz Risø batch 90 [START_REF] Hansen | A new irradiated quartz for beta source calibration[END_REF], which is known to be dominated by the fast decaying luminescence signal. In addition, a linearly modulated OSL (LM-OSL, [START_REF] Bulur | An alternative technique for optically stimulated luminescence (OSL) experiment[END_REF] technique was applied to deconvolute the signal components for this sample to be compared with the calibration quartz. De estimation was carried out using Analyst To be submitted to the journal of Quaternary Science Reviews software v. 4.53 [START_REF] Duller | The Analyst software package for luminescence data: overview and recent improvements[END_REF]. A single saturating exponential function was used for fitting the dose-response curve.

The regenerative dose steps were followed by a null dose delivered to monitor charge transfers. This step was followed by a repetition of the first regenerative dose to verify significant sensitivity changes. In the end, the same dose as the first regenerative was delivered to the samples and the corresponding signal was measured after stimulating with IRSL followed with blue OSL to check any significant contamination by feldspar grains after [START_REF] Duller | Distinguishing quartz and feldspar in single grain luminescence measurements[END_REF].

The luminescence signal was measured at 125°C for 40 s after preheating at 260°C for 10 s.

The preheat plateau test was applied to sample Koz2 to select an appropriate temperature before measuring the signal in SAR protocol [START_REF] Murray | Luminescence dating of quartz using an improved singlealiquot regenerative-dose protocol[END_REF]. In this test, Des were measured at various degrees starting from 200°C and reaching 300°C (ramping up by 20°C at each step). The cut heat was 40°C lower than the corresponding preheat temperature. Doserecovery tests were conducted on some samples to monitor the reliability of the SAR protocol's parameters. Samples were bleached for 100 s using blue light, followed by a pause of 10,000 s. For the second time samples were bleached using blue light and then they received a beta dose inside the OSL reader approximately equivalent to their natural doses (Koz2,4,7,10).

Infrared stimulated luminescence (IRSL50 and pIR50IRSL290)

K-feldspars were used for estimating the Des for some of the samples to obtain independent ages, also revealing the grains' bleaching histories. Particularly, K-feldspar is favourable when the expected temporal range is higher than what the blue-OSL signal of quartz can achieve.

Due to this higher upper age limit, this mineral has been used to unravel the temporal frameworks attributed to the Middle and also Lower Palaeolithic periods of the cave. Having said that, the drawback of using feldspar is athermal anomalous fading (loss of signal over time [START_REF] Wintle | Anomalous fading of thermoluminescence in mineral samples[END_REF]) that impedes the proper estimation of the true De. To overcome this problem, we determined the fading rate [START_REF] Huntley | Ubiquity of anomalous fading in K-feldspars and the measurement and correction for it in optical dating[END_REF]) and corrected its effect on De. In this study, for three samples Koz4, Koz8, and Koz11, IRSL signals were measured (six aliquots for each sample) following fading measurement after a pause of approximately 1 h and 4 h. 'Luminescence' [START_REF] Kreutzer | Introducing an R package for luminescence dating analysis Ancient[END_REF][START_REF] Kreutzer | Luminescence: Comprehensive Luminescence Dating Data Analysis[END_REF] package after methodological work by [START_REF] Huntley | An explanation of the power-law decay of luminescence[END_REF]. These values were then used to correct the De estimates using

Luminescence::calc_Huntley2006() after [START_REF] Kars | A new approach towards anomalous fading correction for feldspar IRSL datingtests on samples in field saturation[END_REF]King et al., 2018 also embedded in the R 'Luminescence' package. A single saturating exponential function was used to generate a dose-response curve.

To be submitted to the journal of Quaternary Science Reviews Furthermore, we measured the IRSL signal at an elevated temperature following the IRSL at 50°C. This method was suggested after [START_REF] Buylaert | IRSL and post-IR IRSL residual doses recorded in modern dust samples from the Chinese Loess Plateau[END_REF][START_REF] Li | Luminescence dating of K-feldspar from sediments: A protocol without anomalous fading correction[END_REF][START_REF] Thiel | Luminescence dating of the Stratzing loess profile (Austria) -Testing the potential of an elevated temperature post-IR IRSL protocol[END_REF][START_REF] Thomsen | Stability of IRSL signals from sedimentary Kfeldspar samples[END_REF][START_REF] Thomsen | Laboratory fading rates of various luminescence signals from feldspar-rich sediment extracts[END_REF] and it is proposed that the signals do not suffer from the fading, or at least the fading rate can be considered as insignificant [START_REF] Buylaert | A robust feldspar luminescence dating method for Middle and Late Pleistocene sediments: Feldspar luminescence dating of Middle and Late Pleistocene sediments[END_REF][START_REF] Buylaert | IRSL and post-IR IRSL residual doses recorded in modern dust samples from the Chinese Loess Plateau[END_REF].

The parameters' setting for measuring this signal was derived from [START_REF] Buylaert | A robust feldspar luminescence dating method for Middle and Late Pleistocene sediments: Feldspar luminescence dating of Middle and Late Pleistocene sediments[END_REF] and was applied to a large number of samples containing 4 to 12 aliquots.

Infrared radiofluorescence (IR-RF)

Aside from the IRSL signal, K-feldspar emits another signal during ionising radiation, which is detected in the infrared wavelength range, known as IR-RF (for a recent review see [START_REF] Murari | Infrared radiofluorescence (IR-RF) dating: A review[END_REF]. The advantage of employing this signal is the larger age limit, which presumably could cover the temporal range attributed to the Lower Palaeolithic sequence of the cave. For IR-RF measurements, we have employed the parameters' setting suggested by [START_REF] Frouin | An improved radiofluorescence single-aliquot regenerative dose protocol for Kfeldspars[END_REF]. The regenerative signals were measured after bleaching of ca 3 h (10,800s). To ensure this time sufficiently reset the natural IR-RF signal, we also carried out a bleaching test.

The Des for samples Koz11, Koz14, Koz16, Koz18, Koz20, and Koz21 were estimated through the approach of horizontal and vertical sliding [START_REF] Kreutzer | Always remain suspicious: a case study on tracking down a technical artefact while measuring IR-RF[END_REF][START_REF] Murari | Further investigations on IR-RF: Dose recovery and correction[END_REF] implemented in Luminescence::analyse_IRSAR.RF().

Violet stimulated luminescence (VSL)

Unlike the relatively early saturation of the fast decaying blue-stimulated luminescence signal, some slow decaying components in quartz are deemed to be of higher saturation dose [START_REF] Jain | Extending the dose range: Probing deep traps in quartz with 3.06eV photons[END_REF]. This makes them attractive for dating older sediment like the Lower Palaeolithic sequences. Thus, we have applied the SAR protocol suggested by [START_REF] Ankjaergaard | Violet stimulated luminescence dating of quartz from Luochuan (Chinese loess plateau): Agreement with independent chronology up to ∼600 ka[END_REF] to estimate the De based on deep traps. This protocol was designed in a way that first the fastdecaying component is depleted using blue light stimulation and violet stimulation is then applied to probe the slow decaying signal with a higher saturation dose. We tested this technique on several samples (Koz8, 14, 16, 18, 21). However, the outcome of this application was not promising (all signals were in apparent saturation), therefore, we proceeded with the multiple aliquot additive dose (MAAD) protocol as an alternative solution to the standard SAR protocol offered in [START_REF] Ankjaergaard | Violet stimulated luminescence dating of quartz from Luochuan (Chinese loess plateau): Agreement with independent chronology up to ∼600 ka[END_REF]. We used one of the youngest samples (Koz4)

for which an independent age using the blue OSL technique was estimated and for which we had enough grains available to prepare a large number of aliquots.

The prepared samples were categorised into several fractions (each group comprising six aliquots) which were exposed to various amounts of beta irradiations on top of the natural dose.

The first irradiative dose was set to 117 Gy and gradually it was increased to the dose of 2760

Gy. In the next step, following the approach of [START_REF] Ankjaergaard | Violet stimulated luminescence dating of quartz from Luochuan (Chinese loess plateau): Agreement with independent chronology up to ∼600 ka[END_REF] we measured only the natural luminescence signal (Ln) which was corrected by the test dose signal (Tn) for each set of irradiated samples to minimise the impact of sensitivity change arising from repeated irritations and heating of the samples. Then the obtained L/T values were subsequently used to construct the MAAD dose-response curve fitted using a Lambert W function [START_REF] Pagonis | A new analytical equation for the dose response of dosimetric materials, based on the Lambert W function[END_REF]. Afterwards, the Ln/Tn ratios measured from the natural samples were used (six aliquots for each natural sample) to estimate the Des (see result section).

Results

Dose rate

Table 2 shows the obtained concentrations of U, Th, and K. The equivalent U concentration was determined from the top (pre-226 Ra) and the bottom of the chain (post-226 Ra). This allows checking potential disequilibria in the U chain [START_REF] Guibert | The importance of U-series disequilibrium of sediments in luminescence dating: A case study at the Roc de Marsal Cave (Dordogne, France)[END_REF]. For samples Koz6, Koz7, Koz18, radioactive disequilibria were observed. To minimise this effect on the dose rates, instead of using the obtained values from the bottom of the chains (which were smaller than the top vaules), the average of the top and bottom concentrations was calculated and employed to estimate the dose rates for those samples (Table 3).

Table 2 The beta-dose rate for all samples ranged from 0.86 0.04 Gy/ka to 1.70 0.08 Gy/ka. The in situ and laboratory (i.e. from the concentrations of K, U and Th) gamma-dose rates ranged from 0.19 (<0.01) Gy/ka to 0.72 0.01 Gy/ka and from 0.49 0.02 Gy/ka to 0.97 0.03 Gy/ka, respectively. We have illustrated the two distinct series of gamma dose rates in Fig. 3. As shown in this figure, most in situ gamma-dose rates were significantly smaller than the corresponding gamma-dose rates obtained through means of gamma-ray spectrometry. It should be noted that we did not have a reliable estimation of the water content at sampling data, thus, it is likely that our estimation of the past water contents (varying between 14 6% to 20 8 %) is not in close alignment with the current water content. However, this does not explain the substantial discrepancies observed between the two dose rate results. It is more likely that the limestone pebbles and gravels that can readily identify in Fig. 2 Particularly, because they cover gamma rays from larger distances and hence provide more

To be submitted to the journal of Quaternary Science Reviews accurate estimates. The internal dose rate for quartz was set to 0.06 0.03 Gy/ka and a value of 0.90 0.14 Gy/ka was obtained for K-feldspar. Cosmic dose rates varied between 0.06 0.01 Gy/ka to 0.09 0.01 Gy/ka. In the end, the final dose rate for quartz and K-feldspar grains vary between 1.21 0.05 Gy/ka and 2.38 0.09 Gy/ka, and between 2.05 0.15 Gy/ka and 3.22 0.16

Gy/ka, respectively.

Table 3 Fig. 3.

Equivalent dose (De) determination

Blue-OSL Fig. S2 in the supplementary material illustrates with pseudo-LM OSL [START_REF] Bos | How to visualize quartz OSL signal components[END_REF] the presence of a fast-decaying signal at the beginning of the curve, but it also highlights at least contamination with another signal, at the end of the curve. Fig. S3 (true LM-OSL) revealed that presumably, the slow decaying signals (slow components) are also contributed to the luminescence signal.

We employed late background subtraction; the signal was measured using channels 1 to 3 (0.3 s) and the background was set to 300 to 399 (10 s). The recycling ratios, (a representative value for expressing the sensitivity change after administrating the cycles of heating and irradiating) were within 20% of the unity for all samples. Furthermore, the recuperation signal ratio following the received zero dose was less than 5%. An IR depletion ratio test [START_REF] Duller | Distinguishing quartz and feldspar in single grain luminescence measurements[END_REF] was carried out at the end of the SAR protocol and the results were within 10% for all samples. Fig. S4 in the supplementary material depicts the outcome of the preheat test for sample Koz2.

We considered 260 ° C and 220° C as acceptable temperatures for the preheat and cut-heat in our SAR protocol; using such elevated temperatures allows subtracting a late background by removing less stable components than the fast component [START_REF] Jain | Characterisation of blue-light stimulated luminescence components in different quartz samples: implications for dose measurement[END_REF]. Dose recovery ratios for samples Koz2, Koz4 (given dose ca 48 Gy), Koz7 and Koz10 (given dose ca 83 Gy) resulted in 1.02 0.04, 1.01 0.03, 0.96 0.08, 0.91 0.07, respectively. These values are consistent with 1 within uncertainties and fulfil the standard criterion (commonly accepted is within 10% of unity), indicating the suitability of the SAR protocol's parameters for our samples. However, the smooth reduction in these values alongside increasing the Des, indicates that reaching the non-linear part of the dose-response curve and may give rise to De underestimation.

Fig. 4a and Fig. 5a show typical dose-response curves, and a visual representation of the De distribution for sample Koz4 as an example. We have estimated the central dose for each To be submitted to the journal of Quaternary Science Reviews sample using the arithmetic average (as suggested for multi-grain OSL by [START_REF] Guérin | Absorbed dose, equivalent dose, measured dose rates, and implications for OSL age estimates: Introducing the Average Dose Model[END_REF]), and its uncertainty was determined through the standard error of the mean. The estimated central doses ranged from 48 2 Gy to 117 9 Gy for samples Koz1 to Koz10 (see Table 4). We also measured blue OSL signals for subsequent samples, however, since we reached the saturation limit of the signal, further results were considered untrustworthy or at least underestimated. The pIR50IRSL290 signals were measured using the first 8 s and the background was subtracted using the last 60 s of the shine-down curves. The recycling ratios were within 10% of unity and recuperation ratios were less than 5%. The dose-response curves and the De distribution for sample Koz4 are depicted in Fig. 4b and Fig. 5b, respectively. The De representing the central equivalent dose for all samples was determined using the arithmetic average of De values and its uncertainty was estimated using the standard error of the mean. The obtained Des for our samples varied between 95 3 (Koz2) to 671 36 Gy (Koz14) (see Table 4). Due to the lack of K-feldspar grains for Koz15, we could not measure any signal, however, for the subsequent samples, the corresponding signals were mainly in or close to saturation and thus the outcome was not considered trustworthy.

Infrared-radiofluorescence (IR-RF)

The result of the bleaching test on sample Koz11 is displayed in Fig. S6. This figure shows that a bleaching time of 10,800 s appears to be sufficient, lying within the plateau. The IR-RF signal was measured using a channel resolution of 10 s. The observed scatter in the points after reaching the plateau is probably induced by negligible movement of grains on the sample To be submitted to the journal of Quaternary Science Reviews carrier. Typical IR-RF signals for samples Koz11, Koz18, Koz21 are shown in Fig. 6. The natural and regenerative signals are illustrated with red and green, respectively. The horizontal and vertical sliding method was successfully applied to natural signal and an excellent match between natural signal and regenerative is observed. The De distributions for samples Koz11, Koz18, and Koz21 are represented in Fig. 7 using kernel density estimate (KDE) plot. The De representing the central equivalent dose for each sample is estimated using the arithmetic average and its uncertainty was determined using the standard deviation. Fig. 6. Fig. 7.

Violet stimulated luminescence (VSL)

A typical dose-response curve alongside the natural and regenerative signals observed after applying the VSL (SAR) protocol is shown in Fig. S7 for one sample. Fig. S8 shows a smooth increase in the obtained Ln/Tn values as samples' depths are increased along with the profile. This gradual rise suggested that VSL has the potential to yield chronological information, as was already observed in [START_REF] Ankjaergaard | Violet stimulated luminescence dating of quartz from Luochuan (Chinese loess plateau): Agreement with independent chronology up to ∼600 ka[END_REF]. Fig. 8 illustrates the L/T distribution obtained for each set of six aliquots prepared originally from Koz4; while the first set corresponds to the natural dose, all other sets were irradiated in the laboratory with increasing doses on top of the natural. We tested several fitting functions on our data set. Finally, as depicted in Fig. 8, we fitted a dose-response curve using a Lambert W function to our dataset, which seems to be more justified than other commonly fitting dose-response curve functions [START_REF] Pagonis | A new analytical equation for the dose response of dosimetric materials, based on the Lambert W function[END_REF]. Subsequently, the Ln/Tn obtained from the natural samples were projected onto the fitted curve (Fig. S9). Because of the nature of the MAAD approach, which adds doses on a natural sample that already has a dose, the dose-response curve necessarily intersects with the y-axis. Theoretically, one could attempt to estimate the equivalent dose for the first sample from the intersection of the curve with the x-axis. Due to the exponential character of the curve, this extrapolation comes at the cost of very large uncertainties; hence it was not used here. However, the Ln/Tn signal of Koz4 project onto the MAAD dose-response curve should result in an equivalent dose around 0. We estimated a value of 11 12 Gy for this sample, which is consistent with the zero dose within uncertainty. To estimate the equivalent dose of all other samples, the natural equivalent dose of sample Koz4 based on quartz blue OSL was added to all estimated results to determine the final De for each sample.

To be submitted to the journal of Quaternary Science Reviews Fig. 8.

Age estimation

The obtained ages based on the various signals are presented in Table 4. The OSL dating for the uppermost sample (Koz1) belonging to layer 5a, exhibited an age of 32 2 ka (one sigma).

Gradually, the OSL ages increase smoothly with depth until layer 10b (sample Koz11), for which two aliquots out of ten were in saturation. Thus, we considered this age (141 12 ka) to be a minimum age. The corresponding ages for samples Koz11 and Koz12 estimated based on pIR50IRSL290 signals led to 161 13 ka and 199 12 ka (one sigma). Subsequently, only the OSL age for layer 10a (59 4 ka) and younger samples above were considered reliable (Fig. 9).

As shown in Fig. 9 (inset) (between layer 5a to layer 10a), the pIR50IRSL290 ages were determined for a number of samples. Generally, at first glance, it can be concluded that the pIR50IRSL290 ages are older than the corresponding OSL ones. However, apart from one sample (Koz2 from layer 5b), for which no overlapping between the OSL and pIR50IRSL290 ages is observed within two sigma (26-33 ka for OSL and 35-43 ka for pIR50IRSL290), the rest overlap and the age outcome from the two approaches can be generally evaluated consistently.

The IRSL ages which were fading-corrected were estimated for three samples. All of them are older compared to both, OSL and pIR50IRSL290 ages. The obtained IRSL ages resulted in 45 3 ka and 66 4 ka, respectively, for samples Koz4 (layer 6/7) and Koz8 (layer 9c). At 95% confidence, these two ages exhibit timeframes of 39-51 ka and 58-75 ka that overall are in alignment with the periods of 33-41 ka and 48-74 ka ascertained from pIR50IRSL290.

Nevertheless, the IRSL age of sample Koz11 (layer 10b) 250 40 ka is significantly larger than the corresponding pIR50IRSL290 age of 161 13 ka. The IR-RF age for this sample resulted in 201 17 ka; thus, it seems that the pIR50IRSL290 age for this sample could be underestimated.

The pIR50IRSL290 ages increase in depth up to layer 10c. Afterwards, we have a clear indication of signal saturation (four aliquots out of six) for sample Koz16 (layer 11b). From this layer down to the lowermost layer, only IR-RF and VSL ages are available.

Besides the VSL ages for the oldest part of the sequence, several ages were also determined for the youngest samples. Although at first glance, the VSL ages for the subsequent samples Koz5, Koz6, and Koz7 do not seem consistent with the outcome of other dating techniques, their corresponding age intervals overlap one another (at two sigma).

Getting back to the older part of the sequence, we can see an excellent match between IR-RFand VSL-based results, 279 29 ka and 324 94 ka, respectively for sample Koz14, layer 10c as an example. This picture does not hold for the older ages, however, still within two sigma, both ages can be considered consistent, given the large error inherent for the VSL ages observed for all samples and related to the shape of the dose-response curve.

The IR-RF ages are more precise (i.e., they show smaller analytical error bars) and are less scattered. The VSL age estimated for the oldest sample Koz21 belonging to layer 13c in the profile was estimated at 704 168 ka. This age seems to be greater than the corresponding age based on IR-RF (387 28 ka), although the large uncertainty obtained for the VSL age, makes a reasonable comparison difficult. Table 4 Fig. 9.

Discussion

How confident are we about the obtained chronologies?

In the previous section, we presented ages based on different luminescence dating techniques.

Since the age discrepancies observed between the various signals were not always negligible, it is worthwhile to critically reflect the ages and discuss reasons for these differences.

The IRSL ages estimated for all three samples led to older ages compared to other results. [START_REF] King | Age determination using feldspar: Evaluating fading-correction model performance[END_REF] argued that employing a single saturating exponential fit can result in overestimated De dependent on the characteristic of the dose-response curve. Nevertheless, it should be noted that a single saturating exponential function perfectly described the doseresponse curve for these three samples in this study. A general-order kinetic (GOK, [START_REF] Guralnik | OSLthermochronometry of feldspar from the KTB borehole, Germany[END_REF] also was tried for fitting sample Koz11 for which the De was not located in the linear part of the dose-response curve, and hence varying fitting functions could potentially lead to different outcomes. Applying this fitting resulted in a larger value compared to the single-saturating exponential function.

It is worth mentioning that the GOK can result in more reliable estimation only if the higher part of the dose-response curve is thoroughly defined by incorporating relatively large regenerative points [START_REF] King | Age determination using feldspar: Evaluating fading-correction model performance[END_REF]. The largest regenerative dose used for generating the dose-response curve for Koz11 was only 600 Gy, while the individual De s of ca 600-700 Gy were obtained for this sample. This could justify why applying GOK resulted in a larger De, and thus it was not a suitable function for our dataset. Having said that, even if the shape of the dose-response curve is considered problematic for Koz11, for younger samples, Koz4 and Koz8 this is not the case. Therefore, we can argue that it is more likely that the obtained gvalues are not displaying the true measures and they are possibly overestimated. The proper estimation of the fading rate in percentage per decade (the g-value) is a challenging task and minute variations in temperature and conditions can result in a very different estimate (cf. [START_REF] Kadereit | A closer look at IRSL SAR fading data and their implication for luminescence dating[END_REF] and the related public discussion). Although the reason behind this is not clear, it is likely that our estimation of the g-value, which is based on maximum pausing of only 4 h, was insufficient to properly determine the fading rate in the laboratory.

Previously in section 3.3, it was stated that pIR50IRSL290 ages for the upper part of the sequence seemed to be slightly older than the OSL ages. However, the pIR50IRSL290 signal may contain a residual dose (hard to bleach signal, [START_REF] Buylaert | IRSL and post-IR IRSL residual doses recorded in modern dust samples from the Chinese Loess Plateau[END_REF] contributing to the measured signal and giving rise to overestimated age. Due to the scarce amount of K-feldspar grains, we did not measure residual dose for our samples, however, in the literature, the amount of 5-20 Gy (e.g., [START_REF] Buylaert | IRSL and post-IR IRSL residual doses recorded in modern dust samples from the Chinese Loess Plateau[END_REF]) was suggested as a reasonable measure. We calculated the residual doses for the upper part of the sequence assuming that the OSL and pIR50IRSL290 ages agree with one another. The values ranged from 9 Gy to 29 Gy, consistent with the values suggested by [START_REF] Buylaert | IRSL and post-IR IRSL residual doses recorded in modern dust samples from the Chinese Loess Plateau[END_REF]. Subtracting this range from the obtained pIR50IRSL290 De resulted in the ages in excellent accordance with the corresponding quartz OSL ages. Thus, the discrepancies observed earlier between these two series of ages might be related to 'unbleachable' residual doses.

Furthermore, taking into account the slow beaching rate of the pIR50IRSL290 signal compared to the OSL signal [START_REF] Murray | Identifying wellbleached quartz using the different bleaching rates of quartz and feldspar luminescence signals[END_REF], we can safely conclude that the sediment samples were well-bleached before burial. This inference is also supported from the sedimentological perspective since the original nature of the sediments in the cave was identified as loessderivates and wind-blown loess [START_REF] Muttoni | An early Brunhes (<0.78 Ma) age for the Lower Paleolithic toolbearing Kozarnika cave sediments, Bulgaria[END_REF]. Although some scholars have argued that these sediments should be categorised as loess-like and not (primary) loess [START_REF] Jipa | The loess-like deposits in the Lower Danube basin. Genetic signif-icance[END_REF].

In addition, the bleaching rate of IR-RF is also considered similar to pIR50IRSL290 (Frouin et al., 2017), hence the IR-RF ages do not suffer from an age overestimation as already observed by an overall concordance between the pIR50IRSL290 and IR-RF ages.

The OSL and pIR50IRSL290 ages increase gradually along the profile until layer 10a.

Afterwards, due to the saturation, no OSL age was reported. [START_REF] Heydari | OSL signal saturation and dose rate variability: Investigating the behaviour of different statistical models[END_REF] demonstrated that the Bayesian models implemented in BayLum [START_REF] Guérin | Towards an improvement of optically stimulated luminescence (OSL) age uncertainties: modelling OSL ages with systematic errors, stratigraphic constraints and radiocarbon ages using the R package BayLum[END_REF][START_REF] Philippe | BayLum -An R package for Bayesian analysis of OSL ages: An introduction[END_REF] can lead to a more accurate dose estimation when the signal is close or in saturation, compared to 'classical analyses' using the Analyst software. Thus, in the future, it might be worth applying Bayesian modelling on single grain measurements for samples older than layer 10a and comparing the outcomes with other independent ages available based on different signals.

The pIR50IRSL290 ages for layer 10b resulted in 161 13 ka (Koz11) and 199 12 ka (Koz12).

The corresponding ages for sample Koz11 using IRSL, VSL and IR-RF yielded 250 40 ka, 321 77 ka, 201 17 ka, respectively, all older than pIR50IRSL290. Thus, the pIR50IRSL290 age for sample Koz11 appears to be underestimated. Albeit the first argument which immediately springs to mind is the possibility of fading, such conjecture was not derived from the rest of the samples, and it does not appear plausible here.

For the next layer (10b), the pIR50IRSL290 age (199 12 ka) with a De of 416 16 Gy is in excellent accordance with the corresponding age using the IR-RF signal (201 17 ka). The same trend is observed for layer 10c where pIR50IRSL290 resulted in 281 19 ka whereas IR-IF exhibited 279 29 ka.

The pIR50IRSL290 ages estimated for layer 10c (De of 659-671 Gy) were the oldest age obtained using this signal. Afterwards, for the lower layer (11b), measurements clearly show evidence for signal saturation (four aliquots out of six). On the other hand, IR-RF measurement for layer 11b resulted in De of 888±23 Gy. Therefore, the interval between 671±36 Gy (the highest De estimated based on pIR50IRSL290) and 888±23 Gy can be considered as a saturation level for the pIR50IRSL290 signal which is consistent with other findings (e.g., [START_REF] Yi | High resolution OSL and post-IR IRSL dating of the last interglacial-glacial cycle at the Sanbahuo loess site (northeastern China)[END_REF]. Although in the literature ca. 1,000 Gy was suggested as a threshold for IR-RF [START_REF] Erfurt | IRSAR -A single-aliquot regenerative-dose dating protocol applied to the infrared radiofluorescence (IR-RF) of coarse-grain K-feldspar[END_REF], comparing this age with the VSL might imply that the signal saturation for IR-RF is around 923±23 Gy (De of oldest sample in layer 13c). By looking at the chronological results from layer 11b to the lowermost layer (13c), it is not inconceivable to speculate that we already gradually reached a saturation level IR-RF signal at even above layer (11b) with the De of 888±23 Gy. Comparing VSL ages determined for these layers with IR-RFs, could even reinforce this possibility although the large errors presented for VSL ages do not allow for a very precise and decisive conclusion.

More importantly, [START_REF] Muttoni | An early Brunhes (<0.78 Ma) age for the Lower Paleolithic toolbearing Kozarnika cave sediments, Bulgaria[END_REF] established a palaeomagnetism sequence for the Kozarnika cave suggesting the layer immediately after 13c lies in the Brunhes-Matuyama boundary with the timing of 780 ka. They also attributed layers 13a-c to layers L6 to L7 in the stratigraphical Danube loess-palaeosol presented in [START_REF] Marković | Danube loess stratigraphy -Towards a pan-European loess stratigraphic model[END_REF] (their Fig. 16) and the link allows them to speculate that a period of 600-750 ka (early Brunhes) might represent the timeframe for layers 13a-c. Indeed, the silt fraction becomes more and more abundant from layers 11a to 13b, and the Kozarnika area is included in Pontic East European loess (domain V) (see [START_REF] Lehmkuhl | Loess landscapes of Europe -Mapping, geomorphology, and zonal differentiation[END_REF]. On the other hand, there is no loess near Kozarnika (see [START_REF] Evlogiev | Evidence for the Aeolian origin of loess in the Danubian Plain[END_REF]. The closest loess deposit is reported in distances from 60 km to 100 km of the cave location. Furthermore, silt remobilized in the karst might be considered loess-like (alluvial and aeolian), which extends more in the west than the loess (see for example Jipa, 2014 or the "lössartige und andere Lehme" of [START_REF] Fotakieva | Der Löß in Bulgarien[END_REF]. The sediment accumulation of the lower layers of Kozarnika is endocarstic and therefore discontinuous, with gaps posing difficulties to any quantification. Therefore, we cannot invalidate the hypothesis of sedimentary lacks under the magnetic "event" at the lower part of layer 11b (around -6.50 m) as identified by Kovacheva et al. (in Guadelli et al., 2005) and [START_REF] Muttoni | An early Brunhes (<0.78 Ma) age for the Lower Paleolithic toolbearing Kozarnika cave sediments, Bulgaria[END_REF]. This might place the Brunhes-Matuyama limit at the base of layer 11b and not in layer 13, which is in agreement with the macro-mammals assemblages of layer 11b (Guadelli et al., 2005, in press;[START_REF] Fernandez | Les Caprinae (Rupicaprini, Ovibovini, Ovini et Caprini) de la séquence pléistocène de Kozarnika (Bulgarie du Nord): morphométrie, biochronologie et implications phylogéniques[END_REF][START_REF] Sirakov | An ancient continuous human presence in the Balkans and the beginnings of human settlement in western Eurasia: A Lower Pleistocene example of the Lower Palaeolithic levels in Kozarnika cave (North-western Bulgaria)[END_REF].

In either way, this independent chronological data (palaeomagnetism sequence) can be seen as a strong argument supporting the saturation assumption of IR-RF for layer 13. The VSL age of 704 168 ka for layer 13c, agrees with the period of 600-750 ka, however, only due to the very large uncertainty. Apart from that, the agreement of the VSL ages with other independent ages obtained for the Kozarnika sequence tends to confirm the overall accuracy of the VSL ages, which is consistent with the observation of [START_REF] Ankjaergaard | Violet stimulated luminescence dating of quartz from Luochuan (Chinese loess plateau): Agreement with independent chronology up to ∼600 ka[END_REF].

These large uncertainty margins, resulting from the leverage effect of the signal uncertainty caused by the flat dose-response curve, are the main drawback of the VSL application.

However, this might be ameliorated by increasing the number of measured samples to specify a more precise probability density distribution after projecting Ln/Tns onto the dose-response curve. This may lead to more confidence in the De estimates in future studies.

Unravelling the timing of humans in the cave

The Kozarnika cave contains three main assemblages associated with the Lower, Middle, and Upper Palaeolithic periods. For the Upper Palaeolithic sequence, some radiocarbon dates are available (Guadelli et al., 2005, in press). Although some of these results were provided nearly two decades ago (new procedures have been suggested since then for sample pre-treatments to improve the accuracy and precision of radiocarbon dates e.g., [START_REF] Wood | From revolution to convention: the past, present and future of radiocarbon dating[END_REF], comparing these dates as independent ages with luminescence-based ages is still advantageous. The uppermost luminescence-age based on OSL technique resulted in 32 2 ka for layer (5a), and 30 2 ka for the layer below (layer 5b). The available radiocarbon age for layer 5b exhibited 31.2 ka cal. BP (original data from [START_REF] Guadelli | Une séquence du Paléolithique inferieur au Paléolithique récent dans les balkans : La grotte kozarnika à oreshets (nord-ouest de la bulgarie)[END_REF] and calibrated from [START_REF] Muttoni | An early Brunhes (<0.78 Ma) age for the Lower Paleolithic toolbearing Kozarnika cave sediments, Bulgaria[END_REF] showing an excellent concordance with OSL ages within uncertainties. New radiocarbon ages To be submitted to the journal of Quaternary Science Reviews fall in the range of 29,030 ka cal. BP to 31,130 ka cal. BP for layer 5a and 36,276 ka cal. BP to 36,875 ka cal. BP for layer 5b (Guadelli et al., in press).

Following that, layer 5c represents the radiocarbon range of 40-44 ka cal. BP (original data from [START_REF] Guadelli | Une séquence du Paléolithique inferieur au Paléolithique récent dans les balkans : La grotte kozarnika à oreshets (nord-ouest de la bulgarie)[END_REF]. These results agree with the new radiocarbon ages suggesting a range of 39,998-41,592 ka cal. BP to 42,085-44,541 ka cal. BP (Guadelli et al., in press).

Although there is no luminescence age available for this layer, the age of the layer immediately beneath (6/7) yielded 33 2 ka and 35 2 ka. Sequentially the radiocarbon dates for this layer ranged from 42,025 ka cal. BP to 43,492 ka cal. BP (Guadelli et al., in press) which is significantly older than the OSL ages.

The OSL ages for the Upper Palaeolithic assemblages (Kozarnikian) uncovered in layer 5a is in good concordance with the radiocarbon age, however the discrepancies between the OSL and radiocarbon ages for layers 5b and 5c (Guadelli et al. in press) are significant. One possible scenario that could explain the discrepancies between the outcome of two distinct dating methods could be the scatter between OSL Des for each sample (see Fig. 5a as an example). In such a case implementing signal grain measurements which allows for properly applying a broad body of statistical models could be preferable.

Layer 6/7 with OSL ages of 33±2 ka and 35±2 ka and radiocarbon age range of 42-54 ka cal.

BP was considered an intermediate level between the Middle and Upper Palaeolithic (Kozarnikian) containing artefacts associated with both affinities. Still, this does not necessarily imply the mixture of Middle and Upper affinities but rather a concatenation of successive and/or alternately rapid occupations (see Guadelli et al., in press for more details).

Although such a situation seems infrequent, it has been already observed for example at Mandrin cave (France) or at Buran Kaya III (Crimea) (Silmak et al., 2017).

The layer after 9ab is attributed to the termination of the Middle Palaeolithic period rendering the Mousterian Levallois elements. The OSL age of 43±3 ka established for layer 9ab within uncertainties agrees with the range of 41,030-39,260 ka cal. BP representing the termination of the Mousterian industry after [START_REF] Higham | The timing and spatiotemporal patterning of Neanderthal disappearance[END_REF]. [START_REF] Heydari | Do Bayesian methods lead to more precise chronologies? 'BayLum' and a first OSL-based chronology for the Palaeolithic open-air site of Mirak (Iran)[END_REF] used Bayesian modelling specifically designed for treating OSL samples to reduce age uncertainties attributed to Palaeolithic assemblages. In particular, they combined OSL ages with radiocarbon dates for one Palaeolithic site and obtained chronologies that were more precise [START_REF] Heydari | Bayesian luminescence dating at Ghāre Boof, Iran, provides a new chronology for Middle and Upper Paleolithic in the southern Zagros[END_REF].

Thus, applying Bayesian modelling to treat luminescence ages in combination with radiocarbon dates seems advantageous in this context to improve the precision for the continuous or discontinuous transition timeframe from the Middle to Upper Palaeolithic period.

A remarkable finding in layer 10b was a radius of an infant attributed to a Neanderthal, which earlier was dated between 128 ± 13 ka and 183 ± 14 ka using OSL and pIR50IRSL290 techniques [START_REF] Tillier | Evidence of Neanderthals in the Balkans: The infant radius from Kozarnika Cave (Bulgaria)[END_REF]. Nevertheless, in the present study, we acknowledge that the OSL age of 128 ± 13 ka is underestimated, and we suggest an age of 201±17 ka, based on IR-RF, revealing older deposit of the Neanderthal infant's skeleton. This result confirms the new biochronological interpretation of the middle sequence of Kozarnika cave (Guadelli et al., in press). Nevertheless, a stark contrast for the oldest layers exists between the biochronological framework according to micromammals [START_REF] Popov | An almost one million year long (Early to Late Pleistocene) small mammal succession from the archaeological layers of Kozarnika Cave in Northern Bulgaria[END_REF] and the one based on large mammal species [START_REF] Guadelli | Une séquence du Paléolithique inferieur au Paléolithique récent dans les balkans : La grotte kozarnika à oreshets (nord-ouest de la bulgarie)[END_REF][START_REF] Fernandez | Les Caprinae (Rupicaprini, Ovibovini, Ovini et Caprini) de la séquence pléistocène de Kozarnika (Bulgarie du Nord): morphométrie, biochronologie et implications phylogéniques[END_REF][START_REF] Sirakov | An ancient continuous human presence in the Balkans and the beginnings of human settlement in western Eurasia: A Lower Pleistocene example of the Lower Palaeolithic levels in Kozarnika cave (North-western Bulgaria)[END_REF].

The Lower Palaeolithic industries were found from layer 11a downwards. The luminescence dating suggested the onset of the Middle Palaeolithic period in the timeframe of 279±29 ka (IR-RF) 324±94 ka (VSL), which is perfectly in line with what we expect from other chronologies established for Middle Palaeolithic assemblages [START_REF] Richter | When Did the Middle Paleolithic Begin?[END_REF].

The youngest ages attributed to the Lower Palaeolithic assemblages in layer 11b were 337±23 ka (IR-RF), and 419±137 ka (VSL) and the oldest ages for layer 13c exhibited 387±28 ka (IR-RF) and 704±168 ka (VSL). Although here luminescence dating is nearly at its age limit, taking into account the smooth increase in the obtained ages along with the profile, does not seem to support a period of 1,400-1,600 ka estimated using large mammal biochronology [START_REF] Sirakov | An ancient continuous human presence in the Balkans and the beginnings of human settlement in western Eurasia: A Lower Pleistocene example of the Lower Palaeolithic levels in Kozarnika cave (North-western Bulgaria)[END_REF] for layers 13a-c or even the revisited biochronological estimation of 1,200/1,000-1,400/1,200 ka for the layer 12 to layers 13a-c suggested by Guadelli et al. (in press).

Surprisingly, the maximum age of 780 ka also questions the biochronology result for the Lower Palaeolithic sequence of the Kozarnika cave (Brunhes-Matuyama boundary) proposed for layer 13c using palaeomagnetic dating [START_REF] Muttoni | An early Brunhes (<0.78 Ma) age for the Lower Paleolithic toolbearing Kozarnika cave sediments, Bulgaria[END_REF]. At this stage, it seems that further chronological work is required to reach a decisive conclusion for the oldest layers of the sequence.

Conclusions

A complete chronology for the Kozarnika cave was established using various luminescencedating techniques. Combining different approaches allowed us to successfully date sediments from ca 30 ka to 700 ka unravelling the timing of human occupations during a large part of the Palaeolithic. Nevertheless, at this stage, we cannot rule out the possibility of underestimated timeframes for the layers containing the Lower Palaeolithic assemblages.

The obtained ages rely on the measurement of different signals for each sample to increase our confidence about the reliability of the outcomes. The overall accordance between the OSL with the pIR50IRSL290 ages for the upper part of the sequence makes it likely that both signals were completely bleached before mineral deposition. The same conclusion was deduced for the IR-RF signal, negating the possibility of age overestimation due to the slow beaching rate that often impedes the application of this method.

This study indicates that the saturation dose for the pIR50IRSL290 signal (at around 700 Gy to 900 Gy) led to ages of 291 22 ka and 281 19 ka as the oldest ages obtained in Kozarnika using this signal. On the other hand, the IR-RF saturated later, which enabled the age limit to be extended until 387 28 ka for this sequence. Moreover, the IR-RF doses measured on samples from layers 11a to 13c are indistinguishable, which probably indicates signal saturation. In that case, a dose of ~900 Gy might correspond to the field saturation limit of IR-RF for our samples.

The IRSL ages were deemed to overestimate other ages presumably due to an inaccurate gvalue estimate. Applying the VSL technique in combination with MAAD protocol resulted in ages that were consistent with other signal-based ages and implied the trustworthiness of outcomes. The drawback of this application was large uncertainty margins, which could be ameliorated by increasing the number of measured samples, and hence improving the precision.

More crucially, [START_REF] Muttoni | An early Brunhes (<0.78 Ma) age for the Lower Paleolithic toolbearing Kozarnika cave sediments, Bulgaria[END_REF] assigned the Brunhes-Matuyama reversal (780 ka) to the layer beneath 13c using Palaeomagnetic dating. Subsequently, a timeframe of 600-750 ka was suggested for layers 13a-c [START_REF] Muttoni | An early Brunhes (<0.78 Ma) age for the Lower Paleolithic toolbearing Kozarnika cave sediments, Bulgaria[END_REF] by correlating layers 13a-c to the Danube loess sequence layers L6-L7 [START_REF] Marković | Danube loess stratigraphy -Towards a pan-European loess stratigraphic model[END_REF]. Thus, the obtained VSL age for layer 13c within error is consistent with that period. However, the large uncertainty on this VSL age renders a precise interpretation difficult. In addition, the independent age available for layers 13a-c (with the range of 600-750 ka) also might indicate that the IR-RF ages represent the minimum age for these layers.

Employing various luminescence signals in this study led to determining the timing of the human inhabitants of Kozarnika cave. The Kozarnikian assemblages discovered in layers 5a-c dates back to 32 2 ka and 33 2 ka. Generally, radiocarbon dates (Guadelli et al., 2005, in press) provided for Kozarnikian assemblages in layers 5a-c and 6/7 are older compared to luminescence dating. This study revealed that the onset of the Middle Palaeolithic period happened at 279±29 ka (IR-RF) 324±94 ka (VSL) and it lasted until between 49±4 ka and 43±3 ka. The lower Palaeolithic assemblages in Kozarnika cave dated back to between 337±23 ka (IR-RF) and To be submitted to the journal of Quaternary Science Reviews 704±168 ka (VSL). However, due to the supposed limit of luminescence dating in this dose range (setting our upper datable age limit), it is likely that they present minimum ages. Therefore, further chronological investigations using dating methods with extended age limit are required to shed light on the reliability of these ages and to further constrain the timing of the humans who lived during the Lower Palaeolithic period in the Kozarnika cave. Fig. 3. Comparison between the in situ and laboratory (gamma-ray spectrometry) gamma dose rate estimates. The in situ results are significantly lower compared to the laboratory measurements due to the heterogeneous environment observed in the sediment layers. Since limestone fragments from the bedrock are the source of radiation heterogeneities (limestone is generally poor in radioelements), this systematic discrepancy is expected. Thus, in situ measurements are more representative of the sediment burial conditions. Fig. 4. Illustrating the dose-response curves for the OSL, pIR50IRSL290, and IRSL measurements derived from sample Koz4. Fig. 5. De distributions generated from OSL (A) pIR50IRSL290 (B) and IRSL50 (C) measurements for sample Koz4 as an example using Abanico plot [START_REF] Dietze | The abanico plot: Visualising chronometric data with individual standard errors[END_REF]. For the fitting, we applied a dose-response curve function based on the Lambert W function. The value of 11 12 Gy was obtained using this fitting for sample Koz4 (not irradiated), which within uncertainties is consistent with zero. Fig. 9. The age-depth plot for all samples in this study covering the period of ca 30 ka to 700 ka using various luminescence signals. To provide better insight into the ages from Koz1 to Koz10 they are depicted in the inset.

Fig

  Fig. 1.

  analyse_FadingMeasurement() implemented in the R ([START_REF] Core | R: A language and environment for statistical computing[END_REF]

  caused a heterogeneous environment in which the in situ gamma dose rate values are smaller than the laboratory measurements. Therefore, we considered the in situ values more representative of the gamma dose rates received by the samples and used those values to determine the final dose rate.

  three samples were measured by integrating the first 6 s and subtracting the last 56 s as the background. The obtained g-values (measured for six aliquots per sample) for samples Koz4, Koz8 and Koz11 resulted in 3.25 0.04 (3.52 0.04, norm. 2 days) %/decade, 3.31 0.04 (3.58 0.04, norm. 2 days) %/decade, 3.35 0.07 (3.63 0.07, norm. 2 days) %/decade (see Fig S5). The fading corrected Des for these three samples were 109 1 Gy, 171 4 Gy, 512 72 Gy, respectively. The typical dose-response curves and De distribution for sample Koz4 are illustrated in Fig. 4c and Fig. 5c, respectively.

Fig

  Fig. 4.

Fig. 6 .

 6 Fig. 6. Typical IR-RF curves were obtained for three samples. Each plot represents results for one multiple grain aliquot. The natural signal is indicated by red-filled circles and the regenerated signal by green crosses. The horizontal and vertical sliding was used to estimate the De for each subsample. Fig. 7. The kernel density estimate (KDE) plots for three samples as examples to show typical Des distributions. Fig. 8. VSL MAAD dose-response curve. The L/T s for each set of irradiated samples are shown in this figure.For the fitting, we applied a dose-response curve function based on the Lambert W function. The value of 11 12 Gy was obtained using this fitting for sample Koz4 (not irradiated), which within uncertainties is consistent with zero. Fig.9. The age-depth plot for all samples in this study covering the period of ca 30 ka to 700 ka using various luminescence signals. To provide better insight into the ages from Koz1 to Koz10 they are depicted in the inset.

Figure

  Figure

Figure

  FigureClick here to access/download;Figure;Figure2.pdf

Table 1 .

 1 Location of the sediment samples taken from the site.

	Sample name	Squares	X	Y	Z	Geological layers
	Koz1	H7	199.000 1498.900 -4.259	5a
	Koz2	H7	199.000 1498.580 -4.570	5b
	Koz3	H7	199.000 1498.333 -4.609	5c
	Koz4	H6	199.000 1497.632 -4.695	6/7
	Koz5	H8	199.000 1499.560 -4.549	6/7
	Koz6	H8	199.000 1499.590 -4.744	9ab
	Koz7	H7	199.000 1498.580 -4.787	9ab
	Koz8	H7	199.000 1498.730 -4.954	9c
	Koz9	H8	199.000 1499.560 -4.896	9c
	Koz10	H7	199.000 1498.610 -5.078	10a
	Koz11	F8	196.000 1499.119 -5.232	10b
	Koz12	F7	196.000 1498.325 -5.310	10b
	Koz13	H7	199.000 1498.770 -5.259	10c
	Koz14	F8	196.000 1499.122 -5.396	10c
	Koz16	F8	196.000 1499.107 -5.880	11b
	Koz18	G11	197.419 1503.000 -7.437	12
	Koz20	G11	197.473 1503.000 -7.551	13a
	Koz21	H11	198.000 1502.000 -8.000	13c

Table 2 .

 2 The concentration of radionuclides obtained through gamma-ray spectrometry. 35±0.02 1.21±0.07 1.17±0.02 4.95±0.05 Koz11 1.13±0.01 0.88±0.06 0.89±0.01 4.08±0.05 Koz12 1.20±0.01 0.92±0.06 0.91±0.01 4.30±0.05 Koz13 1.29±0.02 1.16±0.07 1.11±0.01 5.64±0.05 Koz14 1.34±0.01 1.25±0.05 1.17±0.01 5.75±0.05 Table3. The estimated dose rates for sediment samples in this study.

				U post-
	Sample	K	U pre-226	226 Ra	Th
	name	(%)	Ra (ppm)	(ppm)	(ppm)
	Koz1	1.58±0.02 1.93±0.08 1.97±0.02 6.83±0.07
	Koz2	1.48±0.02 1.17±0.07 1.24±0.02 4.75±0.05
	Koz3	1.54±0.02 1.19±0.07 0.99±0.02 4.95±0.06
	Koz4	1.47±0.02 1.13±0.07 1.33±0.02 5.26±0.06
	Koz5	1.44±0.02 1.11±0.07 1.29±0.02 4.72±0.06
	Koz6	2.15±0.02 2.08±0.08 1.64±0.02 7.92±0.07
	Koz7	2.06±0.02 1.92±0.07 1.55±0.02 7.63±0.06
	Koz8	1.41±0.02 1.14±0.05 1.32±0.01 4.84±0.05
	Koz9	1.37±0.02 1.24±0.07 1.13±0.02 4.78±0.05
	Koz10 1.

Table 4 .

 4 The De and corresponding ages within uncertainties estimated using various luminescence signals.

		OSL		pIR50-		IRSL		VSL		IR-RF
				IRSL290
	Sample	De	Age	De	Age	De	Age	De	Age	De	Age
	names	(Gy)	(ka)	(Gy)	(ka)	(Gy)	(ka)	(Gy)	(ka)	(Gy)	(ka)
	Koz1	57±2	32±2
	Koz2	48±2	30±2	95±3	39±2
	Koz3	56±3	33±2
	Koz4	53±1	33±2	90±3	37±2	109±1	45±3
	Koz5	58±2	35±2
	Koz6	117±9	49±4	176±5	55±3			101±35	42±15
	Koz7	101±6	43±3	162±6	51±3			96±21	41±9
	Koz8	92±4	53±3	158±15	61±6	171±4	66±4	166±119	95±68
	Koz9	82±4	51±3	138±3	57±3
	Koz10	90±5	59±4	168±4	71±4