Chronological control and centennial-scale climatic subdivisions of the Last Glacial Termination in the western Mediterranean region

Chronological control and centennial-scale climatic subdivisions of the Last Glacial Termination in the western Mediterranean region Jon Camuera a , * , Gonzalo Jim enez-Moreno b , María J. Ramos-Rom an a , Antonio García-Alix b , c , Francisco J. Jim enez-Espejo c , Jaime L. Toney d , R. Scott Anderson e a Department of Geosciences and Geography, Faculty of Science, University of Helsinki, Finland b Departamento de Estratigrafía y Paleontología, Universidad de Granada, Spain c Instituto Andaluz de Ciencias de la Tierra (IACT), Consejo Superior de Investigaciones Científicas-Universidad de Granada (CSIC-UGR), Spain d School of Geographical and Earth Sciences, University of Glasgow, UK e School of Earth and Sustainability, Northern Arizona University, USA article info Article history: Received 1 October 2020 Received in revised form 23 December 2020 Accepted 13 January 2021 Available online 3 February 2021 Handling Editor: Donatella Magri Keywords: Paleoclimate Iberian peninsula Pollen analysis Last glacial Heinrich stadial 1 Solar activity abstract The Last Glacial Termination is marked by changing environmental conditions affected by abrupt and rapid climate oscillations, such as Heinrich Stadial 1 (HS1), which is characterized by extremely low sea surface temperatures (SST) and significant changes in northern hemisphere terrestrial landscape (e.g., vegetation) and human dispersion. Previous studies show that overall cold/dry conditions occurred during HS1, but the lack of high-resolution records precludes whether climate was stable or instead characterized by instability. A high-resolution paleoclimatic record from the Padul wetland (southern Iberian Peninsula), supported by a high-resolution chronology and contrasted with other records from southern Europe and the Mediterranean region, shows 1) that the age boundaries of HS1 in this area occurred at ~18.0 kyr BP (median age ¼17,970 cal yr BP; mean age ¼18,030 ±330 cal yr BP) and ~15.2 kyr BP (median age ¼15,210 cal yr BP; mean age ¼15,200 ±420 cal yr BP) and 2) that climate during HS1 was non-stationary and centennial-scale variability in moisture is superimposed on this overall cold climatic period. In this study, we improve the pollen sampling resolution with respect to previous studies on the same Padul-15-05 sedimentary core and suggest a novel subdivision of HS1 in 7 sub-phases, including: i) 3 sub-phases (a.1-a.3) during an arid early phase (HS1a; ~18.4e17.2 kyr BP), ii) a relatively humid middle phase (HS1b; ~17.2e16.9 kyr BP), and iii) 3 sub-phases (c.1-c.3) during an arid late phase (HS1c; ~16.9e15.7 kyr BP). This climatic subdivision is regionally supported by SST oscillations from the Mediterranean Sea, suggesting a strong land-sea coupling. A cyclostratigraphic analysis of pollen data between 20 and 11 kyr BP indicates that the climate variability and the proposed subdivisions characterized by ~2000 and ~800-yr periodicities could be related to solar forcing controlling climate in this area. ©2021 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). 1. Introduction The Last Glacial Termination in southern Europe and Mediterranean areas present one the most inhospitable environmental conditions of the last 130 kyr, reaching one of the lowest SST record for this period (Martrat et al., 2004,2007). During the last deglaciation, HS1 shows abrupt and complex climate signals under the roughest conditions, as observed by the presence of ice-rafted debris during Heinrich events (Hodell et al., 2017), the slowdown of Atlantic and Mediterranean thermohaline circulations (McManus et al., 2004;Sierro et al., 2020) and major genetic bottlenecks in humans (Fern andez-L opez de Pablo et al., 2019). In this respect, deciphering rapid (e.g., millennial-scale) climate changes and environmental impacts due to Dansgaard/Oeschger and Heinrich-like climatic oscillations during the last glacial period and deglaciation have been the aim of ice, marine and terrestrial paleoclimate investigations (Cacho et al., 2006;H€ obig et al., 2012; Panagiotopoulos et al., 2014;S anchez Go~ ni et al., 2008). The different Heinrich nomenclatures, that is Heinrich Events (HEs), Heinrich Layers (HLs) and Heinrich Stadials (HSs), has caused *Corresponding author. E-mail addresses: jon.camuera@helsinki.fi,[email protected] (J. Camuera). Contents lists available at ScienceDirect Quaternary Science Reviews journal homepage: www.elsevier.com/locate/quascirev https://doi.org/10.1016/j.quascirev.2021.106814 0277-3791/©2021 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Quaternary Science Reviews 255 (2021) 106814 much confusion in the scientific literature. HEs refer to massive discharge of ice-rafted debris (IRD) from Laurentide, Fennoscandian and Greenland ice sheets into the North Atlantic, resulting in the deposition of detrital and carbonate-rich sediment layers called HLs (Heinrich, 1988;Hemming, 2004;Hodell et al., 2017). HSs were described as cold temperature intervals revealed in the North Atlantic records during which HEs occurred (Barker et al., 2009). These periods and events represent the most extreme glacial conditions, which resulted from the culmination of the decreasing temperature trends of Bond cycles (Bond et al., 1993). Several marine paleoclimatic records evidenced the effect of especially cold and arid conditions recorded during HSs, also with strong direct influence in terrestrial records (Fletcher and S anchez Go~ ni, 2008;Fletcher et al., 2010a;Hodell et al., 2017;Martrat et al., 2014;Moreno et al., 2010). HS1 is the most recent and one of the coldest Heinrich Stadials of the last glacial cycle and has been described in several marine sedimentary records nearby the study area (Naughton et al., 2016;Salgueiro et al., 2014). However, different records show chronological discrepancies in the onset and the end of this period (Fletcher et al., 2010b;Moreno et al., 2010;S anchez Go~ ni and Harrison, 2010). Causes in age differences in the environmental responses obtained from different records in the same region to HS1 can be a consequence of: a) errors/uncertainties from dating techniques, b) poor age control produced by low sample resolution, c) very low sedimentary rates, precluding an accurate chronological control, d) the reservoir effect in marine sediments deposited during HS1 in the North Atlantic [~500e1300 14 C years according to Stern and Lisiecki (2013)] and Mediterranean Sea [~800 years according to Siani et al. (2001)] and/ or, e) the reworked materials in both marine and continental areas. High-resolution paleoclimatic records show that climate during HS1 was characterized by short-scale internal variability (Dupont et al., 2010;Escobar et al., 2012;Stager et al., 2011;Stríkis et al., 2015;Zhang et al., 2014). However, few studies have focused on short-term internal climate variability of HSs in southern Europe and the western Mediterranean region, and in particular within HS1 (Fletcher and S anchez Go~ ni, 2008). In this regard, a division of HS1 into two and three phases has previously been observed in very few marine records. For example, two phases were described in Iberian margin marine sedimentary records (Naughton et al., 2009;Salgueiro et al., 2014;S anchez Go~ ni et al., 2018) and Nile River Basin (Casta~ neda et al., 2016), characterized by similar wet conditions during the first phase and arid climate during a second phase. Other paleoclimate studies recorded a three-phase division for HS1, including studies from the Alboran Sea (Bazzicalupo et al., 2018;Fletcher and S anchez Go~ ni, 2008), northwestern Mediterranean (Sierro et al., 2005), Iberian margin (Naughton et al., 2016; Sierro et al., 2020;Voelker et al., 2009) and off NW Africa (Bouimetarhan et al., 2012)(Table 1). Nevertheless, the studies showing a three-phase division of HS1 disagree in the paleoenvironmental characterization of each phase and a complete knowledge of the variability within HS1 has yet to be achieved (Hodell et al., 2017). This study aims to improve the resolution of the pollen analysis for the HS1 time-period from a previous palynological study covering the last two glacial-interglacial cycles in the Padul-15-05 terrestrial sedimentary record (southern Iberian Peninsula) (Camuera et al., 2019). As a result, we present high-resolution pollen and sedimentological data between 20 and 11 kyr BP, registering regional vegetation and local lake paleoenvironmental responses to climate changes during the last glacial-Holocene transition, including HS1, Bølling-Allerød (BA) and Younger Dryas (YD). In order to have a good chronological control of HS1 in southern Europe and the Mediterranean region, we have compiled, revised, recalibrated and modelled the HS1 age-boundaries from non-tuned and independently dated high-resolution marine and terrestrial paleoclimatic records from this region. 2. Regional and local settings The Padul wetland (724 m a.s.l.) is located in the western margin of the Sierra Nevada range, 20 km south of Granada city (Andalusia, Spain) and covers an area of 4 km 2 in the Padul-Nigüelas basin (Fig. 1). The NW-SE elongated Padul-Nigüelas endorheic basin developed as a consequence of extensional activity of the main normal fault that delimits the NE edge of the basin (Santanach et al., 1980). The catchment of the basin comprises mainly Triassic limestone/dolomite and Cambrian-Triassic schists, whereas the basin fill sediments are principally Miocene gypsum and detritics, and Quaternary peat and alluvial deposits (Camuera et al., 2018). It bears an estimated sedimentary sequence of about 100 m in the depocenter of the basin (Ortiz et al., 2004). The precipitation in the region is highly controlled by the humidity carried from the westerlies and the North Atlantic Oscillation (Jim enez-Moreno and Anderson, 2012;Lionello, 2012). At present, the Padul area is characterized by a semiarid Mediterranean climate with high temperature and low precipitation during summertime (summer drought), presenting a mean annual temperature of 14.4  C and mean annual precipitation of 445 mm (AEMET, 2016). The present-day vegetation in the Padul wetland is mainly characterized by wetland communities, such as Phragmites australis and Typha domingensis, whereas the surrounding areas are dominated by mesomediterranean vegetation with Quercus rotundifolia, Q. faginea,Q. coccifera,Pistacia terebinthus,Populus alba,Ulmus minor,Fraxinus angustifolia and Celtis australis, among others. 3. Materials and methods 3.1. Padul chronology The 42.64 m-long Padul-15-05 sediment core was drilled in the Padul wetland lakeshore (37  00 0 39"N, 3  36 0 14"W) in July 2015 (Fig. 1). The chronological control of the entire Padul-15-05 core was based on 43 Accelerator Mass Spectrometry (AMS) radiocarbon dates, 4 Amino Acid Racemization (AAR) dates from gastropods (hydrobiid Milesiana schuelei) and two different sediment accumulation rates (SAR) for both peat and carbonate/marl lithologies for the bottom part of the core, showing the record of the last ~200 kyr (Camuera et al., 2018). In addition, six new AMS radiocarbon samples (from ~15.6 to ~19.4 kyr BP) have been analyzed in this study to better delimit the age range of HS1 (Table S1). In this study a new Bayesian age-depth model has been built for the last 30 kyr BP. This new age model, developed using the R-based BACON software (v.2.3.9.1 - July 2019) (Blaauw and Christen, 2011), is based on 40 radiocarbon dates (including specific compound radiocarbon dating), providing a more detailed and accurate chronological control and taking into account age uncertainties (Fig. 2 and Tables S1 and S2). The age model suggests a good fitting (stable output of the log-posterior time-series) and a robust Markov Chain Monte Carlo mixing (Gelman and Rubin Reduction Factor (1.016) under the safety threshold of 1.05 (Brooks and Gelman, 1998)). The parameters of the model are included in Fig. 2. The high-resolution age-depth model shows high sediment accumulation rates (SAR, 0.155 mm/yr) during this time period, and therefore, allows for the development of regional paleoclimate reconstructions from high-resolution multiproxy analyses. Note that all the ages from the Padul record provided throughout the manuscript are calibrated kiloyears before present expressed as kyr BP to simplify. J. Camuera, G. Jim enez-Moreno, M.J. Ramos-Rom an et al. Quaternary Science Reviews 255 (2021) 106814 2 3.2. Chronological analysis of HS1 in southern Europe and the Mediterranean region In order to compare and corroborate the good time delimitation of HS1 in Padul, the age-range (onset and end) for HS1 in southern Europe and the Mediterranean region was statistically calculated using a compilation of available marine and terrestrial (lakes, speleothems) sites recording HS1 in this area. The compilation of these records was done according to the following criteria: 1) Records are located in southern Europe and the Mediterranean region, specifically, between latitude 31  N and 48  N, and between longitude 10  W and 34  E. 2) Records are of high-resolution data only, where HS1 is represented with at least 10 samples (ca. <250-yr mean resolution). 3) In order to avoid circular reasoning, records should not be tuned to ice-core chronologies or to other nearby records, and ages for HS1 should be based on independent absolute dates providing an objective age-depth model. 4) Records should have an exact numerical age-range for HS1 suggested or mentioned in the corresponding original studies. High-resolution marine and terrestrial paleoclimatic records from 17 sites (including Padul-15-05) with non-tuned and independent chronology have been compiled following the above criteria (Figs. S1 and S2). We selected the most sensitive paleoclimatic proxies responding to HS1 from each study site and those chosen by authors as the most representative indicators of Table 1 Marine records presenting a division of HS1 in three phases (Bazzicalupo et al., 2018;Bouimetarhan et al., 2012;Fletcher and S anchez Go~ ni, 2008;Naughton et al., 2016;Sierro et al., 2005,Sierro et al., 2020; Voelker et al., 2009). Marine and continental reconstructions of the early, middle and late HS1 have been schematized. Note the diversity in the interpretations for the three phases identified within HS1. Study Location Record Environment Proxy Early HS1 (HS1a) Middle HS1 (HS1b) Late HS1 (HS1c) Fletcher and S anchez Go~ ni (2008) Alboran Sea MD952043 Marine N. pachyderma-s (cold water indicator) Cool Cold Cool Continental Vegetation Cold/humid Arid Cool/less arid Bazzicalupo et al. (2018) Alboran Sea ODP-976 Marine Calcareous plankton Cold Fresher water Cooler Continental Vegetation Increase aridity Maximum aridity Cold/arid Sierro et al. (2005) Northwestern Mediterranean MD992343 Marine d 18 OG. bulloides (SST, iceberg meltwater) Cold Cool Cold Naughton et al. (2016) Iberian margin MD032697 Marine Alkenone-based SST and % N. pachyderma-s Extreme cooling Warmer (still cool) Cooling (warmer than HS1a) Continental Vegetation Extreme cold/wet Warmer (still cool)/ increase aridity Warming/wet Voelker et al. (2009),Sierro et al. (2020) Iberian margin MD992339 Marine SST and d 18 O seawater Cold Warmer Cold Bouimetarhan et al. (2012) Off NW Africa GeoB95085 Continental Vegetation Dryness Wetter Extreme dry Fig. 1. Geographical location of the Padul-15-05 record in the western margin of Sierra Nevada range and south of Granada city (southern Iberian Peninsula) (modified from Camuera et al., 2018). J. Camuera, G. Jim enez-Moreno, M.J. Ramos-Rom an et al. Quaternary Science Reviews 255 (2021) 106814 3 environmental changes. Even if the compiled paleoclimatic proxies represent different environmental signals (aridity, freshwater pulses, presence of cold water, etc.), the lack of enough records displaying the same signal as in Padul (aridity, temperature) under the criteria mentioned above does not allow us to suggest an age-range for HS1 based on a unique climate signal. The revision of HS1 age boundaries of the compiled records was done by 1) choosing the middle-points of HS1 climate transitions (i.e., Last Glacial Maximum-HS1 transition for the onset of HS1 and HS1-Bølling Allerød transition for the end of HS1), according to previous studies (e.g., Fletcher et al., 2010a,b;Rasmussen et al., 2014), and 2) taking the points presenting significant environmental changes closer to the HS1 boundaries provided by authors in each study (Fig. S2 and Table S3). Uncalibrated ages ( 14 C age ±error) of the revised HS1 age boundaries from each record were calculated by the means of the original 14 C calibration curves used in each study. Afterwards, the uncalibrated ages were independently recalibrated using recent calibration curves (IntCal13, Marine13) (Fig. 3 and Table S3). With respect to the U/Th ages from the CAN speleothem (Pindal cave), we have taken the ages boundaries from the U/Th-based age model ( 230 Th corrected age) of the original study (Moreno et al., 2010). The modelled age-ranges of HS1 for each record were calculated running a Kernel Density Estimation_plot (KDE_plot) function with a Bayesian approach using the Oxcal software (Ramsey, 2017). This methodology was also applied combining all these records (KDE Total) to obtain a Fig. 2. Bayesian age-depth model of the Padul-15-05 record for the last 30 kyr BP. The red square shows the discussed period in this study (20e11 kyr BP). See Tables S1 and S2 for specific information about the radiocarbon samples. Acronyms: acc. shape, accumulation shape; acc. mean, accumulation mean, mem. strength, memory strength; mem. mean, memory mean. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) J. Camuera, G. Jim enez-Moreno, M.J. Ramos-Rom an et al. Quaternary Science Reviews 255 (2021) 106814 4 synthetic regional model for the age-range of HS1 (Fig. 3). 3.3. Palynological analysis In this study, the resolution of the pollen data has been increased with 24 new pollen samples with respect to the previous study from the Padul-15-05 record by Camuera et al. (2019), resulting in a total of 91 pollen samples between 4.82 and 3.47 m depth (20e11 kyr BP). Therefore, the new sample resolution for HS1 (18.4e15.7 kyr BP) is ~61 years, and ~127 years for the BøllingAllerød (BA), Younger Dryas (YD) and the beginning of the Holocene (15.7e11 kyr BP). The pollen extraction was done following a modified methodology of Faegri and Iversen (1989). After the final extraction of the pollen residue, a minimum of 300 terrestrial pollen grains per sample were identified using a transmitted light microscope. Percentages of all pollen taxa were calculated based on the terrestrial pollen sum excluding aquatic plants (Cyperaceae, Typha,Myriophyllum,Utricularia and Potamogeton). The detailed pollen diagram has been represented in Figure S3 using the Tilia software (Grimm, 1987). The zonation of the pollen data was done on the main pollen taxa (i.e., Quercus total, Olea,Pistacia, Cupressaceae, Artemisia and Amaranthaceae) using the constrained cluster analysis (CONISS) in order to identify pollen zones and the climatic subdivision of the HS1, BA and YD periods (Fig. S3). Pinus total was excluded from the cluster analysis because it is overrepresented during some periods, as in other Iberian studies (García-Ant on et al., 2011;Morales-Molino et al., 2011). The Mediterranean forest, xerophytes, Pollen Climate Index (PCI) and Precipitation Index (I p ) have been used as pollen paleoclimatic proxies (Fig. 4cef and Fig. 5a). The Mediterranean forest (sum of Quercus total, Olea,Phillyrea and Pistacia) and xerophytes (sum of Artemisia,Ephedra and Amaranthaceae) have previously been shown to be a good indicator of climate changes in the Mediterranean region (Fletcher and S anchez Go~ ni, 2008;RamosRom an et al., 2018a). The PCI is based on the ratio of mesothermic taxa (sum of Quercus total, Olea,Fraxinus,Phillyrea,Acer, Betula,Alnus,Ulmus,Taxus,Salix,Pistacia, Corylus and Carpinus) divided by steppic taxa (sum of Artemisia,Ephedra,Hippopha€ eand Amaranthaceae) (Combourieu Nebout et al., 1999;Joannin et al., 2011), and has been useful in identifying climate changes mainly related to temperature in this region (Bertini et al., 2015;Camuera et al., 2019). However, some taxa included in the PCI also respond to different precipitation conditions, hence a more precise reconstruction of the precipitation in this area has been obtained using the Ip (Fletcher et al., 2010b). The Ip is expressed as: Ip ¼Quercus deciduous/(Artemisia þEphedra þAmaranthaceae þQuercus deciduous). 3.4. Inorganic geochemistry The inorganic elemental geochemistry data of the last ~200 kyr of the Padul-15-05 record show paleoenvironmental changes related to orbitaland suborbital-scale climate fluctuations Fig. 3. Median ages of HS1 boundaries, obtained by means of a KDE_plot function under a Bayesian modelling. The probability distribution of the Gaussian (light shades) and Bayesian-modelled ages (dark shades) are shown with the 2 s standard deviations and the median values (horizontal and small vertical lines below probability distributions). For the KDE Total, the horizontal lines below the distributions show the 1 s standard deviation. Ages for each non-tuned high-resolution records have been obtained according to our revision of the age boundaries from the paleoclimatic data to HS1 (green dots in Figure S2) and recalibrated to recent age calibrations (IntCal13, Marine13). Red vertical dashed lines indicate the median age distributions of the HS1 boundaries provided by the KDE Total (onset HS1 ¼17,970 cal yr BP; end HS1 ¼15,210 cal yr BP). Codes (C1, C2, C3 …) have been used for an easier identification of records. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) J. Camuera, G. Jim enez-Moreno, M.J. Ramos-Rom an et al. Quaternary Science Reviews 255 (2021) 106814 5 (Camuera et al., 2018). Silicon was used as indicator of siliciclastic input from the Sierra Nevada into the wetland and was taken as the most representative inorganic element for environmental reconstructions (Camuera et al., 2018). In this study, we present Si data with a sampling resolution of ~67 years for the time period between 20 and 11 kyr BP. This element has been used as Si normalized (Si norm ), represented by the silicon data divided by the sum of the total counts (in cps) from the most important elements (Si, K, Ca, Fe, Zr, Br, Sr, Al, Rb, Cl, Zn, Mn, S, Pb, U, Ni and Ti). As calcium presents very high values in carbonate lithologies, normalized Si values excluding Ca from the total counts is also shown in order to better observe changes in Si content in carbonate-rich sediments (Fig. 4b). 3.5. Test of statistical significance (SiZer) In order to observe the internal environmental oscillations occurring during HS1 and the significance level of these changes, the “Significant Zero crossing of derivatives”(SiZer) method was run on the 3-point moving average xerophyte percentages from Padul between 19 and 15 kyr BP. The SiZer method examines the derivatives of a curve for identifying the presence of a threshold and finds where the derivative of a function of an environmental variable changes significantly (Chaudhuri and Marron, 1999). This approach was done using the “SiZer”package from the R software (www.cran.r-project.org), using the locally weighted polynomial regression technique (Fan and Gijbels, 1996) for estimating the two derivatives and threshold. This analysis has also been used for identifying the significance of the environmental and climate changes from several past and present (paleo)environmental data (Giesche et al., 2019; Hald et al., 2004; Sonderegger et al., 2009). The first derivative in the SiZer analysis shows the general trends of statistically significant increasing (red color) or decreasing (blue color) functions, whereas the second derivative gives information about the curvature of the data, pointing into concave up (blue color) or concave down (red color) features. The SiZer map is represented in purple when the derivatives have values 0 (or close to 0) and cannot be concluded to have either decreasing or increasing functions (nor trends), whereas the gray areas indicate regions where the data are too sparse to make statements about significance (Chaudhuri and Marron, 1999)(Fig. 6). 3.6. Spectral analysis A cyclostratigraphic spectral analysis was performed on xerophyte percentages for the age range between 20 and 11 kyr BP from Fig. 4. Paleoclimatic raw data from the Padul-15-05 sediment core for the time period between 20 and 11 kyr BP: (a) Core photograph and age-depth model for the studied period. (b) Normalized silicon values, with calcium excluded (continuous line) and included (dashed line) from total counts (values inverted). (c) Percentage of xerophytes (values inverted). (d) Percentage of Mediterranean forest. (e) Precipitation Index (Ip). (f) Pollen Climate Index (PCI). Yellow shadings show the Younger Dryas (YD) and Heinrich Stadial 1 (HS1). Dark yellow shading within HS1 indicates the slightly warmer/wetter middle phase (HS1b). Blue arrows indicate the moderately warmer/wetter sub-phases within HS1, whereas orange arrows show the colder/more arid sub-phases. Dark and light red shadings show the Bølling-Allerød (BA): dark red shadings correspond to the warmer/wetter periods (similar to Greenland Interstadials 1e, 1c and 1a; Rasmussen et al., 2014), whereas light red shadings correspond to the colder/more arid events (similar to Greenland Interstadials 1 d and 1 b; Rasmussen et al., 2014). See Figure S3 for the detailed pollen diagram of the main pollen taxa. The complete pollen dataset can be found in the PANGAEA data repository (https:// doi.pangaea.de/10.1594/PANGAEA.904053). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) J. Camuera, G. Jim enez-Moreno, M.J. Ramos-Rom an et al. Quaternary Science Reviews 255 (2021) 106814 6 Fig. 5. Paleoclimatic data from Padul and Alboran Sea for the time period between 20 and 11 kyr BP: (a) Xerophyte data from Padul-15-05 with three-point moving average (values inverted). (b) SST (C) from ODP-976 record of Alboran Sea (Martrat et al., 2014). (c) SST (C) from MD95-2043 record of Alboran Sea (Cacho et al., 1999,2006). Within HS1, blue and orange arrows in the Padul record show the relatively humid and arid sub-phases, respectively. In the SSTs from Alboran Sea during HS1, blue arrows marked the warmer temperatures in relation with the relatively more humid sub-phases from Padul (HS1a.2, HS1b and HS1c.2). Vertical dashed lines (from right to left) show transitions between Last Glacial Maximum-HS1, HS1-BA, BA-YD and YD-Holocene for each study. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) Fig. 6. Test of statistical significance based on the SiZer method and run on (A) the 3-point moving average xerophyte percentages from Padul between 19 and 15 kyr BP. The yaxis in the first derivative (figure B) and second derivative (figure C) show the bandwidth parameter hin units log(h). Blue and orange arrows show (as in previous figures) the relatively humid and arid sub-phases, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) J. Camuera, G. Jim enez-Moreno, M.J. Ramos-Rom an et al. Quaternary Science Reviews 255 (2021) 106814 7 the Padul-15-05 record with the purpose of identifying cyclicities related to regional climate oscillations (Fig. 7). Xerophytes abundance has been proven to be a good proxy for regional moisture availability in this area (Pini et al., 2009;Ramos-Rom an et al., 2018b). The spectral analysis was carried out using the PAST 3.19 software (Hammer et al., 2001) with a REDFIT procedure of Schulz and Mudelsee (2002) under the rectangular window function (value of 2 for the segments parameter and value of 3 for the oversample parameter). In order to better observe cyclicities, xerophyte data were filtered under the statistically significant obtained frequencies using the Analyseries 2.0 software (Paillard et al., 1996). 4. Results and discussion 4.1. Chronology of HS1 in southern Europe and the Mediterranean region The obtained KDE ages for the revised, recalibrated and modelled HS1 boundaries of high-resolution records from southern Europe and the Mediterranean region including Padul-15-05 (KDE Total) exhibited ages of ~18.0 kyr BP (median value ¼17,970 cal yr BP; mean value ¼18030 ±330 cal yr BP) for the onset and ~15.2 kyr BP (median value ¼15,210 cal yr BP; mean value ¼15200 ±420 yr) for the end of HS1 (Fig. 3). The terrestrial paleoclimate record from Padul shows overall cold and arid conditions during HS1, deduced by the decrease in mesic forest and abundance of xerophytes between 18.4 and 15.7 kyr BP (Fig. 4c and d). In addition, the centennial-scale variability observed during HS1 allows for the identification of 3 main climatic phases (i.e., HS1a from 18.4 to 17.2 kyr BP, HS1b from 17.2 to 16.9 kyr BP, and HS1c from 16.9 to 15.7 kyr BP) and a further subdivision in 7 smaller-scale phases within them (i.e., HS1a.1, HS1a.2, HS1a.3, HS1b, HS1c.1, HS1c.2 and HS1c.3) (Fig. 4e and f and Fig. 5a), which are statistically supported by the significance test using the SiZer method (see section 4.2 and Fig. 6). 4.1.1. Three main climatic phases within HS1 from Padul The three main climatic phases occurring in Padul (HS1a, HS1b and HS1c) have also been observed in marine sediment records from the southern Iberian margin (e.g., U1389 and MD99-2339 sites), described as early, middle and late HS1 phases (Sierro et al., 2020). These oscillations previously observed in the marine sedimentary records were interpreted as consequence of changes in the physicochemical properties of surface water in the Iberian margin. The first main climatic phase in Padul during HS1 is HS1a (early HS1; 18.4e17.2 kyr BP), characterized by low temperatures with significant variability in precipitation but under generally arid conditions, deduced by high xerophytes and low PCI and Ip values (Fig. 4c, e, f). Especially cold/arid conditions during this early phase are confirmed by high Si norm values, which show that high siliciclastic input from the Sierra Nevada range into the wetland are caused by enhanced erosion during decreased forest cover (Camuera et al., 2019). The general cold/arid conditions shown in Padul during the early HS1a were also documented in nearby marine records presenting the 3 main phases for HS1, such as the pollen records from NW Iberia (Naughton et al., 2016), or the pollen data, SST reconstructions (Fig. 5b and c) and foraminifera/coccolithophore assemblages from Alboran Sea (Bazzicalupo et al., 2018; Cacho et al., 1999,2006;Fletcher and S anchez Go~ ni, 2008;Martrat et al., 2014). The early arid HS1a phase recorded in Padul could also Fig. 7. Cyclostratigraphic analysis of the Padul-15-05 pollen data. (a) Spectral analysis run on Padul xerophyte percentages for the age range between 20 and 11 kyr BP. (b) Percentages of xerophyte taxa from Padul-15-05. (c) Filtered xerophyte data based on the obtained 2055-yr cycle from the spectral analysis (bandwidth parameter of 0.0002). The identified three main phases (HS1a, HS1b and HS1c) within HS1 in Padul are in relation with this ~2000-yr cycle. (d) Filtered xerophyte data based on the 770-yr cycle from the spectral analysis (bandwidth parameter of 0.0005). The internal seven sub-phases (a.1-a.3, b, c.1-c.3) occurring during HS1 in Padul are in relation with this ~800-yr cycle. J. Camuera, G. Jim enez-Moreno, M.J. Ramos-Rom an et al. Quaternary Science Reviews 255 (2021) 106814 8 be related with the enhancement of the cold Portugal current and the low influence of the Azores Current observed in the marine records by the decrease in d 18 O surface seawater and the low SST (Sierro et al., 2020;Voelker et al., 2009). This cold SST could have produced a slowdown of the thermohaline overturning circulation and decreased of the marine heat transport, resulting in more intense Saharan winds and higher aridity (Jim enez-Espejo et al., 2008; Moreno et al., 2002), as observed in southern Iberian Peninsula and in the Padul record. In addition, the North Atlantic SST also affected the north-westerly wind intensity over Europe and the Mediterranean area, resulting in more arid conditions in the northern Mediterranean region (Moreno et al., 2004). HS1b (middle HS1; 17.2e16.9 kyr BP) is characterized in Padul by a moderate increase in temperature and precipitation, deduced from low xerophytes, and higher PCI and Ip values. This is further supported by low Si norm , indicative of low erosion precluding siliciclastic input in the wetland (Fig. 4b, c, e, f). A similar slightly warmer climate during this phase was recorded in the MD95-2043 (Cacho et al.,1999,2006) and ODP-976 (Martrat et al., 2014) records from Alboran Sea (Fig. 5b and c; subjected to age uncertainties for onset/ending of HS1). This warmer/wetter conditions agree with increases in temperate forest recorded in the Iberian margin (Daniau et al., 2007), and in runoff in Lake Estanya (NE Spain) (Morell on et al., 2009). This relatively more humid conditions in Padul, could be a result of the higher influence of the Azores Current, observed by the arrival of the warmer SST waters (and higher seawater d 18 O) to the Iberian margin (Sierro et al., 2020), also affecting the western Mediterranean Sea (Hodell et al., 2017; Martrat et al., 2014). HS1c (late HS1; 16.9e15.7 kyr BP) was climatically similar to HS1a, characterized by cold/dry conditions. This is deduced by the observed increase in xerophytes and Si norm and lowering in Ip between ~17 and 15.9 kyr BP, related with the decreasing moisture (Fig. 4b, c, e). The general cold/arid climate in Padul during this phase is concordant with low SST from Alboran Sea (Martrat et al., 2014)(Fig. 5b), and with increasing salinity and low lake level in Lake Estanya (see Mystery Interval in Fig. 7 from Morell on et al., 2009). This arid HS1c phase in Padul fits well with the late HS1 identified in the Iberian margin marine records, characterized by cooler and less saline marine water conditions (Sierro et al., 2020). The ice-rafted debris (IRD) H1.1 deposition along the Iberian margin (Eynaud et al., 2009; Voelker et al., 2009) at the beginning of this late phase could have been favored by the southward expansion of subpolar waters, which could have generated the coldest conditions recorded in Padul at ~16.5 kyr BP (Fig. 4f). 4.1.2. Subdivision of HS1 in seven centennial-scale climatic subphases in Padul The high-resolution pollen record from Padul-15-05 with the help of the cluster analysis pollen zonation also revealed shorter centennial-scale climatic variability during HS1a and HS1c with a further climatic subdivision of HS1 into 7 sub-phases. The statistical significance of these sub-phases has been validated by the SiZer test (see section 4.2 below and Fig. 6). The climatic subdivisions of HS1 and the correspondent pollen zones are: HS1a.1 (zone 4 g), HS1a.2 (zone 4f), HS1a.3 (zone 4e), HS1b (zone 4 d), HS1c.1 (zone 4c), HS1c.2 (zone 4 b) and HS1c.3 (zone 4a) (Fig. 4c, e, Fig. 5a and Fig. S3). HS1a.1 was characterized by a cold/arid phase between 18.4 and 17.8 kyr BP recorded by high xerophytes, and low PCI and Ip values. Climate changed towards more humidity in HS1a.2 sub-phase at 17.8e17.5 kyr BP, and returned to enhanced aridity during HS1a.3 between 17.5 and 17.2 kyr BP. This arid-humid-arid climatic pattern is further confirmed by oscillations in Si norm (Fig. 4b, c, e, f and Fig. 5a). HS1c also presents a three-phase subdivision, namely HS1c.1, HS1c.2 and HS1c.3. HS1c.1 was characterized by a decrease in precipitation and temperature (low Ip and lowest PCI values), registering the coldest conditions of HS1 at around 16.5 kyr BP (Fig. 4f). Temperature and moisture conditions increased during HS1c.2 at 16.4e16 kyr BP, whereas similar temperatures but under more arid climate conditions are recorded during HS1c.3 at 16e15.7 kyr BP (Fig. 4c, e, f and Fig. 5a). This arid-humid-arid climatic pattern is similar to the earlier HS1a. 4.2. Significance of the seven climatic sub-phases within HS1 from Padul A SiZer statistical method was used in this study to test the significance of the internal centennial-scale climate changes occurring during HS1 (Fig. 6). The first derivative of the 3-point moving average xerophyte data (Fig. 6) shows an increasing trend followed by a decreasing trend between 18.9 and 17.7 kyr BP, coinciding with the concave down function from the second derivative at ~18.2 kyr BP and indicating the first sub-phase of HS1 (HS1a.1). The following concave up (at ~17.7 kyr BP) and the concave down (at ~17.4 kyr BP) functions in the second derivative statistically show the HS1a.2 and HS1a.3 sub-phases, respectively. The mid phase HS1b is observed by the decreasing-increasing trends in the first derivative between ~17.3 and 16.8 kyr BP and in the concave up feature in the second derivative at ~17 kyr BP. The HS1c.1 sub-phase is statistically significant looking at the increasing trend in the first derivative between ~17 and 16.4 kyr BP. Although the HS1c.2 sub-phase exhibit a very weak signal, the small concave up function in the second derivative at ~16.2 kyr BP suggests that it is statistically significant. Finally, the last HS1c.3 is statistically well shown by the increasing-decreasing functions in the first derivative between ~16.1 and 15.5 kyr BP and by the concave down function in the second derivative at ~15.9 kyr BP. 4.3. HS1 record in southern Europe and the Mediterranean region The centennial-scale arid-humid-arid trends recorded during HS1a and HS1c, and the increase in temperature/precipitation during HS1b, are also observed in the SST records from the Alboran Sea in western Mediterranean (Cacho et al., 1999,2006;Martrat et al., 2014), suggesting a similar response in marine and continental environments (Fig. 5aec). However, our study from Padul, supported by the synthetic median ages of the HS1 boundaries obtained with the KDE modelling (Fig. 3) shows an early onset and end for HS1 (onset and end in Padul: ~18.4e15.7 kyr BP; KDE modelling for the onset and end in southern Europe and the Mediterranean region: ~18e15.2 kyr BP) with respect to the SST records from the Alboran Sea. This asynchronicity could be due to several reasons, including uncertainties in radiocarbon dating. Several previous studies suggested that radiocarbon dating of some specific materials in lacustrine environments (e.g., aquatic organisms) could be affected by a reservoir effect due to “old carbon”dissolved in the hard-water and used by algae and aquatic plants as carbon source (Olsson, 1986;Yu et al., 2007). However, d 13 C and C/N values from the radiocarbon samples used in dating this time interval (Table S1) seem to be in agreement with overall vascular C3 land plants (Meyers, 2003;Meyers and Lallier-verg es, 1999), suggesting a prevailing atmospheric organic carbon source and thus, a reduced or negligible reservoir effect. The age offsets between records could also be related with the significant decrease in the atmospheric 14 C between 17.5 and 14.5 kyr, making it difficult to obtain accurate age models based on radiocarbon dating for this time period (Broecker J. Camuera, G. Jim enez-Moreno, M.J. Ramos-Rom an et al. Quaternary Science Reviews 255 (2021) 106814 9