INTRODUCTION
With changes in hydrology imminent as a result of anthropogenic climate change, improved knowledge of past changes in precipitation and freshwater systems is necessary to enhance the accuracy of climate models and better understand the effects of a changing climate on human society. In the Levant, or eastern Mediterranean region, large gradients in climate on small spatial scales in conjunction with predominantly seasonal precipitation generate greater uncertainties in predictions of water availability (Rambeau, Reference Rambeau2010). Furthermore, the region's long history of human habitation makes it uniquely suited to studies on the effects of a changing climate on human society (Frumkin et al., Reference Frumkin, Bar-Yosef and Schwarcz2011; Enzel and Bar-Yosef, Reference Enzel and Bar-Yosef2017). However, few existing paleoclimate records can be correlated with human activity at archaeological sites (Rambeau, Reference Rambeau2010; Enzel and Bar-Yosef, Reference Enzel and Bar-Yosef2017). In addition, very few studies of pre-Holocene paleoseasonality have been performed in this region (Ben Dor et al., Reference Ben Dor, Armon, Ahlborn, Morin, Erel, Brauer, Schwab, Tjallingii and Enzel2018, Reference Ben Dor, Marra, Armon, Enzel and Morin2021; Langgut et al., Reference Langgut, Cheddadi and Sharon2021; Müller et al., Reference Müller, Neugebauer, Dor, Enzel and Schwab2022). During the final Pleistocene, archeologically known as the Epipaleolithic (EP; 25–11.7 ka), human communities in the Levant shifted from small bands of nomadic hunter-gatherers into the earliest sedentary settlements (Belfer-Cohen and Goring-Morris, Reference Belfer-Cohen and Goring-Morris2020). An understanding of weather patterns during this time, when human impact on the environment was minimal, may improve present climate models and enhance insight into ecological and cultural changes during the final Pleistocene (Rambeau, Reference Rambeau2010).
Late Pleistocene paleoclimate in the Levant is largely correlated with glacial fluctuations in the Northern Hemisphere, but while some paleoclimate archives in the Levant, including pollen, lake levels, and speleothems, indicate cold and dry conditions during Heinrich events and warm, wet conditions during the Bølling-Allerød, interpretations of records dated to the last glacial maximum (LGM) are contradictory (Rossignol-Strick, Reference Rossignol-Strick1995; Bar-Matthews et al., Reference Bar-Matthews, Ayalon, Kaufman and Wasserburg1999, Reference Bar-Matthews, Ayalon, Gilmour, Matthews and Hawkesworth2003; Bartov et al., Reference Bartov, Stein, Enzel, Agnon and Reches2002; Robinson et al., Reference Robinson, Black, Sellwood and Valdes2006; Roberts et al., Reference Roberts, Jones, Benkaddour, Eastwood, Filippi, Frogley and Lamb2008; Develle et al., Reference Develle, Gasse, Vidal, Williamson, Demory, Van Campo, Ghaleb and Thouveny2011; Orland et al., Reference Orland, Bar-Matthews, Ayalon, Matthews, Kozdon, Ushikubo and Valley2012; Torfstein et al., Reference Torfstein, Goldstein, Stein and Enzel2013; Miebach et al., Reference Miebach, Chen, Schwab, Stein and Litt2017; Matthews et al., Reference Matthews, Affek, Ayalon, Vonhof and Bar-Matthews2021). In the LGM, levels of Lake Lisan, the precursor to the Dead Sea that stretched northward to the modern Sea of Galilee, were at a very high stand. However, the conditions that led to these high lake levels are not well understood (Bartov et al., Reference Bartov, Stein, Enzel, Agnon and Reches2002; Hazan et al., Reference Hazan, Stein, Agnon, Marco, Nadel, Negendank, Schwab and Neev2005; Torfstein et al., Reference Torfstein, Goldstein, Stein and Enzel2013; Lev et al., Reference Lev, Stein, Ito, Fruchter, Ben-Avraham and Almogi-Labin2019). For example, pollen records indicate a prevalence of steppe vegetation in the region, possibly indicating low water availability (Miebach et al., Reference Miebach, Chen, Schwab, Stein and Litt2017, Reference Miebach, Stolzenberger, Wacker, Hense and Litt2019; Langgut et al., Reference Langgut, Cheddadi and Sharon2021). Meanwhile, speleothem isotope records from caves in the Levant indicate increasing aridity, but speleothems grew in the Negev Desert, where conditions today are too arid for their formation (Bar-Matthews et al., Reference Bar-Matthews, Ayalon, Kaufman and Wasserburg1999, Reference Bartov, Goldstein, Stein and Enzel2003; Vaks et al., Reference Vaks, Bar-Matthews, Ayalon, Matthews, Frumkin, Dayan, Halicz, Almogi-Labin and Schilman2006). Enzel et al. (Reference Enzel, Amit, Dayan, Crouvi, Kahana, Ziv and Sharon2008) suggest that storm path trajectories of the Cyprus Lows were located farther to the south during the LGM than currently, accounting for more moisture in the southern Levant. Oxygen isotope records from speleothems in the region support this hypothesis, with lower δ18Ospeleothem values to the south and north compared with central Israel during the LGM, which may indicate low δ18Oprecipitation values due to Rayleigh distillation processes outside the main storm path (Vaks et al., Reference Vaks, Bar-Matthews, Ayalon, Matthews, Frumkin, Dayan, Halicz, Almogi-Labin and Schilman2006; Bar-Matthews et al., Reference Bar-Matthews, Keinan and Ayalon2019). Evidence of flooding events in Lake Lisan during the LGM support more frequent southerly storm paths (Ben Dor et al., Reference Ben Dor, Armon, Ahlborn, Morin, Erel, Brauer, Schwab, Tjallingii and Enzel2018). Pollen records in the region could also support this inference; models based on pollen records suggest a similar or lower amount of precipitation than today combined with a lower temperature and lower evaporation due to decreased insolation (Miebach et al., Reference Miebach, Stolzenberger, Wacker, Hense and Litt2019; Langgut et al., Reference Langgut, Cheddadi and Sharon2021). Climate modeling efforts also suggest that lower temperatures in combination with lower precipitation and evaporation played a role in the LGM climate of the Levant (Stockhecke et al., Reference Stockhecke, Timmermann, Kipfer, Haug, Kwiecien, Friedrich, Menviel, Litt, Pickarski and Anselmetti2016; Ludwig and Hochman, Reference Ludwig and Hochman2022). Another working hypothesis is that expanded snow cover and frozen ground would result in annual melt with less contribution to soil moisture (Develle et al., Reference Develle, Gasse, Vidal, Williamson, Demory, Van Campo, Ghaleb and Thouveny2011). The lack of speleothem growth on Mt. Hermon during the LGM implies a wider area of freezing temperatures (Ayalon et al., Reference Ayalon, Bar-Matthews, Frumkin and Matthews2013), and dry intervals in a lake record from Yammoûneh, Lebanon, indicate an increase in storage of water as snow and/or frozen ground conditions (Develle et al., Reference Develle, Gasse, Vidal, Williamson, Demory, Van Campo, Ghaleb and Thouveny2011), supported by evidence of glaciers in the Mt. Lebanon area during this interval (Moulin et al., Reference Moulin, Benedetti, Vidal, Hage-hassan, Elias, Woerd, Tapponnier, Schimmelpfennig and Da2022). Zaarur et al. (Reference Zaarur, Affek and Stein2016) reconstructed low δ18Owater values from carbonate clumped isotopes in gastropod shells in the upper Jordan River valley, while speleothems in the more southerly Soreq Cave exhibited higher δ18Ospeleothem values during the LGM (Bar-Matthews et al., Reference Bar-Matthews, Ayalon, Gilmour, Matthews and Hawkesworth2003). Zaarur et al. (Reference Zaarur, Affek and Stein2016) interpret their reconstructed LGM δ18Owater values as being influenced by snowmelt from the surrounding mountains.
High-resolution (subseasonal) records such as those produced from sclerochronologically subsampled organismal remains can help to disentangle moisture sources (Taft et al., Reference Taft, Wiechert, Albrecht, Leipe, Tsukamoto, Wilke, Zhang and Riedel2020; Wiese et al., Reference Wiese, Hartmann, Gummersbach, Shemang, Struck and Riedel2020) and provide a seasonal perspective on geologic-scale climate change (Taft et al., Reference Taft, Mischke, Wiechert, Leipe, Rajabov and Riedel2014; Hartman et al., Reference Hartman, Bar-Yosef, Brittingham, Grosman and Munro2016). Shells from organisms such as gastropods and other mollusks are abundant in the fossil record and can be subsampled to produce stable isotope records along the direction of their growth. The stable isotope ratios in shell carbonate (δ18Oshell and δ13Cshell) vary as the animal grows due to changes in temperature and changes in δ18Owater values and δ13C of dissolved inorganic carbon (DIC) in the freshwater environment (Leng and Marshall, Reference Leng and Marshall2004). Snails of the genus Melanopsis occur in freshwater bodies in the Mediterranean region and are common in the Levant, including the Hula Valley (Tchernov, Reference Tchernov1973; Glaubrecht, Reference Glaubrecht1993; Heller et al., Reference Heller, Mordan, Ben-Ami and Sivan2005). Modern snails of this genus occur as far north as Hungary but are presently limited to thermal springs in that area (Glaubrecht, Reference Glaubrecht1993). Life span estimates for this genus vary widely. Elkarmi and Ismail (Reference Elkarmi and Ismail2006) measured and weighed Melanopsis snails at Azraq Oasis, Jordan, and found five size groups, estimating an age of 5 years for the largest snails. Zaarur et al. (Reference Zaarur, Affek and Stein2016) reported sclerochronological δ18Oshell values from modern specimens in the Jordan River valley that reflect two seasonal cycles. Geary et al. (Reference Geary, Hoffmann, Magyar, Freiheit and Padilla2012) suggest a much longer life span for specimens from paleo-Lake Pannon, with one individual surviving to 16 years based on sclerochronological evidence. Growth rates of these snails can differ based on environmental conditions, with warmer waters generally more favorable to growth and shell growth potentially slowing or pausing during the winter months (Geary et al., Reference Geary, Hoffmann, Magyar, Freiheit and Padilla2012; Bartolini et al., Reference Bartolini, Aquiloni, Nisi, Nuccio, Vaselli and Cianfanelli2017).
The objective of our study is to use high-resolution sclerochronological δ18Oshell and δ13Cshell data from Melanopsis shells in a well-dated geoarchaeological sequence in the northern Jordan River valley to provide insight into hydrology, seasonality, and climate in the eastern Mediterranean during the late Pleistocene.
STUDY AREA
Geologic and hydrological setting
The Hula Basin is a small pull-apart basin associated with the Dead Sea Transform (Horowitz, Reference Horowitz1973; Rybakov et al., Reference Rybakov, Fleischer and ten Brink2003; Schattner and Weinberger, Reference Schattner and Weinberger2008; Heimann et al., Reference Heimann, Zilberman, Amit and Frieslander2009). The basin's area is approximately 150 km2, and it serves as an erosion base level for approximately 1500 km2. The valley floor is currently at an elevation of approximately 70 m above sea level (m asl). The basaltic Golan Heights plateau lies to the east, with Mt. Hermon (2814 m asl) marking the northeastern end of the region and the limestone mountains of the Upper Galilee bordering the west. The Kurazim uplifted basaltic block is the southern boundary of the basin. Subsidence has allowed sediments in the Hula Basin to accumulate since at least the middle to late Miocene (Rybakov et al., Reference Rybakov, Fleischer and ten Brink2003).
The Hula Valley is the uppermost basin in the Jordan River catchment, where the Dan, Banias, and Hasbani Springs contribute the majority of the Jordan River headwaters. Precipitation in the Mt. Hermon area recharges these karst springs (Babad et al., Reference Babad, Burg and Adar2020), with a snowmelt contribution of about 30% (Gil'ad and Bonne, Reference Gil'ad and Bonne1990; Sade et al., Reference Sade, Rimmer, Samuels, Salingar, Denisyuk, Alpert, Borchardt, Bogardi and Ibisch2016). Springs in the basaltic Golan Heights to the east of the valley also contribute to the Jordan River and may have played a larger role in the geologic past (Spiro et al., Reference Spiro, Ashkenazi, Starinsky and Katz2011). Attempts to better understand hydrogeology in the Hula Basin are complicated by its complex geology and a lack of wells and boreholes from which data can be collected. However, recent studies suggest surface discharge of deep groundwater in the southern part of the basin (Babad et al., Reference Babad, Burg and Adar2020).
Before drainage in the 1950s, the Hula Valley held a shallow, up to 3-m-deep lake (~14 km2) at its southern end and extensive swamps to the north (Dimentman et al., Reference Dimentman, Bromley and Por1992). Some wetlands were restored in the Agamon Hula Nature Reserve between 1994 and 2010. The reflooded area contains a small, shallow human-made lake fed by canals (Gophen et al., Reference Gophen, Tsipris, Meron and Bar-Ilan2003). Hydrology of the Agamon Hula is seasonal, with rainy months raising the water level and dry months lowering the water level through evaporative loss (Gophen et al., Reference Gophen, Tsipris, Meron and Bar-Ilan2003; Litaor et al., Reference Litaor, Eshel, Sade, Rimmer and Shenker2008).
Climate
The Hula Basin holds the headwaters of the Jordan River, the primary water source for the entire Dead Sea Rift valley. Modern climate is characterized by moderately cool, wet winters, with mean January temperatures of about 12°C, and hot, dry summers, with mean July temperatures of 28°C. Mean annual temperature is 21°C, and annual rainfall ranges from 400 to 800 mm, with the northern end of the basin receiving significantly more precipitation than the southern end (van Zeist et al., Reference Zeist, Baruch, Bottema, Kaptijn and Petit2009; Sade et al., Reference Sade, Rimmer, Samuels, Salingar, Denisyuk, Alpert, Borchardt, Bogardi and Ibisch2016). Higher altitudes experience a cooler, wetter climate, with mean annual temperatures in the mountains and plateaus surrounding the Hula Basin reaching 16°C and a mean annual precipitation up to 1300 mm (Rimmer and Salingar, Reference Rimmer and Salingar2006; van Zeist et al., Reference Zeist, Baruch, Bottema, Kaptijn and Petit2009; Sade et al., Reference Sade, Rimmer, Samuels, Salingar, Denisyuk, Alpert, Borchardt, Bogardi and Ibisch2016).
Archaeological site
The Jordan River Dureijat excavation site (JRD; Fig. 1) lies on the eastern bank of the Jordan River at the confluence of the small Dureijat (meaning “steps” in Arabic) stream at the southern edge of the Hula Basin (Marder et al., Reference Marder, Biton, Boaretto, Feibel, Melamed, Mienis, Rabinovich, Zohar and Sharon2015; Sharon et al., Reference Sharon, Grosman, Allué, Barash, Bar-Yosef Mayer, Biton and Bunin2020). The site was discovered in 1999 during an archaeological survey before a Jordan River drainage operation. Excavation at the site continued with seven excavation seasons from 2014 to 2021 (Sharon et al., Reference Sharon, Grosman, Allué, Barash, Bar-Yosef Mayer, Biton and Bunin2020). Sediments at the site are composed of silt layers indicating a submerged shallow lake environment and layers rich in gastropod and Unio shells representing a shallow nearshore lake environment (Sharon et al., Reference Sharon, Grosman, Allué, Barash, Bar-Yosef Mayer, Biton and Bunin2020; Fig. 2).
Figure 1. (a) Eastern Mediterranean region and location of Jordan River Dureijat archaeological site (JRD; orange star). (b) Study area showing the location of JRD (orange star), Agamon Hula, Gesher Benot Ya'aqov (GBY), and nearby locations referenced in the text. The approximate outline of Lake Lisan during its LGM high stand is based on Torfstein et al. (Reference Torfstein, Goldstein, Stein and Enzel2013).
Figure 2. Stratigraphy and radiocarbon ages at Jordan River Dureijat archaeological site (JRD; adapted from Sharon et al., Reference Sharon, Grosman, Allué, Barash, Bar-Yosef Mayer, Biton and Bunin2020). Yellow stars indicate the approximate depth of shells sampled for this study.
The long and continuous record of archaeological remains at JRD is rare for the EP of this region. The archaeological sequence of JRD comprises almost the entire Levantine EP era (Goring-Morris and Belfer-Cohen, Reference Goring-Morris, Belfer-Cohen, Enzel and Bar-Yosef2017), including an Early EP lithic tradition, the Middle EP Geometric Kebaran, and the three upper horizons documenting the Late EP Natufian cultural entity. In addition, evidence for Pre-Pottery Neolithic A presence was found in the uppermost archaeological layer of the JRD sequence. A rich assemblage of fishing gear exposed in all archaeological horizons documents more than 10,000 years of intermittent visits to the shore of paleo-Lake Hula (Sharon et al., Reference Sharon, Grosman, Allué, Barash, Bar-Yosef Mayer, Biton and Bunin2020). Similar to other sites on the banks of the upper Jordan River (Goren-Inbar et al., Reference Goren-Inbar, Sharon, Melamed and Kislev2002a, Reference Goren-Inbar, Werker and Feibel2002b; Aharonovich et al., Reference Aharonovich, Sharon and Weinstein-Evron2014; Melamed et al., Reference Melamed, Kislev, Geffen, Lev-Yadun and Goren-Inbar2016), waterlogged layers result in exceptional preservation of botanical remains and other paleoenvironmental proxies, which rarely occurs in EP sites in the region (Belfer-Cohen and Goring-Morris, Reference Belfer-Cohen, Goring-Morris, Bar-Yosef and Valla2013). The long and well-dated archaeological and sedimentological sequence of JRD, the exceptional preservation of paleoenvironmental proxies, its location in the core area of the Mediterranean Levantine EP, and its proximity to other key sites documenting the shift to sedentism and Neolithization together with its geographic location at the uppermost basin of the Jordan River valley, all make JRD an ideal location to study changing paleoenvironments and their impact on key changes in human history (Langgut et al., Reference Langgut, Cheddadi and Sharon2021).
MATERIALS AND METHODS
Shells from the JRD Area B east section (the type section of the site; Sharon et al., Reference Sharon, Grosman, Allué, Barash, Bar-Yosef Mayer, Biton and Bunin2020; Fig. 2) were systematically collected during sediment logging in the 2016 and 2017 excavation seasons. The lowermost layer exposed at the east section (Layer 6) is 14C dated as slightly older than 20 cal ka BP. The uppermost layer (Layer 3-0) is assigned to the beginning of the Holocene at 10.9 cal ka BP (Sharon et al., Reference Sharon, Grosman, Allué, Barash, Bar-Yosef Mayer, Biton and Bunin2020). The sampled layers were chosen from among the sedimentary horizons bearing archaeological material across the chronological sequence and based on the presence of large, well-preserved shells (Fig. 2; Table 1).
Table 1. Shell size and archaeological layer ages and IDs reported by Sharon et al. (Reference Sharon, Grosman, Allué, Barash, Bar-Yosef Mayer, Biton and Bunin2020).a
a JRD, Jordan River Dureijat archaeological site; m asl, meters above sea level.
A comparative control sample of live Melanopsis specimens documenting present-day conditions was collected from Agamon Hula in May 2017 and placed in ethanol. The largest shell in the sample was chosen for analysis. The shell was sandblasted to remove the periostracum, then cleaned with deionized water and a brush and needle to remove stains. It was then dried in an oven at 40°C.
The chosen subfossil shells were placed in ethanol in an ultrasonic bath for 30 minutes. After drying at room temperature, the shells were cleaned with deionized water and a brush and needle and then dried in an oven at 40°C. After cleaning, all shells were photographed and whorls were counted. Locations of growth marks and areas of damage, staining, and poor preservation were noted. Species taxonomy was determined after Heller et al. (Reference Heller, Mordan, Ben-Ami and Sivan2005).
Samples were milled from the shell surface using a modified dental drill placed in a custom holder after Wiese et al. (Reference Wiese, Hartmann, Gummersbach, Shemang, Struck and Riedel2020). The apparatus was placed under a microscope to enable precise milling. A maximum milling depth of 50 μm was set using a gauge on the holder to ensure that only the outermost shell was sampled. Before each sample was milled, 5% HCl was applied to the drill bit and holder using a synthetic swab to remove carbonate contaminants, then deionized water was applied in the same manner to neutralize the surface. Finally, a pressurized air hose was used to dry the apparatus before use. Samples were milled using a 0.5-mm-round diamond bit, except for samples taken from the smallest whorls of Layer 3C Shell 1 and Layer 4 Shell 2. The heights of these whorls were smaller than 0.5 mm and required a tapered diamond drill bit. Shells were milled to collect 30–80 μg from the uppermost (oldest) whorls of all shells and 100–200 μg on larger whorls.
Each shell was individually assessed to determine appropriate sample spacing to collect 50 samples per shell. This number allows for a submonthly sampling resolution for an assumed life span of 2 years. Degraded or stained areas were avoided to the extent possible. Due to size constraints, 45 samples from Layer 5 Shell 2 were collected. Because of shell damage during sampling, 48 samples were collected from Layer 3C Shell 2.
To determine the preservation state of the shells, the microstructures of a shell from Layer 5 and another shell from Layer 3C were examined with a Hitachi TM3000 desktop scanning electron microscope (SEM) at University of Iceland. Shell subsamples of ca. 2 mm2 size were broken off the shells, coated with gold, and visually inspected. Examination with the SEM revealed the presence of crossed-lamellar layers (Iljina and Frolov, Reference Iljina and Frolov2010; Supplementary Fig. 1). Signs of microbial attack or diagenetic shell alteration were not observed; thus, the material was considered to be reliable.
Oxygen and carbon isotope ratios were determined from 30–80 μg aliquots of sample material using a Thermo Scientific KIEL IV Carbonate Device coupled to a MAT 253 isotope ratio mass spectrometer at the GFZ Potsdam, Germany. The system, reference gas, and internal standard C1 were routinely calibrated against the Vienna PeeDee Belemnite (VPDB) standard using the international reference materials NBS19 (δ13C = 1.95‰, δ18O = −2.2‰) and NBS18 (δ13C = −5.014‰, δ18O = −23.2‰). The measured sample data were adjusted using the calibrated internal laboratory standard C1 (δ13C = 2.40‰, δ18O = −1.31‰), which was run six times between 40 samples. Repeated measurements of NBS19 and C1 standard material show a precision better than ±0.08‰ and for NBS18 better than ±0.1‰ for both δ13C and δ18O. To correct for acid fractionation of aragonite, −0.38‰ offset was applied after Kim et al. (Reference Kim, Mucci and Taylor2007a). Isotope ratios of shells are reported in the delta notation relative to VPDB, and the cited and calculated oxygen isotope values of water (δ18Owater) are expressed relative to Vienna Standard Mean Ocean Water (VSMOW).
Due to low intensities in the smallest samples of the upper whorls, δ13C and/or δ18O values for 39 samples could not be determined. An additional seven samples were milled and analyzed where possible. Two split samples were analyzed, Modern Shell (sample 49) and Layer 4 Shell 2 (sample 16). Based on the results, an overall uncertainty of 0.2‰ is assumed for this study. The δ18O and δ13C values are reported to 0.1‰ precision. All shells were photographed after milling to document sample locations (Rice, Reference Rice2018).
Interpretation framework for δ18O and δ13C profiles
As gastropods precipitate their shells, the δ18Oshell values primarily reflect those of the lake water and an offset due to temperature effects (Schmitz and Andreasson, Reference Schmitz and Andreasson2001; Spiro et al., Reference Spiro, Ashkenazi, Mienis, Melamed, Feibel, Delgado and Starinsky2009). Precipitation, evaporation, and inflowing water can alter δ18Owater in the snail's habitat. Evaporation leads to higher δ18Owater values in the lake, while the input of freshwater in the form of precipitation or inflow typically leads to lower values (Leng and Marshall, Reference Leng and Marshall2004). Meltwater from snow and rainwater that falls at high altitudes to recharge springs can contribute to low δ18Owater values in the Jordan River (Gat and Dansgaard, Reference Gat and Dansgaard1972; Zaarur et al., Reference Zaarur, Affek and Stein2016). In addition, some storms can bring high-δ18O rainwater. Sub-cloud evaporation of raindrops and convective rain that originates from evaporated water over land will result in high δ18Owater values in rain (Gat and Dansgaard, Reference Gat and Dansgaard1972). Offsets in δ18Oshell values due to temperature are often smaller than the changes in evaporation and precipitation in a lacustrine environment. Specifically, increases in temperature lead to lower δ18Oshell values but to increased evaporation, which results in higher δ18Owater values. Because δ18Oshell values change by only 0.24‰/°C (Kim et al., Reference Kim, O'Neil, Hillaire-Marcel and Mucci2007b), the changes in δ18Owater values due to evaporation are often more important. Similarly, colder temperatures lead to higher δ18Oshell values but coincide with the rainy season, which leads to lower δ18Owater values in the lake. It is possible that subfossil shells grew during intervals of lower evaporation, in which case temperature effects on δ18Oshell values might be more important. In this case, decreasing δ13Cshell values are still expected with lowering temperatures due to decay of organic matter, and we rely on δ13Cshell values to indicate possible pauses in growth during winter.
Meanwhile, δ13Cshell values will reflect the δ13CDIC values in the water. In lake water, δ13CDIC values are impacted by the δ13C values of inflowing DIC, exchange with the atmosphere, and photosynthesis and respiration of aquatic plants (Leng and Marshall, Reference Leng and Marshall2004). Aquatic productivity is generally the most important factor on short (subannual) timescales and leads to increased values of δ13CDIC in the lake, whereas the decay of organic matter leads to lowered values of δ13CDIC (Leng and Marshall, Reference Leng and Marshall2004). Storms wash soil into the lake, bringing low-δ13C carbon into the lake system.
Growth marks can also aid in interpreting the δ18Oshell and δ13Cshell profiles. Growth marks indicate an interval when the gastropod pauses its shell growth, which may occur for several reasons. Changes in the environment, including changes in temperature or turbidity, can induce growth marks, but reproduction can also result in a pause in shell growth (Schöne et al., Reference Schöne, Rodland, Wehrmann, Heidel, Oschmann, Zhang, Fiebig and Beck2007; Mannino et al., Reference Mannino, Thomas, Leng and Sloane2008; Geary et al., Reference Geary, Hoffmann, Magyar, Freiheit and Padilla2012). Some studies have noted the formation of these features during times of thermal stress, both for cold and hot temperatures (Schöne et al., Reference Schöne, Rodland, Wehrmann, Heidel, Oschmann, Zhang, Fiebig and Beck2007; Mannino et al., Reference Mannino, Thomas, Leng and Sloane2008). Gastropods may respond to thermal stress by thickening their shell rather than lengthening it in order to improve insulation (Mannino et al., Reference Mannino, Thomas, Leng and Sloane2008). The duration of time represented by growth marks varies; while one growth mark could represent an entire wet season, another might represent a pause in growth for a much shorter interval.
When a growth mark coincides with a decrease in δ13Cshell values, this event is interpreted as a pause in growth during winter months, with changing δ13CDIC values due to decay of aquatic organic matter during winter. However, storms can also cause this same pattern, where the water may become turbid, inducing a pause in growth when low-δ13C soil is washed into the lake. For growth marks preceded by several samples with elevated δ13Cshell values or a long trend toward higher δ13Cshell values, the growth mark likely represents winter. Growth marks that coincide with a drop in δ13Cshell values without any long trend beforehand could be attributed to either a storm or to winter. For shells where this occurs early in life, it is particularly difficult to discern due to the small time span represented by each sample. In addition, samples collected from the last part of the shell to precipitate may represent longer intervals of time due to slower growth rates later in the snail's life. This area can also contain several growth marks due to reproduction (Geary et al., Reference Geary, Hoffmann, Magyar, Freiheit and Padilla2012). So, the samples collected from this part of the shell are likely to represent averages of several weeks or even months, and growth marks are less likely to represent winter than in earlier parts of the shell.
After the growth marks likely to represent winter are identified, the δ18Oshell values are interpreted as changes in δ18Owater and temperature. Decreases in δ18Oshell values result from the input of low-δ18O water into the lake, either through direct precipitation, inflow, runoff, or snowmelt. Temperature increases can also result in decreases in δ18Oshell values; however, temperature is likely to be a relatively minor effect, with a temperature increase of 4°C resulting in a decrease in δ18Oshell of about 1‰ (Kim et al., Reference Kim, O'Neil, Hillaire-Marcel and Mucci2007b). Trends toward higher δ18Oshell values generally indicate 18O enrichment of lake water due to evaporation in summer months. High-δ18O rainfall can also raise δ18Olake values. In the modern climate, this occurs due to sub-cloud evaporation of raindrops or due to convective precipitation, where the cloud's moisture source originated on land. In addition, hydrological conditions during the LGM may have led to a greater proportion of snow recharge to the headwater springs, leading to higher δ18Owater values of winter rainfall within the valley compared with summer lake water fed by the springs, which are recharged by low-δ18Owater snowmelt. Elevated values of δ18Oshell immediately following the winter growth mark may indicate high values of δ18Owater during winter compared with summer, likely indicating the effect of snowmelt on δ18O of river recharge.
Although two shells were subsampled per sediment interval, it is unlikely that the shells grew in the same years. Each shell likely represents distinct intervals of time, and due to interannual variability in temperature and precipitation, the conditions of growth differ from one specimen to another. Nevertheless, two shells from the same layer, despite differences in timing of, for example, storms, should display similar, likely typical conditions of the climate in which they grew. Although climate is defined based on multi-decadal averages, it is unlikely that two shells from a sedimentary interval spanning about 23–70 years would both record the same unusual conditions occurring in consecutive years of growth. The likelihood of, for example, a 25 year flood event occurring at least once during a 2 year interval of shell growth is 8%. The likelihood for two independent 25 year events to occur in two sequential years is only 0.16%.
Modeling δ18Oshell patterns
The ShellChron model package (de Winter, Reference Winter2022) includes a script to generate a modeled δ18Oshell record from sinusoidally varying temperature and δ18Owater curves, which was used to validate the data generated from the modern shell assessed in this study and from modern shells in nearby locations analyzed by Zaarur et al. (Reference Zaarur, Affek and Stein2016). Temperature and δ18Owater values reported in the literature were used to generate this record (Gat, Reference Gat1970; Gat and Dansgaard, Reference Gat and Dansgaard1972; Gil'ad and Bonne, Reference Gil'ad and Bonne1990; Litaor et al., Reference Litaor, Eshel, Sade, Rimmer and Shenker2008; Gophen et al., Reference Gophen, Meron, Tsipris, Orlov-Levin and Peres2016; Zaarur et al., Reference Zaarur, Affek and Stein2016; Supplementary Text S1). To mimic weather effects, model temperature and δ18Owater values vary randomly within a normal distribution from the seasonal curve, calculated with standard deviations of 1.5°C and 0.6‰ from the seasonal model temperature and δ18Owater values, respectively. The code was adapted to use the Kim et al. (Reference Kim, O'Neil, Hillaire-Marcel and Mucci2007b) temperature-dependent oxygen isotope relationship between aragonite and the water from which it precipitates:
where α = (δ18Oshell + 1000)/(δ18Owater + 1000). Solving for δ18Oshell:
Values from Eq. 2 yield δ18Oshell values in the VSMOW reference frame, which are then converted to VPDB for comparison with measured δ18Oshell data.
To mimic pauses in growth during winter months, a temperature cutoff value of 14.5°C was implemented under which the shell is assumed not to grow. This cutoff value was determined using the modern shell models and fit to better reflect patterns in the modern shells; it is assumed to be the same for all shells.
The model shell module was also used to assess the plausibility of different temperature and δ18Owater regimes in fossil shells. The temperature and δ18Owater inputs were fit to shell results. Rather than attempting to model each variation in each shell, the overall variability of δ18Oshell data was fit for each paleoclimate interval.
RESULTS
Fossil shells from JRD are identified as Melanopsis buccinoidea (no ribs) or Melanopsis costata (ribbed specimens; Heller et al., Reference Heller, Mordan, Ben-Ami and Sivan2005). However, whether these morphological differences constitute a separate species in Levantine Melanopsis is under dispute (Falniowski et al., Reference Falniowski, Heller, Cameron, Pokryszko, Osikowski, Rysiewska and Hofman2020). The large shells selected for analysis are consistent with observed sizes of adult specimens, ranging in length from 1.43 to 2.56 cm and in diameter from 0.66 to 1.04 cm, and having between 6.00 and 7.75 whorls (Table 1). All shells contained multiple growth marks, some of which are faint and others more readily visible. The duration of interpreted growth cessation based on sclerochronological evidence was not apparent based on the visual appearance of these marks.
Sclerochronological δ18Oshell and δ13Cshell results are provided in Supplementary Table 1. Means of sclerochronological δ18Oshell and δ13Cshell values for each shell range from −7.9 to −6.4‰ and −7.7 to −4.2‰, respectively. The smallest range of δ18Oshell values, 2.3‰, was calculated for Layer 5 Shell 1, and the largest range, 6.0‰, for the modern shell (Table 2). The δ13Cshell values fluctuate by 2.2‰ at their minimum (Layer 5 Shell 2) and by 9.1‰ at their maximum (modern shell). Detailed shell interpretations based on the framework outlined earlier are given in Supplementary Text S2). These are briefly described in the following section and shown in Figure 3.
Figure 3. Sclerochronological results of δ18Oshell (blue circles) and δ13Cshell (red triangles) for the modern shell and subfossil shells annotated with environmental inferences. Growth marks (gray dashed lines) between samples are labeled with letters. Growth direction is from left to right.
Table 2. Summary statistics of δ18O and δ13C values (‰).
Shell sclerochronology interpretations
The modern shell from Agamon Hula has large changes in δ18Oshell values consistent with a large evaporation effect (Fig. 3a). Winter growth marks (G and H) are identified as the largest decreases in δ18Oshell and δ13Cshell values. Two annual cycles are inferred, each of which shows a trend toward high δ18Oshell values with a bifurcated peak that may have come about due to artificial water level adjustment. One inferred storm event occurs in the first year of growth, which may correspond to an April 2015 storm (Israel Meteorological Service, 2018; Supplementary Fig. 2). Although another storm occurred in April 2016, the impact of this storm on Agamon Hula δ18Owater values is not observed.
The shells from Layer 3C have large and concurrent decreases in δ18Oshell and δ13Cshell values that indicate winter (Shell 1 growth marks B and C; Shell 2 growth marks A and D), but show more gradual increases in δ18Oshell values compared with δ13Cshell values (Fig. 3b and c). Layer 3C Shell 1 records two clear annual cycles and has two additional large decreases in δ13Cshell values (growth marks E and I) that could indicate either winter season or storm events. Layer 3C Shell 2 contains a large damaged portion where the surface layer of the shell flaked off upon contact and which was not sampled for stable isotope analysis. In the remainder of the shell, two annual cycles are inferred (winter growth marks A and D), with one major event that lowered δ13Cshell values and increased δ18Oshell values in the second annual cycle (growth mark B). This event may have been another annual cycle and winter growth mark; however, the small area of the shell that these samples represent suggests that this was a storm event.
Shells from Layer 4 exhibit large, rapid changes in δ18Oshell and δ13Cshell values (Fig. 3d and e). The frequency and timing of changes in δ18Oshell and δ13Cshell values and the lack of a clear relationship between these changes make it difficult to reliably discern annual cycles. The complexity in the patterns could be the result of more frequent storms, interannual variability in the source or amount of rain, or variability in δ13CDIC values due to processes other than aquatic productivity and decay, such as input of soil carbon from wind erosion and storm runoff. Growth marks concurrent with decreases in δ13Cshell values occur six times in Layer 4 Shell 1 (growth marks A, B, E, G, I, and J) and four times in Shell 2 (growth marks B, D, E, and F). These are interpreted as winter growth marks. Trends toward lower δ18Oshell values occur after these marks, which could be due to either increased temperature or lower δ18Owater values during summer.
In Layers 5 and 6, minima in δ13Cshell values generally coincide with peaks in δ18Oshell values (Fig. 3f–i). Winter growth marks identified by decreases in δ13Cshell values occur three times in Layer 5 Shells 1 (growth marks A, E, and F) and 2 (growth marks A, C, and D), twice in Layer 6 Shell 1 (growth marks B and C), and five times in Layer 6 Shell 2 (growth marks B, C, D, E, and G). Changes in δ18Oshell values are generally gradual, with some exceptions that may indicate storms: Layer 5 Shell 1 before growth mark D; Layer 6 Shell 1 before growth mark E; and Layer 6 Shell 2 before growth marks B and G.
Shell isotope models
Models of δ18Oshell values for the modern Jordan River and Sea of Galilee reflect the general range of measured δ18Oshell values. However, the modeled Dan Spring shell has a larger deviation from measured values (−8.3 to −6.1‰ modeled vs. ~ −7.2 to −6.3‰ measured; Zaarur et al., Reference Zaarur, Affek and Stein2016), and the modeled Agamon Hula shell also has a different range than measured values (−7.1 to −2.9‰ modeled vs. ~ −8.4 to −3.4‰ measured). It is possible that the modeled temperature or δ18Owater used in these scenarios does not appropriately describe the conditions under which the shells grew. Given these discrepancies, the δ18Owater in fossil scenarios is estimated to vary within ±1‰ of modeled results. When parts of the model shell that formed below 14.5°C are removed, model results closely resemble patterns in measured δ18Oshell values (Fig. 4).
Figure 4. Model results of δ18Oshell for modern (top) and paleo (bottom) scenarios.
Water temperature and δ18Owater parameters that best fit fossil shell patterns (Table 3) suggest that δ18Owater values were seasonally highest at the end of the evaporative summer season in Layer 3C shells, but shells formed in Layers 4–6 exhibit patterns that are better explained by higher δ18Owater values during cold conditions. Furthermore, seasonal amplitudes of modeled δ18Owater values are higher for the Bølling-Allerød (2.6‰) compared with Heinrich 1 and the LGM (1.8‰ and 0.8‰, respectively). Resulting modeled δ18Oshell patterns closely resemble typical patterns in the shells after removing parts of the model shell that formed below 14.5°C (Fig. 4).
Table 3. Model temperature, δ18Owater, and δ18Oshell parameters and comparison with measured δ18Oshell values.
DISCUSSION
Modern shells
Oxygen isotope variations within modern Melanopsis shells from Agamon Hula (this study), Dan Spring, the Jordan River, and the Sea of Galilee (Zaarur et al., Reference Zaarur, Affek and Stein2016) exhibit seasonal fluctuations in magnitudes consistent with their environments. Of these, Dan Spring is the least seasonally varying environment in terms of water temperature (1°C) and δ18Owater values (1‰) and exhibits the narrowest range of seasonal δ18Oshell values (~1.2‰). Variations in δ18Oshell values from the Sea of Galilee, Jordan River, and Agamon Hula shells are increasingly variable due to higher variability in their environments. The general agreement between modeled shells and measured δ18Oshell values also suggests that the oxygen isotope values in the shells faithfully record changes in the environment, as has been found for other freshwater gastropods (Schmitz and Andreasson, Reference Schmitz and Andreasson2001; Spiro et al., Reference Spiro, Ashkenazi, Mienis, Melamed, Feibel, Delgado and Starinsky2009). A strict temperature cutoff of 14.5°C imperfectly describes Melanopsis shell growth. The patterns of modeled δ18Oshell values above this cutoff seem to replicate conditions in Agamon Hula well; however, this cutoff is only rarely reached in the other modern environments. The shells assessed by Zaarur et al. (Reference Zaarur, Affek and Stein2016) exhibit sawtooth-shaped δ18Oshell patterns, suggesting seasonal growth patterns with slow or paused growth in winter. Because temperatures below 14.5°C are rarely reached in these habitats, a strict temperature cutoff might not reflect actual growth conditions. Instead, these gastropods might slow growth during winter in relatively warm waters or might stop growth due to a nonthermal factor, such as storm-induced turbidity at the beginning of the wet season.
Although few Melanopsis shells from modern environments have been examined, model δ18Oshell results suggest that the general variability within the gastropod's environment can be understood through the extent of seasonal variations in δ18Owater values and, to a lesser extent, seasonal temperature variability. Due to model assumptions, short-term changes in δ18Owater values cannot be replicated, and seasonal changes that are not purely sinusoidal have not been assessed. Due to the lack of subseasonal δ18Owater measurements of modern water bodies, it remains open whether the oxygen isotope compositions follow a sinusoidal curve. Because of these limitations in robustly reconstructing growth rates, inferences of the paleoenvironment are focused on general seasonal patterns and the frequency of inferred storms.
Paleohydrological and paleoclimate inferences
With two to five annual growth cycles, an individual specimen provides a subsample of climate that may not be fully representative of the long-term mean, range, and extremes. As a result, interpretation of the data as typical for climate conditions during the interval of shell growth is tentative. However, assessing two shells per interval and comparing trends makes it relatively unlikely that only unusual years are represented by studied materials. Within each layer, the stable isotope patterns generally agree between and within shells, suggesting that the shells and the cycles within them represent conditions typical of the stratigraphic interval. Furthermore, with the exception of Layer 4 shells, annual cycles within a shell usually vary between the same high and low δ18Oshell values, showing seasonality that repeats itself in adjacent years. Higher variability in Layer 4 shells could point to greater interannual variability within the waterbody during the Heinrich 1 event. This would be in agreement with pollen-based biome modeling that suggests more frequent summer rainfall (Langgut et al., Reference Langgut, Cheddadi and Sharon2021), which possibly led to stable isotope patterns with less regularity. In the present climate, small, localized summer storms occur in northern Israel without disrupting the regional synoptic-scale forcing that causes dry summers (Saaroni and Ziv, Reference Saaroni and Ziv2000). Storms can lead to highly variable δ18Owater values in the modern-day Jordan River (Gat and Dansgaard, Reference Gat and Dansgaard1972), with Hula Valley δ18Oprecipitation values reaching as low as −11‰ (Rindsberger et al., Reference Rindsberger, Jaffe, Rahamim and Gat1990), and similar storms may have had analogous effects in paleo-Lake Hula. The longer inferred life span (about 5.5 and 4.5 years) of these specimens increases the likelihood that these specimens experienced more atypical macroweather conditions. Storm events are also inferred for both Bølling-Allerød (Layer 3C) shells, whereas not all LGM shells (Layers 5 and 6) exhibit these features, which could be due to fewer storms during spring and autumn in this time interval or a shorter growing season. Pollen-based modeling suggests greater spring and autumn precipitation during the Bølling-Allerød (Langgut et al., Reference Langgut, Cheddadi and Sharon2021), in agreement with shell-based inferences.
Variations in whole-shell δ18O values can provide further insight into environmental changes. Over the studied time interval, the average δ18Oshell values exhibit an opposite trend to δ18Ospeleothem records from Soreq Cave and Peqiin Cave, in line with previous findings (Zaarur et al. Reference Zaarur, Affek and Stein2016; Fig. 5a). However, while Zaarur et al. (Reference Zaarur, Affek and Stein2016) attribute low Melanopsis δ18Oshell values during the LGM to increased snowfall, the summer growth bias observed in this study suggests that whole-shell averages may more closely reflect changes in summertime evaporation than changes in δ18O values of winter precipitation. In the general trend of increasing δ18Oshell values from the LGM to the Holocene, this seems to fit expectations of greater summer evaporation due to increasing temperatures, which is thought to have controlled the water balance of Lake Lisan during this time interval (Stockhecke et al., Reference Stockhecke, Timmermann, Kipfer, Haug, Kwiecien, Friedrich, Menviel, Litt, Pickarski and Anselmetti2016; Miebach et al., Reference Miebach, Stolzenberger, Wacker, Hense and Litt2019; Langgut et al., Reference Langgut, Cheddadi and Sharon2021; Ludwig and Hochman, Reference Ludwig and Hochman2022). However, the shells from Heinrich 1 (Layer 4), during which lake levels dropped suddenly, have lower δ18Oshell values than those from the LGM. Evidence suggests that this event resulted in cold and dry conditions, in which low precipitation likely caused the lake-level drop (Bar-Matthews et al., Reference Bar-Matthews, Ayalon, Gilmour, Matthews and Hawkesworth2003; Bartov et al., Reference Bartov, Goldstein, Stein and Enzel2003; Langgut et al., Reference Langgut, Cheddadi and Sharon2021). The cold conditions likely also resulted in low summer evaporation, though this has not yet been studied. Nevertheless, the overall pattern of increasing δ18Oshell values indicating greater summer evaporation aligns with sclerochronological evidence from JRD, with highest δ18Oshell values recorded by the shells from the Bølling-Allerød (Layer 3C).
Figure 5. (a) A comparison of Melanopsis δ18Oshell from this study (blue squares) and Zaarur et al. (2016; pink squares) to δ18Ospeleothem records from Soreq (gray), Peqiin (dark green), and Zalmon (light green) Caves (Bar-Matthews et al., Reference Bar-Matthews, Ayalon, Gilmour, Matthews and Hawkesworth2003; Keinan et al., Reference Keinan, Bar-Matthews, Ayalon, Zilberman, Agnon and Frumkin2019) and (b) standard deviation of δ18Oshell values in Melanopsis (blue squares; note the inverted axis) compared with the Lake Lisan lake level curve (Torfstein et al., Reference Torfstein, Goldstein, Stein and Enzel2013). GBY, Gesher Benot Ya'aqov; JRD, Jordan River Dureijat archaeological site.
The results of sclerochronological δ18Oshell analysis suggest that fluctuations in lake level of paleo-Lake Hula mirror those in downstream Lake Lisan, providing evidence that much of the water feeding Lake Lisan flowed from the northern end of the catchment. During the LGM, variations in δ18Oshell values are small, indicating a large, well-buffered waterbody at a time when Lake Lisan stood at a high stand (Bartov et al., Reference Bartov, Stein, Enzel, Agnon and Reches2002; Torfstein et al., Reference Torfstein, Goldstein, Stein and Enzel2013). After the LGM, Lake Lisan was lower, and lake levels dropped dramatically during Heinrich 1 (Bartov et al., Reference Bartov, Goldstein, Stein and Enzel2003), when low water levels are also inferred at paleo-Lake Hula from highly fluctuating δ18Oshell values. Lake Lisan again rose, reaching high levels at about 15 ka, though not as high as LGM levels (Torfstein et al., Reference Torfstein, Goldstein, Stein and Enzel2013). At this time, the waterbody size at JRD was intermediate: neither as well-buffered as it was during the LGM nor as poorly buffered as during Heinrich 1. The qualitative agreement in lake size changes between paleo-Lake Hula at JRD and Lake Lisan provides further evidence that fluctuations in flow from the Hula catchment contributed to changes in lake level at Lake Lisan. Indeed, the standard deviation of δ18Oshell values seems to vary with the Lake Lisan curve from Torfstein et al. (2013; Fig. 5b). When the standard deviation is low, δ18Owater values in the lake exhibit little fluctuation and Lake Lisan levels are high, while greater standard deviations, reflecting large fluctuations in δ18Owater values, coincide with lower levels of Lake Lisan. This is likely due to the buffering capacity of a larger lake, in which storms have lower impacts on δ18Owater values due to mixing a small amount of water from a storm with a large amount of lake water. However, this relationship may not always hold true. Modeled δ18Oshell values show that the timing of maxima in δ18Owater values in relation to maxima in water temperature can impact the magnitude of variability in δ18Oshell values. When summer evaporation leads to the highest δ18Owater values in the lake, the seasonal extent of δ18Oshell values is smaller than when the highest δ18Owater values occur in winter. This is due to the temperature fractionation effect, in which increased temperatures lead to lower δ18Oshell values at a given δ18Owater value. This effect diminishes the peak δ18Oshell values when high temperatures coincide with high δ18Owater values. So, for environments in which evaporation controls δ18Owater values and the highest δ18Owater values occur during the hottest parts of the year, the seasonal variability in δ18Oshell values will be smaller than for environments that experience a different seasonal regime. This complicates a direct comparison of variability in δ18Oshell values from Layer 3C and other fossil layers.
Sclerochronological data of the LGM shells possibly support the inference of a greater proportion of snowfall in the nearby mountains (Develle et al., Reference Develle, Gasse, Vidal, Williamson, Demory, Van Campo, Ghaleb and Thouveny2011; Orland et al., Reference Orland, Bar-Matthews, Ayalon, Matthews, Kozdon, Ushikubo and Valley2012; Ayalon et al., Reference Ayalon, Bar-Matthews, Frumkin and Matthews2013; Zaarur et al., Reference Zaarur, Affek and Stein2016; Moulin et al., Reference Moulin, Benedetti, Vidal, Hage-hassan, Elias, Woerd, Tapponnier, Schimmelpfennig and Da2022). At present, snow accumulates near the summit of Mt. Hermon during the winter season, and melt contributes about 30% of the water in karstic springs that feed the Jordan River (Gil'ad and Bonne, Reference Gil'ad and Bonne1990; Sade et al., Reference Sade, Rimmer, Samuels, Salingar, Denisyuk, Alpert, Borchardt, Bogardi and Ibisch2016). A consistent, strong negative δ18Oshell excursion that could indicate a sudden influx of snowmelt to paleo-Lake Hula was not observed in the examined shells. Only one shell of four dated to the LGM exhibits a possible snowmelt signal (Layer 6 Shell 2 sample 15), which is not repeated in the next interpreted annual cycle. A similar anomalous event appears in the modern shell from Agamon Hula (sample 17), for which snow is not a likely cause. The sampling resolution was possibly not high enough to detect such a signal in other shells, or snowmelt may have occurred during the snail's winter growth hiatus. Another possibility is that snowmelt quickly entered the karstic aquifer near Mt. Hermon, similar to conditions today (Sade et al., Reference Sade, Rimmer, Samuels, Salingar, Denisyuk, Alpert, Borchardt, Bogardi and Ibisch2016). This would contribute water with lower δ18O values to the main springs feeding the lake, resulting in lower δ18Owater values in summer due to recharge from snow-fed springs and higher δ18Owater values during winter, when precipitation falls as rain in the Hula Valley. The patterns in δ18Oshell values from Layers 4–6 closely resemble modeled δ18Oshell patterns in which δ18Owater values reach a maximum during colder temperatures, supporting this interpretation.
Estimated LGM snowfall
Using a linear mixing model, the proportion of snowmelt contribution to the springs that fed paleo-Lake Hula during the LGM can be estimated. During the LGM, δ18Owater values in precipitation were approximately 2.5‰ higher than present (Bar-Matthews et al., Reference Bar-Matthews, Ayalon, Gilmour, Matthews and Hawkesworth2003), and modern monthly δ18Owater values in Mt. Hermon precipitation range from approximately −10.5 to −5.0‰ (Ayalon et al., Reference Ayalon, Bar-Matthews, Frumkin and Matthews2013) and serve as the modern endmember δ18Owater values for snow and rain recharge to the springs, respectively. Thus, the equivalent δ18Owater values of precipitation during the LGM are considered to be −8.0‰ and −2.5‰ for snow and rain, respectively. For a lake with a short residence time of less than 6 months, the highest and lowest values of the modeled δ18Owater values are expected to reflect the endmembers of local precipitation (winter) and spring water (summer), respectively. Based on the model LGM shell, the summer δ18Owater endmember (−7.1‰) is composed of 84% snowmelt. This is much greater than the modern values of 30% (Gil'ad and Bonne, Reference Gil'ad and Bonne1990) and the estimate of at least 40% contribution of snow during the LGM by Zaarur et al. (Reference Zaarur, Affek and Stein2016).
However, these calculations are subject to considerable uncertainty. The model of δ18Owater relies on the modeled temperature during shell formation. Changing the temperature would result in a change in modeled δ18Owater of ~0.24‰/°C. For lower lake temperatures, this would lead to a lower estimate of snow recharge of ~5%/°C. Clumped isotope-based temperatures from nearby LGM Melanopsis specimens fall in a range similar to modern shells or up to 3°C colder (Zaarur et al., Reference Zaarur, Affek and Stein2016). This is in line with our model, which estimates an average lake temperature of ~16°C during the LGM, about 2°C colder than the modern Jordan River. However, clumped isotope evidence from Soreq Cave suggests an ~10−12°C temperature change in the same time interval (Affek et al., Reference Affek, Bar-Matthews, Ayalon, Matthews and Eiler2008; Matthews et al., Reference Matthews, Affek, Ayalon, Vonhof and Bar-Matthews2021). Although the snail-shell temperatures more closely reflect summer-growth temperatures, lake temperatures 10°C colder than modern are unlikely, as this would be colder than modern Melanopsis habitats. The discrepancy could be due in part to temperature buffering in the lake, where water temperatures, especially in cold regions, are warmer than the mean annual temperatures (Hren and Sheldon, Reference Hren and Sheldon2012). For an air temperature of 11°C (~10°C colder than modern), annual average lake temperatures are estimated to be 14.8°C according to the transfer function given by Hren and Sheldon (Reference Hren and Sheldon2012), which would result in a snow-recharge estimate of ~78%. Alternatively, the effect of lower ground temperatures in the mountains could result in snow with lower δ18O values, resulting in an overestimate of the spring-water contribution from snow. Lower sea levels also lead to increased effective altitude, which could also lower δ18Osnow values. On the other hand, evaporation of inflowing water before it reaches the snail's habitat would result in higher δ18Owater values compared with spring discharge, which ultimately results in an underestimate of snow recharge. In addition, if the lake has a longer residence time, which is likely given the larger inferred lake size, then the buffering effect of lake water mixing with the input from local precipitation and spring water would also need to be taken into account.
Archaeological implications
The results presented here are based upon a small sample. Nevertheless, when compared with pollen and archaeological sequences, they demonstrate the great potential of isotopic sclerochronological studies for high-resolution, seasonal reconstruction of the paleoenvironment and its impact on human presence and activity in the region. The origin of the shells from a long, well-dated sequence and defined cultural scheme, combined with subsampling of each shell, makes sclerochronological analysis a powerful tool for comparison and calibration of the different environmental proxies present at JRD. Importantly, the evidence for more snow during the glacial conditions in the Hula Valley can explain the dry signal obtained for this interval from the JRD pollen diagrams (Langgut et al., Reference Langgut, Cheddadi and Sharon2021). Although this signal is normally interpreted as indicating dry conditions, cold climatic conditions with more snow may result in less arboreal pollen.
The JRD archaeological sequence documents 10,000 years of human activity on the shore of paleo-Lake Hula. During this time, the region witnessed the dramatic cultural and economic transition from small bands of nomadic hunter- (and fisher) gatherers to the early sedentary settlements of the Natufian to the agricultural communities of the Neolithic era. The results presented here contribute to a better understanding of the environmental background for these cultural and economic changes. The evidence for warmer and wetter conditions in Layer 3C correlates well with the results of other environmental proxies from JRD (Sharon et al., Reference Sharon, Grosman, Allué, Barash, Bar-Yosef Mayer, Biton and Bunin2020; Langgut et al., Reference Langgut, Cheddadi and Sharon2021), which mark changing conditions at the onset of the Natufian. Interrelations between human cultural changes and paleoenvironment changes are highly complex and debatable; hence, reliable environmental reconstruction is a highly significant foundation for any such discussion.
CONCLUSIONS
Late Pleistocene and modern Melanopsis shells from the Hula Basin were sampled to produce δ18O and δ13C data with submonthly resolution. However, growth marks present in the shells indicate that shell growth slows or stops during winter, inhibiting a full assessment of seasonality.
Applying the results of the modern shell to subfossil isotope data from JRD yielded valuable information about seasonality and the environment of the Hula Basin. Shells from layers dated to the LGM exhibited small fluctuations in δ18Oshell and δ13Cshell values, indicating a large, buffered waterbody. The lack of strong correlation between δ18Oshell and δ13Cshell values indicates that the paleolake was hydrologically open with a relatively short residence time. The data support an increased water contribution from the Hula Valley during the Lake Lisan high stand. Furthermore, relatively high δ18Oshell values immediately after growth marks and a subsequent trend toward lower values may indicate that local precipitation had higher δ18Owater values than spring waters. This observation supports a greater proportion of snowmelt contribution in regional groundwater recharge during this time interval, perhaps more than double the amount of present-day snowmelt.
In contrast, data from shells that likely formed during the drier conditions of Heinrich 1 indicate a significantly smaller waterbody, and data from Bølling-Allerød shells point to an intermediate-sized waterbody. Fluctuations in waterbody size at JRD qualitatively agree with changes in lake level of Lake Lisan, supporting the Hula Valley as a major contributor to changes in water balance in the catchment.
These results agree well with palynological results from JRD (Langgut et al., Reference Langgut, Cheddadi and Sharon2021) and provide additional context for understanding other environmental proxies collected at the site. Furthermore, these data give insight into weather-scale events and seasonality during the final part of the late Pleistocene in the Levant.
Acknowledgments
The authors would like to thank Frank Riedel for providing lab space. We would also like to thank Frank Riedel and Niels de Winter for valuable input on earlier drafts of the article. This research was supported by the University of Iceland Research Fund (Rannsóknasjóður Háskóla Íslands) and the Erasmus+ program of the European Union. The JRD excavation (Israel Antiquity Authority license no. G-75/2018) is supported by the Israel Science Foundation (GS, grant no. 918/17). Constructive comments and suggestions of two anonymous reviewers, associate editor Bob Booth, and the journal's senior editor Nicholas Lancaster helped to improve the original article.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/qua.2023.25.
Open access
