Holocene wet episodes recorded by magnetic minerals in stalagmites from Soreq Cave, Israel

Speleothem paleoenvironmental proxy records focus on stable oxygen and carbon isotope time series (δ¹⁸O and δ¹³C, respectively) (e.g., Fairchild and Baker, 2012; Bar-Matthews et al., 2019). Coupling these data with other climate proxies increases the accuracy and precision of paleoclimate interpretations (Fairchild and Baker, 2012). These other proxies include nontraditional isotopes, trace element geochemistry, and, more recently, mineral magnetism (e.g., Bourne et al., 2015; Zhu et al., 2017). The underlying assumption is that speleothem magnetic properties are controlled by ferromagnetic particles flushed into the cave and encapsulated within calcite laminae (Lascu and Feinberg, 2011). These particles can be detrital in origin (dominated by regionally ubiquitous magnetic minerals such as magnetite, titanomagnetite, and hematite), or authigenic (mainly pedogenic magnetite and goethite) (Lascu and Feinberg, 2011; Bourne et al., 2015; Jaqueto et al., 2016). The concentration, mineralogy, and size of the magnetic particles can be measured indirectly using non-destructive methods, providing a proxy for the physicochemical properties of the overlying soil and the cave hydrology (Lascu and Feinberg, 2011; Feinberg et al., 2020). However, the interpretation of the magnetic data is not straightforward because the link between climate and magnetic parameters is controlled by the local geohydrology. For example, data from mid-latitude stalagmites suggest that more magnetic particles are delivered into a cave as rainfall increases (e.g., Bourne et al., 2015; Zhu et al., 2017), while the opposite relationship is observed in tropical caves, where magnetic particles are delivered primarily during dry periods (Jaqueto et al., 2016; Fu et al., 2021). Thus, adaptation of speleothem magnetism in paleoclimate models requires a site-specific multiproxy analysis that includes both magnetic and geochemical proxies.

Our study examines the coupling between magnetic parameters and δ¹³C values over centennial to millennial timescales in two Holocene stalagmites from Soreq Cave, Israel (31°45′21′′N, 35°01′20′′E). The stalagmite record spans the time interval from ca. 9.7 ka to ca. 5.4 ka, with an overlap between ca. 7.0 and ca. 6.4 ka covering the early-Holocene pluvial conditions associated with the deposition of Sapropel layer 1 (S1) in the eastern Mediterranean Sea (Bar-Matthews et al., 2000), followed by the transition to mid-Holocene wet-dry cycles (Bar-Matthews and Ayalon, 2011). We demonstrate the interplay between magnetic and isotopic proxies during regional climate change derived from speleothems deposited under two drip-sites with different hydrochemical settings. This is part of an effort to further develop speleothem magnetism as a standardized paleoenvironmental proxy.

Soreq Cave (Fig. 1A) is located on the western slopes of the Judean Mountains, 40 km east of the Mediterranean Sea at an elevation of ~400 m above mean sea level. The cave area is ~5000 m², hosted within Cenomanian dolomitic limestone. Bedrock above the cave varies between 10 m to >50 m. The overlying soil occurs as 30–100-cm-thick pockets covering ~50% of the present-day surface (Figure S1 and Section S1 in the Supplemental Material¹). Flora above the cave are intermingled coniferous forest and C₃ shrubland (Danin, 1992). Mean annual rainfall is mm with >95% of rain limited to October through May (Ayalon et al., 1998).

Soreq Cave is situated on the boundary between the Mediterranean semi-arid climate zone to the north and the arid desert environment to the south, thus acting as a sensitive recorder of changes in hydroclimate and seasonality (Bar-Matthews et al., 2019). The composite isotope record from Soreq secondary calcite speleothems spans ~180 ka and is the quintessential proxy time series for the eastern Mediterranean (Bar-Matthews et al., 2019, and references therein; Section S1 in the Supplemental Material). In addition, Soreq data bracket the timing of the terrestrial occurrence of sapropel S1 (Almogi-Labin et al., 2009) and reveal mid-Holocene wet-dry cycles simultaneous with cultural transitions in the Levant (Bar-Matthews and Ayalon, 2011).

Stalagmite 11-24 (from zone 11 in the cave) is a 425-mm-long stalagmite cut along its growth axis. Bar-Matthews and Ayalon (2011) studied one half of the stalagmite for reconstructing paleoclimate, and the second half is used here for paleomagnetic analysis (Fig. 1B; Fig. S2). Stalagmite 2-20 (from zone 2) was cut perpendicular to its growth axis near its base. The ~122-mm-long semi-major axis of the horizontal plane was used for paleoclimate reconstruction (Bar-Matthews et al., 2003; Grant et al., 2012) as well as for paleomagnetic analysis (Fig. 1C, Figs. S2 and S3). The isotopic time series in our study was constructed from 502 drilled samples in stalagmite 11-24 and 162 drilled samples in stalagmite 2-20. The age-depth models for the two stalagmites were recalculated with the StalAge R code (Scholz and Hoffmann, 2011; code is available in Comas-Bru et al. [2020]) using 14 of the published U-Th ages from stalagmite 11-24 (Bar-Matthews and Ayalon, 2011) and 16 U-Th ages from stalagmite 2-20 (12 published ages [Grant et al., 2012] and 4 new ages) (Fig. S3; see Section S2 for the age-depth model).

We performed paleomagnetic experiments in the paleomagnetic laboratory at the Institute of Earth Sciences, The Hebrew University of Jerusalem, using a 2G-RAPID cryogenic sample rock magnetometer equipped with in-line AF demagnetization coils and an ASC pulse magnetizer. The samples are thin slices, 2–5 mm thick with a median thickness of 4.2 mm. Stalagmite 11-24 was cut into 82 slices, and stalagmite 2-20 was cut into 44 slices (Figs. 1B and 1C). The samples were washed with distilled water to remove possible contamination from the saw, and a 1.5 T (tesla) DC field was applied to all specimens to amplify the samples’ weak natural remanent magnetization (NRM) resulting from a low concentration of magnetic minerals. The isothermal remanent magnetization (IRM) was normalized both to mass (IRM_mass) and time (IRM_flux) (Zhu et al., 2017). IRM_mass is a measure of the relative concentration of magnetic minerals, while IRM_flux provides the flux of magnetic particles encapsulated in the calcite per unit time, and compensates for growth rate changes. We used coercivity unmixing to characterize the magnetic components in 11 samples: five samples from stalagmite 11-24, three samples from stalagmite 2-20, two samples from the overlying Terra Rossa Mediterranean soil (at depths of 5 and 30 cm), and one sample from clay collected inside the cave. The procedure included 95 alternating field (AF) demagnetization steps from 1 mT to 250 mT after applying the 1.5 T field. Unmixing data were analyzed using the MAX UnMix program (Egli, 2004; Maxbauer et al., 2016).

A 2 × 3 cm slice cut from the middle of stalagmite 11-24 (Fig. 1B) was used for electron microscopy using an extreme high resolution scanning electron microscope (XHR-SEM; Magellan 400L). The slice was cleaned but not polished, to preserve the magnetic particles occurring between calcite crystals. See Section S3 of the Supplemental Material for the extended methods.

Electron microscopy coupled with energy dispersive spectroscopy (EDS) analysis shows the presence of 0.5 µm to ~2 µm anhedral ferromagnetic grains occurring at the interstices of speleothem calcite grains (Fig. S5). The morphology of the sub-micrometer particles resembles pedogenic magnetite (Strauss et al., 2013), where the larger ones may be weathered detrital particles. Coercivity spectra of the stalagmites samples (Fig. S6a) broadly agree with values reported in previously published speleothem studies (Lascu and Feinberg, 2011; Osete et al., 2012; Bourne et al., 2015; Zhu et al., 2017) and are in the range of pedogenic and detrital magnetite (Egli, 2004). The IRM unmixing data, interpreted with two end-member populations (high-coercivity [HC; median destructive field, MDF, >25 mT] and low-coercivity [LC; MDF <15 mT]), show a general agreement between the stalagmites, the overlying soil, and the clay material found in the cave. The contribution of the HC component in the soil and the clay is lower than that of the stalagmites (Section S2 and Table S2), indicating that the grains incorporated into stalagmites were filtered toward smaller grain sizes when infiltrating through the karst. The IRM unmixing of local dolomite bedrock, characterized by a weak IRM of <10⁻⁷ Am², revealed a significant (30%–50%) component with much higher coercivity (>250 mT), which could not be fully demagnetized with standard AF methods (Fig. S6b). This dolomite-driven component is not incorporated in the stalagmites. Our observations support the underlying assumption that the magnetic minerals in the speleothems are micrometer- to sub-micrometer–sized particles derived from overlying soil and filtered through the karst system.

The interval spanning ca. 9.7 to ca. 5.4 ka (Fig. 2) includes the absolute maxima and minima in the entire Holocene δ¹³C and δ¹⁸O records of Soreq (Bar-Matthews et al., 2019), representing the most extreme climate variations of the interglacial. Inspection of the overlapping time series between ca. 7.0 and ca. 6.4 ka shows that IRM_mass values of stalagmite 2-20 are approximately half an order of magnitude higher than those of stalagmite 11-24 (Fig. 2E). This difference is explained by the drip site–specific, nonlinear nature of karst hydrology (Treble et al., 2013) associated with the stalagmites that formed in different parts of the cave (Fig. S2). Specifically, the geochemical model for Soreq Cave (Burstyn, 2013; Section S1) estimates a higher fraction of soil-derived elements and a significant fraction of vadose flow in the shallower locations, where stalagmite 2-20 was formed. By contrast, IRM_flux values of the two stalagmites are within the same value range, and their mean IRM_flux is indistinguishable within 1σ uncertainty. This suggests that IRM_flux is a potential normalized index for comparing multiple samples from the same cave (Section S4 and Tables S3–S6).

The early-Holocene recorded in stalagmite 2-20 is marked by a rapid decrease in IRM_flux to minimum Holocene values (Fig. 2D). The low IRM_flux coincides with a distinct isotopic event between ca. 9.5 and ca. 7.1 ka, where δ¹³C values are uniquely high (−8‰ to −4‰; Fig. 2C) relative to the Soreq overall mean, approaching bedrock values of 2‰ (Section S1). Conversely, the termination of early-Holocene conditions is observed between ca. 7.2 ka and ca. 7.1 ka as stalagmite 2-20 δ¹³C values decreased back to C₃-like values (−10‰ to −12‰), simultaneous with IRM_flux returning to pre-sapropel values (Figs. 2C and 2D). This excursion in the stalagmite record corresponds with sapropel layer S1 in the eastern Mediterranean sea (striped gray bar in Fig. 2), thought to be warm and humid (Nehme et al., 2015; Bar-Matthews et al., 2019; Goldstein et al., 2020). Bar-Matthews et al. (2003) explained the δ¹³C excursion as meteoric water reaching the cave after little interaction with overlying soil due to either (1) recently eroded soils, as water percolating through minimal soil cover will be undersaturated in low-δ¹³C soil CO₂ (Genty et al., 2001), or (2) high runoff fraction, circumventing the soil pathway altogether. Both hypotheses are supported by a coeval drop in ⁸⁷Sr/⁸⁶Sr and initial (²³⁴U/²³⁸U) ratios toward bedrock ratios (Ayalon et al., 1999). The IRM_flux as a soil contribution proxy may help differentiate between these hypotheses.

Thus, the entire early Holocene is bracketed between the ca. 9.5 ka weakening of the soil-related isotopic component and the large decrease in IRM_flux, followed by the ca. 7.1 ka rehabilitation of the C₃ isotopic signal and IRM_flux values, which is coupled with evidence of soil regeneration in the circum-Mediterranean region (Regattieri et al., 2019). We believe that the structure of the low-IRM_flux event (rapid decrease/increase) and the regional evidence of soil regeneration in the termination of that event supports the hypothesis that a change in soil cover or soil stability resulted in a change in δ¹³C (and IRM_flux) values. Apparently, the post-glacial increase in rainfall destabilized soils, perhaps coupled with a “shut down” of soil production due to a decrease in dust input during the “green Sahara” phase (see Rohling et al., 2015).

The post-sapropel mean IRM_flux values suggest that the soil overlying Soreq Cave was regenerated and stabilized by the mid-Holocene. The four saw-tooth cycles in IRM_flux between ca. 7.0 and ca. 5.4 ka (numbered in Figures 2E and 3) are also evident in the δ¹³C and δ¹⁸O records (Section S4). These cycles were attributed to wet-dry events (Mayewski et al., 2004; Bar-Matthews and Ayalon, 2011). Bar-Matthews and Ayalon (2011) interpreted low δ¹⁸O values as an increase in regional rainfall superimposed with strong annual variations, while the δ¹³C cyclicity was smoothed and amplified by complex soil processes, as they represent changes in soil organic carbon. Thus, the post-sapropel regeneration of soil is linked to the emergence of a positive relationship between IRM_flux and rainfall amount, as also observed in other mid-latitude sites (Zhu et al., 2017; Regattieri et al., 2019; Feinberg et al., 2020).

The mid-Holocene record also highlights an apparent phase shift between the physical flux of magnetic particles into Soreq Cave and the chemical processes affecting the δ¹³C of drip water (Fig. 3). The IRM_flux peaks precede δ¹³C minima by ~60 to 120 yr, where a maximum correlation factor between the two records is achieved by a time-shift of 60 yr (Fig. S7). This lag can be explained by the mechanisms governing the two proxies: (1) An increase in rainfall typically results in increased production of pedogenic magnetic minerals (Maxbauer et al., 2017, and references within) simultaneously enhancing physical flushing of soil particles into the cave. (2) In contrast, the change in soil water content, temperature, vegetation cover, and soil activity that affect biogenic soil-CO₂ concentrations may require several decades to develop and affect δ¹³C values (Genty et al., 2001). Another hypothesis is that a rise in rainfall increases vegetation cover to a certain threshold, inhibiting transfer of fine-grained magnetite (e.g., Chen et al., 2019) and increasing the time lag between peak values of each proxy.

Soreq Cave demonstrates two opposite responses of the magnetic flux parameter to increasing rainfall. During the early-Holocene, low-IRM_flux and high-δ¹³C indicate minimal contribution from soils, adding to the isotopic-based hypothesis of soil removal and regeneration during the prevailing warm-humid climate conditions of sapropel S1 (e.g., Goldstein et al., 2020). During the mid-Holocene, couplets of high-IRM_flux and low-δ¹³C are associated with increased precipitation of wet-dry cycles. A phase shift between the couplets may elucidate soil and vegetation dynamics, crucial to interpreting subannual-resolution climate data (e.g., Orland et al., 2014). We show that IRM_flux in speleothems is a valuable proxy for regional changes in soil-related processes. IRM_flux could be used to resolve controversies in interpreting the complex behavior of isotopes in geological records, and could hold the potential to constrain soil carbon turnover response to climate events.

This research was supported by the U.S.–Israel Bi-National Science Fund grant #2016402, the U.S. National Science Foundation grant EAR-2044535, and the European Research Council Horizon 2020 research and innovation program, agreement No. 804490. We thank V. Gutkin for his help with SEM analysis, and L. Comas-Bru and M. Deininger for their comments. We thank three anonymous reviewers and K. Benison for their constructive reviews.