Purpose and Scope

The purposes of this report are to describe data collection and analysis methods and to infer hydrogeologic characteristics of Hourglass and New Years Lakes within and near Jewel Cave National Monument, southwestern South Dakota. The approach for this study included comparing water-level data collected from Hourglass and New Years Lakes to historical climate data and applying multivariate statistical analyses to evaluate water-chemistry data collected during this study and from previous investigations. Discrete water-level measurements and continuous water-level data from March 2018 to April 2021 and October 2018 to April 2021 at Hourglass and New Years Lakes, respectively, were used in this study. Water-chemistry data from the southern Black Hills and Wind Cave National Monument (fig. 1) were included to build on findings from previous investigations and evaluate Jewel Cave within the larger context that includes both cave systems. Observations from this study were compared to findings from previous investigations and hydrogeologic models to infer hydrogeologic characteristics of Hourglass and New Years Lakes.

Previous Investigations

Previous investigations have described the geologic setting, stratigraphy, and lithologic composition of geologic units in the Black Hills area (Darton and Paige, 1925; Lisenbee and DeWitt, 1993) and have mapped the distribution of geologic structures and geologic units (fig. 1; Martin and others, 2004; Redden and DeWitt, 2008). Investigations from Carter and others (2001) and Carter and others (2003) determined the availability of groundwater in the Black Hills and estimated recharge and hydraulic properties for many of the major (regional) and minor aquifers. Dyer (1961), Deal (1962), Braddock (1963), Fagnan (2009), and Wiles (2013) have discussed the origin of Jewel Cave and the distribution of geologic units within and near Jewel Cave National Monument.

Recharge mechanisms and hydraulic properties of the Madison Limestone and Minnelusa Formation within and near Jewel Cave were investigated by Alexander and others (1989), M.E. Wiles (South Dakota School of Mines and Technology, unpub. data, 1992), and Driscoll and others (2002). Alexander and others (1989) examined recharge to Jewel Cave by conducting dye tracer tests (Wilson and others, 1986). Dye tracer tests consisted of pouring rhodamine water tracer dye in parking lots and septic systems in Jewel Cave National Monument and sampling drip sites within the cave for the presence of the dye. Dye was first recovered about 8 days after insertion, and recovery continued intermittently with variable dye concentrations throughout the yearlong study. Cave lakes were not yet discovered in Jewel Cave when dye tracer tests were conducted by Alexander and others (1989) and the drip sites where dye was recovered were at a higher elevation in the cave than at the cave lakes. M.E. Wiles (South Dakota School of Mines and Technology, unpub. data, 1992) observed that drip sites in Jewel Cave were located where an impermeable unit in the Minnelusa Formation is absent, indicating that groundwater in the Minnelusa aquifer can drain into the underlying Madison aquifer and has sufficient storage capacity to supply drip sites year-round. Driscoll and others (2002) determined the Madison and Minnelusa aquifers in the Black Hills receive appreciable recharge from streamflow losses and precipitation on the outcrops.

Other previous investigations pertaining to the scope of this study were conducted by Long and Valder (2011), Long and others (2012), Anderson and others (2019), and Long and others (2019). Long and Valder (2011) and Long and others (2012, 2019) used multivariate statistical analyses—including principal component analysis (PCA) and cluster analysis—to examine water-chemistry data collected from streams, springs, wells, and cave lakes. Water-quality data collected in various studies established baseline water quality for many sites, which included physical properties (pH, temperature, and specific conductance), major ions, arsenic, strontium, uranium, stable isotopes of oxygen and hydrogen, and radiogenic isotopes of strontium and uranium. The purpose of water-quality data collection was to better understand the newly discovered cave lakes in Jewel Cave and Wind Cave by comparing the cave-lake water to different sources of water in the Black Hills. Five hydrogeologic domains were identified by Long and others (2019) using PCA and cluster analysis and grouped based on similar water-quality characteristics. Samples from cave lakes in Jewel Cave grouped with samples from nearby Madison aquifer wells, which, combined with similar elevations between cave lakes in Jewel Cave and groundwater in nearby observation wells completed in the Madison aquifer, led Long and others (2019) to conclude that cave lakes in Jewel Cave, as well as Wind Cave, were connected to regional groundwater flow in the Madison aquifer.

Anderson and others (2019) prepared a generalized potentiometric map of the Madison aquifer near Jewel Cave National Monument based on water-level data from 24 observation wells completed in the Madison aquifer and 4 cave lakes. Previous investigations focused on preparing potentiometric surfaces of the Madison aquifer in the Black Hills (Strobel and others, 2000; McKaskey, 2013) did not use water-level data from cave lakes because these lakes had not yet been discovered. Before discovery of cave lakes in Jewel Cave, previous investigations had speculated groundwater in the cave was 50 ft below most of the known cave (NPS, 1994). Anderson and others (2019) also evaluated historical and current groundwater recharge to the Madison aquifer near Jewel Cave based on hydrographs of wells and water-level data from Hourglass Lake. While evaluating hydrograph fluctuations, Anderson and others (2019) documented statistical correlation between water levels at Hourglass Lake and cumulative daily precipitation from a nearby climate station. Correlation between water levels at Hourglass Lake and cumulative daily precipitation was attributed to recharge from precipitation on Madison Limestone outcrops. Another observation was that water levels at Hourglass Lake began increasing before cumulative daily precipitation, which was likely from early spring snowmelt contributing to recharge.

Water-Level Data Collection

Water-level data compiled during this study included discrete water levels measured by NPS staff and continuous water-level data measured with pressure transducers. The NPS provided the USGS with discrete water-level elevation measurements for 2 cave lakes that were collected as part of routine trips into the cave (USGS, 2022a). Water-level elevation measurements by NPS staff are performed using in-cave surveys starting from a benchmark of known elevation outside the cave. Discrete water-level measurements have been collected 11 times at Hourglass Lake by NPS staff since 2015, and 4 times at New Years Lake since 2017. Pressure transducers were installed in Hourglass Lake (HGL; USGS site 434258103504201; fig. 1) on March 11, 2018, and in New Years Lake (NYL; USGS site 434238103503501) on October 5, 2018. Water-level data from March 2018 to April 2021 were uploaded to the NWIS database (USGS, 2022a) and will be focus of analysis presented in this report.

Pressure transducers used in this study were unvented Solinst Levelogger LT F15/M5 Model 3001 electronic transducers with a manufacturer accuracy of plus or minus (±) 3 cm (0.01 ft; https://www.solinst.com/). Pressure transducers were unvented, which means that they required separate barometers to convert pressure data measured by the transducer to water-level data by compensating for atmospheric pressure. The barometers were Solinst Barologger LT F15/M5 Model 3001 electronic barometers with a manufacturer accuracy of ± 0.05 kilopascal (kPa; 0.007 pound per square inch; psi; https://www.solinst.com/). Both the pressure transducers and barometers were programmed to record continuously at 1-hour intervals at the start of each hour. NPS staff suspended the pressure transducers from cave ceilings using fishing line and submersed the transducers from 1 to 8 ft below the water surface. Barometers were placed or suspended above the water surface to prevent water damage from rising water levels. Water-level data from pressure transducers were corrected for atmospheric pressure recorded by barologgers using Levelogger Software version 4.6.1 (https://www.solinst.com/).

Continuous water-level data collected at Hourglass and New Years Lakes from 2018 to 2021 were measured using an hourly sampling rate from two transducers. Hourly pressure data measured with the transducers were converted to a height of water above the transducer by compensating for atmospheric pressure recorded by a nearby barometer. The hourly water heights were converted to relative water-level changes by subtracting each hourly measurement from the first measurement when the transducer was installed. The relative hourly water-level changes were added to the water-level elevation recorded by NPS staff when the transducer was installed to calculate water-level elevation for each hourly measurement. The accuracy of the calculated water-level elevations was evaluated by calculating the difference between the calculated water-level elevation on April 16, 2021, and water-level elevation measured by NPS on the same day when the pressure transducer was retrieved. The difference between the calculated water-level elevations on April 16, 2021 (4,694.6 ft at Hourglass and 4,617.0 ft at New Years Lake) and the NPS measured water-level elevations (4,694.5 ft at Hourglass and 4,616.9 ft at New Years Lake) was 0.1 ft. The difference was assumed to be caused by instrument drift, and a time-proportional correction was applied to the recorded water-level differences so that the estimated water elevations were within the accuracy of the discrete measurements (0.1 ft).

Major Ion and Stable Isotope Sample Collection

Water-chemistry data analyzed and presented in this report included samples from previous investigations (Long and Valder, 2011; Long and others, 2012, 2019) and samples collected as part of this study. Water-chemistry data from previous investigations were downloaded from the NWIS database (USGS, 2022a) by selecting only sites with the same chemical constituents and physical property data as the samples collected during this study, including major ions and stable isotopes of oxygen (oxygen-18, δ¹⁸O) and hydrogen (deuterium, δ²H). The list of sites with the physical properties and chemical constituents used in analysis is presented in table 1.2 in appendix 1. Sites with multiple samples were condensed by computing the mean of all available samples in the NWIS database (USGS, 2022a).

Nine water samples were collected from April to August 2021 at 7 sites and analyzed for 12 constituents (table 1). Sites sampled included one stream, three springs, two cave lakes, and one rain collector. Water samples were collected as grab samples using methods described in U.S. Geological Survey (variously dated). Sampling procedures and the constituents analyzed varied by site. Physical properties of sampled water—dissolved oxygen, pH, specific conductance, water temperature, and turbidity—were measured by USGS personnel using a multiparameter sonde (Xylem EXO1) at spring and stream sites in the field before sample collection; however, these physical properties were not measured at rain collector and cave lake sites where NPS personnel collected samples. Water samples from spring, stream, and cave lake sites were analyzed for stable isotopes and major ions. Precipitation samples from the rain collector site were analyzed for stable isotopes and not major ions because the concentration of major ions in precipitation is commonly low and predictable in small geographic areas (Rinella and Miller, 1988). Rain collectors used for precipitation sampling were “ball-in-funnel” type (Michelsen and others, 2018), installed at one site (JECApcp; USGS site 434404103494001) by USGS personnel, and fixed about 2 ft above the land surface in an unobstructed area.

All water samples were analyzed for stable isotopes δ¹⁸O and δ²H. Stable isotope samples were collected using 60-milliliter (mL) glass bottles with Polyseal caps and sent for laboratory analysis at the USGS Reston Stable Isotope Laboratory in Reston, Virginia. The laboratory methods used to determine stable isotope ratios are described by Révész and Coplen (2008a, 2008b). Samples analyzed for major ions were collected using a 125- or 250-mL polyethylene bottle rinsed and filled with a sample passed through a 0.45-micrometer pore size filter. Major ion samples were analyzed by the USGS National Water Quality Laboratory in Lakewood, Colorado, using methods described by Fishman and Friedman (1989), Fishman (1993), and Garbarino and others (2006).

Hydrograph Analysis

Hydrograph analysis involved comparing water-level changes among Hourglass and New Years Lakes and by comparing water-level measurements at each cave lake to climate data. Hydrographs of Hourglass and New Years Lakes were compared using continuous water-level data from October 2018 to April 2021 by calculating the elevation difference between the cave lakes and by qualitatively comparing water-level change differences among the cave lakes. Additionally, annual water-level elevation range and annual water-level change were estimated at Hourglass and New Years Lakes from 2018 to 2021 when pressure transducers were installed. Annual water-level elevation range was calculated by calculating the difference between each year’s maximum water-level elevation and minimum water-level elevation. Annual water-level change was calculated by subtracting water-level elevation on January 1 from water-level elevation on December 31 of the same year. When water-level data were not available for an entire year at Hourglass and New Years Lakes, the earliest and latest available water-level data were used instead to calculate the annual water-level change. The range of available water-level data at Hourglass Lake was from March 11, 2018, to April 16, 2021. The range of available data at New Years Lake was from March 12, 2018, to November 21, 2018, and from April 7, 2019, to April 16, 2021.

Hydrographs of Hourglass and New Years Lakes also were compared to precipitation and snow depth data obtained from the Custer, S.D. (USC00392087) and Custer County Airport, S.D. US (USW00094032) National Oceanic and Atmospheric Administration (NOAA) climate stations (fig. 1; NOAA, 2021a, 2021b). Daily observations from the Custer County Airport, S.D. US (USW00094032) climate station were substituted for missing daily observation data at the Custer, S.D. (USC00392087) climate station to provide a more complete dataset. Daily precipitation observations from Custer, S.D. (USC00392087) were cumulatively summed for each calendar year from 2015 to 2021 and were compared to annual normal precipitation totals from 1991 to 2020 at the Custer, S.D. climate station. Climate normals are means of climatological measurements spanning three decades and include temperature, precipitation, snowfall, and other measurements (NOAA, 2022). Daily snow depth observations were obtained from the Custer, S.D. (USC00392087) and were used to evaluate the effect of winter and spring snowmelt on water levels.

A separate analysis involved obtaining monthly precipitation observations from 2015 to 2021 to evaluate how each month from 2015 to 2021 compared to monthly normal precipitation totals from 1991 to 2020 (NOAA, 2021b). Monthly precipitation observations were obtained from the Custer, S.D. climate station, except where missing data from the Custer, S.D. climate station had to be supplemented with data obtained from the Custer Co Airport, S.D. US (USW00094032) climate station (fig. 1; NOAA, 2021a). Monthly normal precipitation from 1991 to 2020 were obtained from the Custer, S.D. climate station (NOAA, 2021b). The difference between monthly observations and monthly normals (departure from normal value) and the cumulative monthly difference (cumulative departure from normal) are plotted in figure 2.

Stable Isotopic Analysis

Hydrogeological studies commonly use stable isotopes of abundant elements that are present naturally in the environment, such as hydrogen, carbon, or oxygen, to estimate water age, recharge processes, and groundwater-flow paths. Stable isotopes of hydrogen and oxygen are used as flow-path tracers because these isotopes are found naturally in groundwater, and meteoric weather processes can modify their composition (Clark and Fritz, 1997). Stable isotope analysis of a water sample compares ratios of heavier to lighter isotopes to a standard isotope ratio of known composition (Clark and Fritz, 1997). The result is reported in delta (δ) notation in parts per thousand (‰) calculated using equation 1 (Kendall and Caldwell, 1998):

A positive value of δ indicates the sampled isotope ratio is higher than the standard, and a negative value of δ indicates the sampled ratio is lower than the standard (Kendall and Caldwell, 1998). Stable isotope ratios of oxygen (¹⁸O/¹⁶O) and hydrogen (²H/¹H) were measured for water samples collected at one stream, three springs, two cave lakes, and one rain collector for this study. Stable isotope data were converted to δ notation using equation 1 and the Vienna Standard Mean Ocean Water and Standard Light Antarctic Precipitation standards (Révész and Coplen, 2008a, 2008b). The notations δ¹⁸O and δ²H are used in this report to describe the ratios of heavy to light isotopes of oxygen and hydrogen, respectively.

The δ¹⁸O and δ²H data from this study and data from previous investigations in and around Jewel Cave (“West area” sites in table 1.1; fig. 1) were plotted based on a method described by Muir and Coplen (1981). The global meteoric water line (GMWL) from Craig (1961) was plotted with the sample data to compare the stable isotopic composition of samples to the global stable isotopic composition of precipitation. Stable isotope ratios of local precipitation samples deviated from the GMWL and were plotted as the local meteoric water line (LMWL). Linear regression was used to determine a LMWL by relating δ²H to δ¹⁸O for precipitation samples. Stable isotope samples were plotted in stable isotope plots. No samples were enriched (positive δ²H and δ¹⁸O values) relative to the standards used to calculate stable isotopes ratios—indicated by all samples plotting below the origin on stable isotope plots (fig. 3); however, some samples plotted closer to the origin than other samples. Samples plotting closer to the origin in stable isotope plots are heavier (more enriched in heavy isotopes) than samples with more negative δ¹⁸O and δ²H values (lighter—more depleted in heavy isotopes) that plot further below the origin (fig. 3).

Ocean waters are considered to have heavy stable isotope compositions, and as precipitation that originates from ocean evaporation moves inland, it becomes lighter in isotopes as the heavier isotopes preferentially fall to the surface during rainfall (USGS, 2004). Generally, summer rains are heavier (more positive) than lighter winter rains (more negative). Similarly, precipitation from cooler, high latitude, high altitude, and inland sources is lighter than precipitation from warmer, low latitude, low altitude, coastal areas (University of Arizona, 2020). Shallow groundwater stable isotope values are similar to precipitation values, but evaporation, transpiration, and fractionation may alter the stable isotope ratios as the water moves downward toward the saturation zone (USGS, 2004). Additionally, evaporative loss from rivers, reservoirs, and lakes changes the isotopic composition of water, causing it to become heavier (University of Arizona, 2020).

Principal Component Analysis

PCA is a multivariate technique that tests for linear relations among variables in a dataset (Helsel and others, 2020). The PCA technique is performed by linearly transforming the original dataset into new axes, called principal components, that are linear combinations of all the input variables. The principal components responsible for the most variance—principal components one and two—are used to plot multidimensional data in two dimensions so data patterns, clusters, and groupings can be observed and interpreted (Long and Valder, 2011). Additional details and mathematical derivations of PCA are provided in Davis (2002). Python programming language (Rossum and Drake, 2011) was used to implement the PCA technique.

The PCA technique was performed using physical property and chemical constituent data collected for this study and from previous investigations (Long and Valder, 2011; Long and others, 2012, 2019). Data from previous investigations were downloaded from the USGS NWIS database (USGS, 2022a). Additional sites from other investigations assisted in the interpretation of relations among sites within a larger area in the southern Black Hills and improved PCA analysis by increasing the number of observations that have been found to produce a more stable result (MacCallum and others, 1999). Water-chemistry data from other studies were collected from streams, springs, wells, and cave lakes in the southern Black Hills near Jewel Cave National Monument and Wind Cave National Park (fig. 1). Aquifers for groundwater sites were designated in the USGS NWIS database and included the White River, Minnelusa, Madison, and undifferentiated Precambrian igneous and metamorphic aquifers. The naming convention of sites shared between this analysis and that of Long and others (2019) were kept the same for consistency. Sites sampled as part of this study not included (New Years Lake; Hell Canyon Creek) in the analysis by Long and others (2019) adhered to a similar naming convention given in tables 1 and 1.1. Additionally, similar to Long and others (2019), all sites east of the WIL site (south-central area of fig. 1) are referred to as the “East area,” whereas all sites west of the WIL site are referred to as the “West area.”

The number of variables used in PCA analysis was limited to physical property and chemical constituent data shared by all datasets (previously collected data and new data collected during this study) and included specific conductance, pH, carbon dioxide (CO₂), hardness (as CaCO₃), calcium (Ca), magnesium (Mg), sodium (Na), chloride (Cl), sulfate (SO₄), silica (Si), arsenic (As), δ²H, and δ¹⁸O. Sites excluding any of the aforementioned physical properties and chemical constituents were excluded from PCA analysis because the method does not allow for missing values. Physical property and chemical constituent data were normalized, by setting the mean of each variable to zero, and standardized, by setting the standard deviation of each variable to one, before performing PCA to ensure that the distribution of the dataset was independent of measurement units (Davis, 2002).

Water-chemistry data were compiled for 63 total sites in the study area. The mean of each potential variable was computed for sites with more than one sample to reduce the effect of repeat samples on the PCA results. Additionally, the water-chemistry data for the 63 sites were inspected for outliers and certain physical properties and chemical constituents that could affect PCA results. Sites were excluded from PCA if chloride and sulfate concentrations exceeded 100 and 350 milligrams per liter, respectively, to avoid including treated groundwater from wells and sites with elevated ion concentrations from road salts. A total of 8 sites were removed based on the specified criteria, including Prairie Dog Spring (PDGsp) and four other sites from Long and others (2012). Long and others (2019) noticed that PDGsp was elevated by road salts applied during winter. The sodium and percent sodium variables were removed because they explained less than 1 percent of variance in the dataset. In total, 55 observations and 12 variables were used in the PCA analysis (table 1.2).

Principal component biplots were created to visualize the results on principal component axes one and two. Biplots are PCA plots with vectors, called loading lines that relate variance and correlation of the variables. Variance is represented by the magnitude of the vector, with larger vectors indicating larger variance, whereas correlation is represented by the direction of the vector, with opposite direction vectors indicating negative correlation (Jolliffe, 2002). Individual data points that plot close to loading lines have above-mean values for that loading line variable; conversely, data points that plot opposite of loading lines have lower than normal values for that loading line variable (Jolliffe, 2002).

Cluster Analysis

Cluster analysis is a method of assigning data points to groups, called clusters, based on similarities (Long and Valder, 2011). The cluster analysis method selected for this study was the k-means procedure because of its simplicity and frequent use in hydrological studies (Long and Valder, 2011; Masoud, 2014; Marín Celestino and others, 2018). The k-means procedure involves iteratively finding the minimum Euclidean distance between a manually selected number of cluster centroids, k, and an observation, n (Davis, 2002). Iterations continue until the location of each k is optimized and every observation is assigned to a specific cluster.

In this study, the k-means procedure was used to determine clusters for the PCA results of water-chemistry data collected during this study and from Long and others (2012, 2019). The purpose of using the k-means procedure was to statistically group sampling sites without introducing statistical bias, such as site type or spatial location. The number of clusters for the k-means procedure was determined using scree plots (Jolliffe, 2002). Scree plots are developed by plotting the number of cluster groups in numerical order (x-axis) against their respective sum of squared distances of samples to the nearest cluster centroid (y-axis). The optimal number of clusters from the k-means procedure can be determined by observing the “break point” in scree plots. The “break point” is the sharp change in slope (steep to nearly flat) of the scree plot curve that marks where increasing the number of clusters is no longer beneficial to the k-means procedure. After selecting the number of clusters, the k-means procedure was computed using 10 runs at 300 iterations per run (a total of 3,000 iterations) for data plotted in principal component axes one and two. The k-means procedure settings were chosen because the cluster results did not improve with increased runs and iterations. The Python programming language (Rossum and Drake, 2011) was used to perform the k-means procedure in this study.

Hydrograph Comparisons and Annual Water-Level Observations of Hourglass and New Years Lakes

Hydrographs of Hourglass and New Years Lakes are shown in figure 4 and include discrete and continuous water-level data. The water-level elevation of Hourglass Lake was consistently higher than New Years Lake (fig. 4) as expected because Hourglass Lake is in a higher elevation part of Jewel Cave. The elevation difference between Hourglass and New Years Lakes ranged from 61 to 93.5 ft, with mean and median differences of 70 and 71 ft, respectively. The range of elevation difference (32.5 ft) and water-level change dissimilarities between hydrographs of Hourglass and New Years Lakes (fig. 4) were unexpected because the lakes are separated by about 0.4 mi at the land surface and are within the same geologic formation. Hourglass Lake displayed small, gradual water-level changes with elevations ranging from 4,684.3 to 4,695.2 ft (an elevation range of 10.9 ft) from October 2015 to April 2021. New Years Lake displayed relatively large and rapid water-level changes (compared to Hourglass Lake) with elevations ranging from 4,590.9 to 4,629.5 ft (an elevation range of 38.6 ft) from January 2017 to April 2021 (fig. 4). The larger elevation range at New Years Lake—combined with its relatively more dynamic water level increases and decreases—varied the elevation difference between New Years and Hourglass Lakes.

At Hourglass Lake, annual water-level elevation ranges were 4.1 and 4.8 ft in 2018 and 2019, respectively, and decreased to 2.7 and 0.6 ft in 2020 and 2021, respectively (table 2). New Years Lake displayed a similar trend where water-level elevation ranges were largest in 2018 (27.5 ft) and 2019 (12.8 ft) and smallest in 2020 (7.1 ft) and 2021 (2.8 ft; table 2); however, annual water-level elevation ranges at New Years Lake were between 2.6 and 6.7 times greater than Hourglass Lake (table 2). Annual water-level changes at Hourglass Lake showed the overall water level increased from 2018 to 2020 (table 2; fig. 4). In 2021, water levels decreased by 0.6 ft at Hourglass Lake (table 2). Annual water-level changes at New Years Lake were similar to Hourglass Lake where water levels increased in 2018 and 2019 and decreased in 2021, except water levels decreased in 2020 at New Years Lake (table 2). Additionally, annual water-level changes at New Years Lake were between 2.3 and 7.4 times greater than Hourglass Lake (table 2).

Hourglass Lake Hydrograph Compared to Precipitation and Snow Depth

Discrete water-level measurements at Hourglass Lake from mid October 2015 to late December 2016 were compared to precipitation and snow depth datasets to explain water-level changes (fig. 5). Discrete water-level measurements could not be interpreted on a monthly basis with precipitation and snow depth datasets because the measurements were separated by 1 to 22 months; however, some annual interpretations were made using the available data. Annual precipitation in 2015 (21.48 in.) was slightly greater than annual normal precipitation (20.66 in; fig. 5), and cumulative departure from normal, as shown in figure 2, indicates above normal precipitation conditions (wet) until April 2016. The water level of Hourglass Lake increased by 0.8 ft from mid October 2015 to late April 2016 during wet conditions. Annual precipitation in 2016 (17.13 in.) was about 3 in. below annual normal precipitation (fig. 5), and cumulative departure from normal, as shown in figure 2, indicates below normal precipitation conditions (dry) starting in April 2016 and continuing through December 2016. The water level of Hourglass Lake decreased by 0.5 ft from late April to late July 2016 during dry conditions; however, the water level increased by 0.9 ft from late July to late December 2016 while still under dry conditions. The water-level increase from July to December 2016 could result from snowmelt in December; daily snow depth data from December 2016 (fig. 5; NOAA, 2021a) indicated snow accumulations of 6.0 in. and 7.0 in. from two snowfall events that gradually melted in the 2-week period before the December 2016 water-level measurement—indicated by the decreasing snowpack during that time (fig. 5).

In 2017, discrete water-level measurements indicated water levels in Hourglass Lake were decreasing despite above normal annual precipitation (21.83 in.; fig. 5). Six months in 2017 were above normal with the largest departures from normal in July (3.39 in.), June (−1.60 in.), and May (−1.27 in.; fig. 2). Cumulative departure from normal remained negative for the early part of 2017, and the water level at Hourglass Lake correspondingly decreased during the same time. Decreasing water levels continued until March 2018 (fig. 5) while monthly precipitation was above normal for all but 1 month (January 2018) from November 2017 to March 2018 (fig. 2). Cumulative departure from normal (fig. 2) indicates the period of above normal precipitation from November 2017 to March 2018 only slightly increased the cumulative total because the departures from normal were less than 0.5 in. for all but 1 month (0.8 in. in February 2018; fig. 2).

After installing a pressure transducer in March 2018 in Hourglass Lake, the water level was relatively constant until mid-June 2018 when the water level increased by 4.1 ft and peaked in late September 2018 (fig. 5). Increasing water levels at Hourglass Lake likely resulted from a combination of spring snowmelt (fig. 5) followed by above normal spring and early summer precipitation in 2018 (fig. 2). Above normal precipitation fell in February and March 2018 as a mixture of snow and rain (figs. 3 and 5). The longest period of snowpack was from February 4, 2018, to March 8, 2018, where depths ranged from 1.0 in. to 12.0 in. (fig. 5; NOAA, 2021a). Various other snowfall events resulted in March and April 2018, but the snow melted completely within 1 week after each event (fig. 5; NOAA, 2021a). Precipitation in April 2018 was below normal but was followed by above normal precipitation for 3 consecutive months (fig. 2). Cumulative departure from normal shown in figure 2 increased from −1.12 in. below normal in April 2018 to 6.30 in. above normal in July 2018—meaning that from May to July 2018, 7.42 in. more precipitation fell than normal. Additionally, the largest increase in cumulative annual precipitation for 2018 (fig. 5) was from early April 2018 to August 2018; precipitation totaled 20.42 in. from April to August 2018—or about 75 percent of the annual total precipitation of 2018 (27.19 in.; NOAA, 2021a) and 99 percent of the annual normal precipitation (20.66 in.; NOAA, 2021b).

Water-level changes in Hourglass Lake from late September 2018 to late June 2019 were evaluated using precipitation and snow depth datasets. Dry conditions from August 2018 to November 2018—shown by decreasing cumulative departure from normal in figure 2—initially caused the 0.7 ft water-level decline in Hourglass Lake from late September 2018 to late March 2019 (fig. 5); however, water levels continued declining despite above normal precipitation from December 2018 to March 2019 (fig. 2). Continued decline of water levels until March 2019 correlated with cumulative departure from normal indicated in figure 2; the increase in cumulative departure from normal from November 2018 to March 2019 was smaller than the decrease in cumulative departure from normal from July to November 2018—indicating that the overall cumulative departure from normal from July 2018 to March 2019 was below normal. Additionally, precipitation falling in February and March 2018 may not have affected water levels in Hourglass Lake because low temperatures prevented snow melt from February 6, 2019, to March 17, 2019, with snow depths ranging from 1.0 to 11.0 in. (fig. 5). Increasing temperatures in March 2019 began melting snow accumulations from February 2019 and shortly after water levels in Hourglass Lake increased by 0.7 ft until mid-April 2019 (fig. 5). Water levels in Hourglass Lake remained relatively constant from mid-April to late June 2019 (fig. 5) as precipitation was 0.03 and 0.05 in. above normal in March and April 2019, respectively (fig. 2).

Water levels in Hourglass Lake increased by 6.9 ft from late June 2019 to early January 2021 (fig. 5) in response to above normal precipitation from May to December 2019 (fig. 2). Cumulative departure from normal indicated in figure 2 peaked in September 2019 after increasing by 6.2 in. from unusually high summer (May–September) precipitation. Cumulative departure from normal remained above normal and water levels in Hourglass Lake increased despite relatively dry conditions from October to December 2019 (figs. 3 and 5). The effect of dry conditions in fall and winter 2019 was observed when the rate of water-level increase slowed beginning in mid-January 2020 (fig. 5). Annual precipitation in 2020 (17.32 in.) was 3.3 in. lower than annual normal precipitation (20.66 in.); however, water levels in Hourglass Lake increased in 2020 at the slower rate of water-level increase until peaking in January 2021. The rate of water-level increase in 2020 varied slightly (fig. 5). These increases likely resulted from winter or spring snowmelt (increasing rate) and dry summer or fall conditions in 2020 (decreasing rate; fig. 2).

Water levels in Hourglass Lake decreased by 0.7 ft from early January 2021 to mid-April 2021 (fig. 5). Dry conditions from late 2020 continued in 2021 and likely caused water levels in Hourglass Lake to decline. Additionally, the water level at Hourglass Lake did not increase from spring snowmelt in March or April 2021, similar to the increase in previous years. Snow depth data from fall 2020 and winter 2021 show more temporally dispersed snowfall events and a smaller winter snowpack in February and March 2021 that melted completely by early March 2021 (fig. 5). As of November 2021, the annual precipitation of 2021 (13.57 in.) was 7.1 in. below annual normal precipitation. Additionally, cumulative departure from normal indicated in figure 2 decreased by 3.6 in. from July 2020 to February 2021.

New Years Lake Hydrograph Compared to Precipitation and Snow Depth

The New Years Lake hydrograph was compared to the same precipitation and snow depth datasets as the Hourglass Lake hydrograph described in the previous section to evaluate water-level changes (fig. 6). Water levels in New Years Lake declined by 5.0 ft from January 2017 to March 2018 despite above normal precipitation in 2017 (21.83 in.; fig. 6) and increasing cumulative departure from normal from January 2017 to March 2018 (fig. 2). The 5.0 ft water-level decline likely was caused by a deficit in precipitation, as indicated in figure 2, when cumulative departure from normal fell below 0 in April 2016 and remained negative throughout 2017 and early 2018. The precipitation deficit from January 1, 2016, to December 31, 2017, was calculated by first determining the annual departure from normal (difference between cumulative annual precipitation and annual normal precipitation) for 2016 and 2017, and second calculating the difference between annual departure from normal for 2016 (3.53 in. below normal) and 2017 (1.77 in. above normal; fig. 2). The precipitation deficit was 1.76 in. below normal from January 1, 2016, to December 31, 2017.

The water level of New Years Lake increased by 27.5 ft from March 2018 to October 2018 during a period of above normal precipitation (fig. 6). Annual precipitation in 2018 (27.19 in.) was 6.53 in. greater than annual normal precipitation (20.66 in.; fig. 2). Additionally, precipitation from May to July 2018 was 7.42 in. greater than normal (fig. 2). The water level of New Years Lake likely began increasing soon after spring snowmelt in March or April 2018 and continued increasing from May to July 2018 during above normal precipitation (fig. 2). A pressure transducer was installed in New Years Lake on October 7, 2018, shortly before a water-level elevation peak on October 14, 2018 (fig. 6). The largest water-level elevation peak in 2018 may have been during the above normal precipitation period from May to July 2018 rather than on October 14, 2018; however, no water-level data were collected between March 12, 2018, and October 6, 2018, to confirm a larger water-level peak before October 14, 2018. Water levels in New Years Lake declined after October 14, 2018, during a period of below normal precipitation from August to November 2018 (fig. 2) until the water level decreased below the level of the transducer on November 21, 2018, and no data were recorded (fig. 6). The transducer in New Years Lake began recording water levels again beginning on April 7, 2019.

Water levels in New Years Lake began increasing sometime before April 7, 2019, and the increase likely coincided with spring snowmelt. Departure from normal monthly precipitation (fig. 2) indicates that 0.79 in. more precipitation fell in February 2019 than normal in the form of snowfall that accumulated until mid-March 2019 (fig. 6). Snow depth data (fig. 6) indicated a relatively large snowpack starting in February 2019 and fully melted by mid-March 2019. The water level of New Years Lake peaked in early May 2019, and soon after the water level decreased by 2.5 ft until early July 2019 (fig. 6), despite above normal precipitation falling in the area from March to June 2019 (fig. 2). Departure from normal was less than 0.05 in. for March and April 2019 and as a result cumulative departure from normal showed minimal increase from March and April 2019 (fig. 2). Water-level decline at New Years Lake from May to July could result from the return to normal precipitation conditions following above normal conditions; water released by melting of above normal winter snowfall falling over a relatively short time (1–2 months) compared to sporadic near normal precipitation in March and April 2019 (fig. 6).

The water level of New Years Lake fluctuated in response to climate-related events from July 2019 to August 2020. The water level of New Years Lake increased by 11.3 ft from early July 2019 to early October 2019 (fig. 6). Annual precipitation in 2019 was 10.73 in. above annual normal precipitation and cumulative departure from normal (fig. 2) increased from 5.91 in. above normal in April 2019 to 12.11 in. above normal in September 2018, indicating that 6.2 in. more precipitation fell in May to September 2018 than normal. The water level of New Years Lake increased in response to above normal summer precipitation in 2019 until peaking in early October 2019 (fig. 6). The water-level peak in early October 2019 was followed by a water-level decline that persisted until early July 2020 (fig. 6); however, the rate of water-level decline fluctuated in response to precipitation. The first noticeable change in the rate of water-level decline was in January 2020 when the rate of decline increased after 3 months of near or below normal precipitation (fig. 2). The next noticeable water-level change was in March 2020 when the rate of water-level decline decreased—likely in response to spring snowmelt in February and March 2020—and continued decreasing until early July 2020 (fig. 6). The water level of New Years Lake increased by 0.7 ft from early July to late August 2020 (fig. 6) after above normal precipitation in June 2020 and near normal precipitation in July 2020 (fig. 2). Additionally, the start of water-level increase in early July 2020 was less than 1 week after the last week of June 2020 where a total of 4.2 in. of rainfall fell (fig. 6). This total is 0.77 in. greater than the monthly normal precipitation measured in June 2020 (3.43 in.) and about 24 percent of the annual precipitation in 2020 (17.32 in.; fig. 6).

The water level of New Years Lake decreased from August 2020 to April 2021 during dry conditions (fig. 6). Annual precipitation in 2020 and 2021 were below annual normal precipitation (fig. 6), and departure from normal (fig. 2) indicated below normal precipitation for all but one month (March 2021) from July 2020 to April 2021. The period of below normal precipitation from July 2020 to April 2021 resulted in decreasing water levels in New Years Lake and, unlike previous years, the effect of spring snowmelt did not affect water levels in spring 2021 (fig. 6); however, the effects of spring snowmelt on water levels in New Years Lake may be after the last recorded water-level measurement on April 16, 2021.

Stable Isotopic Analysis

Stable isotope results were analyzed by (1) creating a local meteoric water line (LMWL) from precipitation samples collected during this study, (2) comparing stable isotope values from the study area to global meteoric values established by Craig (1961), and (3) comparing stable isotope values among sites to evaluate possible recharge characteristics.

Comparison with Global Meteoric Waters

The δ²H and δ¹⁸O values from samples collected during this study and other investigations (Long and others, 2012, 2019) were plotted with the global meteoric water line (GMWL) defined by Craig (1961; fig. 7). All samples plotted below the GMWL; therefore, a LMWL was calculated for the study area (fig. 7). The LMWL for the study area from samples collected in June 2021 is given in the equation below as

δ

²

H=6.24×δ

¹⁸

O−20.17

(2)

The slope and y-intercept of the LMWL were less than the GMWL and as a result the LMWL plotted below the GMWL (fig. 7). Putman and others (2019) observed the LMWL at arid and snowy sites, similar to the study area, often displayed higher slopes and intercepts relative to the GMWL. The deviation from trends observed by Putman and others (2019) may result from the LMWL not accurately representing the δ²H and δ¹⁸O values of precipitation in the study area. The natural variation of δ²H and δ¹⁸O values in precipitation likely was not captured from the three precipitation samples collected during this study because all three samples were collected during the same month and year (June 2021). Putman and others (2019) also discussed the sensitivity of developing a LMWL and determined the length of record needed to calculate a LMWL depended on the timescale of the processes of interest, such as recharge to an aquifer. Therefore, the LMWL for the study area likely is not an accurate representation of the average δ²H and δ¹⁸O values of precipitation in the study area; however, the LMWL may be useful for comparison of samples from other sites with precipitation samples collected during the same time period.

All stable isotope samples were depleted in heavy isotopes compared with the GMWL (more negative δ²H and δ¹⁸O values) and indicated a greater degree of evaporation from environmental factors, such as latitude, continental position (near coasts versus far inland), and altitude (Dansgaard, 1964; Rozanski and others, 1993). Medler and Eldridge (2021) provide a summary of the effect of latitude, continental position, and altitude on isotopic composition of precipitation. The most notable effects on the δ²H and δ¹⁸O compositions in the study area were from the latitude and continental effects because the study area is at a high latitude (between 43- and 44-degrees north latitude) and near the center of the North American continent. The great distance from low latitudes and coastlines to the study area (hundreds of miles) likely resulted in greater depletion of heavy isotopes for precipitation samples relative to studies at low latitudes or closer to coastlines (Dansgaard, 1964; Gat and Gonfiantini, 1981). The altitude effect was not observed for samples in the study area likely because the altitude difference among all sampled sites is only about 400 ft (table 1).

Streams, Wells, and Cave Lakes

Samples from one stream, two wells, and three cave lakes plotted either on or below the LMWL (fig. 7) and indicated generally greater depletion (more negative values) of heavy isotopes relative to precipitation and some spring samples. The Hell Canyon Creek (HCcr) sample from 2021 displayed the greatest depletion in heavy isotopes of all samples (fig. 7) and was considered anomalous because precipitation, generally enriched in heavy isotopes, was expected to increase the concentration of heavy isotopes in the stream. Springflow from West Hell Canyon Spring (WHCsp) is from the Madison aquifer and groundwater flows directly into HCcr upstream of the HCcr sampling site (fig. 1). Discharge from the Madison aquifer at WHCsp (fig. 7) appears to strongly affect the isotopic composition of HCcr. Samples from wells completed in the Deadwood aquifer (JC1 and JC2) in 2016 plotted close to the LMWL and had varying stable isotopic composition (fig. 7). Jewel Cave well 2 (JC2) was more depleted than Jewel Cave well 1 (JC1) despite the wells being less than 1,000 ft apart and both samples collected on the same day in 2016 (fig. 7). The stable isotopic composition difference between groundwater samples from the two wells completed in the Deadwood aquifer (JC1 and JC2) was not further examined because it was not the focus of this study. Samples from cave lakes—Hourglass, New Years, and Wellspring Lakes—had intermediate stable isotopic composition and were depleted and enriched compared to precipitation samples and samples from HCcr and WHCsp, respectively (fig. 7).

The depleted δ²H and δ¹⁸O values of groundwater samples—combined with the stable isotopic variation among samples from the same aquifer—in the stable isotope plot of the study area (fig. 7) provided insight on recharge mechanisms. Groundwater samples from the Deadwood aquifer (JC1 and JC2) and Madison aquifer (HGL, NYL, WSL, and WHCsp) were depleted in heavy isotopes compared to precipitation samples. Evaporation cannot explain the depleted groundwater samples because evaporation preferentially removes lighter isotopes and leaves behind heavier isotopes (Rozanski and others, 1993). Isotopic compositions typically do not change after precipitation has infiltrated into the ground (Gat, 1971; Naus and others, 2001); therefore, some other source or process was responsible for the depleted groundwater samples.

Gat (1971) provided four explanations that could account for groundwater samples being depleted in heavy isotopes. The first explanation is both aquifers contain large fractions of older water from past climates characterized by colder temperatures and snow containing lighter, more depleted waters, similar to the isotopic composition of snow and ice from modern Arctic and Antarctic environments (Gat, 1971; Rozanski and others, 1993). Age dating of water from the Madison aquifer using tritium (Naus and others, 2001; Rahn, 2018) indicated groundwater near recharge sources is newer and more like modern precipitation, whereas groundwater at greater distances from recharge areas is older and more like past precipitation. The second explanation is heavy rainfall, which generally is depleted in heavy isotopes, provides recharge to aquifers in the study area. The study area typically does not receive steady heavy rainfall; however, Vogel and Van Urk (1975) and Rehm and others (1982) determined that heavy rainfall was an important source of recharge for aquifers in semiarid regions. The third explanation is that snowfall and snowmelt are primary recharge sources for aquifers. Tian and others (2018) observed depleted δ²H and δ¹⁸O values in snowfall relative to spring, summer, and fall precipitation. The fourth explanation is that recharge resulted at a higher elevation where heavy isotopes are often more depleted than at lower elevations. Less than 10 mi from the study area the elevation increases to more than 7,000 ft above NGVD 88.

The first explanation from Gat (1971) likely does not apply to the Madison aquifer in the study area because Jewel Cave is near the recharge area for the Madison aquifer and karstic aquifers, such as the Madison aquifer, usually can be recharged quickly. All samples plot near the LMWL, and recharge areas (outcrops of Madison Limestone shown in figure 8; Redden and DeWitt, 2008) for the aquifer are near the study area—indicating aquifers in the study area are recharged from recent meteoric sources. A key component of the recharge in the study area is from karst features in the Madison Limestone, such as sinkholes and caves, that can recharge the aquifer relatively quickly compared to non-karstic aquifers. Karst features increase the capacity of the Madison aquifer to accept recharge and, in areas of outcrop of the Madison aquifer, water from precipitation and streams can be lost entirely to the Madison aquifer (Greene, 1997; Driscoll and others, 2002). The Hell Canyon Creek watershed delineated at a point downstream (south) of Jewel Cave National Monument is composed of about 80 percent of Madison Limestone outcrop (USGS, 2022b; Redden and DeWitt, 2008; fig. 8), and streamflow in the watershed could be recharging the Madison aquifer, in addition to recharge from direct infiltration of precipitation.

A combination of the latter three explanations from Gat (1971) best account for heavy isotope depletion observed in the Madison aquifer underlying the study area. Recharge to the Madison aquifer in the study area likely is from large precipitation events and spring snowmelt. The study area does not receive steady heavy rainfall; however, the infrequent large precipitation events likely overcome evapotranspiration and would explain the generally depleted isotope values from groundwater samples. Recharge from snowmelt also may contribute to the depletion of heavy isotopes in the Madison aquifer. Climate normal data from 1991 to 2020 shows the study area receives about 57.9 in. of snow per year (NOAA, 2021b). Winter and spring snowmelt could contribute to recharge in the study area when evapotranspiration rates are lower than summer rates (Clark and Fritz, 1997). Large amounts of water released during spring snowmelt can cause normally dry streams, such as Hell Canyon Creek, to flow. Streamflow in Hell Canyon Creek following spring snowmelt could lose all or part of its flow to sinkholes in areas of Madison Limestone outcrop. The location or presence of streamflow loss zones along Hell Canyon Creek is unknown and streamflow data from streamgages upstream and downstream of Jewel Cave National Monument do not overlap to indicate the presence of streamflow losses.

The depletion of heavy isotopes for groundwater samples at Jewel Cave also could be from recharge at higher elevations. The elevation of the Hell Canyon Creek watershed generally decreases from northeast to southwest and ranges from 5,050 ft to 7,170 ft (USGS, 2022b). A generalized potentiometric map from Anderson and others (2019) shows groundwater flow in the Madison aquifer generally follows topography—indicating groundwater in the Madison aquifer flows from higher elevations in the northeastern part of the Hell Canyon Creek watershed toward the study area to the southeast (fig. 8). The effect of the Jewel Cave Fault Zone (fig. 8)—with an estimated vertical offset ranging from 100 (Darton and Paige, 1925) to 300 ft (NPS, 1994)—on groundwater flow in the Madison aquifer is unknown (Dyer, 1961). If the Jewel Cave Fault zone disrupts regional groundwater flow in the Madison aquifer to the cave lakes, then the depleted isotopic signature of the cave lakes could be explained by streamflow losses along Hell Canyon Creek. The isotopic signature of water from Hell Canyon Creek and its source from West Hell Canyon Spring also were depleted in heavy isotopes similar to Hourglass and New Years Lakes (fig. 7).

Springs Discharging from the Minnelusa Aquifer

Samples from springs discharging from the Minnelusa aquifer were qualitatively assessed to evaluate water composition using the stable isotope plot (fig. 7). The δ²H and δ¹⁸O values of springs discharging from the Minnelusa aquifer in the Minnelusa Formation generally were intermediate and plotted close to the LWML between groundwater samples from the Madison aquifer and precipitation samples. The intermediate composition of springs discharging from the Minnelusa aquifer indicated a mixture of groundwater and precipitation; however, the contribution from groundwater (more depleted in heavy isotopes) and precipitation (more enriched in heavy isotopes) varied for each spring. Stable isotopic composition was qualitatively assessed by observing the plotting positions of each site relative to other sites. Samples from A&E Spring (AEsp) and Water Draw Spring (WDWsp) plotted closer to depleted groundwater samples from cave lakes, whereas samples from Lithograph Spring (LTHsp), Chokecherry Spring (CCHsp), and Prairie Dog Spring (PDGsp) were more enriched in heavy isotopes and plotted closer to precipitation samples (fig. 7).

The variability of δ²H and δ¹⁸O values from springs with more than one sample was investigated using the stable isotope plot (fig. 7). Springs with more than one sample were Lithograph Spring (LTHsp) and Prairie Dog Spring (PDGsp). The stable isotopic composition of LTHsp changed from being more depleted in May 2016 to more enriched in August 2021. The stable isotopic composition of PDGsp remained relatively consistent between samples collected in May 2016 and August 2021. The variance or lack thereof in δ²H and δ¹⁸O values at LTHsp and PDGsp likely resulted because of changes in the contribution and type of precipitation falling shortly before sample collection. The contribution of groundwater or precipitation to LTHsp likely was related to (1) the elapsed time between precipitation events and sample collection, and (2) the residence time of infiltrated precipitation before spring discharge. Spring samples collected soon after precipitation events with similar stable isotopic composition as precipitation would indicate a shorter residence time (fast groundwater movement or short flowpaths) of infiltrated precipitation in the aquifer. Conversely, spring samples collected soon after precipitation events with more depleted stable isotopic composition would indicate longer residence time (slow groundwater movement or long flowpaths) of infiltrated precipitation in the aquifer that had not yet reached the spring. Infiltrated precipitation mixed with depleted groundwater, however, could potentially obscure the δ²H and δ¹⁸O values of infiltrated precipitation. A long residence time also may explain variable δ²H and δ¹⁸O values as water discharging from springs on each collection date could be water recharged during different times (spring snowmelt or summer rains). Stable isotope data collected from additional recharge sources (snowmelt and streams close to springs) and more frequently at springs could potentially be used to estimate residence time of water sources at each spring.

The intermediate δ²H and δ¹⁸O values of samples from springs discharging from the Minnelusa aquifer—combined with the stable isotopic variation among spring samples—in the stable isotope plot of the study area (fig. 7) provided insight on recharge mechanisms to the Minnelusa aquifer. The principles of conservation of stable isotopic composition for infiltrated precipitation (Naus and others, 2001; Gat, 1971) and depletion of heavy isotopes in spring samples compared to precipitation (Gat, 1971) also apply to the Minnelusa aquifer. A combination of the second and third explanations provided by Gat (1971) best account for heavy isotope depletion in the Minnelusa aquifer in the study area. The fourth explanation by Gat (1971)—recharge from high-elevation precipitation—is unlikely because higher elevation recharge areas immediately east of the study area are only about 200 ft higher in elevation than the surrounding areas (fig. 9). All spring samples plot near the LMWL and recharge areas for the Minnelusa aquifer are near the study area—indicating recharge to the Minnelusa aquifer is from recent meteoric sources. Recharge to the Minnelusa aquifer is primarily from precipitation on outcrops and streamflow losses; however, recharge from streamflow losses to the Minnelusa aquifer generally is less than to the Madison aquifer (Carter and others, 2003). In the study area, outcrop of the Minnelusa aquifer is mostly south of the Jewel Cave Fault Zone (Redden and DeWitt, 2008; fig. 9). The Hell Canyon Creek watershed is composed of about 20 percent of Minnelusa Formation outcrop (Redden and DeWitt, 2008; fig. 9), and precipitation or streamflow in the watershed could be lost to the Minnelusa aquifer.

A combination of recharge from heavy rainfall, snowmelt, and aquifer residence time could explain depletion of heavy isotopes in samples from springs discharging from the Minnelusa aquifer. Samples from A&E Spring (AEsp) and Water Draw Spring (WDWsp) were collected in May 2016, so effects from spring snowmelt depleted in heavy isotopes could explain why those samples plot closer to groundwater samples. Effects of different precipitation types was apparent for samples from Lithograph Spring (LTHsp) because the sample collected in May 2016 was more depleted in heavy isotopes compared to the sample collected in August 2021 (fig. 7). The LTHsp sample collected in August 2021 likely was affected by enriched summer precipitation, which could indicate shorter aquifer residence time (fast groundwater movement or short flowpaths). Conversely, samples collected in May 2016 and August 2021 at Prairie Dog Spring (PDGsp) did not vary appreciably in stable isotope content. The lack of stable isotope variation at PDGsp could indicate a consistent groundwater source and (or) longer residence time of water in the aquifer resulting from either slower groundwater flow or longer flowpaths. The sample collected in May 2016 at Chokecherry Spring (CCHsp; fig. 7) was considered an anomaly because it was more enriched in heavy isotopes than all other samples. The enriched CCHsp sample could indicate its stable isotopic composition was affected by precipitation falling shortly before sample collection. Additional stable isotope and precipitation data would be required to better understand stable isotope variations and aquifer residence time at springs in the study area.

Another possible explanation for the depletion of heavy isotopes in samples from spring discharging from the Minnelusa aquifer could be from mixing with water from the Madison aquifer. Groundwater from the Madison aquifer was shown to be depleted in heavy isotopes (fig. 7) and mixing between the two aquifers would produce an overall depleted sample; however, hydraulic connection between the Minnelusa and Madison aquifers is not well defined (Carter and others, 2003). It is possible that the hydraulic connection between the Minnelusa and Madison aquifers is present along geologic structures, such as the Jewel Cave Fault Zone, where structure offset has juxtaposed the two units (fig. 9). The stable isotopic composition of samples from HGL—the cave lake sampled closest to the Jewel Cave Fault Zone—is similar to springs discharging from the Minnelusa aquifer, but additional data are needed to determine possible mixing between the two aquifers.

Site Groupings and Water Sources from Principal Component Analysis

The PCA method was used to create a biplot of water-chemistry data from this study and other studies (Naus and others, 2001; Long and others, 2012, 2019; fig. 10). The combined sum of principal components one and two explained 55.3 percent of the total variance in the dataset. Principal component three explained 14.3 percent of the total variance. Multiple variables were responsible for groupings along principal component axes one and two. Separation and grouping along principal component axis one were from a combination of specific conductance, hardness, ions (calcium, magnesium, chloride, and sulfate), and stable isotopes of oxygen and hydrogen. Loading lines for specific conductance, hardness, and ions displayed varying degrees of positive correlation, but all variables had positive loading on principal component one—shown by the loading lines plotting along the positive part of the x-axis (fig. 10). Conversely, loading lines for δ¹⁸O, δ²H, arsenic, silica, and pH had negative loading on principal component one (fig. 10). The δ¹⁸O, δ²H, and silica variables displayed positive correlation to each other and negative correlation to specific conductance, hardness, and ions—indicated by the nearly 180-degree angle between loading lines. The δ¹⁸O and δ²H variables were expected to contribute to groupings because the stable isotope plot (fig. 7) showed strong positive linear correlation. Additionally, the δ¹⁸O, δ²H, and silica variables showed almost no correlation to the pH and carbon dioxide variables—shown by the approximately 90-degree angle between loading lines (fig. 10). Separation and grouping along principal component axis two were from the pH, silica, and carbon dioxide variables.

The plotting positions of sites in the PCA biplot show groupings based on location, aquifer, and the type of site (fig. 10). Sampling sites were identified by color and symbol for each aquifer and site type, respectively, where appropriate to make clear the relations among aquifers in the southern Black Hills (fig. 10). Precambrian aquifer sites (fig. 10) plotted mostly in upper parts of quadrants I and II—indicating lower pH and specific conductance, and higher concentrations of carbon dioxide. Williamson and Carter (2001) also noticed lower pH and specific conductance values, and higher carbon dioxide concentrations in groundwater from Precambrian aquifers compared to the Madison and Minnelusa aquifers. Most Madison aquifer sites (fig. 10) from the East area (fig. 1) plotted in the upper part of quadrant III with other sites plotting in quadrants I, II, and IV. Minnelusa aquifer sites (fig. 10) from the East area were scattered among all four quadrants but plotted mostly in quadrants I and II. Minnelusa aquifer sites were distinguishable from Madison aquifer sites in the East area because of their higher ion concentration, specific conductance, and hardness values. Williamson and Carter (2001) also observed that groundwater from the Minnelusa aquifer generally had higher ion concentrations (especially the calcium and sulfate ions) and greater specific conductance and hardness values than groundwater from the Madison aquifer.

Groundwater and surface-water sites from the West area (fig. 1) displayed distinctive water chemistry and plotted away from most samples from the East area with the exception Hourglass Lake (HGL), Mackenzie Spring (MCKsp), and Jewel Cave Well 1 (JC1; fig. 10). All groundwater samples from aquifer sites from the West area had intermediate ion concentrations and relatively lower δ¹⁸O and δ²H values compared to other samples from the East area; however, these samples from the West area had distinctively high magnesium concentration, and high hardness, pH, and specific conductance values that created separation along the y-axis (fig. 10). HGL and MCKsp both plotted near the y-axis of the PCA biplot between Madison aquifer sites from the East and West areas because of their higher δ¹⁸O and δ²H values and lower ion concentration, specific conductance, and hardness values (fig. 10). Some well and spring sites from the East area (fig. 10) were chemically similar to West area sites and plotted in quadrant I; however, East area sites in quadrant I are not immediately hydrologically connected to the West area because of the large distance (miles) between the areas (fig. 1). The higher ion concentration, and specific conductance, and hardness values for East area sites in quadrant I were likely caused by the evolution of groundwater from dissolution of aquifer material as it travels downgradient in the aquifer (fig. 1). Williamson and Carter (2001) determined ion concentrations, and specific conductance, and hardness values increase in the Minnelusa and Madison aquifers with increasing well depth and with distance from outcrop areas. For example, the COL, CRA, PEK, HBsp, and PAL sites in quadrant I are all downgradient of recharge areas further west (fig. 1) and all had higher ion concentration, and specific conductance and hardness values than sites closer to recharge areas.

The plotting positions of sites displayed in figure 10 from the West area were used to determine differences in water chemistry among the Deadwood, Madison, and Minnelusa aquifers. A similar analysis could be conducted on sites from the East area; however, that analysis is outside the scope of this study. Groundwater samples from Deadwood aquifer sites (JC1 and JC2) had similar water chemistry as Minnelusa aquifer sites (fig. 10); however, the water chemistry of Jewel Cave well 1 (JC1) and Jewel Cave well 2 (JC2) differed appreciably. The sample from JC1 had higher ion concentrations commonly associated with salts (sodium, calcium, magnesium, and chloride). Groundwater from the Madison aquifer generally had the lowest ion concentration, and specific conductance and hardness values in the West area (fig. 10). Additionally, water chemistry of cave lake sites notably differed; HGL plotted closer to Madison aquifer sites from the East area, whereas New Years Lake (NYL) plotted closer to Hell Canyon Creek (HCcr) and Minnelusa aquifer springs (fig. 10). West Hell Canyon Spring (WHCsp)—upgradient from cave lake sites—also differed from cave lake samples because WHCsp had higher ion concentration, and specific conductance and hardness values. The different water chemistry at WHCsp could result from the dissolution of aquifer materials during groundwater flow that increase ion concentrations and specific conductance (Williamson and Carter, 2001). Groundwater samples from Minnelusa aquifer sites displayed the greatest water-chemistry variation among the three aquifers, but generally had the highest ion concentration, and specific conductance and hardness values in the West area (fig. 10).

Cluster Assignments and Water Sources from Cluster Analysis

The k-means cluster procedure was applied to the PCA results shown in figure 10 to statistically group sites (fig. 11). Sampling sites were grouped into five categories based on cluster assignments from the k-means procedure (fig. 11; table 3). Group 1 included samples from Precambrian aquifer springs and wells and one White River aquifer sample. Group 2 mostly consisted of samples from the Madison aquifer in the East area (fig. 1). Group 3 consisted mostly of samples from springs and wells completed in the Minnelusa aquifer in the East area. Group 4 included mostly samples from the West area. Group 5 consisted of springs and wells completed in the Minnelusa and Precambrian aquifers from the East area with high ion concentrations and specific conductance and hardness values.

The k-means cluster procedure grouped samples into clusters sharing similar physical property, chemical constituent, and stable isotope values that were used to identify relations among sites from different aquifers in the East and West areas (fig. 1). Group 1 included only sites from the Precambrian aquifer in the East area (table 3); however, other Precambrian aquifer sites in group 5 (KIR1 and MOR; fig. 11) were not included in group one. Precambrian aquifer sites not included in group 1 were from well sites with water chemistry more similar to the chemistry of the Minnelusa aquifer (group 3). Groundwater from the Minnelusa aquifer could not affect water chemistry for Precambrian aquifer sites excluded from group 1 because all the collected samples were upgradient and isolated from the Minnelusa aquifer outcrop (fig. 1). Possible explanations for Precambrian aquifer sites excluded from group 1 could be from well treatment or elevated ion concentrations from road salts that would increase specific conductance values.

All well sites in group 2 were from the East area (fig. 1) and 12 of 14 of the sites were wells completed in the Madison aquifer (table 3). Separation of sites in the East area from sites in the West area can be explained by the altitude effect (Dansgaard, 1964; Rozanski and others, 1993) because the surface elevation of sites in the East area is about 1,400 ft lower than sites in the West area. As expected, the lower elevation East area well sites had higher δ¹⁸O and δ²H values (enriched with heavy isotopes) than the higher elevation West area sites. The FRA site was the only Minnelusa aquifer site in group 2 (fig. 1; table 3), and the PCA and cluster analysis results indicated that this well could be completed in the Madison aquifer rather than the Minnelusa aquifer as indicated in NWIS (USGS, 2022a). The FRA site did not have a drillers log with lithologic information in the South Dakota Department of Agriculture and Natural Resources (2022) database; however, the well depth of the FRA site was 440 ft, as available from NWIS (USGS, 2022a), and nearby wells with drillers logs indicated the depth to Madison Limestone is about 690 ft. Therefore, the FRA well likely is completed in the Minnelusa aquifer, but groundwater from this well has similar water chemistry as the underlying Madison aquifer. Separation of group 2 from other sites in the East area was from a combination of higher arsenic concentrations and lower ion concentrations, and specific conductance and hardness values than other sites (fig. 11).

Group 3 consisted of 10 sites with 1 from surface water, 1 from the Madison aquifer, and 8 from the Minnelusa aquifer (table 3). Additionally, 9 of the 10 sites in group 3 were from the East area and other was from the West area (table 3). Group 5 was similar to group 3; group 5 consisted of 8 sites from the East area with 2 from the Precambrian aquifer and 6 from the Minnelusa aquifer. Differences between sites completed in the Minnelusa aquifer in groups 3 and 5 likely were related to proximity to recharge areas as compared to other sites. Sites in group 3 often were closer to recharge areas of Minnelusa aquifer outcrop, whereas group 5 included mostly well sites further downgradient and at deeper well depths (fig. 1; table 3).

The KAI site (fig. 1) was the only site in group 3 completed in the Madison aquifer (table 3). The water chemistry of the KAI site was similar to the Minnelusa aquifer—indicating the KAI well could be completed in either the Minnelusa aquifer or both in the Minnelusa and Madison aquifers. The KAI site did not have a drillers log with lithologic information in the South Dakota Department of Agriculture and Natural Resources (2022) database; however, the well depth of the KAI site was 780 ft as obtained from NWIS (USGS, 2022a). Nearby wells with drillers logs indicated that the depth to the Madison Limestone is about 690 ft; therefore, the well probably is completed at least partially in the Madison aquifer and could be partially completed in the Minnelusa aquifer depending on the position and length of the screened (open) interval.

Group 4 contained 16 sites with 2 from the Deadwood aquifer, 5 from the Madison aquifer, 8 from the Minnelusa aquifer, and 1 surface-water site. Additionally, 12 of the 16 sites contained in group 4 were from the West area, and the other 4 were from the East area (table 3). As previously discussed, sites in group 4 differed from other sites because of their low δ¹⁸O and δ²H values (depleted) and high ion concentrations, and specific conductance and hardness values. The higher elevation of West area sites in group 4 compared with sites contained in other groups likely caused the lower δ¹⁸O and δ²H values (Dansgaard, 1964; Rozanski and others, 1993). Anomalies in group 4 from the East area included two wells completed in the Madison aquifer (CON and SVE), one spring discharging from the Minnelusa aquifer (Elk Mountain Spring; EMsp), and one well completed in the Minnelusa aquifer (PAL; fig. 11). Inclusion of the CON and SVE sites in group 4 could result from longer residence times in the Madison aquifer or groundwater mixing with the Minnelusa aquifer. Elk Mountain Spring (EMsp) and the PAL well were included in group 4 because of abnormally high magnesium concentrations and hardness values.

Hourglass Lake

Water-chemistry data—evaluated using PCA and cluster analyses—highlighted differences between Hourglass and New Years Lakes that suggested different recharge mechanisms. PCA and k-means cluster analysis results showed that Hourglass Lake (HGL) grouped with Madison aquifer sites from the West area (fig. 1) but plotted away from other sites in group 4 (figs. 10 and 1) because it was more enriched in heavy isotopes and had lower concentrations of ions than other sites from the West area. Water chemistry of HGL was more similar chemically to precipitation (enriched in heavy isotopes with low ion concentrations) than other West area sites—indicating recharge to the lake results quickly because its water chemistry was less affected by evaporation and (or) dissolution of aquifer materials than other West area sites. Additionally, the water chemistry of HGL was noticeably different than West Hell Canyon Spring (WHCsp; figs. 10 and 11) that likely receives its recharged water from the high elevation Madison Limestone outcrops in the northern part of the study area (fig. 8). The difference in water chemistry between HGL and other Madison aquifer sites in the West area—combined with differences between hydrographs of HGL and NYL—indicates that HGL is recharged at a different rate and from a different or additional source(s) than other Madison aquifer sites in the West area.

Taylor and Greene (2008) provided a conceptual model of karst aquifers showing multiple recharge sources (fig. 14) that were used to hypothesize recharge mechanisms for Hourglass Lake. Recharge sources were differentiated in terms of residence time and the amount of water contributed to the aquifer. Taylor and Greene (2008) provided descriptions of each type of recharge source that will be briefly discussed here. In general, autogenic recharge and allogenic recharge are defined as precipitation falling on karstic and nonkarstic terrane, respectively (Gunn, 1983). The amount and timing of the precipitation determines whether the recharge is considered concentrated or diffuse (Gunn, 1983). Concentrated autogenic recharge (fig. 14) primarily is from filling of sinkhole depressions by surface runoff processes that can be either rapidly drained by swallets or relatively slowly through percolation of soil and other deposits. Diffuse autogenic recharge (fig. 14) from infiltration of water through soils overlying karstic terrane is the most common type of recharge to karst aquifers. Concentrated allogenic recharge (fig. 14) to karst aquifers is water contributed by sinking or losing streams originating as normal gaining streams in non-karstic terrain. Diffuse allogenic recharge (fig. 14) can be groundwater contributed from non-karstic aquifers, but more commonly is water that drains through vertical or near-vertical conduits in the unsaturated (vadose) zone where karstic rocks are overlain by non-soluble caprocks, such as sandstone.

Recharge to HGL likely is a mix of various recharge sources, with the dominant recharge sources likely diffuse allogenic recharge through the Minnelusa Formation (including the Minnelusa aquifer and semiconfining units) and diffuse autogenic recharge in areas of Madison Limestone outcrop (fig. 14). Another possible recharge source could be streamflow losses (concentrated allogenic recharge from stream-sink; fig. 14). The Madison Limestone containing HGL is overlain by about 300–400 ft of the Minnelusa Formation, and its water chemistry was more similar to precipitation than other Madison aquifer sites—indicating the largest fraction of recharged water likely is from quickly infiltrated precipitation or streamflow losses. Dye tracer tests from Alexander and others (1989) indicated precipitation falling on the land surface above the cave infiltrates through the Minnelusa Formation and into Jewel Cave. The relatively quick recovery of dye 8 days after dye insertion may have been through vertical or near-vertical fractures connecting the Madison and Minnelusa aquifers, whereas the intermittent and reoccurring recovery of dye 1 year after insertion may be caused by slow vertical draining of the Minnelusa aquifer into the Madison aquifer.

Diffuse allogenic recharge also can explain water-chemistry similarities between HGL and precipitation. Some fraction of infiltrated precipitation would reach HGL quickly through vertical or near-vertical conduits, and its water chemistry would be similar to precipitation because of its short residence time in the Minnelusa aquifer. The other fraction of infiltrated precipitation would drain through the Minnelusa aquifer slower than water passing through vertical or near-vertical conduits, and its water chemistry would contain higher concentrations of ions and have higher specific conductance and hardness values. Another source of water supplying cave lakes could be groundwater flowing from outcrop recharge areas, such as those north of Jewel Cave (fig. 8), to Jewel Cave, which likely would have higher ion concentrations, and greater specific conductance and hardness values than recently infiltrated precipitation. Additionally, diffuse allogenic recharge can explain the relatively slow response between water levels in HGL and precipitation. Vertical or near-vertical conduits extending through the Minnelusa Formation (including the Minnelusa aquifer) can transport a certain quantity of water that may be small compared to the volume of water in HGL, so water levels would slowly change at HGL. Water-level increases at HGL during dry periods, such as in 2020–21 (figs. 3 and 5), could be related to the relatively slow downward vertical movement of water through the Minnelusa Formation (including the Minnelusa aquifer and the unsaturated zone of the Madison Limestone.

New Years Lake

PCA and k-means cluster analysis results showed New Years Lake (NYL) grouped with sites from the West area and plotted closest to Hell Canyon Creek (HCcr; figs. 10 and 11). Similar water chemistry between NYL and HCcr suggests that NYL may receive recharge from streamflow losses. HCcr is an intermittent stream that flows during wet conditions when springs in Hell Canyon, such as West Hell Canyon Spring (WHCsp), supply enough water to overcome evapotranspiration and streamflow losses to either the Madison or Minnelusa aquifers (NPS, 1994). Cumulative monthly departure from normal (fig. 2) indicates wet conditions throughout the study area from May 2018 to August 2021 and HCcr was flowing during sample collection in August 2021. Streamflow in HCcr is not currently monitored so the duration of streamflow is unknown. Water chemistry of HCcr likely was affected by the water chemistry of springs discharging into it. WHCsp discharges into HCcr; however, PCA analysis shows that WHCsp and HCcr plot away from each other (fig. 10) despite having similar stable isotope compositions (fig. 7) because the carbon dioxide concentration, and specific conductance and hardness values at WHCsp were noticeably higher than HCcr. Water chemistry of HCcr was distinguishable from WHCsp because HCcr received water from a variety of sources, including surface runoff and springs discharging from the Madison and Minnelusa aquifers.

Differences between NYL hydrograph and hydrographs of HGL and CU-95A indicated a different recharge source for NYL. Water-chemistry data indicated similarities between NYL and HCcr; however, comparisons between streamflow from HCcr and water levels at NYL could not be made because streamflow data were not measured during the period of water-level data collection at NYL. Greene and others (1999), Hortness and Driscoll (1998), and Carter and others (2001) compared streamflow data to water levels in wells for other streams in other parts of the Black Hills and noticed correlation between streamflow losses and water levels in the Madison aquifer.

It is unknown if HCcr was flowing during the entire period of water-level data collection at NYL evaluated from 2018 to 2021. Additionally, peak streamflow in HCcr likely would not coincide with peak water levels in NYL because of the time necessary for streamflow losses to recharge the lake. Therefore, a time delay between streamflow losses and water levels would be expected. Another consideration when comparing water-level changes at NYL to streamflow is the difficulty of differentiating recharge from streamflow losses compared to infiltrated precipitation because streamflow increases with increasing precipitation. Water-chemistry and hydrograph similarities between NYL and surface water indicated a possible connection to surface water in the study area, although additional data collection is required to further define the relation between increasing water levels and streamflow losses.

Streamflow loss zones along HCcr have not been mapped but their presence in the Jewel Cave area has been speculated by other studies (Dyer, 1961; NPS, 1994). Dyer (1961) postulated the Jewel Cave Fault Zone (fig. 8) could capture streamflow and (or) shallow groundwater in alluvial deposits in Hell Canyon; however, no evidence was provided to confirm this possibility. NPS (1994) discussed how streamflow in the Jewel Cave area could be affected by exposures of Madison Limestone and, to a lesser extent, the Minnelusa Formation because the porous nature of the substrate could capture streamflow and precipitation. A geologic map of the Jewel Cave area shows outcrops of Madison Limestone along HCcr south of the Jewel Cave Fault Zone (fig. 8). It is possible that the Jewel Cave Fault Zone, some other unnamed structure, or outcrops of Madison Limestone and Minnelusa Formation along HCcr could capture streamflow and (or) shallow groundwater and provide recharge to the Madison aquifer near NYL. The potential relation between streamflow in HCcr and water levels at NYL would be clearer if additional streamflow and water-level data were collected simultaneously.

Water-level changes and water-chemistry data supported the possibility that NYL is recharged by concentrated allogenic recharge (streamflow losses) along HCcr (fig. 14). NYL is overlain by the Minnelusa Formation and its mean water level is about 650 ft below land surface; therefore, a natural conduit connecting HCcr to NYL would have to be fairly large and continuous to quickly transport water from the surface to the cave. A crack in the cave ceiling with water discharging into NYL was observed by NPS staff (Dan Austin, NPS, written commun., December 2021). NPS staff also noted horizontal displacement along the crack that was estimated at about 1 ft (Dan Austin, NPS, written commun., December 2021). The aperture (width of the opening) of the crack was not measured; however, displacement in the cave ceiling indicated the crack was a fault. The fault observed in the cave ceiling above NYL was named “Delmars Fault” by NPS staff and little is known about its extent. It is possible that Delmars Fault is extensive and could provide a natural conduit for direct recharge from HCcr to NYL.

When HCcr is not flowing, the dominant recharge source to NYL probably changes from concentrated allogenic recharge (streamflow losses) to a combination of diffuse allogenic and autogenic recharge (fig. 14). Diffuse allogenic recharge likely results where the Minnelusa Formation overlies the Madison Limestone and would be similar to recharge at Hourglass Lake. Diffuse autogenic recharge likely results either where the Madison Limestone is exposed at the land surface or where alluvial deposits in Hell Canyon directly overlie Madison Limestone south of the Jewel Cave Fault Zone (fig. 8). Diffuse allogenic and autogenic recharge both likely contribute to recharging NYL when HCcr is flowing; however, contributions from diffuse allogenic and autogenic recharge likely are smaller with greater travel time to NYL because of slower vertical flow through the Minnelusa Formation and the unsaturated zone in the upper part of the Madison aquifer.

Susceptibility of Hourglass and New Years Lakes

Susceptibility of Hourglass and New Years Lakes from human-related activities at the land surface can be addressed qualitatively using observations from previous studies and the recharge mechanisms interpreted from the results of this study. Dye tracer tests (Alexander and others, 1989) demonstrated that dyed water inserted at the land surface passed through a septic system and entered Jewel Cave in 8 days. The dye was recovered for over 1 year after the tests began. Additionally, recharge mechanisms discussed in the previous section indicate that Hourglass and New Years Lakes are possibly well-connected to surface processes. Water-level and water-chemistry data from Hourglass and New Years Lakes indicated similarities to atmospheric water sourced from direct infiltration of precipitation and potential streamflow losses along Hell Canyon Creek. Therefore, human-related activities at recharge areas for the Madison and Minnelusa aquifers or in Hell Canyon could potentially affect the water chemistry of cave lakes in Jewel Cave. The theorized connection between Hell Canyon Creek and New Years Lake to recharge sources indicates the water chemistry of New Years Lake likely would be affected before Hourglass Lake; however, the water chemistry of Hourglass Lake could remain affected for a longer time period than New Years Lake because of its complex recharge mechanisms and likely connection to the Minnelusa aquifer. Observations from Alexander and others (1989), M.E. Wiles (South Dakota School of Mines and Technology, unpub. data, 1992), Anderson and others (2019), and from this study indicated recharge to Jewel Cave is complex and takes place on various timescales that are affected temporally by precipitation patterns and spatially by hydrologic connection with the overlying Minnelusa aquifer.