INTRODUCTION
Glacial–interglacial transitions provide precise archives for climate events. This is due to the influence of orbital parameters on the timing of events occurring in these intervals, as well as due to the expansion of the ice sheets and its effect on the oceans and moisture balance. Marine Isotope Stage (MIS) 5 interval (ca. 130–74 ka) is represented by the last interglacial period (LIG) and the onset of the last glacial (Capron et al., Reference Capron, Landais, Chappellaz, Schilt, Buiron, Dahl-Jensen and Johnsen2010). The MIS 5e substage covers the interval of 130–116 ka (130 through ca. 120 ka in the eastern Mediterranean [EM] and eastern Europe based on terrestrial archives [Demény et al., Reference Demény, Kern, Czuppon, Németh, Leél-Ossy, Siklósy and Lin2017]), during which the sea level is estimated to have been 6–9 m higher than the current global average (Kukla et al., Reference Kukla, Bender, de Beaulieu, Bond, Broecker, Cleveringa and Gavin2002; Shackleton et al., Reference Shackleton, Sánchez Goñi, Pailler and Lancelot2003; Dutton et al., Reference Dutton, Carlson, Long, Milne, Clark, DeConto, Horton, Rahmstorf and Raymo2015). The LIG and following glacial inception periods are climatically important periods since they have the potential to allow us to reconstruct the possible responses to climatic variations triggered by current global warming (Govin et al., Reference Govin, Capron, Tzedakis, Verheyden, Ghaleb, Hillaire-Marcel and St-Onge2015).
The continuous ca. 123 ka-long δ18O record of the North Greenland Ice Core Project (NGRIP) ice core is very valuable because it also includes a part of the LIG period (NGRIP members, 2004). High temperatures during the LIG are represented by high values of δ18O (−32‰) (NGRIP members, 2004). During the LIG, the Mediterranean woody taxa reached its highest abundance (Brauer et al., Reference Brauer, Allen, Mingram, Dulski, Wulf and Huntley2007). Later, this warm interval was interrupted by glacial expansion and some cold-water events in the North Atlantic (Kukla et al., Reference Kukla, Bender, de Beaulieu, Bond, Broecker, Cleveringa and Gavin2002). The NGRIP records indicate a gradual depletion of δ18O towards 116 ka with a transition to cold conditions. The end of the cooling phase is defined as the Greenland Stadial 26 (GS 26) interval, which marks the MIS 5e/5d transition.
A cooling of 4°C is suggested to have occurred between the C25 cold-water event and the last interglacial optimum (Oppo et al., Reference Oppo, McManus and Cullen2006). In addition, the planktonic δ18O record indicates the C26 cold-water event during early ice sheet growth, pointing to a 2–3°C decrease in North Atlantic sea surface temperature (Chapman and Shackleton, Reference Chapman and Shackleton1999) and the warm period of MIS 5e lasting ca. 9–10 ka that gradually ended with the C26 and C25 cold-water events (Oppo et al., Reference Oppo, McManus and Cullen2006). The C25 cooling event associated with the end of GS 26 also indicates a decrease in woody taxa, pointing towards a decrease in moisture availability in the eastern Anatolian pollen records (Pickarski et al., Reference Pickarski, Kwiecien, Djamali and Litt2015). Proxies of the LIG and the beginning of the glacial period have been discussed in detail using glacial (Johnsen et al., Reference Johnsen, Clausen, Dansgaard, Gundestrup, Hammer, Andersen and Andersen1997; NGRIP members, 2004), marine (McManus et al., Reference McManus, Bond, Broecker, Johnsen, Labeyrie and Higgins1994; Fronval and Jansen, Reference Fronval and Jansen1997; Chapman and Shackleton, Reference Chapman and Shackleton1999; Sánchez Goñi et al., Reference Sánchez Goñi, Eynaud, Turon and Shackleton1999, Reference Sánchez Goñi, Bakker, Desprat, Carlson, Van Meerbeeck, Peyron and Naughton2012; Oppo et al., Reference Oppo, Keigwin, McManus and Cullen2001, Reference Oppo, McManus and Cullen2006; Knudsen et al., Reference Knudsen, Seidenkrantz and Kristensen2002; Melki et al., Reference Melki, Kallel and Fontugne2010; Galaasen et al., Reference Galaasen, Ninnemann, Irvali, Kleiven, Rosenthal, Kissel and Hodell2014; Irvali et al., Reference Irvali, Ninnemann, Kleiven, Galaasen, Morley and Rosenthal2016), and speleothem records of the North Atlantic and European regions (Drysdale et al., Reference Drysdale, Zanchetta, Hellstrom, Fallick and Zhao2005, Reference Drysdale, Zanchetta, Hellstrom, Fallick, McDonald and Cartwright2007, Reference Drysdale, Hellstrom, Zanchetta, Fallick, Sánchez Goñi, Couchoud, McDonald, Maas, Lohmann and Isola2009; Boch et al., Reference Boch, Cheng, Spötl, Edwards, Wang and Häuselmann2011; Moseley et al., Reference Moseley, Spötl, Cheng, Boch, Min and Edwards2015; Vansteenberge et al., Reference Vansteenberge, Verheyden, Cheng, Edwards, Keppens and Claeys2016, Reference Vansteenberge, Verheyden, Genty, Blamart, Goderis, Van Malderen and Vanhaecke2019; Dumitru et al., Reference Dumitru, Onac, Polyak, Wynn, Asmerom and Fornós2018; Nicholson et al., Reference Nicholson, Jacobson, Hosfield and Fleitmann2021). The effects of climatic changes in the eastern Mediterranean occurring in both the last interglacial and glacial inception periods have been reported in several studies (e.g., Bar-Matthews et al., Reference Bar-Matthews, Ayalon, Gilmour, Matthews and Hawkesworth2003; Nehme et al., Reference Nehme, Verheyden, Noble, Farrant, Sahy, Hellstrom, Delannoy and Claeys2015; Columbu et al., Reference Columbu, Spötl, Waele, Yu, Shen and Gázquez2019), but are not fully known in Türkiye. For this reason, it is important to investigate the potential effects of sub-millennial climate changes during the last interglacial–glacial inception on the Anatolian Peninsula.
Due to its geographic setting, Türkiye is affected by various climatic domains (e.g., Atlantic, Mediterranean, and Asian; Rowe et al., Reference Rowe, Mason, Andrews, Marca, Thomas, van Calsteren, Jex, Vonhof and Al-Omari2012). Regional climate is dominated mainly by humid and cold air masses originating from the north and northwest, as well as by mid-latitude subtropical anticyclones (Akçar and Schlüchter, Reference Akçar and Schlüchter2005; Türkeş, Reference Türkeş, Harmancioglu and Altinbilek2020). The high-pressure systems initiating from Siberia, which are cold and dry in winter, and the extensions of the low-pressure monsoon system originating from Asia in summer are other climate systems affecting the region (Akçar and Schlüchter, Reference Akçar and Schlüchter2005; Türkeş, Reference Türkeş, Harmancioglu and Altinbilek2020; Jacobson et al., Reference Jacobson, Pickett, Gascoigne, Fleitmann and Elton2022). Precipitation in the region is mainly controlled by the North Atlantic Oscillation (NAO) (Türkeş, Reference Türkeş, Harmancioglu and Altinbilek2020). While strong precipitation conditions are observed in winter and autumn months in Türkiye during the NAO negative phase, weaker precipitation conditions prevail during the positive phases (Türkeş and Erlat, Reference Türkeş and Erlat2003, Reference Türkeş and Erlat2005). Therefore, considering its effect on winter precipitation, the Anatolian Peninsula is extremely important in order to explore the eastern extensions of climatic changes in the last glacial inception period.
Stalagmites (secondary cave carbonates) can physically and chemically record climatic and environmental events during the time of their formation in isotopic equilibrium with the dripping water (McDermott et al., Reference McDermott, Schwarcz, Rowe and Leng2006). Precise dating of these deposits up to 600,000 years using the U-series dating technique (White, Reference White, Sasowsky and Mylroie2004) has played an important role in establishing long-term terrestrial climate records. The chronological records of the LIG are generally based on the alignment (correlation) technique of various proxies, which may lead to some discrepancies in ages. High-precision age data obtained from speleothems, however, can provide good age-calibrated climate records (Govin et al., Reference Govin, Capron, Tzedakis, Verheyden, Ghaleb, Hillaire-Marcel and St-Onge2015).
Speleothem-based paleoclimate studies in Türkiye are mostly focused on the last ice age or the Holocene period, some with important linkages with early human settlements (Fleitmann et al., Reference Fleitmann, Cheng, Badertscher, Edwards, Mudelsee, Göktürk and Fankhauser2009; Göktürk et al., Reference Göktürk, Fleitmann, Badertscher, Cheng, Edwards, Leuenberger, Fankhauser, Tüysüz and Kramers2011; Rowe et al., Reference Rowe, Mason, Andrews, Marca, Thomas, van Calsteren, Jex, Vonhof and Al-Omari2012; Ünal-İmer et al., Reference Ünal-İmer, Shulmeister, Zhao, Uysal, Feng, Duc Nguyen and Yüce2015, Reference Ünal-İmer, Shulmeister, Zhao, Uysal and Feng2016a; Koç et al., Reference Koç, Koşun, Cheng, Demirtaş, Edwards and Fleitmann2020; Erkan et al., Reference Erkan, Bayarı, Fleitmann, Cheng, Edwards and Özbakır2021; Jacobson et al., Reference Jacobson, Flohr, Gascoigne, Leng, Sadekov, Cheng, Edwards, Tüysüz and Fleitmann2021, Reference Jacobson, Pickett, Gascoigne, Fleitmann and Elton2022). The oldest records from the region are from Sofular (Badertscher et al., Reference Badertscher, Fleitmann, Cheng, Edwards, Göktürk, Zumbühl, Leuenberger and Tüysüz2011) in NW Türkiye and Dim Cave (Rowe et al., Reference Rowe, Wickens, Sahy, Marca, Peckover, Noble, Özkul, Baykara, Millar and Andrews2020) from SSW Türkiye. The study on Dim Cave (Alanya) by Rowe et al. (Reference Rowe, Wickens, Sahy, Marca, Peckover, Noble, Özkul, Baykara, Millar and Andrews2020) focused on the beginning of the last interglacial period, but the climate record does not extend into the last interglacial–glacial transition. In addition, the Sofular Cave record characterizes different climatic conditions that do not reflect true Mediterranean climate due to significant humidity over the Black Sea. There is a lack of speleothem-based record of the last glacial inception in Türkiye. This study fills this gap by presenting new stable and radiogenic isotope (δ13C and δ18O, 87Sr/86Sr) and trace element (Ca, Mg, Sr, P, Cu, Y, Be, and Zr) data from a MIS 5 stalagmite (YL-1) collected from Yelini Cave (Fig. 1) and provides the first multi-proxy climate record of the MIS 5e/5d transition from the central Anatolia.
Locations of the records: 1 – Yelini Cave (Türkiye); 2 – Sofular Cave (Türkiye); 3 – Van Lake (Türkiye); 4 – Pentadactylos Cave (Cyprus); 5 – Kanaan Cave (Lebanon); 6 – Soreq and Peqiin caves (Israel); 7 – Tenaghi-Phillippon (Greece); 8 – Baradla Cave (Hungary); 9 – Schneckenloch and Hölloch caves (Austria and Switzerland); 10 – Corchia Cave (Italy); 11 – Han-Sur-Lesse Cave (Belgium); 12 – Campanet Cave (Mallorca); 13 – MD95-2040 core (Iberian margin); 14 – MD01-2443 core (Iberian margin); 15 – North Greenland Ice Core Project (NGRIP); 16 – MD03-2664 core (Labrador Basin). (a) Location and geological map of the study area (modified from Demıroğlu, Reference Demiroğlu2008). (b) Map of the Yelini Cave (modified from Nazik et al., Reference Nazik, Törk, Acar, Özel, Mengi, Aksoy, Tuncer, Güner, Ekmekçi and Başal2001).
REGIONAL CHARACTERISTICS
Yelini Cave
The Yelini Cave (also known as the Yılanlı Cave) is located near the town of Kayakent in the Günyüzü district of Eskişehir Province in central Türkiye. The cave is accessed through a narrow entrance along the southwest side of the Sivrihisar Mountains (Fig. 1). Yelini has archaeological importance because artifacts there include remains and cave paintings suggesting prolonged periods of human use in Neolithic times (Nazik et al., Reference Nazik, Törk, Acar, Özel, Mengi, Aksoy, Tuncer, Güner, Ekmekçi and Başal2001).
Yelini Cave (39°14'23.28“N, 31°48'25.33”E) (Fig. 1), located at an altitude of 1105 m above sea level (asl), is a horizontally extending fossil cave with a total length of 271 m (Nazik et al., Reference Nazik, Törk, Acar, Özel, Mengi, Aksoy, Tuncer, Güner, Ekmekçi and Başal2001) (Fig. 1b). Yelini Cave is dry except for the drip waters seeping from the ceilings during wet seasons (Nazik et al., Reference Nazik, Törk, Acar, Özel, Mengi, Aksoy, Tuncer, Güner, Ekmekçi and Başal2001) through a thin cover of limestone in the region. Its deepest point lies ~26 meters below the cave entrance. Relative humidity within the cave varies from 45% in the entrance to 79% in the deepest gallery where the YL-1 stalagmite was collected in 2019 (Fig. 1).
Geology
Yelini Cave lies within the Günyüzü Basin, which is a sub-basin of the Sakarya Basin (Demiroğlu and Örgün, Reference Demiroğlu and Örgün2010). Exposed lithologies include Paleozoic metamorphic rocks (Sivrihisar Metamorphics), Miocene sedimentary rocks, and Quaternary cover units (mainly alluvium) (Demiroğlu, Reference Demiroğlu2008) (Fig. 1). The cave itself formed within Paleozoic marbles, which have been intensely karstified and have been significantly affected by a NW–SE-trending fault (Fig. 1a; Nazik et al., Reference Nazik, Törk, Acar, Özel, Mengi, Aksoy, Tuncer, Güner, Ekmekçi and Başal2001; Şenoğlu, Reference Şenoğlu2006; Demiroğlu, Reference Demiroğlu2008).
Climate
The modern climate of the Anatolian Peninsula is heterogeneous. Located in the middle latitude temperate zone, Türkiye experiences four seasons throughout the year due to its close proximity to the sub-tropical zone (Türkeş, Reference Türkeş, Harmancioglu and Altinbilek2020). While the coastal regions are represented by a Mediterranean climate, the inland Anatolian Plateau experiences a continental climate (Bozkurt and Şen, Reference Bozkurt and Şen2009). Mountain belts extending parallel to the coast in the southern and northern regions prevent the marine climate from reaching the inner regions (Ünal et al., Reference Ünal, Deniz, Toros and İncecik2012). Therefore, precipitation received by the coastal regions is significantly higher than in the inner high plateaus of central Anatolia and the transportation of moisture to the interior is provided through the E–W trending valleys between the mountain ranges extending perpendicular to the Aegean coast in the western part of Anatolia (Bozkurt and Şen, 2009). Because of these topographical barriers, the central Anatolian region, where the Yelini Cave is situated, has semi-arid climate conditions (Altın et al., Reference Altın, Barak and Altın2012). For this reason, winters are cold and rainy, spring months are cool and rainy, and summers are mildly rainy (Türkeş, Reference Türkeş and Bolle2003). According to data from the Turkish General Directorate of Meteorology for the year 2020, the average annual temperature of the Eskişehir region is 11°C, and the average annual precipitation is 365 mm. The highest precipitation occurs in January, March, May, and June, while lowest precipitation is recorded in July and November (Fig. 2). Potential evapotranspiration exceeds precipitation in average conditions throughout the year (Türkeş, Reference Türkeş, Harmancioglu and Altinbilek2020).
The distribution of yearly regional rainfall in Türkiye (5-year average values during 2012–2016; Dilaver et al., Reference Dilaver, Aydın, Özyurt and Bayarı2018).
The predominance of winter precipitation in the region reflects the annual NAO effect. The dominant influence of NAO on Europe and the Mediterranean basin is particularly important in terms of received precipitation in western and central parts of Anatolia (Türkeş and Erlat, Reference Türkeş and Erlat2003). In the negative phase of NAO, weakening of Azores high pressure and shifting of maritime polar air masses to a southward route and evaporation over the western Mediterranean increase the extent of rainy conditions (Bozkurt and Şen, Reference Bozkurt and Şen2009).
MATERIAL AND METHODS
Sample description
The YL-1 stalagmite was collected from the interior of the cave far from the entrance in an area with limited ventilation (Fig. 1b). The total length of the stalagmite is 14.1 cm, and due to the absence of laminae and textural evidence for possible secondary crystallization and deformation, the upper 0.9 cm section was not studied (Fig. 3).
Photograph of the YL-1 stalagmite half showing the growth axes and locations of samples dated by U-series (left). Age–depth model generated on StalAge (Scholz and Hoffmann, Reference Scholz and Hoffmann2011; Scholz et al., Reference Scholz, Hoffmann, Hellstrom and Ramsey2012) using U-series ages of sub-samples at distances of 1, 45, 82, 98, 112, and 131 mm (right). Red lines indicate the age error boundaries.
The laminae are quite prominent along the stalagmite growth. Thinly developed laminae (1 mm thickness or less) of the sample are generally light brown, cream, beige, and white in color. A dark-brown layer is observed at 7.7–8.0 cm from the apex, which divides the stalagmite into two distinct growth periods (Fig. 3). Changing of growth axes is evident for these growth periods. The lower and upper axes are 5.0 and 7.5 cm long, respectively (Fig. 3). The width of the YL-1 stalagmite is ~6 cm, but there are minor variations in thickness.
Mineralogy
Thin-section petrography and XRD studies were carried out for mineralogical investigations. Three thin sections were prepared from one-half of the YL-1 stalagmite at the Ankara University Earth Sciences Application and Research Center (YEBİM) and textural changes of the laminae were examined under a polarizing microscope.
XRD analyzes were performed on laminae from which age data were obtained. The analyses were carried out with the Panalytic-Aeris device in YEBİM by selecting ~500 mg of powder sample for each.
U-Series dating
The U-series dating method was used for dating of the YL-1 stalagmite. Ten sub-samples, each weighing 100–150 mg, were collected from the stalagmite (Fig. 3). Microsampling was performed at the points of 1, 13, 25, 33, 45, 72, 82, 98, 112, and 131 mm from the top of the stalagmite. Sampling was carried out as much as possible from the impurity-free parts of the calcite layers. Drill bits used for microdrilling were cleaned with distilled water and ethyl alcohol after each sampling.
U-Series dating of the samples was performed at the University of Queensland (Australia) Radiogenic Isotope Laboratory (School of the Environment) using the Nu Plasma HR Multi-collector Inductively Coupled Plasma Mass Spectrometer (MC-ICP-MS). The analysis was carried out using the analytical protocols described by Zhao et al. (Reference Zhao, Hu, Collerson and Xu2001) and Ünal-İmer et al. (Reference Ünal-İmer, Shulmeister, Zhao, Uysal, Feng, Duc Nguyen and Yüce2015), and the results are reported in Table 1.
MC-ICP-MS U-Series age results (ka) of the YL-1 stalagmite (Yelini Cave). Ages in bold are out of stratigraphic order and were not utilized in StalAge calculations (Scholz and Hoffmann, Reference Scholz and Hoffmann2011).
Stable isotope analysis
For carbon and oxygen isotopic analyses, 127 calcite sub-samples were collected from the stalagmite at 1 mm intervals using an micro hand drill. In addition, one marble wall-rock sample was also analyzed for its carbon and oxygen isotope composition.
Carbon and oxygen isotope analyses were performed at the Environmental Isotope Laboratory of the University of Arizona (USA). δ18O and δ13C isotope compositions of carbonates were measured with a gas ratio mass spectrometer (KIEL-III) coupled with a carbonate preparation device (Finnigan MAT 252). Powder samples were treated with phosphoric dehydrate under vacuum at 70°C. Isotope ratios were analyzed with an accuracy of ± 0.1‰ for δ18O and ± 0.08‰ (1σ) for δ13C and calibrated based on repeat measurements of NBS-19 and NBS-18. All oxygen and carbon isotope ratios are shown in standard notation and reported relative to Vienna Pee Dee Belemnite (VPDB) (Supplementary Table 1) with calculated StalAge age results and host rock δ18O and δ13C isotope ratios.
In addition, for the Hendy Test (Table 2), carbon and oxygen isotope analyses were carried out on 10 sub-samples collected from 2 prominent layers of 27 and 103 mm from the top of the stalagmite. Microsampling was performed at ~1 mm intervals perpendicular to the growth axis of the stalagmite. Isotope analyses were conducted using the Termofinnigan Delta PlusXP mass spectrometer at the Central Laboratory of the Middle East Technical University (METU), Türkiye.
Results (δ18O and δ13C values, ‰ VPDB) of Hendy Tests.
87Sr/86Sr isotope and trace element concentrations
Nine sub-samples (eight from the stalagmite and one from marble host rock) were selected for trace element and strontium isotope (87Sr/86Sr) analyses. Sampling was performed at the sampling locations selected for dating. Analyses were carried out at the University of Queensland using the Thermo X-series ICP-MS device for trace element analysis and Thermal Ionization Mass Spectrometer (TIMS) for the strontium isotope analysis. Detailed methodology is given in Ünal-İmer et al. (Reference Ünal-İmer, Uysal, Zhao, Işık, Shulmeister, İmer and Feng2016b).
RESULTS
Calcite textures
Different calcite textures have been identified in the YL-1 stalagmite (Fig. 4a–c). Calcite crystals exhibit an irregular structure in the interval of ca. 116.16 to ca. 116.32 ka. Polycrystals of calcite are dispersedly intertwined (Fig. 4c) along with high porosity. A more-arranged calcite fabric appears at ca. 115.49 to ca. 115.83 ka, showing a distinctive layered structure. Layers are characterized by alternations of light- and dark-colored bands. While the light-colored layers vary between ~0.05–0.1 mm in thickness, dark-colored layers have thicknesses changing between 0.2 and 0.9 mm (Fig. 4b). A transition to a coarse-grained compact texture is observed in the layers at around 115.27 ka. Anhedral or subhedral mosaic calcite crystals have variable diameter size, ranging from ~10 μm to 1 mm. There is no significant porosity developed between the crystals (Fig. 4a).
Thin section images of YL-1 stalagmite under crossed nicols. The scales are 500 µm. (a) Transition of fine-grained calcite crystals to mosaic calcite crystals (dashed line), (b) distinctive layered structure of fine-grained calcite, (c) irregular (or dendritic) calcite crystals.
XRD results indicate that the main mineral in the stalagmite is calcite. In addition to calcite, the results also show minor aragonite peaks. However, no residual aragonite was detected petrographically in the thin sections.
U-series dating and age–depth model
Results of U-series dating analyses of 10 samples from the YL-1 stalagmite are presented in Table 1. Accordingly, deposition of the stalagmite started at 117.15 ± 0.60 ka and stopped at 109.05 ± 0.74 ka. The analyses at depths of 13, 25, 33, and 72 mm yielded dates beyond the growth period of the stalagmite. The outlier ages may be related to the transport of clastic thorium (232Th) in groundwater and its inclusion in the stalagmite (White Reference White, Sasowsky and Mylroie2004; Fairchild and Baker, Reference Fairchild and Baker2012). Therefore, the age–depth model was constructed through the StalAge code in R (Scholz and Hoffmann, Reference Scholz and Hoffmann2011; Scholz et al., Reference Scholz, Hoffmann, Hellstrom and Ramsey2012) using the ages (n = 6) in stratigraphic order within the margin of error (Fig. 3).
Stable isotopes
Calculated ages (ka) (StalAge results) and corresponding stable isotope values (in per mil [‰] relative to VPDB) are given in Supplementary Table 1. The δ18O and δ13C values of the host marble sample are 2.27‰ and −0.79‰, respectively (Supplementary Table 1). δ18O values of YL-1 stalagmite calcite samples change from −12.84‰ to −10.56‰ (VPDB) with an average value of −11.63‰ (Fig. 5a). Early layers of YL-1 are represented by relatively high δ18O values including the maximum measured δ18O composition of −10.56‰ at 116.83 (+0.58/−0.42) ka. Then, after a rapid negative shift, an enrichment of 1.33‰ is observed at 116.70 (+0.62/−0.35) through 116.63 (+0.46/−0.40) ka (Fig. 5a). This is followed by a decrease of 1.83‰ at 116.51 (+0.38/−0.50) ka. This period corresponds to the lowest δ18O values throughout the stalagmite (−12.84‰). After this low period, a rapid enrichment of 1.88‰ occurs towards 116.24 (+0.53/−0.86) ka, and a gradual decrease is noticeable towards 115.87 (+0.83/−1.71) ka. From this point onwards, calcite deposition ceases with a gradual enrichment in isotope ratios.
(a) Diagram showing variations in δ18O and δ13C profiles of YL-1 stalagmite with respect to its depth. (b) Correlation between δ18O and δ13C values.
The δ13C values of the YL-1 stalagmite are in the range of −9.63‰ to −7.89‰ (with an average of −8.78‰, VPDB) (Fig. 5a). Variation in δ13C compositions along sample YL-1 is less pronounced compared to variation in δ18O ratios. Earliest stalagmite layers (at 117.05 [+0.48/−0.38] ka) are characterized by relatively enriched δ13C values (−7.90‰). This is followed by a relative depletion in δ13C values (−9.01‰) at 117.01 (+0.51/−0.37) ka and an enrichment of 1.01‰ at 116.65 (+0.51/−0.39) ka. Then δ13C values show a negative shift after a slight enrichment in the range of 116.16 (+0.60/−1.03) to 116.07 (+0.67/−1.24) ka.
There is a weak positive correlation between the δ13C and δ18O values of the YL-1 stalagmite (R 2 = 0.31) (Fig. 5b). The sub-sample from 1 mm (distance from top) has significantly high oxygen and carbon isotope values (−8.41‰ and 0.74‰, respectively). If the dripping is extremely slow or interrupted during the end of calcite deposition, isotopic enrichment may occur as a result of prolonged separation of CO2 from the solution (Frisia et al., Reference Frisia, Borsato, Fairchild and McDermott2000), resulting from low effective infiltration inside the cave.
Hendy Tests
The kinetic fractionation in stalagmites may mislead the paleoclimate data obtained from the calcite (Hendy, Reference Hendy1971). This causes the light isotopes to be preferentially separated from the solution and concomitant enrichment of heavy isotopes in the stalagmite. For this reason, we applied the Hendy Test (Hendy, Reference Hendy1971) to the YL-1 stalagmite. According to the Isotopic Equilibrium Test study, the change in δ18O values along a single lamina is expected not to exceed 0.5‰ (Hendy Reference Hendy1971; Ünal-İmer et al., Reference Ünal-İmer, Uysal, St Pierre, Zhao and Shulmeister2020). δ18O values of YL-1 stalagmite do not show such a variation (Table 2; Fig. 6). In addition, farther from the axis there are no reciprocal changes between δ13C and δ18O that indicate kinetic fractionation (Fig. 6). Therefore, it is suggested that CO2 separation from dripping water was slow during stalagmite formation (Ünal-İmer et al., Reference Ünal-İmer, Uysal, St Pierre, Zhao and Shulmeister2020) and isotopic equilibrium conditions were maintained.
Hendy Test plots for the δ18O and δ13C analyses performed on laminae at 27 and 103 mm of YL-1.
How effectively the Hendy Test applications in stalagmites reflect the isotopic equilibrium is still controversial since the thickness of stalagmite layers at the growth center and edges may vary. Different kinetic fractionation effects may occur in the same layer along with possible problems with the sampling (Dorale and Liu, Reference Dorale and Liu2009). In addition, it is not possible to establish a precise oxygen isotope equilibrium between calcite and dripping water (McDermott et al., Reference McDermott, Schwarcz, Rowe and Leng2006). The Hendy Test is widely used to investigate such kinetic effects, although it does not give as precise results as the comparative isotope records from different stalagmites (Dorale and Liu, Reference Dorale and Liu2009).
87Sr/86Sr and trace element compositions
87Sr/86Sr ratios of YL-1 stalagmite vary in a narrow range from 0.707946 to 0.707987 (Table 3). These values are very close to the Sr isotope ratio of the host rock (0.707851), indicating a close genetic link between the marble host rock and the YL-1 stalagmite. The highest and lowest 87Sr/86Sr ratios in stalagmite are recorded at 71 mm and 82 mm (at the beginning of a dark layer), respectively.
ICP-MS trace elements concentrations and TIMS 87Sr/86Sr ratios of the YL-1 stalagmite and marble host rock.
Magnesium concentrations of the YL-1 stalagmite range between 764.4 ppm and 614.2 ppm (Table 3). Lowest Mg concentrations were observed at layers of 112 mm and 98 mm, whereas highest concentrations were recorded at 82 mm and 33 mm. Sr concentrations, ranging 88.9–110.6 ppm, are lower than Mg. The lowest and highest Sr contents are recorded in laminae at 82–33 mm and 13–25 mm, respectively (Table 3). Mg/Ca and Sr/Ca ratios are only weakly correlated (R 2 = 0.15; Fig. 7a). Phosphorus contents are 50.3–119.8 ppm, and the lowest and highest P concentrations were observed at 82 and 25 mm, respectively. Lower Cu concentrations ranging from 1.8 to 4.6 ppm have a strong positive correlation with P (R 2 = 0.75) (Fig. 7b). Additionally, Y (2–6 ppb), Be (12–43 ppb), and Zr (49–96 ppb) concentrations are very low. Beryllium and Phosphorus are positively correlated (R 2 = 0.73), whereas P weakly relates with Y and Zr (R 2 = 0.12 and R 2 = 0.18, respectively; Fig. 7c).
(a) Weak positive correlation between Sr/Ca and Mg/Ca. (b) Moderate positive correlation of Cu and Be concentrations with P concentration. (c) Weak positive correlation of Zr and Y concentrations with P concentration.
DISCUSSION
Speleothem growth and deposition rate
Speleothem formation is an important indicator of climatic conditions since the growth of speleothem depends crucially on water availability (Genty et al., Reference Genty, Blamart, Ghaleb, Plagnes, Causse, Bakalowicz and Zouari2006). Although growth rate is an important proxy, growth mainly responds to changes in temperature, precipitation, and soil activities (Kaufmann and Dreybrodt, Reference Kaufmann and Dreybrodt2004). Sufficient Ca2+ dissolution from the marble wall-rock and the amount of CO2 in the soil, which controls the solution acidity, are closely correlated with the extent of calcite precipitation in speleothems (Kaufmann, Reference Kaufmann2003; Vansteenberge et al., Reference Vansteenberge, Verheyden, Genty, Blamart, Goderis, Van Malderen and Vanhaecke2019). In addition, water loss through evaporation in the unsaturated region, evapotranspiration, and surface flow also can reduce the effect of dripping water (Vaks et al., Reference Vaks, Bar-Matthews, Ayalon, Schilman, Gilmour, Hawkesworth, Frumkin, Kaufman and Matthews2003). Therefore, lowering of the P/E ratio may cause a decrease in growth rate of speleothems or even hiatuses in speleothem growth. Speleothem growth occurs under humid climatic conditions (Genty et al., Reference Genty, Blamart, Ghaleb, Plagnes, Causse, Bakalowicz and Zouari2006; Fairchild and Baker, Reference Fairchild and Baker2012). Based on the age–depth model, the YL-1 stalagmite was deposited at 117.13 (+0.57/−0.44) to 114.87 (+1.63/−2.89) ka and is characterized dominantly by thin layers of fine-grained calcite (Fig. 4).
The growth rate of YL-1 is found to be ~0.10 mm/yr from the beginning towards ca. 116.94 ka (Fig. 8). After this period, the deposition rate increased to 0.12 mm/yr and slightly dropped to 0.10 mm/year at ca. 116.83 ka. The first notable decrease in growth rate to 0.63 mm/yr is seen at ca. 116.67 ka, which corresponds to the beginning of the detrital-rich brown layer at ca. 116.67 ka (77–80 mm depth; Fig. 3). Age results, however, suggest that the hiatus here was not prolonged. A second noticeable growth rate decrease, down to 0.04 mm/yr, appears at ca. 116.08 ka where the calcite fabric displays interrupted crystallization and dendritic-type porous fabric (Figs. 4 and 8). Generally slow growth rate throughout the YL-1 deposition and even slower growth rates in the upper 1 ka section are linked to the lowering P/E (precipitation/evaporation) ratio in the cave. The onset of stalagmite deposition coincides with the GS 26 cold phase represented in the Greenland ice core records (NGRIP members, 2004). This indicates the presence of a positive moisture balance in the region during the cold phase, which is not surprising considering the semi-arid climatic conditions of the central Anatolian region. Under current climatic conditions, a relatively low annual precipitation rate of 300–500 mm/yr and average summer temperatures reaching 23°C cause sub-cloud and/or surface evaporation in the region (Schemmel et al., Reference Schemmel, Mikes, Rojay and Mulch2013). Therefore, “effective moisture” (precipitation−evaporation) controlling speleothem growth is of great importance (Hodge et al., Reference Hodge, Richards, Smart, Andreo, Hoffmann, Mattey and González-Ramón2008). Increased overall precipitation and lower temperatures in the region indicate higher “effective precipitation” or positive P–E.
Temporal variations of δ18O, δ13C, 87Sr/86Sr, Mg/Ca, Sr/Ca, P (ppm), Cu (ppm), Be (ppb), Y (ppb), Zr (ppb), and growth rate (mm/yr) of the YL-1 stalagmite. The anomalies in all proxies between 116.67 (+0.54/−0.38) ka and 116.65 (+0.51/−0.39) ka are noted (green vertical bar). δ18O at 116.24 (+0.58/−0.86) ka (MIS 5e/5d transition) increases by 1.79‰ (yellow vertical bar). δ13C values record this transition at 116.07 (+0.67/−1.24) ka as a weaker anomaly (yellow vertical bar). The decrease in growth rate is also noticeable during this period.
δ18O and δ13C variations
In order to define the mechanisms controlling changes in δ18Occ, it is important to outline controls on the δ18Owater in that region. There are a number of “isotopic effects” that cause variations in δ18Owater, including the amount of precipitation, latitude differences, ice volume changes, seasonal changes, temperature, and changes in moisture source (Lachniet, Reference Lachniet2009). In addition, other factors such as evaporation at the surface and in the epikarst region, water residence time in the karstic zone, groundwater mixing, and calcite precipitation rate provide important controls on δ18Occ (McDermott et al., Reference McDermott, Schwarcz, Rowe and Leng2006; Rowe et al., Reference Rowe, Mason, Andrews, Marca, Thomas, van Calsteren, Jex, Vonhof and Al-Omari2012). Hence, the factors acting on δ18Occ represent a complex system (Lachniet, Reference Lachniet2009; Rowe et al., Reference Rowe, Mason, Andrews, Marca, Thomas, van Calsteren, Jex, Vonhof and Al-Omari2012).
At Yelini, the rock column (marble/recrystallized limestone) above the cave has low thickness (5 m at most) (Nazik et al., Reference Nazik, Törk, Acar, Özel, Mengi, Aksoy, Tuncer, Güner, Ekmekçi and Başal2001). This might have reduced the possibility of mixing of dripping water with groundwater having different isotopic signatures due to lower residence times (Lachniet, Reference Lachniet2009). Therefore, it is thought that atmospheric processes exerted a major control on δ18Occ.
Previous studies on speleothem δ18Occ in the Mediterranean region examined the “amount effect,” which is based on the inverse correlation between rainfall amount and δ18Oppt, the “source effect” focusing on Atlantic- vs. Mediterranean-sourced humidity that varies during stadial and interstadial transitions, and the “temperature effect” explaining the altitude-based inverse correlation between cold conditions and δ18Owater (Bar-Matthews et al., Reference Bar-Matthews, Ayalon, Kaufman and Wasserburg1999; Jex et al., Reference Jex, Baker, Fairchild, Eastwood, Leng, Sloane, Thomas and Bekaroğlu2010; Boch et al., Reference Boch, Cheng, Spötl, Edwards, Wang and Häuselmann2011; Rowe et al., Reference Rowe, Mason, Andrews, Marca, Thomas, van Calsteren, Jex, Vonhof and Al-Omari2012; Dumitru et al., Reference Dumitru, Onac, Polyak, Wynn, Asmerom and Fornós2018; Columbu et al., Reference Columbu, Spötl, Waele, Yu, Shen and Gázquez2019). In addition to these mechanisms, there are also some studies in Türkiye that have reported a negative correlation between effective winter rainfall and δ18Owater. The speleothem record of Akçakale Cave in NE Türkiye showed a strong and consistent correlation between the October–January rainfall records and δ18Occ values (Jex et al., Reference Jex, Baker, Fairchild, Eastwood, Leng, Sloane, Thomas and Bekaroğlu2010), which is attributed to recharge of the karst aquifer by the autumn–winter precipitation with depleted δ18Owater values. The isotope data of dripping waters, local springs, and streams around the Karaca Cave, which is located in southern Türkiye, also display a strong linkage with contribution from winter rainfall (Rowe et al., Reference Rowe, Mason, Andrews, Marca, Thomas, van Calsteren, Jex, Vonhof and Al-Omari2012). However, in both of these examples, the overburden above the cave is thick with high rainfall amounts during the winter, allowing for longer residence time for infiltrated water, which is less exposed to evaporation and water mixture inside the epikarst. Therefore, the P−E balance is more biased towards winter season, allowing a negative correlation with δ18Occ. In the case of Yelini Cave, the evaporation component in the P−E balance might be important due to atmospheric circulation at the geographical location of the studied cave. The barrier effect of the mountain ranges in the northern and southern parts of Türkiye prevents moisture from moving inland (Erkan et al., Reference Erkan, Bayarı, Fleitmann, Cheng, Edwards and Özbakır2021), resulting in a prevalence of semi-arid climatic conditions around Yelini Cave.
In arid and semi-arid areas, near-surface evaporation processes are reflected as seasonal variations in δ18O (McDermott et al., Reference McDermott, Schwarcz, Rowe and Leng2006, and references therein). Schemmel et al. (Reference Schemmel, Mikes, Rojay and Mulch2013) reported a modern meteoric water study suggesting that the general trend of increase in δD and δ18O from the sheltered sides of the mountain ranges in the southern and northern parts of Türkiye to the interior Central Anatolian Plateau (CAP) is the result of increased sub-cloud and surface evaporation due to drought conditions. Increased evaporation conditions significantly reduce effective precipitation (P–E) in the region. This also indicates significant seasonal changes on δ18Owater values. The δ18O rainfall data for the years between 2012 and 2016 taken from the Turkish State Meteorological Service (MGM) Ankara station, which is the closest meteorological station to Yelini Cave, indicate that winter precipitation is represented by more-negative δ18O values compared to other months (Dilaver et al., Reference Dilaver, Aydın, Özyurt and Bayarı2018). Similarly, stable water isotope data from CAP (Ankara) and Global Network of Isotopes in Precipitation (GNIP) Ankara station data suggest a strongly seasonal (winter dominant) rainfall pattern along the CAP (Schemmel et al., Reference Schemmel, Mikes, Rojay and Mulch2013). This is explained by a positive precipitation–evaporation balance during the cold season and winter-dominated groundwater recharge. Consequently, increased or decreased intensity of winter precipitation is the most likely mechanism reflecting the changes in δ18Occ values of the YL-1 stalagmite (Fig. 8). Lower δ18Occ values of the YL-1 stalagmite are attributed to the increased effect of net effective winter precipitation and higher δ18Occ values are due to decreasing net effective precipitation arising from the increase in atmospheric evaporation. Adding to this is the plausible evaporation within the epikarst of Yelini Cave, due to the location of the cave close to the surface.
Changes in speleothem δ13Ccc are directly associated with variations in isotopic contributions of carbon reservoirs including soil and atmospheric CO2 as well as the carbonate host rock (McDermott et al., Reference McDermott, Schwarcz, Rowe and Leng2006; Dumitru et al., Reference Dumitru, Onac, Polyak, Wynn, Asmerom and Fornós2018). δ13C values are affected by fluctuations in isotopically light soil CO2 (Hodge et al., Reference Hodge, Richards, Smart, Andreo, Hoffmann, Mattey and González-Ramón2008), therefore, speleothem δ13Ccc values strongly depend on carbon isotope composition of vegetation (C3 and C4 plants) growing above the cave (McDermott et al., Reference McDermott, Schwarcz, Rowe and Leng2006). Changes in δ13Ccc values also may be affected by several factors that include stalactite–stalagmite fractionation (evaporation and CO2 degassing), CO2 gas loss of groundwater in the aquifer zone, and short water-residence times in the soil (i.e., inability to establish an equilibrium between soil water and soil CO2) (Baker et al., Reference Baker, Ito, Smart and McEwan1997).
Biogenic CO2 increase due to increase in rainfall during hot and humid periods has a negative-trending effect on δ13Ccc (Genty et al., Reference Genty, Blamart, Ouahdi, Gilmour, Baker, Jouzel and Van-Exter2003). C4-type plants reflecting arid environments and C3-type plants reflecting humid environments cause different δ13C values in speleothems depending on the climate (δ13C for C4 and C3 plants are −6 to +2‰ and −14 to −6‰, respectively) (McDermott et al., Reference McDermott, Schwarcz, Rowe and Leng2006). On the other hand, decreasing CO2 contribution of biogenic carbon during cold and dry periods, increasing dissolved carbon content, PCP, primary carbonate host rock, and atmospheric dissolved carbon concentrations result in high δ13Ccc values in speleothems (Baker et al., Reference Baker, Ito, Smart and McEwan1997; Genty et al., Reference Genty, Blamart, Ghaleb, Plagnes, Causse, Bakalowicz and Zouari2006; Dumitru et al., Reference Dumitru, Onac, Polyak, Wynn, Asmerom and Fornós2018). In general, the δ13C values of a speleothem are products of this equilibrium. In the YL-1 stalagmite, δ13C values varying from −14 to−6‰ indicate C3-type vegetation (McDermott et al., Reference McDermott, Schwarcz, Rowe and Leng2006) similar to the other cave records from Anatolia (e.g., Ünal-İmer et al., 2015). Modern vegetation is quite scarce at the studied cave site, but the δ13C profile of the YL-1 stalagmite is associated with changes in effective winter precipitation (precipitation−evaporation). The changes towards more negative values can be explained by rises in net effective precipitation and growing seasonal moisture in the soil and carbonate rock/atmospheric components during dry seasons (Fig. 8).
87Sr/86Sr ratios
Sr isotope values (87Sr/86Sr) measured in speleothems reflect the relative contributions of Sr dissolved from the marble host rock and exogenic Sr accumulating in the soil (Goede et al., Reference Goede, McCulloch, McDermott and Hawkesworth1998; Bar-Matthews et al., Reference Bar-Matthews, Ayalon, Kaufman and Wasserburg1999; Frumkin and Stein, Reference Frumkin and Stein2004; Li et al., Reference Li, Ku, You, Cheng, Edwards, Ma, Tsai and Li2005). The negligible extent of isotopic fractionation of 87Sr/86Sr during calcite formation enables assessment of environmental factors affecting the 87Sr/86Sr values of dripping water in a cave (Goede et al., Reference Goede, McCulloch, McDermott and Hawkesworth1998). In addition to the Sr isotope composition of the host rock, different radiogenic Sr sources (e.g., dust transport, marine salts, volcanic activity) are reflected in the 87Sr/86Sr values of the speleothems (Banner et al., Reference Banner, Musgrove, Asmerom, Edwards and Hoff1996; Goede et al., Reference Goede, McCulloch, McDermott and Hawkesworth1998; Bar-Matthews et al., Reference Bar-Matthews, Ayalon, Kaufman and Wasserburg1999; Verheyden et al., Reference Verheyden, Keppens, Fairchild, McDermott and Weis2000; Frumkin and Stein, Reference Frumkin and Stein2004; Li et al., Reference Li, Ku, You, Cheng, Edwards, Ma, Tsai and Li2005). This provides insights into the origin of radiogenic Sr, the effect of varying climatic conditions, and changes in the karst.
Previous studies have shown that weathering of dolomitic host rock with an increase in precipitation might cause low 87Sr/86Sr values in speleothems (Ayalon et al., Reference Ayalon, Bar-Matthews and Kaufman1999). Verheyden et al. (Reference Verheyden, Keppens, Fairchild, McDermott and Weis2000) reported the dominant effect of low 87Sr/86Sr bedrock endmembers and the effect of high 87Sr/86Sr clastic-phase endmembers with their time-dependent contributions. In addition to the dominant host rock contribution represented by low 87Sr/86Sr values during humid periods, terrestrial dust can cause high 87Sr/86Sr ratios in speleothems (Frumkin and Stein, Reference Frumkin and Stein2004). The dissolution of silica-rich metamorphic rocks (as cave basement) or quartzo-feldspathic clastics in bedrock also may cause high 87Sr/86Sr in speleothems (Ünal-İmer et al., Reference Ünal-İmer, Shulmeister, Zhao, Uysal and Feng2016a; Carolin et al., Reference Carolin, Walker, Day, Ersek, Sloan, Dee, Talebian and Henderson2019). 87Sr/86Sr ratios of the YL-1 stalagmite mostly indicate genetic signatures that are similar to the host carbonate rock (0.707851) (Fig. 9). Therefore, the contribution of terrestrial dust (0.7082–0.7086; Frumkin and Stein, Reference Frumkin and Stein2004) input (from the Sahara Desert) is thought to be limited (Ünal-İmer et al., Reference Ünal-İmer, Shulmeister, Zhao, Uysal and Feng2016a).
Graph showing the variation of 87Sr/86Sr with respect to stalagmite depth and U-series ages. 87Sr/86Sr ratios of the stalagmite indicate contribution from the marble host rock. 87Sr/86Sr values of the Jerusalem Cave stalagmite (Frumkin and Stein, Reference Frumkin and Stein2004) are taken as the reference for terrestrial dust from the Sahara Desert.
87Sr/86Sr ratios increase from 0.707968 to 0.707985 in the period of 116.94 (+0.49/−0.35) ka to 116.83 (+0.58/−0.42) ka in the YL-1 stalagmite (Fig. 9). A decrease is observed at 116.67 (+0.54/−0.38) ka to the value of 0.707946, which is concurrent with enrichments in the δ18O and δ13C records (Fig. 8). 87Sr/86Sr values initially show a similar increasing trend with δ18O; however, they then follow a different pattern at 116.67 (+0.54/−0.38) ka (Fig. 8). This difference between O–C systematics and 87Sr/86Sr ratios is likely due to the thin soil cover above the cave, increased groundwater residence time, and increase in the extent of water–host rock interaction possibly due to the decrease in groundwater recharge (Rowe et al., Reference Rowe, Wickens, Sahy, Marca, Peckover, Noble, Özkul, Baykara, Millar and Andrews2020). Lowered levels of 87Sr/86Sr at 116.67 (+0.54/−0.38) ka can be attributed to the negative P/E balance due to the decrease in winter precipitation (Fig. 9). This indicates higher contributions from the marble host rock as the water–rock interaction intensified. 87Sr/86Sr ratio rises to 0.707987 at 116.51 (+0.38/−0.50) ka, which might reflect rapid return to wetter conditions after drought (Fig. 9).
Trace element concentrations
Variations in trace element compositions of speleothems result from changes in environmental conditions, which include chemical composition of the dripping water, the amount of precipitation, the types of elements in the soil, water–rock interaction, chemical interactions during the entrance of water to the cave, and trace element exchange reactions during calcite deposition (White, Reference White, Sasowsky and Mylroie2004). Water with a certain PCO2 leaking from the soil zone dissolves calcium carbonate until it attains equilibrium, and absorbs trace elements (White, Reference White, Sasowsky and Mylroie2004). Some elements react rapidly with the infiltrated water while others are leached by slowly flowing seepage waters (Fairchild and Treble, Reference Fairchild and Treble2009).
In the YL-1 stalagmite, both Sr/Ca and Mg/Ca ratios tend to increase from the beginning of deposition up to 116.67(+0.54/-0.38) ka (Fig. 8), despite their low sampling resolution. Sr/Ca ratio attains its maxima at 116.67(+0.54/-0.38) ka. This might be related to the increased water residence times in the epikarst. The soil CO2 increased by the vegetation controls the bedrock dissolution and has a significant effect on the chemical composition of the water entering the cave (Hellstrom and McCulloch, Reference Hellstrom and McCulloch2000). This may cause changes in the amount of calcite or dolomite dissolved, and trace element concentrations in the solution may also change accordingly (Fairchild et al., Reference Fairchild, Borsato, Tooth, Frisia, Hawkesworth, Huang, McDermott and Spiro2000).
Depending on the decrease in the amount of precipitation, residence time of water in the epikarst may increase, which is accompanied by elevated concentrations of elements such as Mg and Sr as a result of extended water–rock interaction (Fairchild et al., Reference Fairchild, Borsato, Tooth, Frisia, Hawkesworth, Huang, McDermott and Spiro2000; Hellstrom and McCulloch, Reference Hellstrom and McCulloch2000). However, linear relationships have been reported between Mg content and increase in the deposition temperature of speleothems (Şenoğlu, Reference Şenoğlu2006; Zhou et al., Reference Zhou, Wang, Huang and Mai2011). Şenoğlu (Reference Şenoğlu2006) discussed this increase using a Holocene-aged stalagmite, also from the Yelini Cave, that is located in a semi-arid climate zone. Changes in Mg/Ca ratios in a specific period of the stalagmite are associated with dry (or humid) conditions. Because the stalagmite studied by Şenoğlu (Reference Şenoğlu2006) was collected close to the entrance of the cave, it probably was affected by changes in air temperature. However, the YL-1 stalagmite was sampled from the distal part of the cave from the entrance (Fig. 1b) so that it should be represented by negligible airflow (Fohlmeister et al., Reference Fohlmeister, Voarintsoa, Lechleitner, Boyd, Brandtstätter, Jacobson and Oster2020).
Concurrent increases of δ13C and Mg/Ca-Sr/Ca ratios in YL-1 at ca. 116.65 ka are attributed to the Prior Calcite Precipitation (PCP) effect during a probable drought event. In addition to the dissolution process, proportional increases in Mg/Ca and Sr/Ca may also occur due to PCP along the flow path (Fairchild et al., Reference Fairchild, Borsato, Tooth, Frisia, Hawkesworth, Huang, McDermott and Spiro2000; Fairchild and Treble, Reference Fairchild and Treble2009). As the PCO2 decreases in the aquifer zone due to drought, CO2 is separated from the solution, which then becomes supersaturated with respect to calcite (Fairchild et al., Reference Fairchild, Borsato, Tooth, Frisia, Hawkesworth, Huang, McDermott and Spiro2000). In this case, Mg/Ca and Sr/Ca ratios may increase with the separation of calcium from the solution (Fairchild and Treble, Reference Fairchild and Treble2009), leading to a simultaneous enrichment of Mg/Ca and Sr/Ca ratios and δ13C compositions of speleothems. This is because the change in speleothem δ13C values is controlled by changes in plant thickness and soil CO2. The simultaneous increases in YL-1 Mg/Ca, Sr/Ca, and δ13C at 116.67 (+0.54/−0.38) ka suggest that calcite was deposited along the flow path and PCP is the dominant mechanism, which could be considered as evidence of a decrease in soil CO2 during that time. A subsequent negative shift in δ13C at 116.51 (+0.38/−0.5) ka is also reflected in variations in Sr/Ca ratios as well as in relatively less–pronounced variations in Mg/Ca ratios (Fig. 8). This can be regarded as a transition to humid conditions during which the soil CO2 increases, possibly due to elevated vegetation density (Fairchild and Treble, Reference Fairchild and Treble2009).
We also examined concentrations of P, Cu, Y, Be, and Zr in the YL-1 stalagmite despite low-resolution sampling for these trace elements (Figs. 8 and 10). Studies on P indicate the existence of a linear relationship between rainfall and P concentration (Huang et al., Reference Huang, Fairchild, Borsato, Frisia, Cassidy, McDermott and Hawkesworth2001; Borsato et al., Reference Borsato, Frisia, Fairchild, Somogyi and Susini2007). The increase in the concentrations of P, Y, and Cu in an Italian stalagmite (Grotta di Ernesto Cave, ER78 stalagmite) was explained by transport of these elements as colloidal particles from the soil zone during high water seepage in autumn (Borsato et al., Reference Borsato, Frisia, Fairchild, Somogyi and Susini2007). P, Y, Cu, and Be concentrations in YL-1 show low values in the stadial period, whereas Zr contents are relatively high during the same interval (< 116.2 ka) (Fig. 8). Simultaneous decreases in P, Y, Cu, Be, and Zr at 116.67 (+0.54/−0.38) ka might be due to the decrease in the winter precipitation. During GI 25, P, Cu, Be, and Zr, as well as δ13C and δ18O, follow a similar (contemporaneous) change (Fig. 10). Relatively high P, Y, Cu, and Be concentrations at the onset of GI 25 (115.87 ka) are associated with increased net effective precipitation and moisture availability in the soil zone (Figs. 8 and 10).
Temporal variations during the period 115–116.2 ka of (a) δ18O, (b) δ13C, (c) Cu (ppm), (d) Zr (ppb), (e) P (ppm), and (f) Be (ppb). Proxies display similar patterns during GI 25, which starts at 115.87 (+0.83/−1.71) ka.
Comparison of the YL-1 stalagmite record with regional records
The climatic changes recognized in the YL-1 stalagmite record during the MIS 5e/5d transition do not match with other regional records (Fig. 11). The lack of a similar climate signal between the YL-1 stalagmite and the δ18O record of the So-17A stalagmite from Sofular Cave (Badertscher et al., Reference Badertscher, Fleitmann, Cheng, Edwards, Göktürk, Zumbühl, Leuenberger and Tüysüz2011) is due to Türkiye's northern mountains (i.e., Pontides) acting as a barrier, preventing moisture from reaching the inland regions of Türkiye via the Black Sea. This can be associated with Sofular Cave in northwestern Türkiye mainly receiving moisture from the Black Sea (Fleitmann et al., Reference Fleitmann, Cheng, Badertscher, Edwards, Mudelsee, Göktürk and Fankhauser2009; Badertscher et al., Reference Badertscher, Fleitmann, Cheng, Edwards, Göktürk, Zumbühl, Leuenberger and Tüysüz2011). The percentages of arboreal pollen from Tenaghi–Phillippon in northern Greece (Koutsodendris et al., Reference Koutsodendris, Dakos, Fletcher, Knipping, Kotthoff, Milner and Müller2023) and from Lake Van in eastern Türkiye (Pickarski et al., Reference Pickarski, Kwiecien, Djamali and Litt2015) also do not show a similar trend with YL-1 stalagmite profiles (Fig. 11). However, an increase in arboreal pollen percentage at 115.6 ka has been labeled as the C25 cool event in the Lake Van pollen record. This event corresponds to GS 26 in the Greenland records (NGRIP members, 2004). With uncertainties in age, a correlation can be argued between this cold event in the Lake Van pollen record and the δ18O increase at 116.2 ka in the YL-1 stalagmite (Fig. 11).
Comparison of the YL-1 record with regional records (from north to south). (a) Sofular Cave So-17A stalagmite δ18O record (Badertscher et al., Reference Badertscher, Fleitmann, Cheng, Edwards, Göktürk, Zumbühl, Leuenberger and Tüysüz2011); (b) Tenaghi–Phillippon cores arboreal pollen (%) record (Koutsodendris et al., Reference Koutsodendris, Dakos, Fletcher, Knipping, Kotthoff, Milner and Müller2023); (c, d) Yelini Cave YL-1 δ13C and δ18O stalagmite records; (e) Van Lake arboreal pollen (%) (Pickarski et al., Reference Pickarski, Kwiecien, Djamali and Litt2015); (f) Pentadactylos Cave speleothem record (Nehme et al., Reference Nehme, Kluge, Verheyden, Nader, Charalambidou, Weissbach and Gucel2020); (g) Kanaan Cave speleothem δ18O record (Nehme et al., Reference Nehme, Verheyden, Noble, Farrant, Sahy, Hellstrom, Delannoy and Claeys2015); (h, i) Soreq and Peqiin caves δ18O speleothem records (Bar-Matthews et al., Reference Bar-Matthews, Ayalon, Gilmour, Matthews and Hawkesworth2003; Grant et al., Reference Grant, Rohling, Bar-Matthews, Ayalon, Medina-Elizalde, Bronk Ramsey, Satow and Roberts2012).
There are no clear synchronous climatic signatures between the YL-1 stalagmite and records from the Levant Basin (Fig. 11). The YL-1 stalagmite does not cover a significant part of the LIG period (130–116 ka); indeed, its growth period (117–114 ka) is considered regionally to fall into the glacial inception period (Nehme et al., Reference Nehme, Kluge, Verheyden, Nader, Charalambidou, Weissbach and Gucel2020). However, the δ18O increase at 116.2 ka in YL-1 stalagmite is synchronous with the end of MIS 5e, which is associated with the LIG. Different ages related to the last interglacial period are present in the EM records. The record from Pentadactylos Cave in Cyprus (Nehme et al., Reference Nehme, Kluge, Verheyden, Nader, Charalambidou, Weissbach and Gucel2020) shows that the last interglacial conditions occurred between ca. 127 and ca. 123 ka, with colder and drier climate conditions prevailing during 123–110 ka (Fig. 11). A decrease in growth rate is also recorded at the end of the LIG. The end of moist conditions during the LIG period is noted at 120.3 ka in the Kanaan Cave speleothem record in Lebanon (Nehme et al., Reference Nehme, Verheyden, Noble, Farrant, Sahy, Hellstrom, Delannoy and Claeys2015), while the caves in Israel (Soreq and Peqiin) (Bar-Matthews et al., Reference Bar-Matthews, Ayalon, Gilmour, Matthews and Hawkesworth2003; Grant et al., Reference Grant, Rohling, Bar-Matthews, Ayalon, Medina-Elizalde, Bronk Ramsey, Satow and Roberts2012) recorded the LIG at ca. 120 ka. Since the age interval of the YL-1 stalagmite does not reach back that far, it is not possible to compare these transitions. The isotopic changes in YL-1 during the glacial inception period (ca. 117–114 ka) also do not correspond to EM records. This discrepancy could be due to topographical effects similar to those observed in the Sofular record in the north. The southern mountain ranges (Taurides) also significantly hinder moisture transport to central Anatolia. Therefore, it is expected that the YL-1 stalagmite would show differences from records primarily reflecting eastern Mediterranean moisture. Another possibility is that these EM records may have lower resolution than the YL-1 stalagmite has during the glacial inception (ca. 114–117 ka) (Fig. 11).
Comparison of YL-1 stalagmite record – global relationships
The glacial inception period at Yelini Cave
According to the age modeling, YL-1 stalagmite growth covers the time interval of 117.13 (+0.57/−0.44) to 114.87 (+1.63/−2.89) ka, with slower growth rates after ca. 116.2 ka (Figs. 8 and 12). The onset of stalagmite formation at 117.13 (+0.57/−0.44) ka corresponds to the GS 26 period during which ice sheet growth was recorded in the Greenland ice cores (NGRIP members, 2004). As a result of reduction in insolation, a number of cooling events were reported in the North Atlantic by the end of MIS 5e during the increase of ice volume at 119–116 ka, which was recognized in the Greenland ice core records at the end of the LIG (Oppo et al., Reference Oppo, McManus and Cullen2006; Galaasen et al., Reference Galaasen, Ninnemann, Irvali, Kleiven, Rosenthal, Kissel and Hodell2014; Irvali et al., Reference Irvali, Ninnemann, Kleiven, Galaasen, Morley and Rosenthal2016). An enrichment of δ18O values in the YL-1 stalagmite is observed at 116.83 (+0.58/−0.42) and 116.63 (+0.46/−0.40) ka, which is close to the MIS 5e boundary. A similar enrichment trend is also noticeable for δ13C at 116.65 (+0.51/−0.39) ka (Fig. 8). In addition, increases in Mg/Ca and Sr/Ca and decreases in 87Sr/86Sr ratios as well as decreases in P, Y, Cu, Be, and Zr concentrations are observed at 116.67 (+0.54/−0.38) ka (Fig. 8). Further, a substantial decrease in growth rate and a significant reduction in the diameter of the stalagmite are also recorded in the same period (Fig. 3). Therefore, we argue that a multi-centennial drought event occurred in central Anatolia at 116.65 (+0.51/−0.39) ka, which is chiefly deduced from the YL-1 stalagmite δ13C profile as well as the anomalies in δ18O and Sr-trace element data (Fig. 8).
Comparison of the YL-1 record with global records. (a–c) MD03-2664 core records (Galaasen et al., Reference Galaasen, Ninnemann, Irvali, Kleiven, Rosenthal, Kissel and Hodell2014; Irvali et al., Reference Irvali, Ninnemann, Kleiven, Galaasen, Morley and Rosenthal2016), IRD (Ice-Rafted Debris), (d) 30°N June insolation curve (Berger and Loutre Reference Berger and Loutre1991), (e) NGRIP δ18O record (NGRIP members, 2004), Vienna Standard Mean Ocean Water (VSMOW), (f) MD01-2443 planktonic foraminifera (G. bulloides) δ18O values (de Abreu et al., Reference de Abreu, Abrantes, Shackleton, Tzedakis, McManus, Oppo and Hall2005; Voelker and de Abreu, Reference Voelker, de Abreu, Rashid, Polyak and Mosley-Thompson2011), (g) MD95-2040 summer SST (°C) record (de Abreu et al., Reference de Abreu, Shackleton, Schönfeld, Hall and Chapman2003), (h) Corchia Cave speleothem record (Drysdale et al., Reference Drysdale, Zanchetta, Hellstrom, Fallick and Zhao2005), (i) Campanet Cave CAM-1 stalagmite δ18O record (Dumitru et al., Reference Dumitru, Onac, Polyak, Wynn, Asmerom and Fornós2018), (j) Han-Sur Lesse Cave Han-9 stalagmite δ13C record (Vansteenberge et al., Reference Vansteenberge, Verheyden, Cheng, Edwards, Keppens and Claeys2016), (k) Schneckenloch Cave speleothem record (Boch et al., Reference Boch, Cheng, Spötl, Edwards, Wang and Häuselmann2011), (l) Hölloch Cave speleothem record (Moseley et al., Reference Moseley, Spötl, Cheng, Boch, Min and Edwards2015), (m) Baradla Cave speleothem record (Demény et al., Reference Demény, Kern, Czuppon, Németh, Leél-Ossy, Siklósy and Lin2017), (n, o) Yelini Cave YL-1 stalagmite δ13C and δ18O records, respectively.
The proposed drought at 116.65 (+0.51/−0.39) ka correlates well with the Eirik Drift MD03-2664 marine record in the Labrador Basin (Galaasen et al., Reference Galaasen, Ninnemann, Irvali, Kleiven, Rosenthal, Kissel and Hodell2014; Irvali et al., Reference Irvali, Ninnemann, Kleiven, Galaasen, Morley and Rosenthal2016; Fig. 12). The decrease in North Atlantic Deep Water (NADW) δ13C values of Cibicidoides wuellerstorfi (an epibenthic foraminifera species; Galaasen et al., Reference Galaasen, Ninnemann, Irvali, Kleiven, Rosenthal, Kissel and Hodell2014) at ca. 116.8 ka fits well with the YL-1 isotopic anomalies (Fig. 12d). Similarly, in the same core, Ice Rafted Debris (IRD) and benthic foraminiferal records (e.g., δ18O of Neogloboquadrina pachyderma and Neogloboquadrina pachyderma %) indicate significant cooling (Irvali et al., Reference Irvali, Ninnemann, Kleiven, Galaasen, Morley and Rosenthal2016) (Fig. 12b, c). The abundance of IRD that starts to increase at ca. 117.1 ka reaches 29% at ca. 116.9 ka, and the relative abundance of N. pachyderma attains 61% (Irvali et al., Reference Irvali, Ninnemann, Kleiven, Galaasen, Morley and Rosenthal2016) (Fig. 12). Massive iceberg discharge and cold-water invasion that are reported at ca. 117 ka agree well with isotopic changes in the YL-1 stalagmite (Irvali et al., Reference Irvali, Ninnemann, Kleiven, Galaasen, Morley and Rosenthal2016). Such reductions in NADW formation and subsequent weakening of the Atlantic Meridional Overturning Circulation (AMOC) caused the southward movement of cold polar waters in the Northeast Atlantic region and decreases in winter cyclone activity (Vansteenberge et al., Reference Vansteenberge, Verheyden, Cheng, Edwards, Keppens and Claeys2016; Tzedakis et al., Reference Tzedakis, Drysdale, Margari, Skinner, Menviel, Rhodes and Taschetto2018). Rowe et al. (Reference Rowe, Mason, Andrews, Marca, Thomas, van Calsteren, Jex, Vonhof and Al-Omari2012) discussed the relation between AMOC and winter cyclones that transported precipitation to the Mediterranean basin during the last glacial period. Therefore, it is likely that the YL-1 record reflects a decrease in winter precipitation at 116.65 (+0.51/−0.39) ka (Fig. 12). Similar changes were observed in speleothem records in this time interval (Fig. 12). Accordingly, δ13C values at ca. 117.7–117.3 ka from the Han-9 stalagmite in Han-Sur-Lesse Cave of Belgium that represent the Eemian–Weichselian Transition (EWT) show a 4‰ increase (Vansteenberge et al., Reference Vansteenberge, Verheyden, Cheng, Edwards, Keppens and Claeys2016, 2019), which, within the margin of error, agrees with the YL-1 stalagmite isotope record (Fig. 12h). Further, decreases in summer Sea Surface Temperature (SST) are reported in the North Atlantic ODP 980 sea core between 117 and 118 ka (Oppo et al., Reference Oppo, McManus and Cullen2006) and in the Alboran Sea ODP 976 record at ca. 118 ka (Martrat et al., Reference Martrat, Jimenez-Amat, Zahn and Grimalt2014).
The YL-1 stalagmite record shows significant correlations with Atlantic and European records. The δ18O increase in the BAR-II stalagmite record from Baradla Cave between ca. 118.5 and 117 ka is interpreted as an arid pulse (Demény et al., Reference Demény, Kern, Czuppon, Németh, Leél-Ossy, Siklósy and Lin2017) (Fig. 12m). This event is recognized with simultaneous increases in δ18O and δ13C values in the BAR-II stalagmite. This period can be compared with the high δ18O values at the beginning of the YL-1 growth. In the MD01-2443 marine record obtained from the Iberian margin (Fig. 12f), a small increase in δ18O values of planktonic foraminifer (G. bulloides) at 116.8 ka corresponds to a decrease in SST (de Abreu et al., Reference de Abreu, Abrantes, Shackleton, Tzedakis, McManus, Oppo and Hall2005; Voelker and de Abreu, Reference Voelker, de Abreu, Rashid, Polyak and Mosley-Thompson2011) (Fig. 12g). The HÖL-10 stalagmite from Hölloch Cave (Moseley et al., Reference Moseley, Spötl, Cheng, Boch, Min and Edwards2015) records a small δ18O decrease at 117.2 ka (Fig. 12l). A similar change is represented by the δ18O increase in the CAM-1 stalagmite from Campanet Cave at 117 ka (Dumitru et al., Reference Dumitru, Onac, Polyak, Wynn, Asmerom and Fornós2018) (Fig. 12i).
MIS 5e/5d transition
The MIS 5e/5d transition in the YL-1 stalagmite is represented by a δ18O enrichment of 1.88‰ at 116.24 (+0.53/−0.86) ka (Fig. 12a). During this transition, the δ13C values show a weaker enrichment than the δ18O profile (Fig. 12c). Similarly, there is no clear trend in other isotopic and elemental proxies. A change in growth rate occurs at 116.08 ka, ca. 200 years after the transition (Fig. 8). However, highly porous (or dendritic), irregularly arranged calcite crystals during ca. 116.32–116.16 ka interval suggest irregular (non-steady) water-drip conditions coinciding with the transition (Fig. 4) (Frisia, Reference Frisia2015).
Elevated δ18O values at 116.24 (+0.58/−0.86) ka in YL-1 coincides with the depletion of δ18O at ca. 115.5 ka, representing the end of the last interglacial period in the Greenland ice record (NGRIP members, 2004) (Fig. 12h). It also matches with the upper limit of the benthic “plateau” of MIS 5e at 116.1 ± 0.9 ka (Shackleton et al., Reference Shackleton, Sánchez Goñi, Pailler and Lancelot2003). Neogloboquadrina pachyderma abundance reaching up to 90–95% in the eastern Nordic seas during the same time interval and the decrease in SST at the MIS 5e/5d boundary at ca. 116 ka indicate an interglacial–glacial transition in northwest Europe (Fronval and Jansen, Reference Fronval and Jansen1997; Knudsen et al., Reference Knudsen, Seidenkrantz and Kristensen2002). The pollen record of Lake Monticchio in Italy also indicates a decrease in Mediterranean woody taxa at ca.115.8 ka (Brauer et al., Reference Brauer, Allen, Mingram, Dulski, Wulf and Huntley2007).
The MIS 5e/5d transition is also recognized in the δ18O records of the Campanet Cave CAM-1 stalagmite in Mallorca (Dumitru et al., Reference Dumitru, Onac, Polyak, Wynn, Asmerom and Fornós2018; Fig. 12i). The CAM-1 stalagmite records the MIS 5e/5d transition through increasing δ18O values towards 115.4 ka. The MIS 5e/5d transition in the YL-1 stalagmite record also overlaps with the transition recognized in the Schneckenloch Cave SCH-7 (Boch et al., Reference Boch, Cheng, Spötl, Edwards, Wang and Häuselmann2011) and the Hölloch Cave HÖL-10 stalagmites (Moseley et al., Reference Moseley, Spötl, Cheng, Boch, Min and Edwards2015) in the northern Alps (Fig. 12). SST records from the Iberian margin MD95-2040 core do not show significant changes during the MIS 5e/5d transition and at 116.65 ka (de Abreu et al., Reference de Abreu, Shackleton, Schönfeld, Hall and Chapman2003). In addition, the CC5 stalagmite record from Corchia Cave in Italy (Drysdale et al., Reference Drysdale, Zanchetta, Hellstrom, Fallick and Zhao2005) does not provide a correlation with the YL-1 stalagmite record during this period (Fig. 12).
Large and rapid climatic improvements, in which stadial phases are temporarily interrupted by the last glacial onset (Govin et al., Reference Govin, Capron, Tzedakis, Verheyden, Ghaleb, Hillaire-Marcel and St-Onge2015) immediately after the MIS 5e/5d transition, are defined as Greenland Interstadial (GI) phases (Capron et al., Reference Capron, Landais, Chappellaz, Schilt, Buiron, Dahl-Jensen and Johnsen2010, Reference Capron, Landais, Chappellaz, Buiron, Fischer, Johnsen, Jouzel, Leuenberger, Masson-Delmotte and Stocker2012). The GI 25 transition at ca. 115.4 ka in the NGRIP record (NGRIP members, 2004) correlates with a 115.87 (+0.83/−1.71) ka anomaly in the YL-1 record (Fig. 12h). Correlation of the Schneckenloch Cave SCH-7 and the Hölloch Cave HÖL-10 records to GI 25 at 115.3 ± 0.5 ka and 114.9 ka, respectively, are contemporaneous (within the error margin) with the YL-1 record (Fig. 12).
The YL-1 δ13C profile at 115.87 (+0.83/−1.71) ka also displays a negative shift during interstadial conditions. Well-laminated, fine-grained crystals in the stalagmite after the transition may reflect regular (or near steady) dripping (and rainfall) conditions (Fig. 4). Mosaic-textured calcite, which replaces fine-grained calcite layers at ca. 115.27 ka, are concurrent with an increase in δ18O profile (Figs. 4 and 5). Relatively coarser calcite crystals suggest a decrease in water availability and a decrease in the dripping regime of cave water, which in turn diminishes (or halts) stalagmite growth (Gonzalez et al., Reference Gonzalez, Carpenter and Lohmann1992; Ayalon et al., Reference Ayalon, Bar-Matthews and Kaufman1999).
CONCLUSIONS
This study presents a detailed investigation of the paleohydroclimatic conditions that prevailed in central Anatolia during the time period of 117.13 (+0.57/−0.44) ka to 114.87 (+1.63/−2.89) ka. δ18O ratio changes recorded in the Yelini Cave YL-1 stalagmite are affected by increases or decreases in effective winter precipitation regulated by annual changes in NAO. The δ13C profile of the YL-1 stalagmite reflects time-dependent changes in soil density and host rock (marble) interactions, which are controlled by variations in net effective precipitation (precipitation−evaporation) in this semi-arid region. In the YL-1 isotope record, the MIS 5e/5d boundary is marked at 116.24 (+0.53/−0.86) ka, which coincides (within the age error margin) with the MIS 5e/5d transition records in several global localities. It is also noted that the YL-1 record shows fewer correlations with regional continental records. We argue that a significant multi-centennial drought event occurred near the end of MIS 5e and ca. 400 years before the MIS 5e/5d transition in central Anatolia. This event at 116.65 (+0.51/−0.39) ka (mainly based on swings of the δ13C profile) is likely due to the decrease in winter precipitation in the Anatolian Peninsula and correlates well with various global paleoclimatic archives. The YL-1 stalagmite also recorded the transition to the GI 25 wet phase at 115.87 (+0.83/−1.71) ka. As shown by trace element and stable isotope data, the transition to the GI 25 phase began with a significant increase in rainfall amount. Stalagmite growth ceased, with a gradual transition towards a colder period. The Yelini Cave YL-1 stalagmite record constitutes the first paleoclimate archive revealing the climatic conditions and timing of the MIS 5e/5d transition in central Anatolia, Türkiye.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/qua.2023.56
Acknowledgments
This study is supported by the Scientific and Technological Research Council of Türkiye (TÜBİTAK) under grant no. 118C275 (awarded to Dr. İ. Tonguç Uysal) and European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No: 842403 (awarded to the corresponding author Dr. Ezgi Ünal İmer). The authors thank to Dr. Ali İmer (METU, Türkiye) for his help during the sampling campaign in the cave. The authors are grateful to two anonymous referees and editors for their critical comments that improved the manuscript significantly.
Open access
