Distinct stores and the routing of water in the deep critical zone of a snow-dominated volcanic catchment

2.1 Site description

The Jemez River Basin Critical Zone Observatory (JRB-CZO) within the Valles Caldera National Preserve (VCNP) is situated in the Jemez Mountains in northern New Mexico, northwest of Albuquerque (Fig. 1a). This region is located in the transition zone between the snow-dominated Rocky Mountains and the deserts of the southwestern United States dominated by the North American monsoon (NAM) (Broxton et al., 2009). The JRB-CZO is in a montane, continental, sub-humid-to-semiarid climate characterized by a bimodal precipitation pattern (Zapata-Rios et al., 2015b). The VCNP is in a 21 km wide caldera that formed approximately 1.25 Myr ago (Bailey et al., 1969; Self et al., 1986). Ongoing volcanic activity as recent as 40 kyr ago caused the emplacement of several resurgent domes throughout the caldera (Wolff et al., 2011) of which Redondo Peak (3432 m above sea level, m a.s.l.) is the largest. The JRB-CZO comprises several headwater catchments that drain different aspects of Redondo Peak. The geology of Redondo Peak is characterized by several faults and is dominated by Pleistocene Bandelier Tuff, rhyolite, and andesitic rocks (Broxton et al., 2009; Vazquez-Ortega et al., 2015) that were intermixed by collapse breccias in some locations (Hulen and Nielson, 1991), which created a highly heterogeneous and complex geology (Fig. 1b).

La Jara catchment and the mixed-conifer zero-order basin, which is nested within the headwaters of La Jara, drain the eastern side of Redondo Peak (Fig. 1b). La Jara catchment ranges in elevation from 2702 to 3429 m a.s.l. with a mean slope of 15.7^∘ and drains an area of 3.66 km² (Perdrial et al., 2014). The ZOB consists of SE- and SW-facing slopes that are separated by a convergent zone and drains 0.15 km² (Vazquez-Ortega et al., 2016). The convergent zone of the ZOB, just above the ZOB flume, is characterized by boggy land in which standing water is present during the wettest parts of the year. Further evidence of near-surface saturation in this area is the presence of marshy plants in standing-water areas (e.g., broad-leaf cattails (Typha latifolia), Rocky Mountain irises (Iris missouriensis), and skunk cabbage (Veratrum californicum)), which contrast with those in the surrounding upland area of the ZOB (mostly bunch grasses (e.g., Deschampsia, Festuca, Koeleria, Muhlenbergia, and Poa) and aspen (Populus tremuloides) regrowth following a wildfire in 2013). Surface water from the ZOB flows through a Parshall flume before exiting the ZOB and then drains into La Jara creek.

Daily average precipitation and daily average temperatures were recorded at the Redondo Peak Weather Station (35.8839^∘ N, 106.5536^∘ W, 3231 m a.s.l.) maintained by the Desert Research Institute's Weather Regional Climate Center. Snow water equivalent (SWE) values were measured at the nearby Quemazon snow telemetry (SNOTEL) site (35.9167^∘ N, 106.4000^∘ W) located at a similar site elevation (2896 m a.s.l.) in the Santa Fe National Forest about 5 km east of Redondo Peak. Cumulative precipitation depths are summed over 1 water year (from 1 October to 30 September of the next year).

2.2 Groundwater well completions

Nested groundwater monitoring wells were drilled to total depths ranging from 6.7 to 47.2 m below the ground surface (m b.g.s.) at each of three locations within the ZOB in June 2016 (Figs. 1b and 2). At each location, wells are separated laterally by no more than 2 m from the next well so that the well casings stand in a line. A LiBr tracer was mixed with fluids injected during the drilling and well development process. Br⁻ concentrations were measured in initial samples to ensure that drilling fluids were flushed from monitoring wells prior to analysis of groundwater samples. All screened sections of the well casing have 0.051 cm slotted intervals. All wells were drilled to depth with a diamond-impregnated TT coring drill bit at an HQ3 diameter (9.6 cm annulus diameter; Table S1 in the Supplement). The annuli of all wells were packed with 8–12 Colorado silica sand (CSS) sand surrounding screened intervals, and the sand packs were sealed with hydrated bentonite pellets.

Well 2D is situated above the deeper site 2 wells with approximately 20 m separation between the base of well 2D and the maximum water table of the next deepest well (well 2C) with a hydraulic head gradient between wells 2D and 2C slightly greater than unit gradient at 1.1 m. Without a fully screened well between the depths of 2D and 2C (Table S1), we do not have direct evidence of an unsaturated zone between those wells, and we did not measure the saturation of cores collected during drilling because of the use of drilling fluids. However, cores collected during drilling provide evidence of redoximorphic features between the two wells that suggest the shallow well 2D is likely a perched aquifer.

2.3 Water quality and age tracers

Differences in well casing diameter (Table S1), the depth of water column, sampling frequency, and seasonal site accessibility necessitated different sampling collection among wells. Groundwater samples were collected from site 1 wells (2.54 cm diameter casing) using a Waterra inertial pump (Waterra USA Inc., Peshatin, WA, USA). Groundwater samples from site 2 wells (5.08 cm diameter casing) were collected with a Geotech SS Geosub Pump (Geotech Environmental Equipment, Inc., Denver, CO, USA) except during times when snow accumulated leaving the site inaccessible by vehicle during which times samples were collected with a 42.1 mm stainless-steel bailer (Geotech Environmental Equipment, Inc., Denver, CO, USA). Approximately three borehole volumes were discarded before collecting samples from each well to ensure that formation water was retrieved. Surface water samples of La Jara and ZOB stream water were collected at flume sites as grab samples.

Groundwater and surface water samples were collected in acid-washed polypropylene bottles (for cation analysis) and deionized-water-washed (DI-washed), combusted amber glass bottles (for carbon content and δ¹⁸O and δ²H analysis). Bottles were triple rinsed with sample water prior to sample collection, and samples were kept cold until they were promptly filtered at the University of Arizona. Sample splits for cations and trace elements were filtered through 0.45 µm nylon filters and acidified with trace metal grade nitric acid, while all other splits were filtered through 0.7 µm glass fiber filters. Surface water samples were collected biweekly during dry seasons and twice daily at the highest and lowest daily discharge using an automatic sampler (Teledyne ISCO, NE, USA) during the spring 2017 snowmelt. Groundwater samples were collected from all wells producing water on a biweekly basis during the dry seasons and daily from the shallowest well (well 2D; total depth 6.7 m b.g.s.) during the spring snowmelt using an autosampler (Teledyne ISCO, NE, USA).

All water samples were analyzed for cations, dissolved inorganic carbon (DIC), and δ¹⁸O and δ²H. Cations were measured by inductively coupled plasma mass spectrometry (ICP-MS) at the University of Arizona Laboratory for Emerging Contaminants (ALEC). Dissolved inorganic carbon was measured with a Shimadzu TOC-VCSH carbon analyzer in ALEC. δ¹⁸O and δ²H were measured with a DLT-100 laser spectrometer at the University of Arizona. Analytical precision (1σ) for all samples was 0.4 ‰ for δ²H and 0.1 ‰ for δ¹⁸O. Stable isotope data from several previous studies at the JRB-CZO and surrounding locations were also incorporated (Vuataz and Goff, 1986; Longmire et al., 2007; Gustafson, 2008; Broxton et al., 2009; Zapata-Rios et al., 2015b).

The tritium analysis was completed at the University of Arizona Environmental Isotope Laboratory on select groundwater and surface water samples filtered through 0.45 µm nylon filters. Mean residence time (t) was calculated using the two most common lumped parameter models, the piston flow and exponential models, according to Eqs. (1) and (2), respectively (Zuber and Maloszewski, 2000; Manga, 2001; Suckow, 2014; Zapata-Rios et al., 2015b).

where C is the measured tritium concentration (reported as tritium units, TU) in groundwater, C₀ is the measured TU in local precipitation, and λ is the tritium decay constant ($\mathrm{5.576}\times {\mathrm{10}}^{-\mathrm{2}}$ yr⁻¹) based on a tritium half-life of 12.43 years. The tritium analysis of precipitation in Albuquerque, NM (70 km from the study site), was measured monthly from 1962 to 2005 and is available as part of the Global Network of Isotopes in Precipitation (GNIP) through its Water Isotope System for data analysis, visualization, and Electronic Retrieval (WISER) database. Eastoe et al. (2012) found that volume-weighted mean (VWM) tritium concentrations in Albuquerque precipitation have remained stable since the early 1990s and are similar to local prebomb atmospheric tritium concentrations. Zapata-Rios et al. (2015b) used SWE data from the previously described Quemazon station to calculate a VWM of background tritium (8.6 TU from 1992 to 2005) that is used as C₀ in mean residence time calculations herein. The uncertainty of age calculations was computed according to Eq. (3) following Bevington and Robinson (1992), Scanlon (2000), and Zapata-Rios et al. (2015b).

where σ_t is the uncertainty of the age calculation, ${\mathit{\sigma }}_{c_{\mathrm{0}}}$ is the uncertainty of the background tritium concentration calculated as the mean of all lab-reported errors of tritium measurements reported for Albuquerque, and σ_c is the uncertainty (lab-reported error) in the measurement of tritium in samples in this study. It is important to note that the calculated residence times are biased towards shorter times because tritium in groundwater with a residence time greater than 25 years would have entered the groundwater system before the time frame of constant tritium levels over which the VWM tritium input was calculated. Therefore, the actual tritium input was likely higher than the VWM tritium input used in residence time calculations, and calculated residence times are likely underestimates.

Groundwater samples from wells 2D and 2C were collected with HydraSleeve^™ (GeoInsight, Las Cruces, NM, USA) passive collector bags in February 2018 for ¹⁴C analysis. Self-sealing bags were deployed in the wells and left to equilibrate for 24 h before sample retrieval. The radiocarbon analysis was completed at the University of Arizona Accelerator Mass Spectrometry (AMS) Laboratory on DIC in groundwater samples. Measured ¹⁴C activities are expressed as percent modern carbon (pmC) that were calculated as weighted averages of combined machine runs to reduce overall error, and δ¹³C_DIC measured values are expressed as per mil.

Uncorrected radiocarbon ages were computed following Eq. (4) for radioactive decay using a ¹⁴C half-life of 5730 years and an initial ¹⁴C activity (A₀) of 100 pmC, according to data from the neighboring Española Basin bound to the west by the Jemez Mountains (Manning, 2009). In order to calculate corrected radiocarbon ages, the δ¹³C_DIC value of the recharge water (−16.13 ‰) was first calculated based on an assumed temperature of recharge of 9.6 ^∘C (the average temperature of soil water from the ZOB over the last year of consecutive measurements, 2011), a pH of 7.51 (the average pH of soil water from the ZOB over the last year of consecutive measurements, 2011), and a ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{{\mathrm{CO}}_{\mathrm{2}}}$ of the ZOB vegetation of −24.7 ‰ (unpublished data from Tjasa Kanduc) using equilibrium constants of the carbonate system (Drever, 1997), alpha values between CO_2(g) and CO_2(aq) according to Vogel (1970), and alpha values between ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ and CO_2(aq) and between ${\mathrm{HCO}}_{\mathrm{3}}^{-}$ and ${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ according to Mook et al. (1974).

Corrected ages were calculated using the δ¹³C mixing model of Eqs. (5) and (6) following Pearson (1965), Pearson and Hanshaw (1970), and Clark and Fritz (1997) using an assumed δ¹³C_DIC of calcite of 0 ‰, a calculated value of the δ¹³C_DIC value of the recharge water, and A₀ of 100 pmC.

2.4 Streamflow and groundwater depths

Streamflow measurements were recorded at 15 min intervals from pressure transducers inside the stilling well of a Parshall flume at the base of the ZOB at 2996 m a.s.l. and La Jara catchment at 2739 m a.s.l. (flume locations shown in Fig. 1b). Transducer data were calibrated by depth measurements taken by hand at the time of each sampling event.

All monitoring wells were installed with vibrating wire piezometers (VWP) during the drilling campaign to provide nearly continuous (15 min) measurements of the hydraulic head in each well (Table S1). In this study, time series of groundwater depths are only shown from two wells (well 1A and well 2D) with continuous monitoring via VWPs. Electronic sounder measurements of the depth to water were taken before every groundwater sampling event, converted into the depth of the water column above the VWP, and used to calibrate the VWP and transducer measurements.

2.5 Saturated hydraulic conductivity estimates

Nearly continuous monitoring of the hydraulic head in groundwater wells enabled individual sampling events to be treated as rising-head bail-down slug tests (Butler, 2015). The AQTESOLV software program (Duffield, 2007) was used to curve fit the aquifer response (VWP measurements of depth of water converted into a normalized change in water table height against the logarithm of elapsed time) using the Kansas Geological Survey (KGS) model (Hyder et al., 1994) for three sampling events in wells 1A and 2D to estimate saturated hydraulic conductivity (K_sat) values using the AQTESOLV automatic estimation procedure. Early time recovery (for the first 15 min after a disturbance) was missed between time intervals; so the largest normalized head displacement is approximately 0.8 (m m⁻¹). The loss of early time data (beginning H∕H₀ of 0.8) from the Pratt 4-2 slug test data (Butler, 2015; data available from http://www.aqtesolv.com/examples/uncslug1.htm, last access: 1 November 2018) produced a 1.07 % error in K_sat values. The slug test analysis assumed isotropy.

2.6 Volumetric water content in soil pedons

Decagon EC-5 soil moisture sensors measured volumetric water content (VWC) in six instrumented soil pedons within the ZOB (Fig. 1) ranging from depths of 9.5 to 65 cm below the ground surface. Decagon soil moisture sensors measure the dielectric constant of soil, which is a sensitive measure of water content because of the much higher dielectric constant of water compared to that of soil and air. Dielectric constant measurements in mV are converted to VWC in m³ m⁻³ according to the Decagon EC-5 factory calibration equation. The VWC was measured at multiple depths in four of the six pedons. The VWC was recorded every 10 min.

2.7 Downhole neutron probe surveys

Water content profiles with depth were determined from neutron probe (model 503DR, Campbell Pacific Nuclear, Concord, CA, USA) surveys within the vadose zone in the top 18.3 m b.g.s. in wells 2B and 2C and the total depth of well 3B (12.9 m b.g.s.) over four different events. Raw neutron counts were recorded by the detector using a 32 s interval. Measurements were recorded every 0.3 m and a minimum of three readings were taken at each depth. Standard counts were measured in an acrylic sleeve before and after measurements of each well and were used to correct for radioactive decay of the Am-241 : Be source. Wet and dry calibrations were completed by measuring neutron counts within a 55-gallon (208 L) drum filled with Viton sand surrounding a PVC casing (same material as well casing). The calibration between neutron counts and water content is generally linear up to water contents of 0.4 (Rempe, 2016), which is greater than the maximum water contents measured here; however, neutron probe counts are sensitive to changes in bulk density and variable solid-phase chemical composition (Gardner et al., 2000; Rempe, 2016). Because of the highly heterogeneous geology and mineralogy in the JRB-CZO ZOB (Moravec et al., 2018), one calibration material (Viton sand) was used for all wells and all depths, which limits the interpretation of neutron probe surveys to comparison over time at the same site.

2.8 Fracture traces and density estimates

Downhole images captured with an optical televiewer at the wells of sites 1 and 2 were used to identify fracture traces and calculate fracture density using the fracture pattern quantification (FracPaQ) MATLAB^™ toolbox available as open-source code via MathWorks^™ File Exchange (Healy et al., 2017). The near-surface borehole instability required steel casing from the surface down to approximately 15 m b.g.s., which necessitated that downhole images, and thus fracture density characterization, started at approximately 15 m b.g.s. Because of the variable mineralogy and corresponding color transitions in the boreholes, it was necessary to trace fractures seen in optical televiewer images of wells 1A and 2A onto a new layer in Adobe Illustrator^™. The scalable vector graphics (SVG) file of the polyline fracture traces without the underlying downhole image was used as the input file for FracPaQ (Healy et al., 2017), which extracted the x and y fracture-trace coordinates from the file. The fracture density is defined as the number of fractures per unit area and has the unit of m⁻² (Dershowitz and Herda, 1992). FracPaQ (Healy et al., 2017) calculates fracture density from the fracture segment network as m∕2πr² according to Mauldon et al. (2001), where m is the number of trace endpoints in a circle with a radius r by resampling the exposed traces using circular scan windows that eliminate orientation, censoring, and length biases.

3.1 Temporal analysis of hydrographs and soil moisture

The temporal analyses of the climatic parameters (daily precipitation and temperature, Fig. 3a, and snow water equivalent, Fig. 3b) that drive streamflow response in ZOB and La Jara surface waters (Fig. 3c and d, respectively) and groundwater from the fractured tuff (well 1A, Fig. 3e) and shallow aquifers (well 2D, Fig. 3f) are used to examine seasonal hydrologic responses. The snow water equivalent peaked at 213 mm during WY 2017, while the annual cumulative precipitation was 673.5 mm for water year 2017. The average temperatures during the summer and winter months of WY 2017 were 13.4 and −4.2 ^∘C, respectively. As temperatures rose above freezing, the snowpack began to melt, and SWE values dropped rapidly (Fig. 3a and b). In response to snowmelt, daily averaged discharge of the ZOB flume peaked, while La Jara streamflow peaked. Temperatures dropped below freezing again causing the ZOB streamflow to freeze, producing no measured discharge, while La Jara streamflow reached a local minimum. During freezing temperatures, a second short-lived snowpack accumulated 18 mm of SWE. As temperatures increased above freezing again, maximum flows were reached in La Jara flume because streamflow remained high above baseflow conditions between snowmelt peaks, and a local streamflow maximum was reached in the ZOB flume. While streamflow peaks were greatest in response to the spring snowmelt, there were also obvious, smaller peaks in ZOB and La Jara surface waters following summer monsoons and fall storms. Both snowmelt and rainfall events influenced La Jara creek and the ZOB streamflow leading to increases above baseflow (Fig. 3c and d, respectively). However, the major driver of streamflow and groundwater depth response was spring snowmelt.

Figure 3(a) Daily averaged precipitation falls down from the top axis. Average daily temperature is shown in red when temperatures are above 0 and blue at or below 0 ^∘C. (b) Snow water equivalent (SWE) for WY 2017 reaches a maximum of 213 mm. (c) Discharge of ZOB flume in 15 min intervals for WY 2017. (d) Discharge of La Jara flume in 15 min intervals for WY 2017. (e) Well 1A depth of water below the ground surface from a vibrating wire piezometer placed 41.5 m below the ground surface in 15 min intervals. Drops in water depth correspond to sample collection. (f) Well 2D depth of water below the ground surface from a vibrating wire piezometer placed 6.4 m below the ground surface in 15 min intervals. Drops in water depth correspond to sample collection. The shaded region marks the timing of the North American monsoon (NAM) season from 15 July through 15 September.

Download

The depth of water in wells 1A and 2D (Fig. 3e and f, respectively) also responded to spring snowmelt, summer monsoons, and fall rainfall. Well 1A water depths peaked just 2 d after SWE values dropped to zero and on the same day that La Jara streamflow reached a local maximum. As temperatures froze again and a second snowpack developed, the depth of water from well 1A receded until 6 April 2017, after which it quickly rose to a second local maximum on 11 April 2017 4 d before ZOB streamflow reached a second max and 6 d before La Jara streamflow reached its second peak. The water depths from well 1A also peaked sharply on 2 October 2017 in response to fall precipitation. Conversely, the water depths from well 1A reached one gradual peak during NAM season on 16 August 2017, despite several smaller La Jara streamflow peaks in response to summer monsoon storms.

Water depths in well 2D (Fig. 3f) did not have pronounced, sharp maxima and minima in response to snowmelt dynamics like those seen in surface waters and site 1 groundwater. Instead, water depths in well 2D continued rising, with slight changes in rate on 21 March 2017 and 4 April 2017, while temperatures dropped below freezing again and a second snowpack accumulated and melted before well 2D water depths peaked on 22 April 2017. Well 2D water depths did not respond to summer monsoons and, instead, continuously decreased after their spring snowmelt peak until reaching a minimum after summer monsoons on 4 October 2017. However, the well 2D water table gradually increased in response to fall storms until reaching peak water depths on 25 October 2017. It appears that water slowly infiltrated into the shallow aquifer and recharged the near-surface groundwater store.

Figure 4Volumetric water content (VWC) of six soil pedons in the ZOB. Pedons are grouped spatially with pedons 4 and 3 on the western side of the ZOB (near the site 2 wells), pedons 5 and 2 in the convergent zone (near the site 3 wells), and pedons 6 and 1 on the eastern side of the ZOB (near the site 1 wells). The shaded region marks the timing of the North American monsoon (NAM) season from 15 July through 15 September.

Download

Table 1Estimated saturated hydraulic conductivity in meters per day for three sampling events from wells 2D and 1A and their mean. Estimates were made by curve fitting 15 min VWP data with the KGS model in AQTESOLV.

Download Print Version | Download XLSX

Volumetric water content (VWC) in ZOB soils from six pedons ranging from 9.5 to 65 cm b.g.s. also varied seasonally with generally higher VWC during spring snowmelt, lower VWC at the onset of NAM season, and intermediate VWC during fall storms (Fig. 4). VWC changes across the water year at the shallowest depth in pedons 4 and 3, located in the western area of the ZOB nearest site 2 wells, were small (0.20 m³ m⁻³ range in pedon 4 and 0.16 m³ m⁻³ range in pedon 3), while the VWC at a greater depth in pedon 3 increased drastically and remained elevated during spring snowmelt and fall storms. Changes in the VWC are likely most pronounced at depth due to subsurface lateral flow through more-saturated, more-conductive media, as suggested by previous studies in the JRB-CZO (Liu et al., 2008b; McIntosh et al., 2017; Olshansky et al., 2018). Changes in the VWC at pedons 5 and 2, located in the convergent zone of the ZOB nearest site 3 wells, were most pronounced during spring snowmelt. Again, an increased VWC in pedon 2 persisted longer at depth than at the near surface. While located in the same tuff rock type as pedon 5 situated in the ZOB's convergent zone, soils in the eastern area of the ZOB near site 1 wells at pedons 6 and 1 were located upslope of the convergent zone and did not remain wet for extended periods, but rather increases in the VWC of pedons 6 and 1 were flashy, while a decreased VWC during NAM season was sustained. Estimates of the saturated hydraulic conductivity of wells 1A and 2D were computed to explore differences in hydrologic properties of the two groundwater reservoirs following three sampling events that served as bail-down slug-out tests. Averaging K_sat from the three events produced a mean for well 2D of $\mathrm{7.22}\times {\mathrm{10}}^{-\mathrm{3}}$ and $\mathrm{1.22}\times {\mathrm{10}}^{-\mathrm{4}}$ m d⁻¹ for well 1A (Table 1), despite the fact that the fracture density of the deeper site 2 wells are estimated to be approximately 5 times less than the fracture density of site 1 wells (Fig. 5). However, it is important to note that fracture density could not be calculated across the water table depths of well 2D or site 1 wells because of surface instability and the presence of water inhibiting downhole televiewer images.

Figure 5(a) Fracture density of sites 1 (left) and 2 (right) in m⁻². Notice that fracture density is approximately 5 times greater in site 1 than 2 according to the color scale. (b) Fracture traces of sites 1 (left) and 2 (right) show where fractures exist in meters below the ground surface.

Download

3.2 Neutron probe profiles

Soil moisture content values in two site 2 boreholes (wells 2B and 2C) and one site 3 borehole (well 3B) were estimated from a neutron probe soil moisture gauge that was run downhole on four dates. Due to the textural shifts and complexity of the mineral composition as a function of depth at each site and across sites (Moravec et al., 2018), water content estimates are used to qualitatively examine changes in water content with depth and over time within each respective borehole.

Figure 6Daily precipitation in mm from the Redondo Weather Station above profiles of the water content with the depth in meters below the ground surface for wells 2B, 2C, and 3B. Colors of the profiles correspond to the timing shown on the precipitation figure above. It is important to note that the elevation above sea level for the wells of sites 2 (3024 m a.s.l.) and 3 (2989 m a.s.l.) are different. The dotted lines outline the boundary of the coarse gravel-like layer noted in the site 2 well logs.

Download

The profiles of water content with depth below the ground surface in borehole 2B (Fig. 6b), reach a local maximum water content of approximately 0.32 cm³ cm⁻³ at 1 m b.g.s. during all events. At an increased depth, the water content recedes to an overall minimum of 0.15 cm³ cm⁻³ at 2.4 m b.g.s. in the June, August, and February measurements, while the October measurement recedes to only 0.23 cm³ cm⁻³ at 1.8 m b.g.s. before spiking to 0.36 cm³ cm⁻³ at 2.4 m b.g.s. and receding slightly to 0.35 cm³ cm⁻³ at 4 m b.g.s., below which depth changes in water content are consistent across all time series.

In borehole 2C (Fig. 6c), just 2 m away from borehole 2B, differences in water content profiles are also seen over time. In October, the water content increases from the surface until reaching an overall maximum of 0.4 cm³ cm⁻³ at 3.7 m b.g.s., whereas profiles from the other three events peak at 0.29 cm³ cm⁻³ at only 0.6 m b.g.s., recede to 0.08 cm³ cm⁻³ at 1.2 m b.g.s., increase gradually to 0.12 cm³ cm⁻³ at 3.4 m b.g.s,. and drastically jump up to a range of water content from 0.34 to 0.39 cm³ cm⁻³ (July to February, respectively) at 4 m b.g.s. From 4 to 15.5 m b.g.s., the water content profiles for all events are relatively consistent between 0.28 and 0.30 cm³ cm⁻³ (June and August are identical; October is 0.1 greater than the summer profiles; and February is 0.1 greater than October).

Finally, water content profiles with depth in borehole 3B (Fig. 6d) are nearly identical across all measurement events and show three peaks in the water content of 0.36, 0.33, and 0.25 cm³ cm⁻³ at 1.2, 8.5, and 11.6 m b.g.s., respectively.

3.3 Geochemistry

In a volcanic setting such as the Valles Caldera, silicate mineral weathering is the primary driver of stream water chemical fluxes (McIntosh et al., 2017), and larger concentrations of base cations have been found in waters with longer flow paths as mineral dissolution fluxes increase with increasing water transit times (Zapata-Rios et al., 2015b). The analysis of the quantitative mineralogy of cores collected during the June 2016 drilling campaign in the ZOB (Moravec et al., 2018) found that quartz, potassium feldspar, plagioclase, volcanic glass, and smaller percentages of mica are the primary minerals ubiquitous in the cores of sites 1 and 2. Smectite, iron oxides, illite, and magnesite, as well as diagenetically altered minerals like Ca-zeolites (clinoptilolite and mordenite), are present in smaller percentages within the top 15 m b.g.s. of site 1, along with some 2 : 1 clays from 15 to 17 and 20 to 26 m b.g.s. Ca-zeolites, smectite, illite, iron oxides, and trace talc and tremolite are also found throughout site 2 cores, as well as secondary minerals like calcite and illite in the top 15 m b.g.s. Greater percentages of 2 : 1 clays are also found throughout site 2 cores, especially from 12 to 16 and 22 to 30 m b.g.s. A previous analysis of saturation indices of ZOB groundwater found that well 2D shallow groundwater fluctuated between near equilibrium and oversaturation with respect to calcite (Olshanksy et al., 2018), while previous work by Zapata-Rios et al. (2015b) found that springs across the JRB-CZO were undersaturated with respect to calcite, albite, and sanidine; therefore, interaction with those minerals is expected to influence groundwater and surface chemistry in the ZOB and La Jara catchment.

Table 2Number of samples (n) and average concentrations of major ions and pH of surface waters and groundwaters. Standard deviations are shown italicized and in parentheses.

Download Print Version | Download XLSX

The primary groundwater cations from each monitoring well are Ca²⁺, Na⁺, and Mg²⁺ (Table 2); however, the percentages (in terms of  µeq L⁻¹) of each cation differ between sites (sites 1 and 2). Major ion concentrations also differ with depth between site 2 wells, while site 1 groundwaters (wells 1A and 1B) are geochemically similar to one another (Fig. S1). These differences in cation concentrations are used in further analyses to distinguish streamflow sources.

Figure 7Time series of Na⁺, Ca²⁺, and DIC over WY 2017 for surface waters from La Jara flume and ZOB flume and groundwaters from sites 1 and 2. Concentrations are plotted over La Jara discharge in grey to highlight temporal changes.

Download

A temporal analysis of major cation concentrations in groundwater and surface water over WY 2017 (Fig. 7) again shows clear separation of groundwater concentrations between the two sites. Ca²⁺ and DIC concentrations are the highest in the shallow aquifer (well 2D), where calcite is known to be present, and all site 2 groundwater Ca²⁺ and DIC concentrations are considerably greater than surface and site 1 water year-round (Fig. 7). Furthermore, Ca²⁺ and DIC concentrations are most variable in wells 2D and 2C, which change together in time, suggesting a connection between the two water stores. Finally, both Ca²⁺ and DIC concentrations of shallow groundwater increase simultaneously during the spring snowmelt at the time that La Jara streamflow increases above baseflow, which is likely from calcite dissolution. In contrast, Ca²⁺ and DIC concentrations of La Jara surface waters do not change markedly as streamflow increases, while the Ca²⁺ and DIC concentrations of ZOB surface waters decrease slightly during this time. Mg²⁺ concentrations are greatest in the deepest well (well 2A) and decrease with the decreasing depth at site 2, are lower still in site 1 groundwater, and are lowest in La Jara stream water (Table 2). Na⁺ concentrations of the shallow aquifer increase steadily at the onset of the snowmelt, while La Jara and ZOB surface water Na⁺ concentrations decrease slightly (Fig. 7). As the site 1 water table slowly recedes following the spring snowmelt, Na⁺ concentrations in well 1A increase, likely because of the weathering of feldspars. Unfortunately, the lack of site 1 groundwater data during the snowmelt makes it impossible to determine if the variability of site 1 groundwater chemistry resembles that of the surface waters; however, surface water concentrations of major ions are generally closer in magnitude to those of site 1 groundwaters throughout the remainder of the water year.

Figure 8(a) ${\mathrm{Ca}}^{\mathrm{2}+}/{\mathrm{Mg}}^{\mathrm{2}+}$ molar ratios compared to Na⁺ concentrations show that groundwater from well 2D has ${\mathrm{Ca}}^{\mathrm{2}+}/{\mathrm{Mg}}^{\mathrm{2}+}$ ratios that are chemically distinct from the other water stores, which suggests that deeper groundwater is more representative of streamflow than shallow groundwater. (b) ${\mathrm{Ca}}^{\mathrm{2}+}/{\mathrm{Mg}}^{\mathrm{2}+}$ ratios compared to DIC concentrations show that site 1 groundwaters are more representative of the contribution to streamflow.

Download

A comparison of ${\mathrm{Ca}}^{\mathrm{2}+}/{\mathrm{Mg}}^{\mathrm{2}+}$ molar ratios with Na⁺ concentrations (Fig. 8a) and DIC concentrations (Fig. 8b), which can differentiate between the weathering of silicate minerals rich in Ca²⁺ and Mg²⁺, also show distinct groupings of groundwater between sites and depths. The ${\mathrm{Ca}}^{\mathrm{2}+}/{\mathrm{Mg}}^{\mathrm{2}+}$ molar ratios of the shallow aquifer are greater than those of the deeper waters. The ${\mathrm{Ca}}^{\mathrm{2}+}/{\mathrm{Mg}}^{\mathrm{2}+}$ molar ratios of surface waters overlap in space with those of deeper groundwater. The DIC concentrations clearly differentiate between sites 1 and 2, and the surface waters plot with similar DIC concentrations like site 1 groundwater, which indicates that deeper groundwater from site 1 is more representative of streamflow in La Jara catchment.

3.4 Age tracer analysis

Oxygen (δ¹⁸O) and stable-hydrogen (δ²H) isotope values of groundwaters and surface waters are comparable to one another. In addition, δ¹⁸O and δ²H values of both groundwaters and surface waters are more similar to the lower isotope values of snow than to those of summer precipitation (Fig. 9), indicating that snowmelt is the dominant source of recharge to groundwaters and surface waters in the ZOB and La Jara catchment. The majority of samples plot along the local meteoric water line (LMWL; from Broxton et al., 2009), showing consistency in δ¹⁸O and δ²H values over time. However, several surface waters and a few samples of deep groundwater from site 1 wells and shallow groundwater from well 2D have higher δ¹⁸O and δ²H values that plot to the right of the LMWL, along an enrichment line with a slope of 3.4 (Fig. 9).

Figure 9Oxygen (δ¹⁸O) and stable hydrogen (δ²H) isotope values of surface waters from La Jara and ZOB flumes, groundwaters, and volume-weighted mean ranges of summer precipitation (precip, taken from Zapata-Rios et al., 2015b) and snow (taken from Gustafson, 2008). The global meteoric water line (GMWL, solid line, slope of 8.0) and local meteoric water line (LMWL, slope of 6.8, dashed line from Broxton et al., 2009) are plotted for reference. Surface waters and site 1 groundwaters plot along an evaporation trend (dotted line, slope of 3.4).

Download

Tritium was detected in groundwater from each sampled well, which indicates that there is a component of modern recharge to all groundwater stores (Manning, 2009). Tritium content from wells analyzed in June 2017 and February 2018 are within 2 standard errors of one another, indicating little difference between the tritium content of groundwater stores between the summer dry season and the winter. The highest tritium content (4.4 TU in February 2018) and therefore the shortest residence time waters (12 and 17 years according to the piston flow and exponential models, respectively) are those of the shallow aquifer, while the lowest tritium content (0.7 TU in February 2018) and longest residence times (45 and 202 years, according to the piston flow and exponential models, respectively) are from wells 2B and 2C (Table 3). As expected, there is more tritium present in the shallowest groundwater compared to the deeper waters from site 2 wells; however, the deepest site 2 groundwater from well 2A has more tritium than the shallower wells 2B and 2C. The differences in tritium content (Table 3) across similar depths from sites 1 and 2 (Fig. 2) indicate the presence of separate groundwater stores of water within the ZOB. Radiocarbon age calculations were computed for shallow groundwater from well 2D and groundwater from well 2C beneath the shallow aquifer based on ¹⁴C activity and δ¹³C of the DIC of the two wells. The ¹⁴C activity of the DIC from well 2D was 75.34±0.19 pmC, while the δ¹³C_DIC was −13.1 ‰, which corresponds to a corrected radiocarbon age of 621 years. The ¹⁴C activity of the DIC from well 2C was 60.02±0.17 pmC, and the δ¹³C_DIC was −12.4 ‰, which corresponds to a corrected radiocarbon age of 2050 years (Table 3). As expected, there is less-modern ¹⁴C-DIC in the deeper groundwater (well 2C) indicating longer residence times at greater depth.

Table 3Tritium content ± lab estimated error and mean residence time calculated with the piston flow and exponential models ± age calculation uncertainty of groundwater from wells 2D, 2C, 2B, 2A, and 1B from two low-flow sampling events (June 2017 and February 2018). Radiocarbon analysis (pmC of ¹⁴C±lab estimated error and δ¹³C) of dissolved inorganic carbon from wells 2D and 2C with uncorrected ¹⁴C age calculated with A₀ of 100 pmC and corrected.

Download Print Version | Download XLSX

4.1 Hydrologic response

Site 1 wells are situated in highly fractured welded tuff with a maximum fracture density of approximately 5000 m⁻². This contrasts with highly weathered volcanic breccia in the top 15 m b.g.s. at site 2 wells that corresponds with a caldera collapse breccia deposit overlying more-consolidated, less-fractured welded tuff with a maximum fracture density of approximately 1000 m⁻² at depth (Fig. 5). We hypothesize that the differences in subsurface structure, presence of a confining layer, and greater fracture density of site 1 wells compared to site 2 wells (Fig. 5) influences the different hydrologic response of the two sites (Fig. 3). Well 1A responds rapidly to snowmelt, reaching its first peak on the same day as La Jara streamflow's first peak, and it responds gradually to summer monsoon events (Fig. 3e). Well 2D has a muted response to snowmelt with the first inflection point occurring on the same day as the first ZOB discharge peak, while the first well 2D water table maximum occurs 45 d after the onset of snowmelt and does not respond to summer monsoon events at all (Fig. 3f).

Mean K_sat values from both wells (Table 1) are lower than K_sat estimates of the same Tshirege member of Bandelier Tuff from the nearby Los Alamos National Laboratory that range from $\mathrm{7.6}\times {\mathrm{10}}^{-\mathrm{2}}$ to 1.12 m d⁻¹ (Rogers and Gallaher, 1995; Smyth and Sharp, 2006) and $\mathrm{9.3}\times {\mathrm{10}}^{-\mathrm{1}}$ to 1.6×10¹ m d⁻¹ (Kearl et al., 1990; Smyth and Sharp, 2006). This may be, in part, because slug tests provide localized estimates of the hydraulic conductivity directly surrounding the screened interval of the well in contrast to pumping tests that provide larger scale volumetric averages of hydraulic properties. Butler (2015) notes that hydraulic conductivity estimates from slug tests should be viewed as lower bounds of the conductivity in the vicinity of the well. Furthermore, the degree of welding and presence of alteration may account for the discrepancies in conductivity estimates as Rogers and Gallaher (1995) noted that the degree of welding of the Bandelier Tuff was greatest closer to the Valles Caldera, the tuff's volcanic source and site of this study.

Surprisingly, the mean K_sat of the shallow aquifer (well 2D) is more than 1 order of magnitude greater than that of the deep well 1A ($\mathrm{7.22}\times {\mathrm{10}}^{-\mathrm{3}}$ and $\mathrm{1.22}\times {\mathrm{10}}^{-\mathrm{4}}$ m d⁻¹, respectively, Table 1). Many more fractures are present at site 1 relative to site 2 (Fig. 5a), which generates a fracture density approximately 5 times greater at site 1 wells than site 2 wells (Fig. 5b). Unfortunately, it is not possible to directly compare the fracture density across the screened intervals from which hydraulic conductivity estimates were made because downhole images could not be captured within either well at those depths, but we expect that the general trend of dense fractures at site 1 and few fractures at site 2 would persist. The lower mean K_sat of well 1A could suggest that fractures may be backfilled by mineral precipitates and weathering rinds decreasing the ability of fractures to act as preferential flow paths. It is noteworthy, however, that such a mechanism is not clearly supported by the downhole images or core samples; these images and cores do not show obvious evidence of backfilled fractures.

Despite unexpected K_sat differences between the two wells, the very similar shape and timing of the hydrograph at well 1A compared to those at La Jara stream and the ZOB surface water (Fig. 3) is a function of pressure propagation and indicates a hydraulic connection between the fractured aquifer system and streamflow (Sophocleous, 1991a; Welch et al., 2013). The rapid pressure pulse transfer between site 1 groundwater and the stream suggests that the fractured welded tuff aquifer system has low specific storage and high transmissivity. Worthington (2015) noted much more-rapid changes in head in low storage fractured bedrock aquifers compared to granular aquifers, and Sophocleous (1991a) found that hydraulic diffusivity is the major controlling factor of the extent and speed of pressure pulse propagation.

The rapid response of the welded tuff aquifer system (well 1A, Fig. 3e) contrasts with the muted response of the shallow aquifer (well 2D, Fig. 3f), which does not show evidence of pressure pulse propagation between shallow groundwater and the ZOB and La Jara surface waters. The slower response of groundwater levels in well 2D suggests that the shallow aquifer is not directly hydraulically connected to the stream likely because it appears to be separated from deeper groundwater by a low permeability confining layer and the significantly lower fracture density of site 2 wells (Fig. 5). The comparison of site 1 and 2 groundwater hydrographs to surface water hydrographs bears a striking resemblance to the juxtaposition of stream-flood wave propagation in monitoring wells drilled into buried river channels (similar to well 1A) and the absence of pressure pulses in monitoring wells not associated with buried channels (similar to well 2D) seen in the Great Bend alluvial aquifer in Kansas (Sophocleous, 1991a). Furthermore, the cubic shape of the rising water table in the shallow aquifer identified by the slow initial rise of groundwater levels followed by an inflection point and subsequent rapid rise in well 2D groundwater levels (Fig. 3f) is indicative of groundwater recharge (Sophocleous et al., 1988; Sophocleous, 1991a, b). Further analysis of the rates of increase before and after the well 2D hydrograph inflection point can be found in Olshansky et al. (2018).

The hydrologic response to incoming precipitation exhibits strong dependence on season, as indicated by differences following spring snowmelt, summer monsoons, and fall precipitation. Specifically, well 1A groundwater responds rapidly to spring snowmelt and fall precipitation, in contrast to the welded tuff aquifer's response to summer monsoon rain that is smaller and much more gradual, suggesting that spring snowmelt and fall precipitation induce a different hydrologic flow regime than that of summer monsoons (Fig. 3e). Different hydrologic flow regimes across seasons also exist in the shallow aquifer (well 2D). While all changes in the shallow water table are gradual, spring snowmelt and fall precipitation produce water table peaks in late April and October, while the shallow water table steadily decreased during summer monsoons indicating no water table changes induced by summer storms (Fig. 3f).

The less-pronounced hydrologic response to summer monsoons, compared to snowmelt and fall rain, in shallow groundwater of site 2 and deeper groundwater at site 1 is likely a function of increased evapotranspiration and thus, drier antecedent soil moisture conditions at the onset of the NAM season, as indicated by decreased VWC in shallow soils from six instrumented pedons in the ZOB immediately before monsoon storms began compared to wetter soils during the spring snowmelt and fall storms (Fig. 4). This agrees with previous studies in the VCNP, which found that soil moisture was lowest in early summer after soil moisture from snowmelt had receded, and it increased after the arrival of monsoon storms (Vivoni et al., 2008; Molotch et al., 2009). Furthermore, Zapata-Rios et al. (2015a) found that the NDVI (normalized difference vegetation index) in the VCNP increased during the NAM season suggesting that precipitation was partitioned to plant use during monsoon rains and was not available to recharge groundwater stores. Moreover, smaller, sporadic precipitation during summer monsoon storms increased evapotranspiration due to higher temperatures (Fig. 3a), and it increased plant use create a wetting and drying effect in shallow soils that can be seen as small fluctuations in the VWC (Fig. 4). This effect likely inhibits the infiltration of water into the deep subsurface and agrees with the model of Langston et al. (2015) of unsaturated zone flow in two seasonally snow-covered hillslopes in Colorado, which found that dry periods between low-intensity precipitation and infiltration events inhibited recharge because of the need for shallow soil to re-wet after each precipitation event.

4.2 Water routing in the unsaturated zone

Precipitation inputs differ before each neutron probe measurement (cumulative precipitation of 163.18 mm between June and August neutron probe measurements, 67.55 mm between August and October measurements, and 162.58 mm between October and February measurements; Fig. 6a); however, water content is nearly identical in the top 4 m b.g.s. in three of the four profiles for site 2 wells and in all site 3 data collection events (Fig. 6b–d). We hypothesize that the response seen at site 2 in the October survey is a function of both the slightly larger precipitation depth prior to this survey and the wetter shallow soil conditions preceding the survey as compared to the conditions preceding the other three surveys. Higher-frequency sampling during the wet season are needed to determine the impact of precipitation depth and the potential for precipitation thresholds to induce vertical flow. Perhaps temporal changes in water content with depth were missed because of the sporadic timing of neutron probe surveys due to the arduous transportation and permitting issues involved with the use of this instrument. While it does not appear that water content profiles with depth captured progressive enrichments in rock moisture as seen by Salve et al. (2012) and Rempe and Dietrich (2018), they do indicate that the minimum dry-season rock moisture storage is consistent across dry seasons, suggest differences in water routing, and identify lithologic discontinuities in the subsurface.

Neutron probe surveys show small shifts in water content with depth that are likely associated with small-scale heterogeneities in bulk density created by the lithologic discontinuities in volcanic collapse breccia deposits, variable degrees of ash consolidation, welding, and secondary mineral precipitation, which are evident in quantitative mineralogical analyses (Moravec et al., 2018). For instance, well logs showed that a thin layer of coarse gravel-like material underlies the soil at site 2 wells, which is expected to drain quickly (high permeability) and retain little water (low porosity). At the depth corresponding to that gravel-like layer from 1.5 to 2.3 m b.g.s., water content across all times recedes quickly (Fig. 6b) or remains constant (Fig. 6c) before increasing rapidly just above the water table of the shallow aquifer. Salve et al. (2012) also found that moisture content variation from neutron probe surveys in weathered argillite were strongly linked to changes in material properties with depth, which suggested different flow processes through the unsaturated zone. Here, layers of increased water content (like that from 9 to 11 m b.g.s. in well 2B) above layers of relatively lower water content (like that of 12 to 13.5 m b.g.s. in well 2B) are indicative of subsurface lateral flow through more-saturated, more-conductive media that can be seen in Fig. 6.

Evidence of vertical infiltration is also seen in the site 2 wells. The marked change in shape (increased water content from 1 to 4 m b.g.s.) of the October water content profile suggests vertical infiltration and subsequent recharge to the shallow aquifer. The analysis of the well 2D hydrograph (Fig. 3f) confirms that this October neutron probe survey was completed while the shallow water table was rising. Despite similar shallow water table depths for all four surveys (2.7 m b.g.s. on 27 June 2017, 2.9 m b.g.s. on 15 August 2017, 2.9 m b.g.s. on 12 October 2017, and 3.2 m b.g.s. on 6 February 2018), the October survey was the only survey that corresponds with a rising rather than receding water table (Fig. 3f) and is the only water content profile that captured the vertical infiltration of recharge to the shallow aquifer. Furthermore, the wet October profile returns to dry conditions within the top 4 m in February, and February water content beneath 4 m exceeds that of previous surveys, which suggests that water drains vertically at depths greater than 4 m b.g.s. in February.

Geologic maps of the Valles Caldera (Goff et al., 2011) indicate that a concealed fault bisects the ZOB, and it appears that site 3 wells were drilled immediately next to (possibly within) the fault zone, which coincides with the convergent outlet of the ZOB (Fig. 1). This is likely why site 3 wells did not produce water in the time period of this study (Fig. 2). Furthermore, water content is nearly identical across all four measurement events in borehole 3B (Fig. 6). Despite being located in a convergent zone subjected to seasonal wetland-saturated conditions at the land surface, neutron probe surveys do not indicate that water infiltrates vertically in site 3. Rather, data suggest that water moves laterally in the subsurface as indicated by the three zones of increased water content seen in Fig. 6. This is further supported by high clay content (up to 50 %) observed below 1.5 m depth at site 3 (Moravec et al., 2018), which likely impedes vertical infiltration and induces lateral flow in the subsurface.

4.3 Contribution of distinct groundwater stores to streamflow

In the same headwater catchment in the JRB-CZO, Olshansky et al. (2018) observed temporal changes in major ion concentrations of soil, surface, and groundwater during the spring snowmelt of 2017. They found pulses of high Si concentrations during the falling limb of surface water hydrographs that were hypothesized to result from increasing groundwater contributions during receding surface flows because surface water Si concentrations were similar to those of shallow groundwater (well 2D). However, surface water concentrations of other major ions (Ca²⁺, Na⁺, Mg²⁺, Mn, and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$) did not coincidentally rise despite higher shallow groundwater concentrations of those ions compared to Si, which suggests, instead, a deeper source of groundwater to La Jara streams. Herein, we present further evidence that the deep groundwater source to La Jara and the ZOB surface water is the deep welded tuff aquifer system found throughout the greater La Jara catchment and represented at site 1 wells.

High ${\mathrm{Ca}}^{\mathrm{2}+}/{\mathrm{Mg}}^{\mathrm{2}+}$ molar ratios of shallow groundwater (well 2D; Fig. 8) are attributed to calcite dissolution, which is only present in the shallow subsurface at site 2 and leads to increased Ca²⁺ concentrations but does not impact Mg²⁺ concentrations. The clear separation of ${\mathrm{Ca}}^{\mathrm{2}+}/{\mathrm{Mg}}^{\mathrm{2}+}$ molar ratios in shallow (well 2D) and deep groundwater (wells 1A, 1B, 2A, 2B, and 2C) and the overlap of surface water ${\mathrm{Ca}}^{\mathrm{2}+}/{\mathrm{Mg}}^{\mathrm{2}+}$ molar ratios with those of deeper groundwater (Fig. 8a) suggests that the shallow aquifer does not contribute substantial volumes to surface flow, but instead deep groundwater is the dominant source of streamflow.

The dissolution of calcite in the shallow aquifer also leads to higher DIC concentrations in waters in the ZOB (Table 2; Appelo and Postma, 2005). Higher DIC concentrations in all site 2 waters is consistent with some degree of vertical connection between site 2 wells that, in turn, suggests that the suspected confining layer beneath the shallow aquifer acts as more of an aquitard than aquiclude. The ${\mathrm{Ca}}^{\mathrm{2}+}/{\mathrm{Mg}}^{\mathrm{2}+}$ molar ratios and DIC concentrations of surface waters are very similar to those of site 1 groundwaters (Fig. 8b), which again indicates that deep groundwater from site 1 wells is more representative of groundwaters that contribute to La Jara stream. The welded tuff rock type (Fig. 1) of site 1 wells is also more representative of the geology, by volume, throughout the greater La Jara catchment.

Temporal analysis of major ion chemistry indicates that deep groundwater from the fractured tuff (site 1) sustains stream baseflow, as streamflow concentrations and trends in concentrations over time are consistent with site 1 groundwater concentrations. In contrast, pronounced changes in shallow site 2 groundwater (wells 2D and 2C) major ion chemistry are not reflected in streamflow concentrations over time, which suggests that shallow groundwater represents only a small volumetric contribution to streamflow. Furthermore, recent work by Olshansky et al. (2018) found that soil water was an important contributor to surface water during the spring snowmelt of 2017 and may explain the subtle trends in Ca²⁺ and DIC concentrations in surface waters, particularly ZOB surface water which was correlated with soil $P_{{\mathrm{CO}}_{\mathrm{2}}}$ concentrations at the onset of snowmelt.

4.4 Mix of old and young snowmelt-dominated waters

Snowmelt is the dominant source of recharge to all groundwater stores in this study, as the δ¹⁸O and δ²H values of surface water and groundwater plot much closer to the volume-weighted mean of snow (Gustafson, 2008) than to the VWM of summer precipitation (Zapata-Rios et al., 2015b) (Fig. 9). Furthermore, the detection of tritium in groundwater from each sampled well indicates a component of modern recharge to each groundwater store, which agrees with previous work that found springs surrounding Redondo Peak are composed of modern water (Zapata-Rios et al., 2015b). The presence of tritium in groundwater suggests that snowmelt slowly infiltrates into all groundwater stores (Fig. 10). However, radiocarbon analysis of wells 2D and 2C also indicates the presence of much older waters (corrected ages of 621 and 2050, respectively; Table 3) in these shallow groundwater reservoirs (Fig. 10). The detection of tritium and less than 100 pmC in these wells suggests a mixture of old and young waters (Bethke and Johnson, 2008; Jasechko et al., 2017), which is expected to persist in each groundwater store.

Figure 10Schematic of intricacies of the hydrologic structure and function across two contrasting sites within the ZOB (sites 1 and 2). Site 2 has multiple, separate stores of groundwater across depth that are distinct from each other and distinct from deep groundwater stores at site 1. All groundwater is recharged via snowmelt and seasonal differences in hydrologic response to precipitation inputs exist at both sites with a less-pronounced response to summer monsoons. There is evidence of vertical infiltration and subsurface lateral flow at site 2 and a mix of young and older waters, which are expected to persist across all groundwater stores. The fracture density at site 1 is approximately 5 times greater than at site 2, and the CZ structure and architecture of site 1 is most representative of the greater catchment. Deep groundwater from the fractured aquifer system at site 1 is hydrologically connected to streamflow, and site 1 deep groundwater chemistry is most representative of water contributing to streamflow, while the distinct chemical signature of shallow groundwater from site 2 is not evident in streamflow.

Download

Decreasing tritium content with depths from 2D to 2C to 2B (Table 3) agrees with previous studies that have suggested residence times increase with depth (Zapata-Rios et al., 2015b). This trend along with the distinct ${\mathrm{Ca}}^{\mathrm{2}+}/{\mathrm{Mg}}^{\mathrm{2}+}$ molar ratios and DIC concentrations from the presence of shallow calcite deposits in site 2 wells is consistent with a vertical connection between the three shallowest site 2 wells. However, the increased tritium content in well 2A (closer to that of site 1 than other site 2 wells) suggests that 2A is not vertically connected to the above site 2 wells, and, instead, it is likely laterally connected to younger waters upgradient.

The few samples that plot to the right of the meteoric water line may be related to evaporation and/or the presence of geothermal waters. The slope of the evaporation trend (m=3.4) is consistent with the slope of evaporation trends seen in arid environments like southern Arizona and northern Mexico (Gray, 2018; Zamora, 2018). The timing of the lower surface water δ¹⁸O and δ²H values is indicative of evaporative enrichment; however, the sporadic timing during surface water baseflow of the lower δ¹⁸O and δ²H values for all groundwaters is surprising (Fig. S2).

The second possible explanation of this enrichment trend is the presence of geothermal waters in the ZOB. The observed enrichment trend is consistent with that seen by Vuataz and Goff (1986) from geothermal waters sampled from lower elevations in the Valles Caldera and surrounding area. However, the analysis of other geothermal water indicators (pCO₂ and temperature) does not suggest the presence of geothermal waters in the ZOB. Furthermore, ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ and Cl⁻ concentrations (data not shown) are several orders of magnitude lower than those of known geothermal waters from Vuataz and Goff (1986).

In summary, we propose a schematic model (Fig. 10) to conceptualize the details of the hydrologic structure and hydrologic function at two contrasting hillslopes within the ZOB (sites 1 and 2). Site 2 has multiple, separate stores of groundwater across depth that are distinct from each other and distinct from deep groundwater stores at site 1. All groundwater is recharged via snowmelt and seasonal differences in hydrologic response to precipitation inputs exist at both sites with a less-pronounced response to summer monsoons at both sites linked to drier antecedent shallow soil moisture at the onset of NAM season. There is evidence of vertical infiltration and subsurface lateral flow at site 2 and a mix of young and older waters, which are expected to persist across all groundwater stores. The fracture density at site 1 is approximately 5 times greater than at site 2, and the CZ structure and architecture of site 1 is most representative of the greater La Jara catchment. Deep groundwater from the fractured aquifer system at site 1 is hydrologically connected to streamflow, and site 1 deep groundwater chemistry is most representative of water contributing to streamflow, while the distinct chemical signature of shallow groundwater from site 2 is not seen in streamflow. This study is not able to discern where the site 2 shallow groundwater drains, but it is possible that it transmits to streamflow in a quantity too low to detect because of dilution. The lateral extent of the shallow aquifer is not known and therefore, calculations of the volume of water in that aquifer cannot be made at this time.