Position Determination

Survey grade positions were determined at several locations from 2019 to 2021 around the Zuni Salt Lake using the National Geodetic Survey’s Online Positioning User Service (National Geodetic Survey, 2015) at a Geodetic Survey marker (base station [BS]) on the maar’s rim. A real time kinematic (RTK) survey was made from the BS to establish control points and positions at other sites. The quality indicators for the rapid static occupation at the BS were 87 percent of the observations used, a quality indicator of 8.67/27.57, and a normalized RMS of 0.262 relative to the common datums North American Datum of 1983 (NAD 83) and North American Vertical Datum of 1988 (NAVD 88). However, the BS station position was verified to within 1.6 in. horizontally and 3.5 in. vertically by comparing prior reported elevations of the BS survey mark and secondary Online Positioning User Service occupations on positions determined by RTK survey from the BS. In addition, all of the locations and elevations relative to each other can be considered to have a relative accuracy of within about 0.4 in. horizontally and about 0.6 in. vertically, because all of the locations were determined by RTK survey less than 1.25 mi from that BS.

Digital Elevation Model Construction

The surrounding topography of the Zuni Salt Lake was mapped using UAS based imagery and analyzed with Agisoft Metashape Professional version 1.6.1 build 10009 (Agisoft, 2023). The UAS survey was conducted in the fall of 2019 and covered 1.1 mi with 1,404 images. The structure from motion (SfM) process was used to create the resulting DEM, which was verified with five control points established by the RTK survey at various locations around the lake. The root mean square error (RMSE) of the elevations from the SfM compared with surveyed point elevations was 0.1 ft. The bathymetry of the lake was estimated by interpolating water depths measured with staff rods from kayaks at 133 locations and recorded to the nearest 0.01 ft and UAS imagery shoreline estimates at 403 locations (total observation data count = 536). Measured water depths ranged from 0.22 to 5.50 ft (median = 0.83 ft, average = 1.01 ft); however, 75 percent of the depths were between 0.25 and 1.00 ft, indicating a relatively flat bottom. A small depression near a former pump house, and the current location of a water-level gage in the lake, was where most depths deeper than 1 ft were recorded.

Lake bottom elevation was interpolated by fitting a thin-plate spline to the point measurements. The spline model was fit by using the Tps function of the fields package (version 9.8-6; Nychka and others, 2021) in R (version 3.5.3; R Core Team, 2023). Spline complexity (in other words, effective degrees of freedom) was determined by using five-fold cross validation (CV; James and others, 2013). This CV methodology randomly reserves 20 percent of the data as test data while using the remaining 80 percent of the data to train the spline model. The fitted spline model is then used to predict the test data to assess its ability to accurately make predictions. This process is repeated four additional times in one five-fold CV run, so that eventually all of the data are used as test data in the analysis. The five-fold CV procedure was run 500 times for spline complexities ranging from 3 to 200 degrees of freedom for a total of 495,000 model fits. The RMSE was calculated for all model runs as a function of spline complexity.

Mean CV RMSE was minimized by the spline model with 121 degrees of freedom. A “one-standard error” rule was used to choose the most parsimonious model whose mean CV RMSE was no more than one standard error above the minimum mean RMSE; this is commonly done in CV in favor of simplicity and to avoid over-fitting the data (Hastie and others, 2009, p. 244). This approach reduced the final model complexity to 98 degrees of freedom (−19 percent) while increasing mean CV RMSE by just 0.0005 ft (0.15 percent).

Leave-one-out cross validation (LOOCV) was also performed on the final model to further evaluate its predictive abilities and spatial sensitivities. LOOCV is similar to five-fold cross validation, except that only one observation is withheld as test data for each model fit (James and others, 2013). LOOCV errors were generally small (mean error = 0.0002 ft) and only increased appreciably where there were abrupt changes in the Zuni Salt Lake bottom elevation. Increased model sensitivity near these types of features is inherently expected because of practical observation data density limitations and suggests that some of these more localized features are smeared (that is, smoothed over) where observation data are sparse. As a whole, LOOCV RMSE was 0.3450 ft, and LOOCV predictions correlated well with observed data (Pearson correlation coefficient = 0.82). Model predictions also lacked notable correlation with LOOCV residuals (Pearson correlation coefficient = 0.03), thereby indicating minimal systematic bias in model predictions.

The final thin-plate spline model was used to make predictions at rasterized locations throughout the Zuni Salt Lake bottom to estimate continuous lake bathymetry. A uniform raster cell size of 1 ft was selected as an adequate balance between computational expense and minimizing interpolation errors introduced by rasterizing the domain. The model parameters, predictive performance criteria, and LOOCV results are summarized in table 3, whereas the estimated bathymetric surface is presented in figure 4 (Bosch and others, 2025). Overall, the spline model represents the observed bathymetry data well, despite smearing some of the more localized elevation anomalies, and demonstrates adequate predictive capabilities for the purposes of this study.

The interpolated bathymetric surface was combined with local UAS-based topography in ArcGIS Pro (Esri, Redlands, Calif.) to permit Zuni Salt Lake volume estimation as a function of measured lake stage. The topography within the adjacent cinder cone was removed from the combined surface to facilitate hydrologically accurate lake volume estimates. The final raster is shown in figure 5.

Volume Determination

The lake volume was estimated using the resulting DEM to develop a rating curve between the Zuni Salt Lake stage (elevation) and volume using the U.S. Army Corps of Engineers Reservoir Inundation Calculator tool in ArcGIS. The greatest lake depth, 5.5 ft at the depression near current location of the lake gage, corresponds to an elevation of 6,216.1 ft above NAVD 88. However, smoothing in the spline interpolation and rasterization of the surface, along with the small, localized nature of the deeper parts of the lake, yield negligible storage estimates (less than 1 acre-foot [acre-ft]) below an elevation of about 6,220.0 ft above NAVD 88 (1.7 ft depth). The analysis steps from 6,219.0 to 6,227.0 ft above NAVD 88 at 0.1 ft increments. The resulting rating curve may be approximated by two linear segments, from elevations of 6,219.9 to 6,220.9 ft above NAVD 88 and from 6,221.0 to 6,227.0 ft above NAVD 88 (fig. 6).

Based on the elevation uncertainty estimates of the bathymetric survey (LOOCV elevation errors = 0.0002 ft) and the RMSE in the DEM created from the aerial survey (0.01 ft), the uncertainty in the Zuni Salt Lake storage was estimated as a function of the lake elevation from the rating curve. Storage changes corresponding to elevation changes of 0.01 ft as computed by the Reservoir Inundation Calculator tool, are between 0 and 2 acre-ft, with a mean change of 1.5 acre-ft and a median change of 1.7 acre-ft. The percent errors associated with this relationship decrease with increasing storage; for example, the percent errors in storage above a lake elevation of 6,221.5 ft above NAVD 88 are just under 20 percent, and the percent error is under 10 percent for lake elevations higher than 6,222.2 ft above NAVD 88. Approximately 80 percent of the daily lake elevations were higher than 6,221.5 ft above NAVD 88, so assuming a 20-percent error in the storage would represent a conservative estimate.

Precipitation

Precipitation was measured at the Zuni Salt Lake using a tipping-bucket rain gage. The tipping-bucket rain gage measures precipitation events using two identical chambers that alternately fill and drain. As each chamber fills, it tips, simultaneously draining it, bringing up the second bucket under the collector funnel and recording a known amount of precipitation, 0.01 in. The precipitation total for each day is computed by summing the number of tips that occur on that date. Data from tipping-bucket instruments are reported to the nearest 0.01 in. and are available at the USGS NWIS website (USGS, 2023).

Evaporation

The Collison Floating Evaporation Pan (CFEP; Collison, 2019) was used to estimate evaporation from the Zuni Salt Lake from September 2020 to January 2023. The automated design of the CFEP for measuring atmospheric conditions and water-level changes in the pan with native water was considered to be an important feature for this remote location. The CFEP is designed so that a decrease in water level within the evaporation pan per unit of time equals the evaporation rate (Collison, 2019). The water level is measured with a linear potentiometer and two pumps are used to fill, or drain, the pan to a set water level each night. Errors associated with this technique include thermal expansion of the pan and (or) the water inside the pan, pan oscillations from wave and wind action, and other types of instrument error. Based on the results from Collison (2019) and Reclamation (2021, 2023), the CFEP measurements are in reasonable agreement with measurements obtained from other evaporation techniques, including very low error (average percent difference less than 2 percent) when compared with three discrete hemispherical evaporation chamber (Stannard, 1988) measurements and seasonally variable errors (cumulative average errors on the order of 10 percent) when compared to Hamon equations and adjusted Penman equations. Comparisons of the estimated evaporation measured with the CFEP with eddy covariance and aerodynamic mass-transfer estimates indicated the CFEP estimates generally agree with the eddy covariance and aerodynamic estimates in the winter and greatly exceed them in the summer. Comparing daily and cumulative evaporation estimates from the CFEP with modified Penman and Hargreaves-Samani (Hargreaves and Samani, 1985) estimates show large variability day to day, with only about 50 percent of the daily data having differences of less than 0.04 inch per day [in/d], but good agreement on seasonal trends and excellent agreement over the period of record (RPD less than 2 percent; Reclamation, 2024). Without a clear determination of the inherent error in this method, a 10-percent error was assumed for the evaporation data and considered to be a conservative estimate, based on previous studies (Collison, 2019; Reclamation, 2021, 2023).

During the period, there were 159 days where the CFEP’s water-level change measurement was unusable because of logger issues, major rain events, and sensor malfunctions because of salt and (or) salt encrusting the fill pump. To gap-fill these days, a modified Penman evaporation estimation equation, equation 3.8 in Collison (2019), was fit to the CFEP’s estimated evaporation on a 15-minute time step:

Atmospheric data collected by the CFEP (air temperature, relative humidity, wind speed, barometric pressure, surface-water temperature, and net radiation) were used in the modified Penman equation. This equation was then fit to the CFEP’s estimated evaporation using an average monthly fitting coefficient (equaling for example, the average of three September monthly averages), with a total of twelve individual monthly fitting coefficients representing the entire year. The fitting coefficients are 0.33, 0.37, 0.43, 0.45, 0.57, 0.41, 0.55, 0.67, 0.47, 0.41, 0.39, and 0.34 for January to December, respectively. The daily evaporation estimates produced by the fitted modified Penman equation were used as the evaporation dataset in other analyses throughout this report. Estimated evaporation was lowest (average daily rates of 0.04 in/d) in the winter months of December and January and increased to an average daily rate of over 0.31 in/d in May until the summer monsoons reduced the vapor pressure deficit (fig. 7).

A review of the daily evaporation estimates with atmospheric data displays predicted responses in which increasing evaporation is associated with increasing temperatures and wind speed and decreasing evaporation is associated with decreasing temperatures and precipitation events (Reclamation, 2024; USGS, 2023). The daily estimated evaporation volume was computed as the product of the daily estimated evaporation rate and the average daily Zuni Salt Lake surface area.

Groundwater Temperature

In order to test whether vertical groundwater flow is occurring in the area, temperatures and depths in groundwater wells and springs were compiled and evaluated. Temperature as a function of depth can aid in identifying areas and structural features where groundwater has a downward or upward flow component. Temperatures measured as a function of depth in wells can be affected by regional heat flow, changes in surface temperature, and vertical movement of groundwater (Bredehoeft and Papadopulos, 1965). In the absence of groundwater flow, the geothermal gradient usually is represented by an increase of 25–50 degrees Celsius per kilometer of depth (°C/km) (Anderson, 2005). However, temperatures can also be disturbed by advective heat transfer through fluid flow.

Although no temperature profiles were collected for this project because of the lack of open wells, several groundwater temperatures associated with groundwater wells having reported depths were compiled and compared with data from a geothermal study of the area (Levitte and Gambill, 1980). The Levitte and Gambill (1980) study explored the geothermal potential in the area and estimated the thermal gradient near the Zuni Salt Lake to be about 35 °C/km. That study also identified areas with high thermal gradients to the north (60 °C/km) and to the west, in Arizona (>70 °C/km). The authors concluded that the high temperatures measured in Arizona, about 35 miles southwest of the study area, were due to warm water circulating up from depth.

Groundwater temperatures were compiled from NWIS (USGS, 2023) and are presented in relation to well or borehole depth in figure 9. Because the maar represents a depression of about 300 ft, depths of wells and springs within the maar were normalized to the surrounding landscape by computing the difference between the feature elevation and the rim elevation (6,500 ft above NAVD 88). This assumption does not fully account for the seasonal conductive heating at the land surface in Smith Spring, East Shore seeps, Spring Box, and MW1, but does represent a depth in the absence of the maar. The temperatures in the surficial zone (depths less than 50 ft) are influenced by seasonal heating and cooling of the land surface (Anderson, 2005). The average annual temperature reported for the nearby town of Quemado, N. Mex. is 8.9 °C (the period of record is 1915–2005; WRCC, 2024). The higher temperatures relative to the average annual air temperature in the two shallow samples may reflect the season (spring) in which they were collected or indicate that more recharge is occurring in the warmer months of the monsoon. The groundwater temperatures associated with deeper wells are generally warmer than predicted by the estimated thermal gradient of 35 °C/km (Levitte and Gambill, 1980). It should be noted that the temperature of groundwater being discharged from a well or spring is not equivalent to a bottom hole temperature or temperature/depth profile, because the temperature of the discharging groundwater may represent fluid temperatures from depths above or below the reported well depth or multiple depths. That said, it is often the case that groundwater wells are completed in the most productive zones, and therefore these temperatures can provide an indication of subsurface temperature profile. The least-squares regression of the groundwater temperature as a function of the reported well or hole depth yields a temperature gradient of about 79 °C/km. The trend is largely driven by the two San Andres Limestone/Glorieta Sandstone aquifer wells, but one can be reasonably confident in the relative temperature to depth relationship because of the reported thickness of the formation and overlying confining Chinle Formation. This high temperature to depth relationship indicates a substantially higher thermal gradient than previously reported. It is unclear if this finding represents a localized shallow high temperature source or advective heat flux through fluid flow from depth. It is notable that the groundwater temperature in DKOT-17 is near the predicted temperature for the 35 °C/km gradient for its reported depth, whereas the two San Andres Limestone/Glorieta Sandstone aquifer wells are substantially warmer. This may be explained by the ability for upflow to the San Andres Limestone/Glorieta Sandstone aquifer but not to the Dakota Sandstone because of the presence of the Chinle Formation between them. Upflow may also be elevating the temperatures in groundwater samples collected from in and near the maar because of fractures associated with the volcanic intrusion that created the maar.

Noble Gas

Vertical groundwater flow can also be indicated by the noble gas content of groundwater (Ballentine and others, 2002). Deviations from atmospheric concentrations can aid in understanding the age, source, and flow conditions of groundwater. Recharge water, having atmospherically derived noble gases, percolates through the vadose zone, and eventually mixes with groundwater. Being in equilibrium with the atmosphere, the noble gas concentrations in the recharge water have a defined composition and can be assumed to be air equilibrated water (AEW). In the groundwater system, additional noble gasses may be produced from radioactive decay of uranium, thorium, and potassium in the Earth’s crust (known as radiogenic noble gases) or from the migration of primordial noble gases from the mantle in areas of continental extension or magmatic activity (Ballentine and others, 2002). In this study, we focused on concentrations and isotopic compositions of helium in groundwater. Additional helium in the groundwater is normalized to the atmospheric component by comparing helium-4 (⁴He) to neon-20 (²⁰Ne), the dominant isotopes in the critical zone (Ballentine and others, 2002), and expressed as the He/Ne ratio where the AEW He/Ne is 0.290. The process by which the additional helium is added to the groundwater can be delineated by the isotopic composition of helium. Although both processes (radiogenic helium production and addition of mantle-derived gases) increase the concentration of helium, ⁴He is produced by radioactive decay, whereas ³He is remanent primordial helium from the mantle. The ³He/⁴He is expressed as R and is compared with the AEW value (Ra; Ballentine and others, 2002). The mantle end member was assumed to have an R/Ra of 8±1 (Graham, 2002) and a He/Ne of 10,000 (Morikawa and others, 2008). The radiogenic end member was assumed to have an R/Ra of 0.02 (Oxburgh and others, 1986) and He/Ne of 10,000 (Morikawa and others, 2008).

Based on the work of Pinti and others (2019), a mixing line was developed between AEW and the radiogenic end member (fig. 10). This line, with increasing He/Ne ratio, may be thought of as representing increasing time since recharge; in other words, a sample with higher radiogenic ⁴He and therefore higher He/Ne is likely to be older than one with lower He/Ne. Although the radiogenic ⁴He production rate for the formations in the area are not known, the extremely low rates reported in other formations (Solomon, 2000; Ballentine and others, 2002) indicated that groundwater with elevated He has been in contact with aquifer material for at least 1,000 years. Hypothetical “old water” end members (fig. 10) were selected along the mixing line between AEW and the radiogenic end member and a new mixing line was developed between the “old water” end members and the mantle end member.

The East Shore seeps, the Cinder Cone well, DKOT-17, and the two San Andres Limestone/Glorieta Sandstone aquifer wells have a He/Ne ratio greater than 150, indicating that less than 1 percent of the helium content is from modern recharge (AEW; fig. 10). Additionally, the East Shore seeps, Cinder Cone well, and the two San Andres Limestone/Glorieta Sandstone aquifer wells have R/Ra values greater than 1, indicating that more than 12 percent of the dissolved helium is derived from a mantle source, whereas DKOT-17, with an R/Ra value of 0.46 only has about 4 percent mantle contribution. The difference in the amount of primordial helium between the DKOT-17 sample and the two San Andres Limestone/Glorieta Sandstone aquifer wells supports the previous interpretation of restricted vertical flow in Dakota Sandstone at DKOT-17 caused by the presence of overlying and underlying aquitards and adds evidence that there may be some deep exchange of fluid or gas in the Permian units beneath the Chinle Formation. The relatively high R/Ra values in the Cinder Cone well and East Shore seeps also support the concept that the intrusion that formed the maar may be acting as a pathway for deeper fluids (or gases) to make their way to the surface of the maar. Wells MVRD-24 and MW2 had about 34 percent helium from recent recharge and only 5 percent primordial helium, whereas the Spring Box had 8 percent helium from a mantle source, possibly indicating mixing between nearby Mesaverde Group aquifer water and deeper mantle-rich groundwater like the San Andres Limestone/Glorieta Sandstone aquifer or the East Shore seeps. The sample collected from DKOT-35 has a helium content very similar to AEW, indicating relatively recent recharge, whereas the helium content in MVRD-58 has only about 1 percent mantle helium, further supporting the interpretation that vertical groundwater flow is restricted through confining units beyond the maar.

Water Budget

Between August 2019 and March 2023, daily stage at the Zuni Salt Lake fluctuated between 6,221.11 and 6,222.72 ft above NAVD 88, which corresponds to lake storage ranging from 22.7 to 261.3 acre-ft. Lake storage typically began declining in March because of dry conditions and increasing temperatures and began recovering in June–July with the onset of monsoon moisture (fig. 11). Change in the lake stage was recorded on more than 60 percent of the days but did not exceed 0.02 ft/d, except during large precipitation events. The annual volume of direct precipitation falling on the lake was calculated to be 135 acre-ft in calendar year 2021 and 132 acre-ft in calendar year 2022, which is very close to the average annual volume of 135 acre-ft calculated from average annual precipitation at a local meteorological station and the area of the lake. The estimated annual evaporation from the lake was determined to be 627 acre-ft in 2021 and 742 acre-ft in 2022. Based on pan evaporation data from several nearby stations, Bradbury (1971) estimated annual evaporation losses of about 872 acre-ft.

Groundwater contributions were estimated using a water balance approach by comparing daily storage changes with inflow (precipitation) and outflow (evaporation). The resulting imbalance is assumed to be groundwater inflows. The daily balance is highly variable, ranging from −2.1 to 18.2 cubic feet per second (ft³/s) with an average of 0.7 and a median of 0.5 ft³/s indicating a normal distribution of data with a slight bias toward positive imbalances (fig. 12). In order to minimize the error associated with unaccounted inflows from surface runoff following precipitation events, a second water balance calculation was conducted, whereby the inflows and outflows around precipitation events were not included. The data were reduced by removing daily budget data on the day prior to recorded precipitation (because of daily averaging) and 2 days following the event to avoid response lags. This data reduction resulted in the elimination of 398 days (from 829 to 431) from the analysis. The resulting imbalance was far less variable, ranging between −1.0 and 3.1 ft³/s, but with a comparable mean and median of 0.6 and 0.5 ft³/s respectively, indicating some stability in the method. The uncertainty associated with this imbalance was estimated to be 22 percent by combining errors in the components of the water budget using the root sum of squares (20 percent for storage changes and 10 percent for evaporation).

The groundwater contributions of 0.61 ±0.13 ft³/s (441 ±94 acre-feet per year [acre-ft/yr]) to the Zuni Salt Lake is inferred from the mean budget imbalance. By visually inspecting the balance over time (fig. 12), the groundwater contributions appear to be relatively stable (30-day moving average generally being within estimated uncertainty) over the period of data collection (September 2020–January 2023) and over variable lake volumes. The increasing imbalance during the periods of falling lake levels may indicate higher groundwater discharge or may represent increasing uncertainty in the evaporation estimates.

Salt Budget

The disparity of dissolved solids concentrations in groundwater samples necessitates an understanding of the salt contributions to the Zuni Salt Lake in order to understand the end-member chemistry and resulting lake composition. Several authors report the existence of salt beds beneath the lake in Paleozoic units in the subsurface (Levitte and Gambill, 1980; Maxwell and Nonini,1977; Zillmer, 1973). Based on the concentration of sodium and chloride measured as a part of this study and the groundwater seepage estimated above, the amount of potential halite entering the maar from different groundwater sources can be estimated.

Near-Lake Water Levels

Within the maar, water levels between the Zuni Salt Lake, Cinder Cone pool, Cinder Cone well, Spring Box, and East Shore seeps are very similar and require survey-grade positions to determine the gradients. Assessment of these gradients can enhance understanding of whether the lake is the source of salinity to the surrounding features or these features are a source of water and salt to the lake.

Water-level comparisons between the Zuni Salt Lake and the Spring Box show that between summer 2019 and fall 2021 the groundwater elevation of the Spring Box was consistently about 3 ft above the lake. The variability in the Spring Box groundwater elevations could either be related to seasonal changes in the source groundwater or be the result of transpiration demands of nearby phreatophytes (fig. 13). During the period of investigation, the Cinder Cone pool elevation only varied by about 7 in. and was more than 6 in. higher than the lake elevation until the lake elevation rose sharply in July 2021; thereafter, the difference in water-surface elevation between the lake and pool remained within 3 in. (fig. 13). The Cinder Cone pool elevations likely respond similarly to the lake elevations, with increases resulting from precipitation events and decreases resulting from evaporation. However, because of the smaller surface area of the pool and the sheltering provided by high, steep walls, those responses are minimized. The consistently higher elevations observed in the Cinder Cone pool indicate additional groundwater seepage discharging directly into the Cinder Cone pool. The groundwater elevation in the Cinder Cone well varied by about 8 in. during the investigation and averaged about 5 in. above the Zuni Salt Lake elevation. That difference, however, decreased to less than 1 in. during the high lake levels in late 2019 and early 2020, and was lower than the lake elevation recorded in September 2021 (fig. 13). Lastly, the elevations of the East Shore seeps are not well constrained; however, observations of the seeps along the shoreline indicated flow into the lake.

On the basis of the observations just presented, each of the features discussed is inferred to be a potential source of water to the Zuni Salt Lake. Likewise, the evaluation of hydraulic potential and the possible evidence of advective flux does not exclude any of the regional aquifers as a possible source of water and salt to the lake.

Geochemistry and Isotopic Composition

Because of the numerous inflows to the Zuni Salt Lake and the successive precipitation and redissolution of minerals, simple mixing models with a few end members are not adequate to identify a source of water or salt to the lake. Numerous aqueous chemical concentrations and isotopic compositions were analyzed from the lake and surrounding aquifer samples in order to better characterize the potential inflow sources.

Major Ions and Bromide

The major ion composition of natural water in the study area can offer insights into groundwater mixing and geochemical evolution. The major ion composition of surface-water and groundwater samples collected from and near the Zuni Salt Lake are presented on a Piper diagram (fig. 14), which is a trilinear diagram designed to illustrate the hydrochemical facies of a water sample (Hem, 1985). The percentages on the axes of the diagram represent the relative abundance of ions in percent milliequivalents.

Examination of the major ions (fig. 14) illustrates the different hydrochemical facies between the Zuni Salt Lake and surrounding groundwater. Both the lake and Cinder Cone pool are dominated by sodium (Na⁺) and chloride (Cl⁻) ions, as are the samples collected the Cinder Cone well and East Shore seeps. Each of these samples contains over 100,000 mg/L of dissolved solids and is classified as hypersaline. There is relatively minimal abundance of Cl⁻ in all other samples, with most samples having proportionally more bicarbonate (HCO₃⁻). The anion abundance in the two San Andres Limestone/Glorieta Sandstone aquifer samples is approximately 50 percent HCO₃⁻ and almost 40 percent sulfate, likely because of the abundance of limestone and gypsum in the formation (Baldwin and Anderholm, 1992). The sample from MVRD-58 has a similar anion composition as the San Andres Limestone/Glorieta Sandstone aquifer samples, but with slightly more sulfate and slightly less Cl⁻, proportionally. Other than MVRD-58, anion composition in Cretaceous groundwater samples is dominated by HCO₃⁻. The Spring Box sample contains similar amounts (about 40 percent) of HCO₃⁻ and Cl⁻, possibly indicating mixing between the Cl⁻-dominated hypersaline waters and most waters in Cretaceous aquifers. Cations are dominated by Na⁺ in all the samples except for the two San Andres Limestone/Glorieta Sandstone aquifer wells and the DKOT-35 well. Given the similar anion compositions in many of the groundwater samples and the reported presence of clay, cation exchange may be responsible for the distinct facies. The elevation of Smith Spring coincides with the Mesaverde Group aquifer; therefore, its major ion composition would be expected to resemble samples from the Mesaverde Group aquifer, which is supported by its plotting near samples from Mesaverde Group wells on the Piper diagram.

In Na⁺-Cl⁻ type brines, like those present in the Zuni Salt Lake and Cinder Cone pool, the enrichment of Na⁺ and Cl⁻ in solution through mineral precipitation processes can be difficult to distinguish from the addition of Na⁺ and Cl^- by halite dissolution. Relations between bromide (Br⁻), Na⁺, and Cl⁻ plotted as combinations of isometric log ratios can be helpful in determining which of these processes is occurring (fig. 15; Davis and others, 1998; Engle and Rowan, 2012). Br⁻ is more soluble than both Na⁺ and Cl⁻; therefore, halite precipitation has the effect of increasing the amount of Br⁻ in solution relative to Na⁺ and Cl⁻. As Br⁻ is excluded from the crystal lattice of halite, halite dissolution increases the content of Na⁺ and Cl⁻ in solution relative to Br⁻ while driving the Na⁺/Cl⁻ ratio of the solution towards unity (zero for isometric log ratios). The Na⁺/Cl⁻ ratio and Na⁺ and Cl⁻ enrichment relative to Br⁻ in hypersaline groundwater at the site (Cinder Cone well and East Shore seeps) indicates that the surrounding groundwater at the site could not have geochemically evolved to the composition of hypersaline groundwater without dissolving halite. Furthermore, the Zuni Salt Lake and Cinder Cone pool compositions have a similar Na⁺/Cl⁻ ratio (near unity) as hypersaline groundwater, although the surface waters exhibit lower ratios of Na⁺ and Cl⁻ to Br⁻, supporting the interpretation that the lake and Cinder Cone pool are sufficiently evaporated to precipitate halite, and that the source is hypersaline groundwater.

Stable Water Isotopes

Samples of the Zuni Salt Lake and possible groundwater sources were collected and analyzed for the isotopic compositions of the stable isotopes of water (δ²H and δ¹⁸O). The compositions were analyzed to constrain groundwater recharge and flow spatially and temporally. The variation in the isotopic composition of various potential recharge waters may help identify the source of recharge because of the mass-dependent fractionation of water isotopes resulting from temperature changes during the formation of precipitation and during evaporation prior to infiltration into the aquifer (Genereux and Hooper, 1998; Ingraham, 1998). This mass-dependent fractionation results because the various isotopic forms of water have different vapor pressures. Therefore, warmer-temperature precipitation contains a greater number of heavier isotopes (less negative δ¹⁸O and δ²H values) than cooler-temperature precipitation (Ingraham, 1998). Similarly, the water left after evaporation will contain a larger number of heavier isotopes than water that has not been subjected to evaporation. The stable-isotope ratios, reported in per mil (‰) relative to Vienna Standard Mean Ocean Water, of water from various sources are plotted in figure 16.

Groundwater δ²H and δ¹⁸O values in the study area plot below (to the right of) the Global Meteoric Water Line (GMWL) (Craig, 1961) (fig. 16). The overall position of the δ²H and δ¹⁸O values compared to the GMWL is typical of the arid and semiarid southwestern U.S. climate, where there is a large evaporation component to the water budget (Friedman and others, 1992). For reference, the isotopic composition of water collected from high mountain springs in the Zuni Mountains (Frus and others, 2020) are also plotted in figure 16. The spring samples are far more depleted in δ²H and δ¹⁸O than the water samples collected in and around the Zuni Salt Lake and indicate that recharge to the groundwater surrounding the lake is from warmer, lower elevation precipitation. The study samples are further classified in figure 16 using the noble gas concentrations reported in the “Noble Gas” section herein as a proxy for relative groundwater age. Lower He/Ne ratios indicate limited addition of helium from mantle or radiogenic sources and therefore shorter residence times relative to groundwater with higher He/Ne ratios. Groundwater samples were classified in terms of their He/Ne ratio being either greater than 2 or less than 2. The least-squares regression of groundwater isotopic compositions with low He/Ne ratios near the Zuni Salt Lake yields a strong correlation (R² = 0.99) and, for the purpose of comparison, is considered a local water line. The position of the line relative to the GMWL and the smaller slope indicates considerable evaporation following precipitation events and before this water recharged the aquifer. Groundwater samples with higher He/Ne ratios plot both above and below this local water line, indicating large variations in recharge temperature and subsequent evaporation.

The Zuni Salt Lake and Cinder Cone pool samples are far more enriched in heavy isotopes than any other samples and are indicative of evaporative effects. The isotopic composition of these samples relative to the nearby groundwater and the spring samples in the Zuni Mountains (Frus and others, 2020) are within expected range of evaporative effects (slopes between 3 and 6) according to Coplen and others (2000) and close to a slope of 5 reported by Phillips and others (2003) for samples representing increasing evaporation from the Rio Grande, N. Mex. The isotopic compositions of the Zuni Salt Lake (and Cinder Cone pool) water indicate that evaporation of many of the nearby samples could result in an isotopic composition similar to the lake, and therefore each aquifer may be a plausible source of the lake’s water.

Another process that may affect the isotopic compositions of the groundwater samples near the Zuni Salt Lake is the so-called oxygen isotope shift. This shift was first interpreted by Craig (1961) to result from isotopic equilibration of water with rock minerals rich in oxygen with high δ¹⁸O values. Increasing temperature allows for increasing exchange between groundwater and the higher δ¹⁸O rock content, resulting in groundwater having higher δ¹⁸O without a corresponding change in δ²H. That appears to be the condition of the DKOT-17 sample. If this shift did occur, it would also indicate the recharge water was derived from a cooler source, such as higher elevation precipitation (similar to the Zuni Mountain springs) or much older recharge when global temperatures were cooler, as during the Pleistocene (Dam, 1995). Although the relatively cool groundwater temperature measured at DKOT-17 does not support the high temperature oxygen isotope shift, the enriched δ¹⁸O may simply result from longer residence times. The oxygen isotope shift could also be interpreted to be affecting the Cinder Cone well and East Shore seeps samples; however, it is also entirely plausible that evaporation is the dominant process. The isotopic compositions that, through evaporation, could yield the isotopic composition of the hypersaline samples (East Shore seeps and Cinder Cone well) include those similar to Spring Box, MVRD-58, and DKOT-35, or depleted compositions associated with cooler recharge like the springs in the Zuni Mountains. The strong evaporation signal of the Zuni Salt Lake and Cinder Cone pool and the variability of nearby groundwater samples do not indicate any one source of water to the lake nor to the hypersaline groundwaters.

Age Tracers

Radioactive isotopes are frequently used to approximate the age of groundwater. Tritium (³H) and carbon-14 (¹⁴C) are used to determine the amount of time that groundwater has been disconnected from the atmosphere. ³H is a radioactive isotope of hydrogen that is short-lived, with a half-life of 12.32 years (Lucas and Unterweger, 2000). During the late 1950s to the mid-1960s, testing of nuclear weapons raised the atmospheric concentrations of ³H hundreds of times above normal background concentrations (Plummer and Friedman, 1999). ¹⁴C has a half-life of 5,730 years and is produced in the upper atmosphere by the interaction of cosmic rays and atmospheric nitrogen. The cosmogenic ¹⁴C is rapidly oxidized and becomes incorporated into biological and hydrologic processes as carbon dioxide (¹⁴CO₂) (Plummer and others, 2012). Because of its long half-life, the amount of ¹⁴C in dissolved inorganic carbon (DIC) may be used to estimate groundwater ages ranging from 1 ka to about 30 ka (Clark and Fritz, 1997). ³H is incorporated as part of the water molecule itself, whereas the ¹⁴C is measured in DIC and can have multiple sources.

DIC in natural waters is most often the result of the dissolution of calcite and (or) dolomite, which is greatly enhanced by dissolving gaseous carbon dioxide (CO₂), usually from soil gas, to create a weak acid (Han and others, 2012). The stable carbon isotope composition of delta carbon-13 (δ¹³C) is a useful indicator of the sources of dissolved inorganic carbon in groundwater (Chapelle and Knobel, 1985) because of the different isotopic composition of CO₂, calcite, and other organic constituents that may become part of the dissolved phase in groundwater.

A graphical analysis of the δ¹³C and ¹⁴C data from samples collected at and near the Zuni Salt Lake (fig. 17; Han and others, 2012; Han and Plummer, 2016) was performed by comparing the sample compositions to the probable isotopic compositions of dissolved CO₂ equilibrated with soil gas (point A1, where δ¹³C is −18‰, from Truini and Longsworth [2003]), DIC in an open system in equilibrium with soil gas (point A3), DIC in a closed system in equilibrium with carbonates (point M”), and DIC when water has equilibrated with soil CO₂ and carbonates in a closed system (point O). Prior to carbonate dissolution, the DIC of the water (made up mostly of dissolved CO₂) will be near point A1 in figure 17.

The Zuni Salt Lake tritium content is close to the atmospheric level (between 5 and 9 TU; Eastoe and others, 2012) indicating that the lake water is in equilibrium with atmospheric levels through diffusive exchange and mixing with air-equilibrated precipitation and surface runoff entering the lake (fig. 17). The Cinder Cone pool tritium content is also considered modern but is below the atmospheric level. Because the sample was collected below the halocline, the lower tritium content of the Cinder Cone pool sample compared with the Zuni Salt Lake sample may indicate that there is minimal mixing between water above and below the halocline and that the diffusive processes at depth are slower. This hypothesis is also consistent with the inferred existence of seepage from deep groundwater directly to the Cinder Cone pool. Except for Smith Spring, all groundwater is premodern based on tritium content. The presence of some tritium in Smith Spring likely results from mixing of older groundwater and recent recharge from precipitation in and adjacent to the ephemeral stream channel near the spring.

The ¹⁴C and δ¹³C content of the Zuni Salt Lake samples is unique and has an uncertain origin. The δ¹³C in the lake is substantially enriched (5.9 ‰) above any potential groundwater source or calcite or CO₂ source (fig. 17) and indicates that the enrichment of δ¹³C is occurring within the lake. Processes that may be responsible for the enrichment of δ¹³C include continuous dissolution and re-precipitation of calcitic materials (Smith and others, 1975) and methanogenesis (Han and others, 2012). Both processes are viable at the lake and may be co-occurring. Calcite dissolution and re-precipitation can occur as freshwater inflows reduce the saturation indices and evaporation drives dissolved solid concentrations higher. Methanogenesis would have to occur at depth where oxygen has been consumed in the lake. Although observations of methane are not known, hydrogen sulfide observations made during this study and by Bradbury (1971) when lake sediments are disturbed confirm the presence of reducing conditions at depth.

The Cinder Cone pool, DKOT-35, and the two San Andres Limestone/Glorieta Sandstone aquifer samples plot near the equilibration and exchange pathway between dissolved CO₂ in equilibrium with soil gas (point A1) and the composition of DIC in equilibrium with solid carbonates in a closed system (point M”; δ¹³C = −1.8‰ and ¹⁴C = 0 percent modern carbon) (fig. 17). The San Andres Limestone/Glorieta Sandstone aquifer and DKOT-35 samples represent different degrees of carbon exchange between DIC and solid carbonates in the aquifer, indicating that the exchange reaction may be limited at DKOT-35 (for example, if the availability of carbonate solids is limited) or that DKOT-35 is less geochemically evolved (and likely younger) than water at SAGA-9 and SAGA-10. The samples collected from the Cinder Cone well, the East Shore seeps, and the DKOT-17 well plot below the carbon exchange pathway, indicating additional radioactive decay at various stages of exchange between the soil gas and carbonates. Samples from Mesaverde Group aquifer groundwater wells, Smith Spring, and the Spring Box plot within, or very close to, the lower left quadrant of figure 17, which Han and others (2012) suggest may be due to additional carbon introduced by oxidation of old organic matter. Oxidation of organic fossil material would increase the DIC, lower the ¹⁴C, and lower the ¹³C value.

Sulfate Isotopes

Variations in the sulfur-34 (δ³⁴S) values of sulfate are caused by two kinds of processes: reduction of sulfate to sulfide by anerobic bacteria, which results in an increase in the δ³⁴S of the residual sulfate, and various kinds of exchange reactions, which result in δ³⁴S being concentrated in the compound with the highest oxidation state of sulfur (Mitchell and others, 1998). Depending on the reaction responsible for sulfate formation, between 12.5 and 100 percent of the oxygen in sulfate is derived from the oxygen in the environmental water, whereas the remaining oxygen comes from oxygen gas, that is in addition to any fractionation effects (Clark and Fritz, 1997; Taylor and others, 1984).

The Zuni Salt Lake sulfate isotopic composition is very similar to the isotopic composition of the East Shore seeps and the Cinder Cone well, with δ³⁴S around 15‰ and δ¹⁸O above 14‰ (fig. 18). The Cinder Cone pool is slightly more enriched in both δ³⁴S and δ¹⁸O, whereas the San Andres Limestone/Glorieta Sandstone aquifer samples, MW1 and MW2, and the DKOT-17 are slightly more depleted (around 11‰ and 10‰, respectively), but all plot as marine sediments or evaporites (Clark and Fritz, 1997) and close to Paleozoic marine Yeso Formation and San Andres Limestone water measured in southern New Mexico (Szynkiewicz and others, 2012). The δ³⁴S in groundwater samples from the Mesaverde Group aquifer range from 4.5 to 9.8‰. The Spring Box, Smith Spring, and a MVRD-24 plot within atmospheric values, as do MW1 and MW2. DKOT-35 and MVRD-58 plotted in the range typical of terrestrial evaporites, partially because the δ¹⁸O was below 2‰. (It is unclear why MVRD-58 is depleted in δ¹⁸O.) DKOT-35 had the only negative δ³⁴S composition. Based on the older age interpreted for the deep groundwater wells (SAGA-09, SAGA-10, and DKOT-17) and the more recent recharge interpreted in DKOT-35, the variability in the sulfate isotopic composition may be related to the equilibrium between the sulfate content at recharge, or at the surface, and contact time with the Mesozoic marine sediments. The range of sulfate isotope values in groundwater samples collected from the Mesaverde Group aquifer indicates differences in dominant sulfate source. However, enrichment of the sulfate isotopes by the reduction of sulfate to sulfide cannot be ruled out, given the presence of organic carbon and reducing conditions.

Strontium

Differences in the isotopic composition of strontium (Sr) (⁸⁷Sr/⁸⁶Sr) arise from the parent mineral, and given sufficient time, the radiogenic production of strontium-87 (⁸⁷Sr) from the decay of rubidium-87 (⁸⁷Rb) (Bullen and Kendall, 1998). The half-life of ⁸⁷Rb is such that considerable variations exist in the present-day ⁸⁷Sr/⁸⁶Sr ratio in minerals. ⁸⁷Sr/⁸⁶Sr in minerals and rocks ranges from about 0.7 to greater than 4.0 (Faure, 1986). Aqueous Sr, derived from any mineral through weathering reactions, will have the same ⁸⁷Sr/⁸⁶Sr as the mineral largely because Sr isotopes are not measurably fractionated by geological processes, unlike the isotopes of the light elements. Therefore, differences in ⁸⁷Sr/⁸⁶Sr among waters require either differences in mineralogy along the contrasting flowpaths, or differences in the relative amounts of Sr derived from the different minerals.

Because of the ionic and size similarities between calcium and strontium, there is an expected linear relation between the concentrations of each ion in the samples collected in and near the Zuni Salt Lake (fig. 19). The range in Ca/Sr values results from relatively high calcium concentrations in the San Andres Limestone/Glorieta Sandstone aquifer, a very low strontium concentration in the Spring Box sample, and a relatively high strontium concentration in the lake sample. The ratios of ⁸⁷Sr/⁸⁶Sr of the samples at and near the lake are considerably depleted in the heavy isotope and have little variability, indicating the isotopic composition of Cretaceous rocks and the groundwater moving through them are indistinguishable in the area. The notable exceptions are the two samples from the San Andres Limestone/Glorieta Sandstone aquifer wells, which have ⁸⁷Sr/⁸⁶Sr values at the higher end of most marine materials and are more commonly associated with weathering of continental rocks (Capo and others, 1998). The distribution of ⁸⁷Sr/⁸⁶Sr values in the lake and the groundwater samples indicates large contributions of strontium from nearby Cretaceous aquifers and minimal contributions from groundwater resembling the two San Andres Limestone/Glorieta Sandstone aquifer samples.