Implications of submonthly oxygen and carbon isotope variations in late Pleistocene Melanopsis shells for regional and local hydroclimate in the upper Jordan River valley

Abstract

Many water-stressed regions of the globe have a highly seasonal precipitation regime. However, seasonality in the past and under changing climates is little studied. Submonthly records of sclerochronological δ¹⁸O and δ¹³C values of Melanopsis shells from the Jordan River Dureijat archaeological site (JRD) in the upper Jordan River valley presented here document the hydrology of paleo-Lake Hula. These records were assessed for changes in seasonal hydrology in the lake and compared with modern shells collected from present-day waterbodies in northern Israel and with models of δ¹⁸O_shell. Results from shells in sediments dating from the last glacial maximum (LGM) to the Bølling-Allerød imply changes in waterbody size that qualitatively parallel changes in the late Pleistocene Lake Lisan levels; Hula Lake was well buffered when Lake Lisan stood at a high stand and poorly buffered when water levels were lower. Furthermore, data from shells dated to the LGM suggest inflowing water with lower δ¹⁸O values than local rainfall, providing evidence for a greater proportion of snow in the catchment than today. Reconstruction of water δ¹⁸O and mixing-model calculations suggest that snowmelt contribution to spring water during the LGM may have been more than twice the amount in the modern-day catchment.

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, 2010). 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., 2011; Enzel and Bar-Yosef, 2017). However, few existing paleoclimate records can be correlated with human activity at archaeological sites (Rambeau, 2010; Enzel and Bar-Yosef, 2017). In addition, very few studies of pre-Holocene paleoseasonality have been performed in this region (Ben Dor et al., 2018, 2021; Langgut et al., 2021; Müller et al., 2022). 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, 2020). 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, 2010).

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, 1995; Bar-Matthews et al., 1999, 2003; Bartov et al., 2002; Robinson et al., 2006; Roberts et al., 2008; Develle et al., 2011; Orland et al., 2012; Torfstein et al., 2013; Miebach et al., 2017; Matthews et al., 2021). 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., 2002; Hazan et al., 2005; Torfstein et al., 2013; Lev et al., 2019). For example, pollen records indicate a prevalence of steppe vegetation in the region, possibly indicating low water availability (Miebach et al., 2017, 2019; Langgut et al., 2021). 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., 1999, 2003; Vaks et al., 2006). Enzel et al. (2008) 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 δ¹⁸O_speleothem values to the south and north compared with central Israel during the LGM, which may indicate low δ¹⁸O_{precipitation} values due to Rayleigh distillation processes outside the main storm path (Vaks et al., 2006; Bar-Matthews et al., 2019). Evidence of flooding events in Lake Lisan during the LGM support more frequent southerly storm paths (Ben Dor et al., 2018). 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., 2019; Langgut et al., 2021). 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., 2016; Ludwig and Hochman, 2022). 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., 2011). The lack of speleothem growth on Mt. Hermon during the LGM implies a wider area of freezing temperatures (Ayalon et al., 2013), 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., 2011), supported by evidence of glaciers in the Mt. Lebanon area during this interval (Moulin et al., 2022). Zaarur et al. (2016) reconstructed low δ¹⁸O_water values from carbonate clumped isotopes in gastropod shells in the upper Jordan River valley, while speleothems in the more southerly Soreq Cave exhibited higher δ¹⁸O_speleothem values during the LGM (Bar-Matthews et al., 2003). Zaarur et al. (2016) interpret their reconstructed LGM δ¹⁸O_water 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., 2020; Wiese et al., 2020) and provide a seasonal perspective on geologic-scale climate change (Taft et al., 2014; Hartman et al., 2016). 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 (δ¹⁸O_shell and δ¹³C_shell) vary as the animal grows due to changes in temperature and changes in δ¹⁸O_water values and δ¹³C of dissolved inorganic carbon (DIC) in the freshwater environment (Leng and Marshall, 2004). Snails of the genus Melanopsis occur in freshwater bodies in the Mediterranean region and are common in the Levant, including the Hula Valley (Tchernov, 1973; Glaubrecht, 1993; Heller et al., 2005). Modern snails of this genus occur as far north as Hungary but are presently limited to thermal springs in that area (Glaubrecht, 1993). Life span estimates for this genus vary widely. Elkarmi and Ismail (2006) 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. (2016) reported sclerochronological δ¹⁸O_shell values from modern specimens in the Jordan River valley that reflect two seasonal cycles. Geary et al. (2012) 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., 2012; Bartolini et al., 2017).

The objective of our study is to use high-resolution sclerochronological δ¹⁸O_shell and δ¹³C_shell 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, 1973; Rybakov et al., 2003; Schattner and Weinberger, 2008; Heimann et al., 2009). The basin's area is approximately 150 km², and it serves as an erosion base level for approximately 1500 km². 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., 2003).

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., 2020), with a snowmelt contribution of about 30% (Gil'ad and Bonne, 1990; Sade et al., 2016). 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., 2011). 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., 2020).

Before drainage in the 1950s, the Hula Valley held a shallow, up to 3-m-deep lake (~14 km²) at its southern end and extensive swamps to the north (Dimentman et al., 1992). 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., 2003). 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., 2003; Litaor et al., 2008).

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., 2009; Sade et al., 2016). 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, 2006; van Zeist et al., 2009; Sade et al., 2016).

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., 2015; Sharon et al., 2020). 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., 2020). 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., 2020; 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. (2013).

Figure 2. Stratigraphy and radiocarbon ages at Jordan River Dureijat archaeological site (JRD; adapted from Sharon et al., 2020). 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, 2017), 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., 2020). Similar to other sites on the banks of the upper Jordan River (Goren-Inbar et al., 2002a, 2002b; Aharonovich et al., 2014; Melamed et al., 2016), 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, 2013). 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., 2021).

MATERIALS AND METHODS

Shells from the JRD Area B east section (the type section of the site; Sharon et al., 2020; 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 ¹⁴C 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., 2020). 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. (2020).^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. (2005).

Samples were milled from the shell surface using a modified dental drill placed in a custom holder after Wiese et al. (2020). 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 mm² 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, 2010; 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 (δ¹³C = 1.95‰, δ¹⁸O = −2.2‰) and NBS18 (δ¹³C = −5.014‰, δ¹⁸O = −23.2‰). The measured sample data were adjusted using the calibrated internal laboratory standard C1 (δ¹³C = 2.40‰, δ¹⁸O = −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 δ¹³C and δ¹⁸O. To correct for acid fractionation of aragonite, −0.38‰ offset was applied after Kim et al. (2007a). Isotope ratios of shells are reported in the delta notation relative to VPDB, and the cited and calculated oxygen isotope values of water (δ¹⁸O_water) are expressed relative to Vienna Standard Mean Ocean Water (VSMOW).

Due to low intensities in the smallest samples of the upper whorls, δ¹³C and/or δ¹⁸O 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 δ¹⁸O and δ¹³C values are reported to 0.1‰ precision. All shells were photographed after milling to document sample locations (Rice, 2018).

Interpretation framework for δ¹⁸O and δ¹³C profiles

As gastropods precipitate their shells, the δ¹⁸O_shell values primarily reflect those of the lake water and an offset due to temperature effects (Schmitz and Andreasson, 2001; Spiro et al., 2009). Precipitation, evaporation, and inflowing water can alter δ¹⁸O_water in the snail's habitat. Evaporation leads to higher δ¹⁸O_water values in the lake, while the input of freshwater in the form of precipitation or inflow typically leads to lower values (Leng and Marshall, 2004). Meltwater from snow and rainwater that falls at high altitudes to recharge springs can contribute to low δ¹⁸O_water values in the Jordan River (Gat and Dansgaard, 1972; Zaarur et al., 2016). In addition, some storms can bring high-δ¹⁸O rainwater. Sub-cloud evaporation of raindrops and convective rain that originates from evaporated water over land will result in high δ¹⁸O_water values in rain (Gat and Dansgaard, 1972). Offsets in δ¹⁸O_shell 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 δ¹⁸O_shell values but to increased evaporation, which results in higher δ¹⁸O_water values. Because δ¹⁸O_shell values change by only 0.24‰/°C (Kim et al., 2007b), the changes in δ¹⁸O_water values due to evaporation are often more important. Similarly, colder temperatures lead to higher δ¹⁸O_shell values but coincide with the rainy season, which leads to lower δ¹⁸O_water values in the lake. It is possible that subfossil shells grew during intervals of lower evaporation, in which case temperature effects on δ¹⁸O_shell values might be more important. In this case, decreasing δ¹³C_shell values are still expected with lowering temperatures due to decay of organic matter, and we rely on δ¹³C_shell values to indicate possible pauses in growth during winter.

Meanwhile, δ¹³C_shell values will reflect the δ¹³C_DIC values in the water. In lake water, δ¹³C_DIC values are impacted by the δ¹³C values of inflowing DIC, exchange with the atmosphere, and photosynthesis and respiration of aquatic plants (Leng and Marshall, 2004). Aquatic productivity is generally the most important factor on short (subannual) timescales and leads to increased values of δ¹³C_DIC in the lake, whereas the decay of organic matter leads to lowered values of δ¹³C_DIC (Leng and Marshall, 2004). Storms wash soil into the lake, bringing low-δ¹³C carbon into the lake system.

Growth marks can also aid in interpreting the δ¹⁸O_shell and δ¹³C_shell 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., 2007; Mannino et al., 2008; Geary et al., 2012). Some studies have noted the formation of these features during times of thermal stress, both for cold and hot temperatures (Schöne et al., 2007; Mannino et al., 2008). Gastropods may respond to thermal stress by thickening their shell rather than lengthening it in order to improve insulation (Mannino et al., 2008). 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 δ¹³C_shell values, this event is interpreted as a pause in growth during winter months, with changing δ¹³C_DIC 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-δ¹³C soil is washed into the lake. For growth marks preceded by several samples with elevated δ¹³C_shell values or a long trend toward higher δ¹³C_shell values, the growth mark likely represents winter. Growth marks that coincide with a drop in δ¹³C_shell 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., 2012). 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 δ¹⁸O_shell values are interpreted as changes in δ¹⁸O_water and temperature. Decreases in δ¹⁸O_shell values result from the input of low-δ¹⁸O water into the lake, either through direct precipitation, inflow, runoff, or snowmelt. Temperature increases can also result in decreases in δ¹⁸O_shell values; however, temperature is likely to be a relatively minor effect, with a temperature increase of 4°C resulting in a decrease in δ¹⁸O_shell of about 1‰ (Kim et al., 2007b). Trends toward higher δ¹⁸O_shell values generally indicate ¹⁸O enrichment of lake water due to evaporation in summer months. High-δ¹⁸O rainfall can also raise δ¹⁸O_lake 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 δ¹⁸O_water values of winter rainfall within the valley compared with summer lake water fed by the springs, which are recharged by low-δ¹⁸O_water snowmelt. Elevated values of δ¹⁸O_shell immediately following the winter growth mark may indicate high values of δ¹⁸O_water during winter compared with summer, likely indicating the effect of snowmelt on δ¹⁸O 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 δ¹⁸O_shell patterns

The ShellChron model package (de Winter, 2022) includes a script to generate a modeled δ¹⁸O_shell record from sinusoidally varying temperature and δ¹⁸O_water 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. (2016). Temperature and δ¹⁸O_water values reported in the literature were used to generate this record (Gat, 1970; Gat and Dansgaard, 1972; Gil'ad and Bonne, 1990; Litaor et al., 2008; Gophen et al., 2016; Zaarur et al., 2016; Supplementary Text S1). To mimic weather effects, model temperature and δ¹⁸O_water 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 δ¹⁸O_water values, respectively. The code was adapted to use the Kim et al. (2007b) temperature-dependent oxygen isotope relationship between aragonite and the water from which it precipitates:

(Eq. 1)$$1000 \; {\rm ln} \; \alpha = 17.88 \times 10^3/{\rm T} - 31.14$$

where α = (δ¹⁸O_shell + 1000)/(δ¹⁸O_water + 1000). Solving for δ¹⁸O_shell:

(Eq. 2)$$\delta ^{18}{\rm O}_{\rm shell} = ( \delta ^{18}{\rm O}_{\rm water} + 1000) \times e^{( 17.88 \times {10}^3/T-31.14) /1000}-1000$$

Values from Eq. 2 yield δ¹⁸O_shell values in the VSMOW reference frame, which are then converted to VPDB for comparison with measured δ¹⁸O_shell 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 δ¹⁸O_water regimes in fossil shells. The temperature and δ¹⁸O_water inputs were fit to shell results. Rather than attempting to model each variation in each shell, the overall variability of δ¹⁸O_shell 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., 2005). However, whether these morphological differences constitute a separate species in Levantine Melanopsis is under dispute (Falniowski et al., 2020). 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 δ¹⁸O_shell and δ¹³C_shell results are provided in Supplementary Table 1. Means of sclerochronological δ¹⁸O_shell and δ¹³C_shell values for each shell range from −7.9 to −6.4‰ and −7.7 to −4.2‰, respectively. The smallest range of δ¹⁸O_shell values, 2.3‰, was calculated for Layer 5 Shell 1, and the largest range, 6.0‰, for the modern shell (Table 2). The δ¹³C_shell 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 δ¹⁸O_shell (blue circles) and δ¹³C_shell (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 δ¹⁸O and δ¹³C values (‰).

Shell sclerochronology interpretations

The modern shell from Agamon Hula has large changes in δ¹⁸O_shell values consistent with a large evaporation effect (Fig. 3a). Winter growth marks (G and H) are identified as the largest decreases in δ¹⁸O_shell and δ¹³C_shell values. Two annual cycles are inferred, each of which shows a trend toward high δ¹⁸O_shell 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 δ¹⁸O_water values is not observed.

The shells from Layer 3C have large and concurrent decreases in δ¹⁸O_shell and δ¹³C_shell values that indicate winter (Shell 1 growth marks B and C; Shell 2 growth marks A and D), but show more gradual increases in δ¹⁸O_shell values compared with δ¹³C_shell values (Fig. 3b and c). Layer 3C Shell 1 records two clear annual cycles and has two additional large decreases in δ¹³C_shell 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 δ¹³C_shell values and increased δ¹⁸O_shell 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 δ¹⁸O_shell and δ¹³C_shell values (Fig. 3d and e). The frequency and timing of changes in δ¹⁸O_shell and δ¹³C_shell 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 δ¹³C_DIC 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 δ¹³C_shell 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 δ¹⁸O_shell values occur after these marks, which could be due to either increased temperature or lower δ¹⁸O_water values during summer.

In Layers 5 and 6, minima in δ¹³C_shell values generally coincide with peaks in δ¹⁸O_shell values (Fig. 3f–i). Winter growth marks identified by decreases in δ¹³C_shell 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 δ¹⁸O_shell 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 δ¹⁸O_shell values for the modern Jordan River and Sea of Galilee reflect the general range of measured δ¹⁸O_shell 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., 2016), 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 δ¹⁸O_water used in these scenarios does not appropriately describe the conditions under which the shells grew. Given these discrepancies, the δ¹⁸O_water 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 δ¹⁸O_shell values (Fig. 4).

Figure 4. Model results of δ¹⁸O_shell for modern (top) and paleo (bottom) scenarios.

Water temperature and δ¹⁸O_water parameters that best fit fossil shell patterns (Table 3) suggest that δ¹⁸O_water 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 δ¹⁸O_water values during cold conditions. Furthermore, seasonal amplitudes of modeled δ¹⁸O_water 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 δ¹⁸O_shell 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, δ¹⁸O_water, and δ¹⁸O_shell parameters and comparison with measured δ¹⁸O_shell 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., 2016) 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 δ¹⁸O_water values (1‰) and exhibits the narrowest range of seasonal δ¹⁸O_shell values (~1.2‰). Variations in δ¹⁸O_shell 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 δ¹⁸O_shell 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, 2001; Spiro et al., 2009). A strict temperature cutoff of 14.5°C imperfectly describes Melanopsis shell growth. The patterns of modeled δ¹⁸O_shell 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. (2016) exhibit sawtooth-shaped δ¹⁸O_shell 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 δ¹⁸O_shell results suggest that the general variability within the gastropod's environment can be understood through the extent of seasonal variations in δ¹⁸O_water values and, to a lesser extent, seasonal temperature variability. Due to model assumptions, short-term changes in δ¹⁸O_water values cannot be replicated, and seasonal changes that are not purely sinusoidal have not been assessed. Due to the lack of subseasonal δ¹⁸O_water 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 δ¹⁸O_shell 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., 2021), 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, 2000). Storms can lead to highly variable δ¹⁸O_water values in the modern-day Jordan River (Gat and Dansgaard, 1972), with Hula Valley δ¹⁸O_{precipitation} values reaching as low as −11‰ (Rindsberger et al., 1990), 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., 2021), in agreement with shell-based inferences.

Variations in whole-shell δ¹⁸O values can provide further insight into environmental changes. Over the studied time interval, the average δ¹⁸O_shell values exhibit an opposite trend to δ¹⁸O_speleothem records from Soreq Cave and Peqiin Cave, in line with previous findings (Zaarur et al. 2016; Fig. 5a). However, while Zaarur et al. (2016) attribute low Melanopsis δ¹⁸O_shell 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 δ¹⁸O values of winter precipitation. In the general trend of increasing δ¹⁸O_shell 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., 2016; Miebach et al., 2019; Langgut et al., 2021; Ludwig and Hochman, 2022). However, the shells from Heinrich 1 (Layer 4), during which lake levels dropped suddenly, have lower δ¹⁸O_shell 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., 2003; Bartov et al., 2003; Langgut et al., 2021). The cold conditions likely also resulted in low summer evaporation, though this has not yet been studied. Nevertheless, the overall pattern of increasing δ¹⁸O_shell values indicating greater summer evaporation aligns with sclerochronological evidence from JRD, with highest δ¹⁸O_shell values recorded by the shells from the Bølling-Allerød (Layer 3C).

Figure 5. (a) A comparison of Melanopsis δ¹⁸O_shell from this study (blue squares) and Zaarur et al. (2016; pink squares) to δ¹⁸O_speleothem records from Soreq (gray), Peqiin (dark green), and Zalmon (light green) Caves (Bar-Matthews et al., 2003; Keinan et al., 2019) and (b) standard deviation of δ¹⁸O_shell values in Melanopsis (blue squares; note the inverted axis) compared with the Lake Lisan lake level curve (Torfstein et al., 2013). GBY, Gesher Benot Ya'aqov; JRD, Jordan River Dureijat archaeological site.

The results of sclerochronological δ¹⁸O_shell 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 δ¹⁸O_shell values are small, indicating a large, well-buffered waterbody at a time when Lake Lisan stood at a high stand (Bartov et al., 2002; Torfstein et al., 2013). After the LGM, Lake Lisan was lower, and lake levels dropped dramatically during Heinrich 1 (Bartov et al., 2003), when low water levels are also inferred at paleo-Lake Hula from highly fluctuating δ¹⁸O_shell values. Lake Lisan again rose, reaching high levels at about 15 ka, though not as high as LGM levels (Torfstein et al., 2013). 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 δ¹⁸O_shell values seems to vary with the Lake Lisan curve from Torfstein et al. (2013; Fig. 5b). When the standard deviation is low, δ¹⁸O_water values in the lake exhibit little fluctuation and Lake Lisan levels are high, while greater standard deviations, reflecting large fluctuations in δ¹⁸O_water 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 δ¹⁸O_water 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 δ¹⁸O_shell values show that the timing of maxima in δ¹⁸O_water values in relation to maxima in water temperature can impact the magnitude of variability in δ¹⁸O_shell values. When summer evaporation leads to the highest δ¹⁸O_water values in the lake, the seasonal extent of δ¹⁸O_shell values is smaller than when the highest δ¹⁸O_water values occur in winter. This is due to the temperature fractionation effect, in which increased temperatures lead to lower δ¹⁸O_shell values at a given δ¹⁸O_water value. This effect diminishes the peak δ¹⁸O_shell values when high temperatures coincide with high δ¹⁸O_water values. So, for environments in which evaporation controls δ¹⁸O_water values and the highest δ¹⁸O_water values occur during the hottest parts of the year, the seasonal variability in δ¹⁸O_shell values will be smaller than for environments that experience a different seasonal regime. This complicates a direct comparison of variability in δ¹⁸O_shell 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., 2011; Orland et al., 2012; Ayalon et al., 2013; Zaarur et al., 2016; Moulin et al., 2022). 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, 1990; Sade et al., 2016). A consistent, strong negative δ¹⁸O_shell 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., 2016). This would contribute water with lower δ¹⁸O values to the main springs feeding the lake, resulting in lower δ¹⁸O_water values in summer due to recharge from snow-fed springs and higher δ¹⁸O_water values during winter, when precipitation falls as rain in the Hula Valley. The patterns in δ¹⁸O_shell values from Layers 4–6 closely resemble modeled δ¹⁸O_shell patterns in which δ¹⁸O_water 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, δ¹⁸O_water values in precipitation were approximately 2.5‰ higher than present (Bar-Matthews et al., 2003), and modern monthly δ¹⁸O_water values in Mt. Hermon precipitation range from approximately −10.5 to −5.0‰ (Ayalon et al., 2013) and serve as the modern endmember δ¹⁸O_water values for snow and rain recharge to the springs, respectively. Thus, the equivalent δ¹⁸O_water 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 δ¹⁸O_water values are expected to reflect the endmembers of local precipitation (winter) and spring water (summer), respectively. Based on the model LGM shell, the summer δ¹⁸O_water endmember (−7.1‰) is composed of 84% snowmelt. This is much greater than the modern values of 30% (Gil'ad and Bonne, 1990) and the estimate of at least 40% contribution of snow during the LGM by Zaarur et al. (2016).

However, these calculations are subject to considerable uncertainty. The model of δ¹⁸O_water relies on the modeled temperature during shell formation. Changing the temperature would result in a change in modeled δ¹⁸O_water 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., 2016). 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., 2008; Matthews et al., 2021). 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, 2012). 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 (2012), 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 δ¹⁸O 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 δ¹⁸O_snow values. On the other hand, evaporation of inflowing water before it reaches the snail's habitat would result in higher δ¹⁸O_water 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., 2021). 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., 2020; Langgut et al., 2021), 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 δ¹⁸O and δ¹³C 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 δ¹⁸O_shell and δ¹³C_shell values, indicating a large, buffered waterbody. The lack of strong correlation between δ¹⁸O_shell and δ¹³C_shell 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 δ¹⁸O_shell values immediately after growth marks and a subsequent trend toward lower values may indicate that local precipitation had higher δ¹⁸O_water 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., 2021) 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.