Cave ice can contain a wealth of paleo-climatic and geochemical information that is rapidly being lost with the melting of the temperate zone cryosphere. The karst areas of Slovenia host over 260 ice caves. We collected samples for stable water isotope, major ion, and nutrient analyses from two Slovenian ice caves. Samples included two shallow ice cores in Snežna Cave, collected ∼5 m apart, and an ice face profile in Ivačičeva Cave. All ice isotopic ratios reflected modern precipitation that could be described by high-elevation meteoric water lines. An offset suggested that fractionation and mixing processes of melted ice affected the isotopic signals. Cation concentrations of ice in both caves showed Ca ≫ Mg > Na > K. The high Ca2+ and Mg2+ contents and elevated HCO3 concentrations indicate that CaCO3 dissolution within the local karst landscape is a primary control on ice chemistry. Low concentrations and inconsistent profile patterns of other major ions and nutrients suggest atmospheric deposition and vadose zone leaching were also primary sources of ions to the ice. Differences in Cl and SO42− profile concentrations at similar depths in Snežna Cave imply that ice melting, water mixing, and re-freezing processes can affect the primary climatic signal stored in the ice. While temperate ice caves can be repositories of climatic information, secondary diagenetic processes that affect ice chemical composition alter the original signal. In addition to chemical analysis, physical processes within the caves must be studied at a small spatial scale to understand and interpret ice chemistry.

Ice caves occur in bedrock with perennial ice accumulation (Lipar et al., 2021b) that is hoar, stalactitic, massive (floored and stratified), and extrusive (Yonge and Macdonald, 1999). The preservation of ice in these caves is dependent on cave climate, which is closely tied to the climate and environment outside of the cave. Ice caves are a widely distributed portion of the cryosphere, but the ice in many caves is rapidly disappearing (Perşoiu and Onac, 2019). Ice caves are found throughout the Northern Hemisphere, often occurring at lower elevations and in locations where surface ice is absent, such as in Central Europe (Kern and Perşoiu, 2013). Permafrost conditions at high latitudes or high elevations enhance ice cave preservation, but ice caves have also been identified in locations as warm as the Mediterranean climate of Turkey (Perşoiu and Onac, 2019). During the past 20 to 25 years, studies have investigated ice caves around the world (Citterio et al., 2004; Fórizs et al., 2004; Clausen et al., 2007; Fang et al., 2010; Kern et al., 2011b; Perşoiu et al., 2011; Sancho Marcén et al., 2012; Higham and Palmer, 2017; Yonge et al., 2017; Carey et al., 2019, 2020; and Smirnov and Sokolov, 2022).

Ice may be created by snow compaction, freezing of karst drip water, or condensation of water vapor. Climatic factors that are modulated by seasonality, such as temperature and precipitation, control the rate of ice accumulation and stability within ice caves (Perşoiu and Pazdur, 2011). The preservation and accumulation of the ice formed are also facilitated by the geographic location and geomorphologic setting of a particular cave. When caves are located at lower latitudes, mechanisms such as cold air trapping, unidirectional air flow, and evaporative cooling preserve ice masses (Perșoiu and Onac, 2019). In order to understand the formation processes and preservation of ice in caves, the subsequent diagenetic processes (e.g., melting, re-freezing, mixing, sublimation, evaporation) must be analyzed for each individual cave, or even for a specific location within a cave (Clausen et al., 2007; Perşoiu et al., 2011). Water can enter a cave quickly through the vadose zone via fissures in the rock or slowly through water drip. Snow can accumulate through vertical and inclined cave entrances, and once in the cave, snow may melt and then freeze as ice or directly form ice through densification (Lipar et al., 2021b). Differences in the fast inflow and slow inflow of water delivery change the mechanisms that drive ice formation and, consequently, the chemical composition of ice. Mixing of incoming water and meltwater from older precipitation may further complicate interpretations (Kern et al., 2011a), as these physical hydrologic alterations can alter ice chemical compositions. In addition, evaporation can concentrate a chemical compound or cause chemical precipitation. Progressive freezing in the top ice layers may also concentrate solutes in residual water, which is then frozen as a solute-enriched bottom layer of ice (Citterio et al., 2004).

Nevertheless, ice caves can provide an archive of chemical solute and stable isotope composition in recent precipitation (Carey et al., 2019, 2020). There is discussion on the extent to which the ice found in caves may also serve as paleo-climate archives and indicators of climate change (Perşoiu et al., 2011; Zorn et al., 2020). Recent studies have demonstrated the utility of ice caves as climate archives when paired with analysis of ice formation, diagenesis, and dating methods (Perşoiu et al., 2017; Sancho et al., 2018; and Racine et al., 2022a). Cave ice δ18O and δ2H compositions and other climate proxies, such as pollen and microfossils stored in the ice, can be used to reconstruct climate conditions during the time of ice formation (Bǎdǎlutǎ et al., 2020). Many of the ice caves now under investigation have entrances that have only recently been uncovered as surface snow and ice cover have melted (Zorn et al., 2020). The possible future rapid disappearance of ice caves as paleo-climate archives (Kern and Perşoiu, 2013) has ironically been made possible by the current warming climate, which has melted surface ice and snow and allowed access to the ice caves.

The present and the future of these cave ice deposits are closely tied to climate. Ablation of ice within these caves has been related to the increases in percolation of warm water within the karst system (Holmlund et al., 2005), to ice melt through warming temperatures (Fórizs et al., 2004), and to enhanced summer rainfall (Perşoiu et al., 2021). Increased warming, atmospheric perturbations, and enhanced precipitation events in Central Europe may well lead to a potentially rapid reduction of cave ice deposits in the region over the next decades (Colucci et al., 2016). This cave ice loss would coincide with the rapid loss of small glaciers that has occurred in the karst regions of northern Italy and Slovenia over the past century (Triglav-Čekada et al., 2014; Čekada et al., 2016; and Lipar et al., 2021a).

More than 260 caves with permanent ice features have been found and mapped in Slovenia (Mihevc, 2018; Speleological Association of Slovenia, 2022; and Blatnik et al., 2023). We have previously reported on results of coarsely sampled vertical profiles from Paradana Cave, Snežna Cave, and Ivačičeva Cave (Carey et al., 2019, 2020). Here, we present results of the geochemical analysis of detailed profiles from Ivačičeva Cave and Snežna Cave, wherein we measured stable isotopic compositions along with major ions and nutrient elemental concentrations. Ivačičeva Cave is the highest cave in Slovenia with known stratified layers of ice and is geographically located in the northwest region of Slovenia, which hosts the majority of ice caves. These data were compared to Snežna Cave, which is a tourist cave with stratified layers of ice and is located further east. To our knowledge, this is the first attempt to drill ice cores from ice caves in Slovenia. The goal of this research was to investigate the utility of cave ice as a paleo-climatic proxy through the quantification of stable water isotopic compositions and chemical concentrations of the ice to identify patterns and discern environmental variables that control ice chemistry.

We hypothesized that the ice would reflect precipitation chemistry that has been altered by interactions with the vadose zone and epikarst as it flows into the caves. We tested this hypothesis by collecting samples for chemical analysis from a 2 m exposure in Ivačičeva Cave at the end of the 2019 melt season during the brief time when the cave is safely accessible for sampling and from two shorter (36 cm) ice cores drilled in Snežna Cave to explore differences on a small spatial scale. To investigate the meteorological influence on ice chemistry, we compared the data to pristine snow and firn chemistry of Alto dell’Ortles (3,830 m above sea level [a.s.l.]) in Italy. Mt. Ortles is among the highest summits in the Eastern European Alps, and the snow/firn layers demonstrate evidence of preserving the climate record (Gabrielli et al., 2010).

Stable water isotopes oxygen (18O/16O; δ18O) and hydrogen (2H/1H; δ2H) are useful natural tracers that can be used as a climate proxy to reconstruct paleo-climate. Previous studies have linked the isotopic composition of ice in caves to winter ambient temperature changes and changes in moisture sources (Perșoiu, 2018). The utility of δ2H and δ18O as climate proxies is limited by the fractionation effects that control the meteoric source water and occur during ice formation and by secondary diagenetic processes. Stable water isotope data are reported in delta notation (δ):
formula
in reference to the Vienna Standard Mean Ocean Water (VSMOW) standard. The global meteoric water line (GMWL) linear relationship between δ2H and δ18O is well defined in precipitation (Craig, 1961):
formula

This equation describes the average δ2H–δ18O relationship for global precipitation. On a local or regional scale, the δ2H–δ18O relationship can be characterized by local meteoric water lines, where the δ2H–δ18O divergence from the GMWL is caused by geographic, topographic, and/or meteorologic variables that are driven by kinetic and equilibrium liquid–vapor fractionation (Dansgaard, 1964). The effects that control local precipitation include the (1) latitude effect, where increasing latitude corresponds to a decrease in temperature and greater fractionation; (2) continental effect, where heavier isotopologues preferentially condense as an air mass travels inland from the ocean; (3) altitude effect, where heavier isotopologues preferentially condense as an air mass cools with increasing altitude; and (4) amount effect, where heavier isotopologues are preferentially condensed from cloud vapor, so precipitation incorporates lighter isotopologues with continuous rainfall (Dansgaard, 1964). All of these effects cause continual loss of the heavier isotopes of 2H and 18O in the precipitation and depletion of δ2H and δ18O signatures. The combination of effects can cause the localized isotopic signatures of meteoric waters that feed ice caves. The deuterium excess (d-excess = δ2H − 8 × δ18O) values preserved in ice can also provide information on the original moisture source of the water (Racine et al., 2022b).

When meteoric waters enter a cave as water or snow, further fractionation of source water occurs during ice formation. Under equilibrium conditions, the heavier isotopes (2H and 18O) are preferentially incorporated into ice due to differences in bond strength between heavy and light isotopes (Clark and Fritz, 1997). However, the degree of fractionation is dependent on ambient variables that may speed up or slow down the freezing process, creating kinetic fractionation conditions. Rapid freezing may result in a lesser degree of fractionation between heavy and light isotopes in the ice and residual water (Jouzel et al., 1999; Souchez and Tison, 1987). Partial melting and re-freezing can cause additional isotopic fractionation through the same enrichment process. Seasonal temperature fluctuations within a cave can result in ice melting, mixing of ice layers, and partial re-freezing of ice (Citterio et al., 2004; Clausen et al., 2007; and Ersek et al., 2018). The mixing, melting, and re-freezing of ice cause the isotopic composition of ice to deviate from the source water isotopic composition. Meltwater will likely be influenced by evaporation prior to re-freezing, and the evaporative influence will result in enrichment of 2H and 18O in the residual water (Craig and Gordon, 1965). Similar to isotopic fractionation during ice formation, the degree of residual water fractionation from evaporation is determined by ambient variables. These processes further alter the initial δ2H and δ18O values of the source water.

Sample Collection

Snežna Cave (Cave Register No. 1254) is situated at 1,504 m elevation, below the tree line in northern Slovenia (Figure 1). The cave ceiling thickness, from the surface to the roof of the chamber, is 40–50 m in the Entrance Chamber, where samples were collected from massive floor ice. The source of water to the cave is drip from water that has traveled through the vadose zone and epikarst. Samples from Snežna Cave (Figure 2) were collected from the 10-m-deep floor ice, which represents approximately 400 m3 of ice (Mihevc, 2018). We sampled at two locations ∼5 m apart using a SIPRE ice auger that had been rinsed in deionized water and wrapped in plastic wrap prior to use. Cores were drilled to the depth of a rock layer incorporated into the ice, which we assumed to be widespread throughout the ice block since the inflow of debris was constant from the entrance slopes and walls of the cave. Core A was 36 cm in length, and core B was 24 cm in length. The ice cores were returned to the laboratory, sectioned, and wrapped in new Al-foil for subsequent storage under freezing conditions prior to sampling for analysis. Both cores were then cut into 2 cm sections that were melted into pre-cleaned high-density polyethylene (HDPE) bottles. Deionized water was placed in contact with the Al-foil, and multiple samples of this water were collected and analyzed as “blanks” along with the melted ice samples (see below).

Ivačičeva Cave (Cave Register No. 2399) is situated at 2,471 m elevation, above the tree line in northwest Slovenia (Figure 1). The cave ceiling thickness is 60–70 m, and samples were taken from a stratified section of wall ice within the cave that is fed by drip water and condensation within the cave (Novak and Kustor, 1983). There was no massive floor ice in the section of the cave where samples were collected. Wall ice samples from Ivačičeva Cave were collected from the ice face described in Carey et al. (2020); however, this data series is much more detailed (40 samples collected every 5 cm over 200 cm) (Figure 3) than the six samples published in Carey et al. (2020) over about the same distance. An ∼1–2 cm layer of ice was removed from the ice face in Ivačičeva Cave before samples were drilled with Petzl Laser ice screws and chipped directly into polyethylene bags, with nitrile gloves worn during sample collection. The ice screws had a diameter of 1.4 cm. Therefore, the edges of consecutive samples were 3.6 cm apart, and the center of each sample was 5 cm apart. After melting, aliquots for ion and nutrient analyses were filtered through 0.45 μM pore-size Millipore filters into HDPE bottles that had been clean with 18.2 MΩ deionized water. Aliquots for isotope analyses were poured directly into 20 mL vials with no headspace. These water samples were shipped to The Ohio State University (Columbus, OH, USA) and refrigerated at 4°C prior to chemical and isotopic analyses within 1 week of arrival.

δ2H and δ18O Analysis

The 2H/1H and 18O/16O ratios of melted ice samples, here reported as δ2H and δ18O values relative to the VSMOW standards (δ2H = 0‰, δ18O = 0‰), were measured on a Picarro Wavelength Scanned-Cavity Ring Down Spectroscopy Analyzer Model L2130-I. Seven injections of 2.3 μL were analyzed per sample. The first three injections were discarded to avoid memory effects, while the last four injections were averaged to produce uncorrected δ2H and δ18O values. Samples were corrected by internal laboratory standards that had been calibrated to VSMOW at the Institute of Arctic and Alpine Research (INSTAAR) at the University of Colorado at Boulder through analysis by a dual-inlet mass spectrometer. Internal laboratory standards were Colorado (δ2H = −126.3‰, δ18O = −16.53‰), Nevada (δ2H = 104.80‰, δ18O = −14.20‰), Ohio (δ2H = −61.80‰, δ18O = −8.99‰), and Florida (δ2H = −9.69‰, δ18O = −2.09‰). Instrumental precision was 0.34‰ and 0.028‰ for δ2H and δ18O, respectively. In-run accuracy was 2.0‰ and 0.35‰ for δ2H and δ18O, respectively. Ice δ2H and δ18O data were compared to the GMWL (Craig, 1961), Slovenia regional meteoric water line (RMWL), Kredarica meteoric water line (MWL), Rateče MWL, and Zgornja Radovna MWL. The Slovenia RMWL was determined from samples collected at Ljubljana (n = 336, 1981–2010) at 282 m a.s.l., Kozina (n = 39, 2000–2003) at 484 m a.s.l., and Portorož (n = 84, 2000–2006) at 2 m a.s.l., reported by the International Atomic Energy Agency (IAEA) as part of the Global Network of Isotopes in Precipitation (GNIP) program (IAEA/WMO, 2022; Vreča et al., 2006, 2011). The unweighted Slovenian RMWL was calculated though an ordinary least square regression method (Hughes and Crawford, 2012). The Kredarica (2514 m a.s.l), Rateče (864 m a.s.l), and Zgornja Radovna (750 m a.s.l) MWLs were calculated from samples collected as part of the Slovenian Network of Isotopes in Precipitation (SLONIP) (Vreča et al., 2022; SLONIP, 2023). SLONIP is a platform that provides detailed information on the isotopic composition of precipitation from eight locations in the country to capture differences driven by climate, geography, and precipitation sources (Vreča et al., 2022). The reported precipitation weighted MWLs were calculated with a reduced major access (RMA) regression (Crawford et al., 2014).

Chemical Analyses

Un-acidified melted ice samples were analyzed for major cations (Na2+, Mg2+, K+, Ca2+) and anions (Cl, SO42−) using ion chromatography following previously described methods (Welch et al., 2010). Nutrients (N-NO3 + NO2, N-NH4+, and H4SiO4 [dissolved silica]) were measured calorimetrically using a Skalar SAN++ nutrient auto-analyzer system according to methods recommended by the manufacturer. N-NO3 + NO2 and N-NH4+ are hereafter referred to as NO3 + NO2 and NH4+. Bicarbonate (HCO3) values were calculated using the major cation and anion charge imbalance approach (Welch et al., 2010). Propagated error from instrument precision and accuracy was ≤10 percent for Ca2+, Mg2+, Cl, and N-NO2 + NO3. Propagated error for SO4 was 22 percent. Na, K, and H4SiO4 had errors <5 percent. Analytical blanks (N = 5) were run on deionized water in contact with Al-foil used in the field. Blanks for ice collection bags were also analyzed (Carey et al., 2020). All analytes in the blanks were below the limits of detection, with the exception of Ca2+ (0.9 micromole, μM) and Mg2+ (3.5 μM). All chemical analyses from this study are reported in Supplemental Material Table S1. Correlations between chemical components were conducted in JMP Statistical Software Version Pro 15 using a pairwise method.

Ice δ2H and δ18O Values

Snežna Cave δ2H and δ18O profiles for cores A and B ranged from −99.6‰ to −67.6‰ and from −14.1‰ to −9.99‰, respectively. Deuterium excess (d-excess) ranges for cores A and B were +11.1‰ to +13.4‰ and +11.3‰ to +12.8‰, respectively (Supplemental Material Table S1). The ice profiles had mean values for δ2H, δ18O, and d-excess of −81.2‰, −11.69‰, and +12.3‰ in core A and −78.8‰, −11.36‰, and +12‰ in core B. These mean δ2H, δ18O, and d-excess values were not significantly different (α = 0.05), but core A showed a larger range of variation in δ2H and δ18O (−31.9‰, −4.13‰) than core B (−13.6‰, −1.63‰). The mean δ2H (−80.3‰) and δ18O (−11.57‰) values of all the Snežna Cave ice analyses were significantly different than the mean values of the Ivačičeva Cave wall ice profile of −60.7‰ and −9.29‰, which were enriched in 2H and 18O in comparison to the Snežna Cave samples. Ivačičeva Cave δ2H and δ18O values ranged from −75.2‰ to −50.0‰ and from −11.1‰ to −7.93‰, respectively (Supplemental Material Table S1). The d-excess values for Ivačičeva Cave profile ranged from +12.1‰ to +15.7‰ (mean +13.6‰).

As illustrated in Figure 4, the δ18O and δ2H data from all profiles fall above the GMWL and the Slovenia RMWL (Craig, 1961; IAEA/WMO, 2022). Snežna Cave ice samples plot on or above the Rateče and Zgornja Radovna MWLs and below the Kredarica MWL. Ivačičeva Cave samples plot above all other MWLs, and either on or below the Kredarica MWL (Figure 4). The regression line of all Snežna Cave ice samples is δ2H = 7.82 × δ18O + 7.6, and the regression line of Ivačičeva Cave ice samples is δ2H = 7.80 × δ18O + 11.7. The Ivačičeva Cave ice intercept matches that of the Rateče MWL, but the cave ice slope is lower than the MWL.

The most enriched isotopic values occurred at profile depths of 20, 70, and 150–155 cm for the Ivačičeva Cave samples (Figure 5). The most depleted isotopic compositions for Ivačičeva Cave occurred at depths of 5, 30–35, 55, 105–110, 125, 170, and 200 cm, in the profile; however, these variations are considered minor. The most enriched values for Snežna Cave occurred at depths of 12 cm and 36 cm (core A) and 4 cm and 12 cm (core B). The most depleted Snežna Cave values were at depths of 8 cm (core A) and 24 cm (core B).

Ice Inorganic Chemistry

The ice profiles for Snežna Cave had higher concentrations of both sulfate (SO42−) and chloride (Cl) than concentrations observed in Ivačičeva Cave, with maximum values of ∼18 μM and ∼20 μM at ∼10 cm depth in the Snežna Cave cores A and B, respectively. All the SO42− and Cl values measured for Ivačičeva Cave were <5 μM. The concentrations of Ca2+ were similar for ice from both caves, with the highest values of ∼300 μM at 5 cm and ∼100 cm depth in the Ivačičeva Cave. NO3 + NO2 concentrations were the most variable, ranging between ∼9 and 1 μM, with both the highest and lowest values observed in the Ivačičeva Cave ice. For Snežna Cave, ice NO3 + NO2 concentrations ranged only from 1 to 5.5 μM (Figure 5 and Supplemental Material Table S1). As with the stable isotope profiles, there was no consistent pattern for the ion concentrations with depth in the Snežna Cave cores. Except for the bottom 20 cm, Cl and SO42− concentrations in the Ivačičeva Cave profile were universally low, whereas Ca2+ and NO3 + NO2 concentrations showed much more compositional variation.

Pearson correlation coefficients (r) of stable water isotopes, major ions, and nutrient species for all sampled ice profiles are presented in Table 1. Ice δ2H and δ18O values were negatively correlated with chemical ion concentrations at the 95 percent confidence interval, except for Ca2+. The depleted isotopic values (lower H and O isotope ratios) corresponded to higher chemical concentrations of Cl, SO42−, Na+, K+, and Mg2+. However, these correlations were not as robust or consistent when analyzing relationships for ice from each cave individually (Supplemental Material Tables S2 and S3). These same grouping of ions all had significant positive relationships with each other (Table 1). In contrast, Ca2+ had weak, positive, or negative relationships with these major constituents, with the exception of Mg2+ (r = 0.36, p < 0.05). Ca2+ and H4SiO2 also had a significant positive relationship (Table 1). NO3 + NO2 only showed weak positive and negative relationships with other constituents. NH4+ had significant positive relationships with most major cations, except Ca2+. Concentrations of NH4+ and Cl were also significantly correlated. We interpreted these relationships with caution because they are driven by 22 percent of the samples from Snežna Cave having NH4+ concentrations >1 μM. NH4+ concentrations in Snežna Cave core profiles ranged from <1 to 4.6 μM, whereas NH4+ concentrations at Ivačičeva Cave were <1.2 μM (Supplemental Material Table S1).

Correlations among data from both Snežna Cave and Ivačičeva Cave ice were not consistent between the two sites (Supplemental Material Tables S2 and S3). Many of the correlations observed in Snežna Cave ice were not observed for Ivačičeva Cave ice samples. Snežna Cave cores A and B had significant negative correlations between Cl and stable water isotopes, and among Cl and major ions Na+, K+, Ca2+, and nutrients NO3 + NO2, NH4+, and H4SiO4. Dissolved silica (H4SiO4) also had significant positive correlations with Mg2+ and Ca2+ and significant negative correlations with stable water isotopes (Supplemental Material Table S2). NO3 + NO2 had significant positive correlations with Na+ and K+. Concentrations of Cl, SO42−, Na+, and K+ for Ivačičeva Cave ice had significant positive correlations with each other, and NO3 + NO2 concentration in Ivačičeva Cave ice also had a significant positive relationship with K+ (Supplementary Table 3).

Cave Ice δ2H and δ18O

The δ2H and δ18O values for ice in Snežna Cave and Ivačičeva Cave reflect modern precipitation, but additional environmental processes played roles in controlling their isotopic compositions. The δ18O profile of Ivačičeva Cave wall ice shows greater variation than the Snežna Cave floor ice profiles (Figure 5). Due to differences in geographic location and geomorphological setting, we did not expect these data sets to reflect the same isotopic enrichment or depletion patterns. Furthermore, we cannot constrain the ages of ice in either cave, and the ice in these two locations could represent different time frames. The mean Ivačičeva Cave ice δ2H and δ18O values were similar to the mean of congelation ice (δ2H = −61.99‰, δ18O = −9.25‰) measured in a cave ice sample (M-17; Cave Register No. 5878) also in the Julian Alps in northwest Slovenia (Racine et al., 2022b). Snežna Cave δ2H and δ18O values were isotopically depleted in comparison to Cave M-17. However, the high (>10‰) d-excess values of both Ivačičeva Cave and Snežna Cave ice samples resembled those reported by Racine et al. (2022b), which are characteristic of a west and central Mediterranean moisture source (Gõmez-Hernández et al., 2013; Kern et al., 2020).

As shown in the δ2H-δ18O plot (Figure 4), the Ivačičeva Cave samples were isotopically enriched in comparison to Snežna Cave samples, and while the slopes of the regression lines were similar, the intercepts differed by ∼4‰. Both cave ice data sets plotted above the Slovenia RMWL (IAEA/WMO, 2022), indicating that the RMWL from lower elevations does not adequately describe the precipitation source of ice to the cave. This finding supports the purpose of the SLONIP program, which seeks to characterize the localized variation of Slovenian precipitation isotopic compositions (Vreča et al., 2022). It is important to measure the stable isotopic composition of localized precipitation to identify the initial isotopic composition of ice source water.

Ivačičeva Cave δ2H and δ18O

The mean Ivačičeva Cave ice δ2H and δ18O values are similar to the mean Rateče (−65.3‰, −9.55‰) and Zgornja Radovna (−60.6‰, −8.94‰) values, but the data plot above those MWLs (Figure 4). At 2,471 m a.s.l., the isotopic composition of water entering Ivačičeva Cave is probably better described by the Kredarica MWL (Figure 4). The weighted mean precipitation isotopic composition at Kredarica (n = 34) is lower than the Ivačičeva Cave samples, but this mean precipitation value likely accurately reflects the source water to ice in Ivačičeva Cave. Freeze-thaw cycles alter the δ2H and δ18O values of the ice via fractionation and mixing processes through the preferential incorporation of heavy isotopes into ice under non-equilibrium conditions, partial ice melting and freezing at different rates, and evaporation of melt water that enriched the residual water in heavy isotopes (Craig and Gordon, 1965; Souchez and Tison, 1987; Jouzel et al., 1999; Citterio et al., 2004; Clausen et al., 2007; and Ersek et al., 2018).

The Ivačičeva Cave ice samples plot below the Kredarica MWL and are described by a regression line slope <8. The regression line of the ice sample offset from the MWL is termed the freezing line (Jouzel and Souchez, 1982; Souchez and Jouzel, 1984). The freezing line is idealized for a “lake ice model,” where an isotopically uniform water body freezes from the top down. The initial frozen layer of surface ice is enriched in 18O and 2H because of the preferential incorporation of heavy isotopes into ice. The frozen surface creates a closed system in which the underlying water freezes from the top, and the resulting ice is increasingly depleted due to the continual incorporation of heavy isotopes into newly created ice (Citterio et al., 2004). Therefore, the freezing line can be extrapolated to identify the isotopic composition of the initial water. This “lake ice model” is not ideal for Ivačičeva Cave wall ice because the system was not sealed off by an ice boundary layer. However, this same fractionation process has also been observed for open systems (Jouzel and Souchez, 1982; Souchez and Jouzel, 1984), where ice formed by thin layers of water freezing can be described by an isotopic freezing line (Souchez et al., 2000). This idea of open-system fractionation is supported by the lack of an inverse d-excess–δ2H (or δ18O) relationship, which would arise under equilibrium conditions while freezing (Supplemental Material Figure S1). The absence of the inverse d-excess–δ2H relationship indicates that incomplete freezing of the water deposited and/or rapid freezing due to quick temperature changes resulted in non-equilibrium processes within the cave (Souchez et al., 2000; Perşoiu et al., 2011).

The Kredarica precipitation weighted mean isotopic composition plots on the Ivačičeva Cave ice data regression line (Supplemental Material Figure S2), which implies that it represents the initial source of ice in Ivačičeva Cave. In addition to the subsequent freezing of thin water layers along the ice wall, it is possible that the regression line describing the Ivačičeva Cave ice data may also be caused by ice melt and the subsequent evaporation of lighter isotope water prior to re-freezing. Seasonal temperature fluctuations can cause melting and mixing of layers, where the meltwater undergoes an evaporative influence. Mean annual air temperature at the Mt. Kredarica meteorological station (2,514 m a.s.l.) above the Ivačičeva Cave is −0.7°C (1991–2020; ARSO, 2023). While the assumed average annual temperature is below freezing in the portion of the cave covered in ice, the maximum temperature measured in the cave where ice is located is >0°C, and there is a seasonal freeze-thaw cycle that occurs in the caves (Novak and Kustor, 1983). Therefore, the preferential evaporation of 1H and 16O from meltwater likely causes further enrichment of the residual water, which is eventually frozen again due to seasonal temperature fluctuations in the cave (Clausen et al., 2007).

Snežna Cave δ2H and δ18O

Seasonal freeze-thaw events likely altered the stable isotopic compositions in Snežna Cave (Mihevc, 2018). Papi and Pipan (2011) reported freeze-thaw activity at the site in Snežna Cave where our samples were collected, and elevated temperatures may have caused ice to melt and evaporation of standing water. The mean annual air temperature at the Krvavec meteorological station (1,740 m a.s.l.), located approximately 20 km SW from the Snežna Cave, is 3.8°C (1991–2020; ARSO, 2023). The average annual temperature of the Entrance Chamber in Snežna Cave is −5°C (±2°C), whereas in other parts of the cave, the average annual temperature is ∼4.5°C (Mihevc, 2018). There can also be interaction between ice and meltwater when the ice is close to the melting point, further altering the initial stable isotopic composition (Clausen et al., 2007).

The Snežna Cave samples plot between the Kredarica and Rateče and Zgornja Radovna MWLs (Figure 4). The Snežna Cave ice is likely not described by the MWLs presented in Figure 4. Work reported in the Adriatic-Pannonian region identified the altitude effect as the main driver of precipitation isotopic depletion of −1.2‰/km for δ18O and −7.9‰/km for δ2H, which is essentially what we observed from differences in Kredarica and Rateče and Zgornja Radovna weighted precipitation means (Kern et al., 2020). Situated at an elevation between these collection sites, and in a different geographic location (Figure 1), the combination of continental and altitude effects on precipitation in Snežna Cave is not captured in the Kredarica and Rateče and Zgornja Radovna MWLs.

We observed that Snežna Cave floor ice values were more isotopically depleted compared to the wall ice in Ivačičeva Cave (Figure 3). While the Snežna Cave samples were floor ice samples, the lack of depletion with depth indicates that the “lake ice model” closed-system freezing process, discussed above, is also not applicable to describe these data. This is further supported by the lack of a d-excess–δ2H relationship in these samples (Souchez et al., 2000; Perşoiu et al., 2011). We explain the observed depletion of Snežna Cave floor ice in two potential ways: (1) the cold air trap nature of Snežna Cave, where air descends into the cave from a single entrance, compared to continual wind flow through Ivačičeva Cave, which is reported at a quarter of the external wind velocity (Novak and Kustor, 1983; Papi and Pipan, 2011; and Perșoiu, 2017), and (2) an enhanced seasonal bias through an ice chimney fed by vadose zone drip water (Mihevc, 2018). Ice formation in Snežna Cave may be biased towards early spring, when isotopically depleted snowmelt increases water percolation through the vadose zone, and this drip water freezes within the cave (Perşoiu et al., 2021). An enhanced seasonal bias is also likely at Snežna Cave because it is below the tree line, which can result in a greater proportion of water percolating through the vadose zone during early spring months, when evapotranspiration is low, and ice and snow are melting (Jasechko et al., 2014). Kern et al. (2011a) also observed a seasonal bias of ice formation due to snowmelt in an Austrian cave, and because of fractionation processes during freezing, they found samples that were enriched in 2H and 18O compared to average local precipitation. This has been observed elsewhere; for example, Fórizs et al. (2004) observed δ18O values in the uppermost ice in Focul Viu Cave, Romania, that were more positive than the annual mean of local precipitation due to partial melting and re-freezing of ice. It is likely that we observed enrichment fractionation of residual ice for both δ2H and δ18O from depleted cool-month (October–March) precipitation. The monthly weighted mean isotopic composition of cool-month precipitation at Kredarica are lower than −75‰ and −10‰, for δ2H and δ18O (Vreča et al., 2022).

Other work has shown that partial freezing in the vadose zone, where the heavier isotopes are concentrated in ice in the surface, results in isotopically depleted drip water delivery to karst ice caves (Ersek et al., 2018). Enriched ice near the surface will melt and contribute to water flow into the cave in the early spring (ARSO, 2023). However, Ersek et al. (2018) observed that the water availability from partial melting is low, and so this is likely only a small contribution to ice formation in the cave. This water may mix and affect the original primary isotopic signature (Ersek et al., 2018). Regardless, the primary source of this water would be from depleted winter precipitation. These processes could explain the 2H and 18O depletion that we observed in Snežna Cave ice core profiles compared to Ivačičeva Cave, where the geomorphologic setting of Ivačičeva Cave and consequent air flow may also cause a greater range of seasonal isotopic signatures there.

We clearly observed that neither Snežna Cave nor Ivačičeva Cave exhibited a consistent pattern of enrichment or depletion with depth in floor or wall ice. Previous studies by Perşoiu et al. (2011, 2017), and Ersek et al. (2018) have suggested that the isotopic composition of cave ice should reflect a dampened local meteoric water signature that fluctuates both over the short term (seasonally) and long term (decadal temperature changes). The isotopic composition of the two Snežna Cave core profiles had the same δ2H and δ18O values at the ice surface (top of profile), indicating the same source of water. However, the profile patterns differed until approximately 12 cm depth, where they displayed similar isotopic compositions until 22 cm depth (Figure 5). Although these cores were collected in close proximity (∼5 m apart), the data suggest differences in freezing-thaw processes on a small horizontal scale, which could also be influenced by the influx of water from ice masses thawing on the entrance slope. Differences in isotopic composition in the upper portion of the two Snežna Cave ice core profiles imply that fractionation and mixing processes caused by differences in the freezing rate, re-freezing, meltwater enrichment, and mixing of ice layers during formation played a role in the observed ice isotopic compositions.

Cave Ice δ2H and δ18O in Context

The age of the Snežna Cave and Ivačičeva Cave ice cores cannot be constrained by the isotopic composition of the ice due to mixing and fractionation processes. However, the ice data in both caves can be contextualized for the region and climatic variations. The cave ice δ18O results from this study are compared in Figure 6 with meteoric precipitation, glacier ice, snow, and ice δ18O values from Gabrielli et al. (2010), Vreča and Malenšek (2016),Carey et al. (2020), and Vreča et al. (2022) for Slovenia. Locations near the coast receive the most enriched 18O precipitation, and the average δ18O values for each location decrease with increasing elevation (Kern et al., 2020). Analyses from the Ivačičeva Cave and ice and meltwater near Triglav Glacier (Figure 1) were published in Carey et al. (2020). These previously published values closely resemble the average Ivačičeva Cave ice profile δ18O values in this study, but a difference in mean values between Ivačičeva Cave δ18O data from this study and data presented by Carey et al. (2020) is likely due to the different sampling resolution on the same ice wall. In this study, 40 samples were taken over 2 m, and Carey et al. (2020) presented six samples from the same 2 m (Carey et al., 2020). Overall, Ivačičeva Cave samples from both studies differ from the overarching pattern of decreasing 18O/16O ratios with elevation (Figure 6). The Snežna Cave ice core δ18O results exhibit a better fit with the pattern of δ18O depletion with increased elevation. The rainfall sample at a similar elevation (Zavižan, Croatia, 1,594 m) is near the coast of the Adriatic Sea, which likely is the cause for 18O-enriched precipitation there (Vreča and Malenšek, 2016).

Ice Cave Geochemistry

Although caves associated with permanent ice are found throughout Central Europe (Perșoiu and Onac, 2019), with more than 260 in Slovenia alone (Speleological Association of Slovenia, 2022), there are few published studies on the chemistry of the ice, with a much greater focus on ice δ18O and δ2H data as climate proxies. The work that does exist suggests that the primary inputs of ions to cave ice are derived from three main sources: (1) the original precipitation, (2) water that has percolated through the vadose zone, and (3) water derived from the karst system (Citterio et al., 2004; Clausen et al., 2007; Kern et al., 2011b; and Carey et al., 2019). The cation concentrations in ice from both the Snežna Cave and Ivačičeva Cave varied as follows: Ca ≫ Mg > Na > K, strongly suggesting that the primary source of Ca2+ and probably Mg2+ was from the dissolution of CaCO3 within the karst system or carbonate bedrock itself, as previously described by Citterio et al. (2004), Clausen et al. (2007), and Kern et al. (2011b) for caves in Italy, Slovakia, and Croatia, respectively. This supposition is supported by our data from these Slovenian caves on a Ca +Mg versus HCO3 + SO4 diagram, where data plot on a 1:1 equivalence line, indicating the ice chemistry in these two caves is dominated by the dissolution of the karstified limestone within the cave system (Figure 7). The Mt. Ortles snow pit mean Ca2+ concentration is at least 20 times lower than cave ice mean values, and the mean Mt. Ortles snow pit Mg2+ concentration is 2.5–5 times lower than the cave ice concentrations (Table 2; Gabrielli et al., 2010; Carey et al., 2019). Na+ and K+ have been determined in previous studies to be derived from both primary precipitation and soil water (Kern et al., 2011b). The very low concentrations of H4SiO4 in our cave ice samples suggest that the primary source of these cations in our samples is from precipitation. However, in our previous works, we have suggested that at least a portion of the K+ could be derived from organic matter decomposition (Carey et al., 2019).

Citterio et al. (2004) measured what they considered were lower than anticipated NO3 and NH4+ concentrations in cave ice and speculated that the loss of fixed nitrogen may have occurred prior to the freezing of the ice within the cave environment. Clausen et al. (2007) argued that the major source of NO3 in their cave system was derived from precipitation, while NH4+ was derived from bat droppings. Papi and Pipan (2011) reported NO3 concentrations of 24–50 μM in epikarst drip water in Snežna Cave. Concentrations at five drip sites in different locations within Snežna Cave had a similar range. This indicates that the source of NO3 + NO2, where NO3 is assumed to be the dominant species due to oxic conditions, is from epikarst drip waters that have been leached from within the vadose zone, as also suggested by Kern and Perşoiu (2013) and Carey et al. (2019). The elevated NO3 concentrations measured in drip water combined with diverse (N = 10) and abundant (N = 185) invertebrate taxa measured by Papi and Pipan (2011) support the hypothesis proposed by Citterio et al. (2004) that fixed nitrogen could be lost prior to freezing. The NH4+ concentrations measured in the Snežna Cave ice cores (<1–4.6 μM) had a greater range than the Ivačičeva Cave wall ice (<1–1.2 μM). All Snežna Cave samples contained < 2 μM NH4+, except for samples at 10–12 cm in core A and 20–24 cm in core B (Supplemental Material Table S1). The three samples with higher NH4+ concentrations may have been derived from biological material and mineralization of organic matter within the vadose zone (Clausen et al., 2007).

Differences in Ivačičeva Cave and Snežna Cave Ice Chemistry

Geochemical data from our previous work in both Snežna Cave and Ivačičeva Cave are compared to our newest data in Table 2. We cannot determine the age of the ice in either cave from previous observations or observations in this study. Our original work in Snežna Cave (Carey et al., 2019) was done by sampling a wall of ice close to the major entrance of the cave that had been formed by a large doline collapse into a large subsurface gallery. The data presented herein are from ice cores obtained by auguring into the massive floor ice of Snežna Cave in the large collapsed region, which has a mean annual temperature of −5°C (Mihevc, 2018). Both locations are within close proximity and ∼75 m from the entrance (Mihevc, 2018). The mean Cl and Na+ concentrations are similar in both sample sets (Table 2). If these ions are considered as conservative tracers, then the concentration differences in other ions may suggest variation in dissolution of the karst system (Ca2+), organic material decomposition in the vadose zone (K+), and/or variation in atmospheric aerosols (SO42−). The Ivačičeva Cave wall ice samples have lower mean Cl, SO42−, Na+, Mg2+, and K+ concentrations than the lower elevation Snežna Cave ice (Table 2). In addition, the Ivačičeva Cave ice is more similar in composition to higher elevation snow and ice from the Italian Alps, as previously noted in Carey et al. (2020).

Even though Ivačičeva Cave ice samples have Cl, SO42−, Na+, and K+ concentrations that are similar to high-elevation precipitation (snow from Mt. Ortles; Table 2), Ivačičeva Cave δ18O signatures resemble those of precipitation from lower elevations (Figure 6). Figure 8 shows the similarity in Mt. Ortles snow pit and Ivačičeva Cave ice ion chemistry, using Cl as a conservative tracer, but it also shows a difference in δ18O signature. The isotopically enriched δ18O values of the Ivačičeva Cave wall ice and most Snežna Cave floor ice samples in comparison to the Mt. Ortles average snow pit may be explained by a mixture of snow and rainfall contributing to epikarst drip waters at lower elevations, which results in 2H- and 18O-enriched signatures when compared to the average at Mt. Ortles. The 3‰ standard deviation for Mt. Ortles δ18O values is caused by the seasonality captured in the ice core, particularly from warm years that caused ablation and melt percolation (Gabrielli et al., 2010). Another explanation may be incomplete freezing of epikarst drip waters where non-equilibrium fractionation resulted in the preferential incorporation of heavier isotopes into ice.

Snežna Cave ice had greater ion concentrations compared to Ivačičeva Cave ice and Mt. Ortles snow pit (Table 2 and Figure 8), but the mean Snežna Cave δ18O signatures plot between the two locations (Figure 8). These differences in ice chemistry and isotopic composition between Snežna Cave and Ivačičeva Cave drive the significant inverse relationships between stable water isotopes and ions that we observed in the entire data set (Table 1). However, the individual cave δ2H and δ18O relationships with ions are not as robust (Supplemental Material Tables S2 and S3). Figure 8 shows an example of this visually. The apparent correlation coefficient (r) for δ18O and Cl in the entire data set is −0.82 (Table 1); however, Ivačičeva Cave ice samples have an insignificant δ18O and Cl relationship (r = −0.22; Supplemental Material Table S3), and Snežna Cave samples demonstrate a significant inverse relationship for δ18O and Cl (r = −0.34; Supplemental Material Table S2). Snežna Cave ice δ2H and δ18O values also have significant inverse relationships with H4SiO4 (Supplemental Material Table S2). This may be an indication of changing H4SiO4 concentrations in precipitation by season or year; however, this relationship seems to be driven by a few samples. Ivačičeva Cave ice δ2H and δ18O values have significant inverse relationships with Ca2+ concentrations, and δ18O values are inversely related with K+ concentrations (Supplemental Material Table S3). These significant relationships, where depleted stable isotopic values correlate with higher ion concentrations, suggest cold-season precipitation had an elongated residence time in the vadose zone and epikarst.

However, Ivačičeva Cave sits above the tree line with scarce patches of soil and grasses, so the soil depth and vegetation density are low relative to the location of Snežna Cave, which is below the tree line. Infiltrating waters derived from local meteoric precipitation likely have an overall shorter residence time in the shallow vadose zone above the tree line. This may limit weathering and the significance of vegetation decomposition, which results in lower ion concentrations in epikarst drip waters in Ivačičeva Cave compared to Snežna Cave, except for Ca2+ and HCO3. Ice samples in Snežna Cave had higher ion concentrations, which is indicative of longer water residence times and weathering/dissolution within the vadose zone and supports our hypothesis.

Given the strong chemical signal of karst water input (i.e., high concentrations of Ca2+ and HCO3), it is clear that these ice samples include the major components derived from CaCO3 dissolution. The complex topography at these locations in a karst landscape has previously allowed the preservation of snow and ice through the warmer portions of the year (Securo et al., 2022). However, these glaciokarst environments in the Dinaric and southeastern Alps, with their collapsed dolines, “schachtdolines,” and shafts, have prompted high rates of chemical weathering of the underlying limestones and dolomites (Securo et al., 2022). Clearly, primary processes such as direct precipitation input as well as hydrologically associated processes related to infiltration and water transport within the triple-porosity karst system control the chemical character of cave ice. Secondary processes such as melting, mixing, and re-freezing also play a role, depending upon the elevation and aspect of the cave location (Citterio et al., 2004). In Slovenia, most ice caves are located where the mean annual temperature is above 0°C (Mihevc, 2018). Because of this, conditions within caves have to be favorable for the ice to survive (Mihevc, 2018). Thus, the freezing of epikarst waters must play a role in permanent ice formation and therefore influences ice chemistry.

Cave Ice as a Paleo-Climate Proxy

Previous work on the chemistry of cave ice in Central Europe has explored the use of ice cave profiles as a potential proxy for understanding paleo-climate variations and past air-pollution history (Clausen et al., 2007; Kern et al., 2011b). In order to do this, a knowledge of the depth-age relationship is needed, and the assumption is made that the ice profiles represent unperturbed primary precipitation, and that melting and mixing of water within the profile have not occurred to a great extent. The horizontal variations in the Snežna Cave ice cores, especially between the Cl and SO42− concentrations at the same depths, suggest either a different age-depth relationship between the two cores or a potential pattern of differential melting and re-freezing within the floor ice at Snežna Cave, paralleling the interpretations made from stable water isotope data. The age of these cores is not known, and based on the geochemical data, it is probable that these ice profiles do not represent a true historic record of the chemical composition of precipitation. It is highly likely that layers of different precipitation ages have been melted, mixed together, and re-frozen to produce chemical and isotopic profiles that integrate periods of time. The depth to which this may occur is unknown at this time. The difference in fast inflow and slow inflow water delivery to the cave can also affect the chemical composition of ice. Previous work in a Romanian ice cave has demonstrated that annual ice accumulation can vary from 0.84 cm over the long term to 2 cm over the short term (Fórizs et al., 2004). Yet, many of these ice caves in Central Europe have rather thick deposits of ice associated with them, and thus many contain longer records that could reflect pre–late 20th century warming (Perșoiu, 2018).

The chemistry of ice formations of two caves in Slovenia reflects their precipitation source waters that have interacted with epikarst and have been altered by a set of complex processes (e.g., melting, re-freezing, mixing, sublimation, evaporation) that probably varied temporally. The δ2H and δ18O values of Snežna Cave and Ivačičeva Cave ice samples plot above the Slovenia RMWL and GMWL in isospace, but within the range of MWLs from higher elevations (Vreča et al., 2022). The analysis of samples with MWLs in a δ2H-δ18O biplot shows that precipitation from lower elevations is not appropriate in relating cave ice to source water and deducing the environmental mechanisms that cause isotopic fractionation (i.e., elevation, temperature). The Ivačičeva Cave ice regression line demonstrates the relationship between ice and weighted mean source water isotopic composition. This relationship is likely derived from (1) fractionation caused by freezing of thin layers of ice in an open system where ice is formed in cool months, (2) seasonal freeze-thaw activity that drives mixing and partial freezing, and (3) partial evaporation of meltwater, which is caused by fluctuating cave temperatures and relative humidities. Compounding changes in air temperature and hydroclimate are needed to extract patterns of ice accumulation and melt in ice caves in the region (Racine et al., 2022b). Localized differences in freeze-thaw activity were observed on a small (∼5 m) spatial scale, as ice core profiles in Snežna Cave had the same isotopic composition at the top of the core but did not follow the same isotopic pattern with depth.

Ice chemistry data support meltwater or precipitation mixing based on stable isotopic character, in that Cl and SO42− concentrations are also inconsistent at similar depths in Snežna Cave core A and core B. Our data show that the chemical profiles of ice in caves differ between caves and within a cave, despite similar (or identical) source waters. Ice chemistry is dominated by high Ca2+ and HCO3 concentrations from carbonate material dissolution in the epikarst domain above both caves. Variations in concentrations of NO3 + NO2 suggested biological utilization of nitrogen prior to freezing, and Na+ and K+ in ice were derived from vadose zone input. While we cannot determine the age of the ice, we do conclude that the ice reflects a precipitation signature that has been modified isotopically due to freeze-thaw processes that induce fractionation, and chemically from interaction with the vadose zone and epikarst before entering the cave. Regardless of these modifications in Snežna Cave and Ivačičeva Cave, previous studies have demonstrated that the ice in caves can be used to reconstruct climate (Fórizs et al., 2004; Perşoiu et al., 2011, 2017). However, as noted above, cryosphere loss in Central Europe, and particularly in Slovenia, is rapidly occurring, and it has been predicted that increased warming through the 21st century will lead to the demise of these features over the next decades (Čekada et al., 2016; Zorn et al., 2020; and Racine et al., 2022a).

Ice caves can be important repositories of historical climate information and atmospheric chemistry despite ice formation and freeze-thaw processes that can alter the initial chemistry. Furthermore, cave ice masses are currently being lost due to warming conditions. Warming may also increase the degree of diagenesis, further altering the original chemistry recorded during ice formation. To extract historical climate data from these ice caves, the ice formation and freeze-thaw patterns must be examined at a small scale. Further chemical research and age dating of cave ice cores and ice wall profiles in the southeastern Alps will enhance the utility of these caves as climatic proxies.

We would like to thank The Ohio State University’s Slovene Research Initiative for partial support to A.E.C. and W.B.L. to visit Slovenia and collect samples. We are extremely grateful to Professor Lonnie Thompson at Ohio State University for the use of his SIPRE auger. Additional support for fieldwork came from research program no. P6-0101, research project no. J6-3141 and J6-50214, and bilateral project BI-US/22-24-023 of the Slovenian Research and Innovation Agency. Support for laboratory instrumentation used for geochemical analyses came from National Science Foundation Directorate for Geosciences (GEO) grants EAR-0744166, EAR-0930016, and EAR-1342632. We would also like to thank Russell Harmon for the suggestion that we to participate in this special issue.

We acknowledge that The Ohio State University has long served as a site of meeting and exchange of Indigenous peoples, including the Shawnee, Miami, Wyandot, and Delaware Nations. We honor and respect the diverse Indigenous peoples connected to the territory where data analysis and writing of this work took place.

Supplemental Material associated with this article can be found online at https://doi.org/10.2113/EEG-D-23-00001 and https://www.aegweb.org/e-eg-supplements.