The unexpectedly short Holocene Humid Period in Northern Arabia

The early to middle Holocene Humid Period led to a greening of today’s arid Saharo-Arabian desert belt. While this phase is well defined in North Africa and the Southern Arabian Peninsula, robust evidence from Northern Arabia is lacking. Here we fill this gap with unprecedented annually to sub-decadally resolved proxy data from Tayma, the only known varved lake sediments in Northern Arabia. Based on stable isotopes, micro-facies analyses and varve and radiocarbon dating, we distinguish five phases of lake development and show that the wet phase in Northern Arabia from 8800–7900 years BP is considerably shorter than the commonly defined Holocene Humid Period (~11,000–5500 years BP). Moreover, we find a two century-long peak humidity at times when a centennial-scale dry anomaly around 8200 years BP interrupted the Holocene Humid Period in adjacent regions. The short humid phase possibly favoured Neolithic migrations into Northern Arabia representing a strong human response to environmental changes. The Holocene Humid Period in Northern Arabia may only have lasted for less than 1000 years and was characterised by substantial regional climatic variability, according to a high-resolution multi-proxy varved lake record from Tayma, Saudi Arabia.

P ast millennial-scale pluvial periods driven by precessionforced intensification of summer monsoons and northward migration of associated rainfall are thought to have facilitated human dispersal out of Africa 1-3 by providing 'green corridors' through today's arid Saharo-Arabian desert belt [4][5][6] . Research on human-climate interaction on the Arabian Peninsula has intensified only recently, even though the region demonstrates high ecological sensitivity to climatic changes and represents the geographic nexus between Africa and Asia [1][2][3]7,8 . This recent wave of research in Arabia has already fundamentally transformed our perception of Arabian prehistory, including discoveries of Middle Palaeolithic (MIS 5 or even older) sites in Central Arabia 9 or traces of Homo sapiens in the Nefud desert at approximately 87 ka 8 , i.e. phases associated with conditions more humid than today 2,6 . One of the emerging topics is the role of the early to middle Holocene Humid Period (HHP) in Neolithic migrations and cultural progress 10,11 .
Climate models suggest that the African Summer Monsoon (ASM) was the dominant moisture source on the Arabian Peninsula during pluvials 1,12 . Yet, this remains a matter of debate for the Northern Arabian desert [13][14][15] , as stronger insolation intensified and extended both ASM 1,2,12,[16][17][18] and, possibly, Mediterranean winter rains 3,19 . The latter system is the main source of moisture in this region today. In addition, tropical plumes (TPs), i.e. synoptic disturbances conveying water vapour as continuous mid-upper tropospheric cloud bands from the Intertropical Convergence Zone (ITCZ) to >15°N, are known to affect Northern Arabia during winter and spring 14,20,21 . Higher frequency of such patterns during past pluvials was suggested to have contributed increased rainfall to the Saharo-Arabian desert 14,22,23 , even though their past role as a moisture source remains poorly understood.
The rich archaeological heritage of Arabia is currently being investigated by major research initiatives 7,11,24 . "Potentially thousands of water bodies" have been reconstructed for past pluvials 25 , but it is still unknown how these water bodies and human habitats exactly looked like and for how long they existed 14,15 . Only a few climate records are available from speleothems in the wider region, i.e. the Levant [26][27][28] and Southern Arabia 29,30 . The entire lack of high-resolution palaeoclimate data from Northern Arabia leads to an inconsistent picture about the timing and magnitude of the HHP for this culturally important corridor to the Middle East, where some lowerresolution lacustrine records have pointed to more humid conditions during the last interglacial (MIS 5) and the early to mid-Holocene 10,17,[31][32][33] .
The Tayma palaeolake record 4,34 is the only known highresolution archive of the HHP in Northern Arabia providing insights into the early to mid-Holocene hydroclimate variability in unprecedented detail. The 20 km 2 -sized inland sabkha of Tayma has a 660 km 2 hydrological catchment ( Fig. 1; Supplementary Figs. 1 and 2) and is located at 27°40′N within the arid desert's interior. It receives only scarce rains (on average 45 mm a −1 ) from Mediterranean winter storms, occasional cross-Saharan TPs or the Red Sea cyclones between autumn and spring 14 (Fig. 1). Previous investigations of shoreline deposits  and sediment cores from the sabkha basin indicate the existence of a >17 m deep, perennial groundwater-supported lake 15 and the spread of grassland 4 during the early Holocene. The catchmentlake ratio ( Fig. 1b; Supplementary Fig. 1b) and the short duration of this peak lake phase 4,15,34 exclude the influence of tectonics on lake-level changes, emphasising the significance of the lake as a palaeoclimate archive that is mainly controlled by rainfall, groundwater inflow and evaporation (see also Supplementary Note 1 for details on controlling processes of the sedimentary archive). Yet, a precise determination of the timing of the lake phases was still missing, preventing a detailed view of the evolution of the palaeolake and the palaeoclimatological implications.
In this paper, we decipher five phases of lake development at Tayma and show that the HHP in Northern Arabia from 8800 to 7900 years BP is shorter than commonly defined (~11,000-5500 years BP). Interestingly, we identify a two century-long peak humidity at Tayma overlapping with a centennial-scale dry anomaly around 8200 years BP that interrupted the HHP in adjacent regions. The results clearly demonstrate that regional patterns in palaeoclimate were more complex and spatially smallscale than previously assumed.

Results
Chronology of the Tayma palaeolake record. The Tayma palaeolake record partly contains annually laminated sediments that were counted under the microscope (see Methods, Fig. 2, Supplementary Fig. 9). The new high-resolution age-depth model integrates AMS radiocarbon ages of pollen concentrates, microscopic varve counting and the independent age of a cryptotephra 35 in a Bayesian model (see Methods, Supplementary Data 1, Supplementary Fig. 7). The floating varve chronology comprising 650 ± 40 couplets is anchored to the radiocarbon age scale and constrains the varved lake phase at Tayma to 8550-7900 ± 40 cal varve yr BP (±90 cal yr BP including the 14 C calibration error). A robust time marker is provided by the identification of the central Anatolian 'S1' tephra in the lower part of the record, dated to 8983 ± 83 cal yr BP in the Dead Sea record 34 . The lacustrine and wetland sediments in the Tayma basin were deposited from ca. 9250 to ca. 4200 cal yr BP (Supplementary Fig. 7).
Changes in precipitation and evaporation signals in the Tayma record. Compound-specific hydrogen isotope compositions of plant-wax n-alkanes (δD nC29, nC31 ), as well as porewater, rainwater and groundwater isotopes (δ 18 O water and δD water ) trace variations in moisture supply and rainfall amount (Fig. 2, Methods and Supplementary Fig. 8). Stable oxygen and carbon isotope compositions of single primary aragonite laminae from the varved core section (δ 18 O arag and δ 13 C arag ) and of bulk carbonates from the entire core (δ 18 O carb and δ 13 C carb ) indicate changing groundwater and surface-water inflow, lake-water evaporation and lake-internal productivity (Fig. 2). The most striking finding was terrestrial plant wax δD nC29, nC31 values of about −100‰ during the shallow lake or wetland phase and much lighter values down to −150‰ during the palaeolake phase. These data reflect higher rainfall between 8800 and 7950 cal yr BP due to increased precipitation and a probable amount effect (Methods and Supplementary Fig. 8).
Evolution of the Tayma palaeolake. The evolution of the lake can be separated into phases I-V, followed by phases VI (wetland) and VII (sabkha). A basal zone of brownish-grey carbonate mud with irregular, coarse carbonate lamination from 9250 to 8800 cal yr BP (lake phase I) ( Supplementary Fig. 6) represents a shallow lake initiated by increasing rainfall and recharge of the local Saq aquifer when clastic sediments were deposited in a deflated endorheic basin from a prevailing desert environment 4 . Carbonates precipitated with very high δ 18 O carb values of around +11‰ and low δ 13 C carb values of around -8‰.
At ca. 8800 cal yr BP (lake phase II), a sharp decrease of δ 18 O carb to +8‰, increasing δ 13 C carb (Fig. 2) and the in situ deposition of the brackish-water ostracod Cyprideis torosa (Fig. 2e) indicate reduced lake-water evaporation and the initial establishment of a shallow, but perennial and increasingly productive water body 36 as a response to wetter conditions.  [3][4][5] and stratified permanent lake, similar and coeval to, e.g. the Awafi palaeolake at Ras Al-Khaimah, United Arab Emirates 37 . From ca. 8550 to 8250 cal yr BP (lake phase III), variable but continuously decreasing plant wax δD nC29, nC31 values between −150 and −100‰ indicate a humid period with enhanced seasonality 38,39 . The alternating deposition of dark clastic, clay-and organic-rich laminae and white, primary aragonite laminae reflects pronounced wet and dry seasons. Suspension grading observed in some of the clastic varves indicates wadi activation and fluvial sediment input into a standing water body during the wet season. The aragonite laminae represent the arid season of enhanced evaporation, reduced or absent surface-water input, lake contraction and lakelevel fall, as well as increased concentration of soluble matter 40 . The highest concentrations of brine occur near the water surface, where aragonite crystals form and settle through the water column as pelagic rain 41,42 . The δ 18 O carb values generally decrease from +8‰ to +6‰ simultaneously with progressively increasing δ 13 C carb values from about −6 up to +2‰ toward enhanced lake productivity 43,44 . The positive excursion of δ 18 O carb to >+10‰ centred at ca. 8400 cal yr BP reflects a decadal-to centennial-scale drawback to even stronger dryseason evaporation. This intensified dry-season evaporation was compensated by enhanced humidity during the rainy season and groundwater inflow, indicated by the decreasing trend of δD nC29, nC31 at that time, data which mostly reflects the wet season of leaf growth 39 . This rainy-season moisture surplus was sufficient to sustain a high lake level and varve formation.
The phase between ca. 8250 and 7950 cal yr BP (lake phase IV) shows the highest production rate of organic matter in the lake and annual blooms of planktonic diatoms (mainly Cyclotella cf. choctawhatcheeana) (Supplementary Figs. 6 and 9). In combination with the greatest abundances of foraminifera, the lowest δD nC29, nC31 values down to −155‰ and weakest dryseason evaporation with lowest δ 18 O carb values of +4‰, they indicate the highest lake stand and most humid period at Tayma during the Holocene. This is supported by a distinct change in varve composition from evaporation-driven aragonite varves to pronounced productivity-driven diatom-aragonite varves and total organic carbon (TOC) contents of up to 5%. The highest ratio between δ 13 C carb and δ 18 O carb in this part of the core is related to phytoplankton bloom controlled by 12 C depletion due to photosynthetic uptake of CO 2 , and higher precipitation-evaporation balance 43 . From about 8200 cal yr BP the δD nC29, nC31 values again start to vary between −140‰ and −100‰, and the δ 18 O carb values increase from +4 to +8‰ ( Fig. 2; Supplementary Fig. 6).
At ca. 7950 cal yr BP, ceasing diatom and aragonite laminae and more abundant clastic quartz grains of aeolian origin, as well as the first appearance of gypsum (Supplementary Figs. 6 and 9), show a rapidly declining lake level (lake phase V). The increased mobility of sand grains indicates a regional decline of vegetation cover, the retreat of grasslands and the establishment of droughtresistant steppe vegetation 4 . The occurrence of gypsum is related to the enrichment of Ca 2+ ions in the contracting water body and sulphate dissolved in the reduced groundwater and surface water inflow 45 . This led to the disappearance of varves within a few decades, accompanied by a sharp reduction in TOC content and a decline of δ 13 C carb back to a level comparable to the early shallow-lake phase II, prior to 8550 cal yr BP. Progressively enriched δD nC29, nC31 values of up to −60‰ and δ 18 O carb values toward +12‰ reflect the decrease in surface-water and groundwater inflow and a strongly increasing evaporation, indicating a gradual end of the humid phase over 100-150 years until ca. 7800 cal yr BP.
Further increasing gypsum and Mg-carbonate precipitation and δ 18 O carb rising to +12‰ ( Fig. 2; Supplementary Fig. 6) point to a shrinking lake and concentration of brine under an increasingly arid climate between 7800 and ca. 6800 cal yr BP.  14 . Isohyetes of 200 mm and 500 mm are adapted from the wettest simulation of the COSMOS climate model at 6000 BCE in ref. 12 . b Reconstructed extent of the Tayma palaeolake during its peak phase, today's surface catchment and groundwater catchment 79 . The topography is based on GTOPO30 data 80 .
At that time, wetland conditions set in with TOC levels close to 0 and further increasing aeolian influx (phase VI). Around ca. 4200 cal yr BP, greyish mud is replaced by reddish-brown, oxidised, mainly aeolian clastics mixed with gypsum (phase VII) (Supplementary Figs. 6 and 7). This pattern reflects temporary desiccation and a further aridification pulse correlating with a dry event recorded at several sites in the Eastern Mediterranean/Middle East, e.g. the Northern Red Sea 46 .

Discussion
Our data support existing low-resolution Northern Arabian palaeoenvironmental records 10,11,31 , but they also show that the HHP in Northern Arabia was remarkably short, its peak lasting only ca. 650 years from 8550 to 7900 cal yr BP. A temporary intensification and northeastward shift of the ASM to provide increased rainfall in northwestern Arabia around Tayma during the short HHP is supported by the δ 18 O water and δD water data. The evaporation line calculated from modern water samples around Tayma and the palaeolake data best fits the isotope data from δD-depleted ASM precipitation at Khartoum, Sudan 47 ( Supplementary Fig. 8). The ASM fuelling the HHP in Northern Arabia has also been suggested by several proxy-based 17,32 and climate modelling 1,12,16 studies. We propose a moderate influence of the ASM during the HHP, as Tayma is located at the fringes of the ASM influence at that time 12 .
In addition, we observe an intriguing regional hydroclimatic diversity, since the aforementioned peak humidity in Northern Arabia from 8550 to 7900 cal yr BP coincides with a widespread, centennial-scale dry and cool anomaly centred around the 8.2 ka cold event at other low-latitude sites in the Northern Hemisphere such as the Eastern Mediterranean or Southern Arabia 48 (Fig. 3). A low-latitude dry period between ca. 8500 and 7800 cal yr BP was the most pronounced hydroclimatic drawback of the HHP (Fig. 3a). It is evidenced, e.g. in the desiccation or temporary lowstands of North African lakes 49 (Fig. 3f), slightly heavier δ 18 O values in Northern Red Sea planktonic foraminifera 13 (Fig. 3e) and diminished runoff of the Nile River 50 , leading to reoxygenation of the Eastern Mediterranean Sea and interruption of sapropel S1 51,52 (Fig. 3b). Drier conditions mostly resulted O carb/arag and δ 13 C carb/arag measured on bulk carbonates (solid lines) and on single aragonite layers (blue and green diamonds); c δD wax of n-alkanes nC 31 and nC 29 (δD nC29, nC31 ); d Tayma lake phases I-V and wetland phase VI; e varve sublayers expressed as % of varve thickness for the varve chronology 8550-7900 ± 40 cal varve yr BP, and microscope photographs of thin sections highlighting different micro-facies of lake phases II-V. from reduced summer-monsoon rainfall 29,30 . However, speleothem records from the Levantine region 26 (Jeita and Soreq caves; Fig. 3b, c) and marine records from the Eastern Mediterranean suggest that Mediterranean winter rains might have been reduced as well, as a result of temporary, meltwater-related deceleration of the North Atlantic thermohaline circulation 48,53 (Fig. 3). The inference of moderately higher than average rainfall during the HHP from the Soreq speleothem record, mostly during the phase of excess winter rains between 8000 and 7000 cal yr BP, is based on the recalculation of the δ 18 O data in ref. 51 eliminating the source effect, i.e. the bias through changes of isotopic composition in the Eastern Mediterranean surface waters 28 .
These smaller-scale, regional discrepancies may be explained by the additional contribution of synoptic-scale patterns, which are scarce today. They may have played a more dominant role in delivering moisture to Northern Arabia between 8250 and 8000 cal yr BP resulting in a humidity peak at Tayma. In particular, more frequent TPs may have led to a regional moisture surplus in Northern Arabia. In contrast to short and localised convective cells of the Active Red Sea Trough pattern triggering flash floods in the Southern Levant, TPs promote long-lasting moderate rains and thus more effective moisture over a slightly larger region 20 . TP formation was likely favoured by oceanatmosphere feedbacks during the 'cool poles-dry tropics' anomaly around 8.2 ka: lower sea-surface temperatures in the North Atlantic and the Mediterranean Sea promote deeper, southwards penetrating mid-latitude troughs and stronger subtropical anticyclones (i.e. drier air masses). This leads to an intensification of tropical moisture advection and the sub-tropical jet stream, inducing jet streaks that reach as far as northern tropical West Africa and convey moist air to Northern Arabia at mid-to-upper tropospheric levels 21 . The regionally and chronologically confined enhanced contribution of TPs may have compensated the weakened ASM around 8.2 ka and contributed to the strong seasonal pattern reflected by high-level, highamplitude δ 18 O carb and δD nC29, nC31 data. Mediterranean westerlies as a source of additional winter rainfall, however, seem unlikely, as climate models indicate an even lower intensity during the HHP compared to today 12 .
There is multiple evidence that the observed moisture surplus of 8550-7900 cal yr BP identified in the Tayma record, in combination with charged aquifers, had distinct short-term impacts on the local environment and probably also on human migration. Vegetation resources 4 and the abundance of prey animals 12 increased and stimulated Neolithic migrations into Northern Arabia as indicated by abundant Levant-type Pre-Pottery Neolithic A and B assemblages identified in the Northern branch of the Nefud desert 10,11 .

Methods
Tayma sediment cores. Drilling on today's sabkha of the Tayma palaeolake basin was performed in 2011 and 2013 using an Atlas Copco vibracoring device (Cobra mk1) fitted with closed steel auger heads and PVC liners with a diameter of 5 cm. Two series of ca. 6 m long sediment cores (Tay 220/221 and Tay 253/254/255/256) capturing the entire Holocene sequence and reaching down to Ordovician sandstone (Qasim Formation) were obtained in close vicinity ('mastercores' in Supplementary Fig. 1a). They each consist of two parallel, overlapping core sequences A and B with 1 m-long core sections. The cores were opened and photographically documented at the University of Cologne (Laboratory for Physical Geography) and GFZ Potsdam, Germany. The construction of composite profiles and correlation of the sediment cores is based on 24 macroscopic lithological marker layers (fixed marker horizons, FMH). Archive cores of the master site are stored and accessible in the core storage facility of the Institute of Geography, Heidelberg University, refrigerated at 4°C.
Tayma sediment cores were analysed for their sedimentology (XRF [X-ray fluorescence] core scanning, quantitative XRF on discrete samples, semiquantitative X-ray diffraction, micro-facies analyses on thin sections), geochemistry (elemental analyses, stable isotopes, lipid biomarkers), palynology (vegetation reconstruction through pollen analysis) and micropalaeontology (assemblages of foraminifera, ostracods and diatoms). Here, we used stable isotopes of oxygen and carbon measured on aragonite laminae of the annually laminated (varved) core section (δ 18 O arag , δ 13 C arag ) as well as bulk carbonates from the entire core (δ 18 O carb , δ 13 C carb ) in combination with micro-facies analyses of the varved sediments to trace the evolution of the early to mid-Holocene palaeolake at Tayma. Further proxy data have partially been published 4,36,54 or will be presented in forthcoming publications.
In Supplementary Fig. 6, the lithological profile of the ca. 6.5 m-long composite core, TOC content 54 , δ 18 O carb and δ 13 C carb (see Methodological details further down), and statistical clustering results of the XRF core-scanning record are shown. The elemental composition of the sediment core was determined by nondestructive XRF core scanning on the split-core sediment surface using an ITRAX elemental scanner at GFZ Potsdam. Measurements were obtained every 0.2 mm using a Cr X-ray source, operated at 30 kV, 30 mA and 10 s, to capture intensities of the elements Si, S, Cl, K, Ca, Ti, Fe, Sr and Zr. A centred log-ratio (clr = ln [element intensity/geometric mean of all nine elements]) transform was performed for all elements of each measurement to eliminate the influences of physical properties, sample geometry and matrix effects 55,56 and to enable robust statistical analyses 57 . Statistical clustering (Ward's method) of XRF core-scanning results indicates four main sediment groups (Supplementary Fig. 6): Cluster 1 (light grey) is dominated by Si, Ti and Fe and describes the siliciclastic sediments and occurs predominantly in the upper part of the Tayma profile (VII-sabkha phase). Cluster 2 (green) does not show a clear preference, but is rather a mixture of all considered elements, describing clastic, carbonate and evaporitic 'background' sediments. Cluster 3 (blue) is dominated by Sr and Ca and describes aragonite, which occurs exclusively in the varved sediments of the Tayma core representing the deep-lake phases III and IV (Fig. 2e). Cluster 4 (orange) is dominated by the elements S and Ca and mainly describes gypsum that was deposited during the terminal lake phase (V) and thereafter, when wetlands occupied the Tayma basin (phase VI).
TOC indicates the production and preservation of organic matter in the lake and was measured on in situ decalcified samples using an elemental analyser (NC2500 Carlo Erba) at GFZ Potsdam, Germany. Ca. 3 mg of freeze-dried and powdered sample material was weighed into Ag-capsules, treated with 20% HCl, heated and dried for 3 h at 75°C. The calibration of the data was carried out applying element standards (Acetanilide, Urea). It was verified using a soil reference sample (Boden3, HEKATECH). The reproducibility for replicate samples was~0.2% 54 . The sediments deposited in the Tayma basin are mainly composed of clay, silt and sand, evaporites (sulphates), authigenic carbonates and in parts high amounts of diatoms (at the mastercore site in the centre of the basin), gastropod shells, ostracods and foraminifera (along the margins of the basin representing the palaeo-shoreline) (Supplementary Fig. 6). Clay-and silt-sized detritus is dominant in the core and was deposited as dark grey, mm-to cm-thick, occasionally graded layers. Coarser silt-to sand-sized minerals (mainly quartz) are scattered in the sediments or are concentrated in the uppermost part of the Tayma profile. Evaporites were mainly identified in the form of whitish-beige, finer-grained laminae of gypsum and other sulphates. Carbonates are present in the form of white, sub-mm thick primary aragonite laminae, and, in the upper part, Mg-carbonate layers. Biogenic carbonates (ostracods, foraminifera, and barnacle and gastropod shell fragments) are dominant along the shorelines of the palaeolake. Their concentration is much lower in the centre of the basin where the mastercore was taken.
Varve micro-facies analysis. We used changes in varve micro-facies, i.e. the composition of seasonal sublayers of the annual laminations, to infer changing seasonality and the interannual variability of lake-internal evaporation and productivity. The thickness and composition of varve sublayers were analysed under the microscope along with varve counting on petrographic thin sections. A total of eleven different sublayer types were grouped into five main sediment components (carbonate, organic, clastic, diatoms and gypsum). Data are given as relative contribution (in %) to the varve thickness (Fig. 2e). Raw data of micro-facies sublayer thicknesses are presented in Supplementary Fig. 9.
Age model construction. Due to the absence of datable terrestrial macroscopic plant remains in Tayma cores and reported hard-water effects altering radiocarbon ages from gastropods, ostracods and Ruppia seeds for up to 1500 years 4,15 , preliminary age models 4,54 were based on AMS radiocarbon dating of pollen grains, as these are unsusceptible to incorporating old carbon 58,59 . Pollen extraction from a total of 33 samples of 1-13 cm long sediment sections followed a combination of physical and chemical separation protocols [58][59][60][61] . Sample preparation included sieving (at 6, 20, 40 and 70 µm), treatment with heated HCl, KOH and H 2 SO 4 , and heavy-liquid density separation using CsCl and sodium polytungstate.
Varve counting was performed on 14 large-scale (10 cm × 1.5 cm) petrographic thin sections using a Leica DMLP petrographic microscope under semi-/fully polarised light and with 50× magnification. Thin sections were prepared following standard procedures for soft sediments 62 including freeze-drying and impregnation with epoxy resin (Araldite 2020). Sawing and polishing were performed manually under dry conditions to avoid salt crystallisation. Multiple counting and the definition of correlation marker layers ensured a negligible subjective counting error. Counting uncertainty due to poor sublayer quality is ±40 varves (6.2%).
The age-depth model was constructed with Bacon v2.2 using flexible Bayesian modelling 63 including implemented outlier analysis and the IntCal13 atmospheric calibration curve 64 . All 38 radiocarbon ages of pollen concentrate, other plant remains (Ruppia seeds, non-pollen palynomorphs, charred plant particles), two mollusc samples, as well as a tephrochronological anchor identified close to the base of the Tayma sediment record (the central Anatolian 'S1' tephra dated at 8983 ± 83 cal yr BP in the Dead Sea) 35 , were considered for age modelling (Supplementary Data 1). The floating varve chronology of 650 ± 40 varve years served to refine the Bayesian model within the varved section. The start of varve formation is defined by 14 C dating to 8549 cal yr BP (8470-8605 cal yr BP for the 95.4% probability range). Based on this fixed point (F1), the varve age of a turbidite layer at 8138 ± 40 varve years BP and the end of varve formation at 7900 ± 40 varve years BP were used as further fix points (F2 and F3) in the adjusted Bayesian model (Supplementary Fig. 7).
Outlier analysis reliably discarded samples (n = 6) containing ≤50% pollen or hard-water-affected material (gastropod shells, Ruppia seeds), and 13 samples with ≥50% pollen unsuitable for the Bayesian age-depth model. All remaining 18 14 C ages of pollen concentrate included in the final model contained high pollen concentrations of at least 50% (Supplementary Data 1, Supplementary Fig. 7).
Reconstruction of hydroclimatic conditions. The stable isotope composition of lake water (δ 18 O water and δD water ) in closed lakes is mainly controlled by precipitation and evaporation and reflects hydrological changes and moisture sources 43 . The δ 18 O carb (δ 13 C carb ) values of bulk carbonates and δD nC29, nC31 from fossil leaf waxes in lake sediments are proxies for hydroclimatic conditions and were used to reconstruct the precipitation, lake-water evaporation and temperature during the early to mid-Holocene. To assess the hydrological balance of the Tayma palaeolake (8800-7950 cal yr BP), the wetland (7800-6800 cal yr BP), and the potential moisture sources during the HHP, we compared calculated δD p (precipitation) and δ 18 O water (lake water) values with the isotopic characterisation of the main regional atmospheric systems, recent precipitation, as well as surface and groundwater isotope compositions ( Supplementary Fig. 8).
Stable oxygen and carbon isotopes. Stable oxygen and carbon isotope measurements (δ 13 C carb and δ 18 O carb ) were performed on the carbonate fraction of a total of 262 freeze-dried and ground samples taken in cm slices from core Tay 220. Bulk samples of~0.4 mg were loaded into 10 ml Labco Exetainer vials automatically flushed with He and reacted in phosphoric acid (100%) at 75°C for 60 min 65 . The stable isotope compositions were determined at GFZ Potsdam using a Finnigan GasBenchII with carbonate option coupled to a DELTAplusXL IRMS (isotope ratio mass spectrometer) (Thermo Fisher Scientific). For finer resolution and to exclude biogenic carbonates, we sampled pure primary aragonite of laminae from dried and impregnated sediment blocks by drilling. The aragonite samples of about 0.06 mg were measured for δ 18 O arag and δ 13 C arag at GFZ Potsdam with an automated carbonate device (KIEL IV) coupled to a Finnigan MAT253 IRMS (Thermo Fisher Scientific) on cryogenically purified CO 2 released by dissolution with 103% H 3 PO 4 at 72°C for 10 min. Oxygen and carbon isotope compositions are given relative to the VPDB (Vienna Peedee Belemnite) standard in conventional delta notation δ (‰). Calibration was performed using international reference standards (NBS18 and NBS19). For both methods, standard deviations (1σ) for reference and replicate analyses are better than 0.08‰ for both δ 18 O and δ 13 C.
In closed lakes, δ 18 O carb values mainly reflect hydrological changes and are used as a proxy for precipitation, groundwater influx and lake evaporation because: (i) seasonality and temperature have little effect on oxygen isotope fractionation of precipitation in low-latitude regions 43,66 ; (ii) the lake-water oxygen isotopic composition in an endorheic basin is governed by evaporation under arid climate conditions resulting in increased δ 18 O water ; (iii) equilibrium oxygen isotope fractionation is assumed for inorganic carbonates; and (iv) primary inorganic carbonates precipitate during the spring-summer season induced by evaporation and/or phytoplankton bloom in the epilimnion. The latter is consistent with increasing δ 13 C carb values, indicating 12 C depletion of the total dissolved inorganic carbon due to atmospheric release and/or preferential use of aquatic plants.
The calculation of δ 18 O VSMOW palaeolake water from δ 18 O VPDB values using the re-expressed relationship of ref. 67 in the simplified Eq. (1) according to ref. 43 T ¼ 13:8 À 4:58ðδc À δwÞ ð 1Þ under equilibrium water (w)-calcite (c) precipitation at a temperature (T in°C) of 21°C (as the average temperature in spring) and an offset of +0.6‰ for aragonite and magnesium bearing calcite reveals comparable δ 18 O values between precipitated carbonates and hot water. The modelled mean δ 18 O water for the palaeolake water is high with +8.4‰ and, thus, much lighter due to freshwater inflow of surface and groundwater than for the wetland with a calculated mean δ 18 O water of +13.1‰ due to lower precipitation and high evaporation.
Stable hydrogen isotopes of leaf-wax n-alkanes (δD wax ). Stable hydrogen isotopes of leaf-wax n-alkanes (δD wax ) were measured on 64 samples from core Tay 255. Samples were taken in 1 cm-slices, freeze-dried and ground for lipid biomarkers extraction at the Max Planck Institute for Biogeochemistry in Jena. 5-15 g of the sample was extracted using a 40 ml dichloromethane:methanol (9:1) mixture at 100°C and 120 bar for 15 min in two consecutive cycles using a BÜCHI Spee-dExtractor. The total lipid extract was separated into aliphatic, aromatic and alcohol/fatty acid fractions using solid-phase extraction on silica gel according to the method presented in ref. 68 . The aliphatic hydrocarbon fraction was desulfurized using HCl-activated copper (15% HCl). Identification and quantification of n-alkanes were accomplished using a GC-MS (Agilent Technologies, 7890A GC-System; 220 Ion trap MS) by comparing peak areas and retention times with an external n-alkane standard mixture (nC16 to nC36). Compound-specific hydrogen isotope ratios (expressed as δD) of the n-alkanes were measured on a DELTA V plus IRMS (Thermo Fisher Scientific) coupled to an Agilent 7890 GC (Agilent Technologies) at GFZ Potsdam. Every sample was measured in triplicates. The mean standard deviation of all measured samples was 3‰. The δD values were normalised to the Vienna Standard Mean Ocean Water (VSMOW).
The changes in δD wax of the lake records are interpreted as indicators for the variability in precipitation, humidity and vegetation type 66,68,69 . Hydrogen isotopes δD nC29, nC31 of terrestrial leaf wax nC 29 and nC 31 n-alkanes were used to calculate δD p between precipitation (p) according to refs. 38,39,69 . The negative isotopic fractionation from δD p to δD wax due to the incorporation of hydrogen in leaf waxes has been calculated using Eq. (2) with ε = −130 for nC 29 and nC 31 n-alkanes representing the mixture of C3/C4 plant waxes of grasses, shrubs and trees that can be expected for the region 4 ( Supplementary Fig. 8). Following ref. 5 , we inferred the precipitation rate from δD p values. The relationship of precipitation and rainfall amount for the Sahara region described a non-linear dependence with a steep slope in δD p values below 100 mm yr −1 and a strong influence of the amount effect.
δ 18 O water and δD water isotopes of water samples. Filtered water samples of groundwater from the historical Bir Haddaj well of the Tayma oasis, from the Tay 255 borehole in the palaeolake, and evaporated rainwaters from small water pools south of the palaeolake taken shortly after a rain event in December 2015, were triple-measured at the GFZ for δ 18 O water and δD water relative to VSMOW using Cavity Ring-Down spectrometers (PICARRO L2120-i and L2130-i). Analytical precision of VSMOW and SLAP calibrated analyses was <1‰ for both δ 18 O water and δD water . The isotopic fractionation of lake-water evaporation was calculated for the remaining lake water sampled in December 2015 (δ 18 O rw ) using initial groundwater-supported lake water with δ 18 O iw of −3‰ with simple Rayleigh distillation after Eq. (3) where α = fractionation factor between water and vapour at 21°C 70 and ƒ = fraction of remaining lake water.
Reconstruction of palaeo-moisture source and lake-water evaporation. The few meteoric water samples from Tabuk (IAEA) plot with δ 18 O p~− 1‰ closely to the global meteoric water line (GMWL), except for one lighter sample tending more to the Eastern Mediterranean meteoric water line (EMMWL). Recent (12/ 2015) evaporated rainwater samples collected in water pools in a wadi SW of the Tayma palaeolake show slightly enriched δ 18 O water and δD water values of around −0.5‰. Using δD wax to estimate past precipitation rates as explained above reveals the lowest values of about −2‰ δD p for the wetland phase and values as low as −28‰ for the palaeolake (Supplementary Fig. 8), indicating higher rainfall amounts between 8800 and 7950 cal yr BP due to increased precipitation and a probable amount effect. The stable isotopes of the historical Bir Haddaj well in Tayma, one of the largest and most prominent historical wells on the Arabian Peninsula 71,72 that taps the uppermost level of the Saq aquifer, reflect subsurface groundwater with −3.5‰ δ 18 O water and −24.6‰ δD water , similar to the middle of the Saq aquifer further south 73,74 . The water from the palaeolake sampled in 1.5 m depth of the well Tay 255 in 2015 with +7.4‰ δ 18 O water and +16.2‰ δD water probably reflects porewater isotope composition. The portion of surface-water evaporation calculated using Rayleigh fractionation between −3‰ δ 18 O water for groundwater and +7.4‰ δ 18 O water for the Tay 255 well reaches about 70% for the deep lake phase and >80% for the wetland phase. The isotopic difference between the deep palaeolake and the wetland water is mainly influenced by decreases in precipitation and increasing evaporation.
Although the three main atmospheric systems affecting the northwestern Arabian Peninsula (Indian Monsoon, Mediterranean Westerlies and African Monsoon) show isotopic fingerprints with more or less variation and deuterium excess, it is unreasonable to decipher the moisture sources during the time of palaeolake formation due to the determination of precipitation δD p using δD wax and δ 18 O carb being indirect, as well as associated fractionation effects. In general, the moisture source-related isotope fingerprints over the Arabian Peninsula are masked by strong evaporation, continental and altitude effects, sub-cloud evaporation, moisture recycling and the amount effect 75 .

Code availability
This study does not use custom code or mathematical algorithm that is deemed central to the conclusions.