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

This study combines major ion and isotope chemistry, age tracers, fracture density characterizations, and physical hydrology measurements to understand how the structure of the critical zone (CZ) influences its function, including water routing, storage, mean water residence times, and hydrologic response. In a high elevation rhyolitic tuff catchment in the Jemez River Basin Critical Zone Observatory (JRB-CZO) within the Valles Caldera National Preserve (VCNP) of northern New Mexico, a periodic precipitation pattern creates different hydrologic flow regimes during spring snowmelt, summer monsoon rain, and fall storms. Hydrometric, geochemical, and isotopic analyses of surface water and groundwater from distinct stores, most notably shallow groundwater that is likely a perched aquifer in consolidated collapse breccia and deeper groundwater in a fractured tuff aquifer system, enabled us to untangle the interactions of these groundwater stores and their contribution to streamflow across 1 complete water year (WY).

Despite seasonal differences in groundwater response due to water partitioning, major ion chemistry indicates that deep groundwater from the highly fractured site is more representative of groundwater contributing to streamflow across the entire water year. Additionally, the comparison of streamflow and groundwater hydrographs indicates a hydraulic connection between the fractured welded tuff aquifer system and streamflow, while the shallow aquifer within the collapse breccia deposit does not show this same connection. Furthermore, analysis of age tracers and oxygen (δ¹⁸O) and stable hydrogen (δ²H) isotopes of water indicates that groundwater is a mix of modern and older waters recharged from snowmelt, and downhole neutron probe surveys suggest that water moves through the vadose zone both by vertical infiltration and subsurface lateral flow, depending on the lithology. We find that in complex geologic terrain like that of the JRB-CZO, differences in the CZ architecture of two hillslopes within a headwater catchment control water stores and routing through the subsurface and suggest that shallow groundwater does not contribute significantly to streams, while deep fractured aquifer systems contribute most to streamflow.

1 Introduction

Understanding the interconnections of groundwater and surface water is fundamental to water resource management as groundwater and surface water should be considered a single resource (Winter, 1998); however, their interactions in different hydrogeologic settings are varied and complex (Winter, 1999). Discerning stream water sources and groundwater dynamics are even more important in the context of a changing climate, especially in the semiarid, mountainous environment of the western United States, where warming trends are expected to threaten water supply (Barnett et al., 2005). Furthermore, identifying compartmentalized groundwater stores is necessary to sufficiently account for all components of the water balance (McDonnell, 2017). Therefore, characterizing localized water stores and the hydrologic connection of those aquifers to streams in mountainous environments that act as water towers (Viviroli et al., 2007) has important implications for water resource availability of large population centers downstream.

The influence of the hydrogeologic environment (i.e., geology, topography, and climate) on the groundwater flow system of a given geographic region has long been accepted as the theoretical framework used to conceptualize groundwater flow. Building on the conceptual model of Toth (1970) and the understanding that one part of the framework informs our knowledge of the other, several studies have focused on topography (Beven and Kirkby, 1979; Woods et al., 1997; Hutchinson and Moore, 2000; Kirchner et al., 2001; McGuire et al., 2005) and rock type (Farvolden, 1963; Freeze, 1972; Kelson and Wells, 1989; Mwakalila et al., 2002) as controls of groundwater flow systems. However, subsurface heterogeneities, which can be abundant and are challenging to identify, can give rise to complex, localized groundwater stores, whose contribution to streamflow can be very difficult to discern. There is still much to learn about the extent to which structural heterogeneities exist and how, specifically, they control groundwater storage, routing, and contributions to streamflow. For instance, evidence of perched aquifers transmitting shallow subsurface flow has been shown across variable rock types (Salve et al., 2012; Kim et al., 2014, 2017; Brantley et al., 2017; McIntosh et al., 2017). Furthermore, Brooks et al. (2015) highlighted the need to understand the influence of subsurface structure on water routing and residence time as they concluded that surface water, across several catchments and flow regimes, substantially interacted with or spent time within various soil and groundwater reservoirs. In fact, the lack of subsurface characterization in high-elevation groundwater systems often forces hydrologic modeling studies to focus solely on shallow groundwater systems and make the generalizing assumption that hydraulic conductivity decreases with depth and deep fractured aquifers support little flow (Manning and Caine, 2007; Welch and Allen, 2014; Markovich et al., 2019).

Heavily instrumented and intensively studied sites, such as Critical Zone Observatories (CZOs), which are part of a network of field-based laboratories arrayed across a variety of rock types, land uses, elevations, and climates (Anderson et al., 2008), are ideal locations to examine the interplay between subsurface structure and function. Moreover, recent focus on characterizing the deep subsurface structure or architecture is beginning to elicit a deeper understanding of the role of weathering, lithology, and hydrology in the overall function of the critical zone (CZ) and its landscape evolution (Riebe et al., 2017). The CZ is the near-surface terrestrial layer of the Earth that spans from the tops of trees down to unweathered bedrock where water, rock, air, and life meet and interact (Brantley et al., 2006, 2007; Anderson et al., 2007; Chorover et al., 2007; Küsel et al., 2016). Herein, we use the term subsurface structure or architecture to refer to physical properties of the subsurface such as its lithology, fracture density, and location and the extent of geologic heterogeneities that may impact the movement of water through the subsurface. Understanding the coupling between CZ architecture, developed over geologic timescales, and CZ function on short event timescales is a primary goal of CZ science (Chorover et al., 2011; Brooks et al., 2015). In particular, there is limited knowledge about the structure of the deep CZ and its direct influence on water storage (Holbrook et al., 2014; Dralle et al., 2018) and routing, mean residence times (McGuire et al., 2005), and streamflow sources. Integrated studies that simultaneously examine both CZ architecture and CZ hydrology through hydrometric, geophysical, geochemical, and residence time analyses are needed to understand the distribution of groundwater stores, their connection to streamflow, and the underlying impact of CZ architecture on hydrologic response to climatic drivers.

A current focus of hydrology is identifying and quantifying groundwater stores (Holbrook et al., 2014; McDonnell, 2017; Rempe and Dietrich, 2018; Dralle et al., 2018; Bhanja et al., 2018), and geophysics is an important tool for examining CZ architecture and its influence on water storage and movement. For example, McGuffy (2017) used seismic refraction surveys to estimate porosity and found that initial porosity plays a significant role in bedrock weathering in granitic and rhyolitic tuff CZs. Flinchum et al. (2018) took those porosity calculations a step further by using geophysics to estimate the water holding capacity of another granitic CZ; however, both studies noted the strong influence of, and uncertainty associated with, the degree of saturation of the media. Rempe and Dietrich (2018) used downhole surveys with a neutron probe to estimate rock moisture in the CZ, and Dralle et al. (2018) used geophysics-based storage estimates from Rempe (2016) and Rempe and Dietrich (2018) to suggest differences in direct and indirect storage within the CZ from a coupled mass balance and storage–discharge function. The complexity of these estimates and their interactions highlights the need to couple geophysical approaches with subsurface interrogation, such as drilling and field characterization of hydraulic properties, to resolve this complexity, particularly in fractured heterogeneous environments.

In a headwater catchment and nested zero-order basin (ZOB) within the complex volcanic Jemez River Basin Critical Zone Observatory (JRB-CZO), a considerable amount of research has been done to characterize the hydrology of the system. For instance, previous studies have explored energy limitations and topographic controls on hydrologic partitioning and water transit times (Zapata-Rios et al., 2015a, b). Other studies used carbon pool and rare earth elements and yttrium (REY) as biogeochemical tracers of streamflow generation (Perdrial et al., 2014; Vazquez-Ortega et al., 2015, 2016) and estimated groundwater contributions using end-member mixing analyses (Liu et al., 2008a, b). The most recent JRB-CZO studies explored concentration–discharge relationships to study seasonal shifts of hydrologic flow paths (McIntosh et al., 2017) and identify the hydrochemical processes governing the transport behavior of five distinct groups of solutes (Olshansky et al., 2018). Furthermore, studies agree that there is little overland flow contribution to streamflow in headwater catchments (Liu et al., 2008b; Perdrial et al., 2014; Zapata-Rios et al., 2015a) and subsurface flow is the primary contributor to streamflow (Liu et al., 2008a, b; Perdrial et al., 2014; Vazquez-Ortega et al., 2015; Zapata-Rios et al., 2015a; McIntosh et al., 2017; Olshansky et al., 2018).

Studies spanning several water years (WYs) have shown that spring snowmelt and summer monsoons induce different surface water flow regimes. More specifically, groundwater recharge appears to be restricted to winter snowmelt (McIntosh et al., 2017) and large evaporative fluxes diminish streamflow in summer months (Zapata-Rios et al., 2015a). However, seasonal groundwater changes have not been previously observed here, and the interaction of different stores of water within the subsurface and the timing of their connection to streamflow are not understood. This gap motivated the current study, which sought to relate groundwater response, geochemistry, and age tracers across a full water year to the characterization of subsurface structure, mineralogy, and hydraulic properties. We hypothesized that there will be a more-dramatic hydrologic response of shallow groundwater to spring snowmelt and a more gradual, small change after summer monsoon events. This study also aimed to elucidate how multiple groundwater stores within the CZ contribute to streamflow during different seasonal hydrologic flow regimes.

Figure 1 (a) The Jemez River Basin Critical Zone Observatory (JRB-CZO) is located on the Valles Caldera National Preserve (VCNP) in northern New Mexico, north of Albuquerque. (b) Notice the geologic complexity of La Jara catchment and the zero-order basin (ZOB) and the location of surface water flumes and springs. (c) Note the different geology of the three well sites, the location of the soil pedons in relation to the wells, and the fault at the site 3 well. Rock types beginning with the letters “Qx” are associated with breccia deposits (Goff et al., 2011).

Most upland catchment studies to date have used springs as proxies for groundwater in the JRB-CZO; however, recent work by Frisbee et al. (2013) showed that while groundwater is a significant component of most springs, no springs are consistently composed entirely of groundwater, and they may also include soil water, unsaturated flow, and precipitation. With the recent drilling of a set of nested monitoring wells in a headwater catchment at the JRB-CZO (Fig. 1b), we can now directly access groundwater from several depths within the CZ (Fig. 2). This enabled the geochemical and isotopic analysis of groundwater and surface water from the JRB-CZO to answer the following research questions:

1. What is the seasonal hydrologic response of groundwater as a function of depth below the ground surface in two hillslopes with a contrasting lithology and CZ architecture?

2. How does the CZ architecture, such as the fracture density, the lithology, and subsurface heterogeneities, influence water routing, mean residence times, and the seasonal contribution of distinct groundwater stores to streams?

Figure 2 Profile of the nested groundwater wells within the ZOB showing the surface topography, depth of the wells, and seasonal changes in water column heights between their maximum (dotted) and minimum (solid) levels in meters above sea level. Lines capped with square ends represent the base of each well. It is important to note the different rock type at each site, the suspected presence of a perched aquifer at well 2D, the disconnection between the wells of sites 1 and 2 across the same elevation, as well as the absence of water in site 3 wells.

To address these questions, we integrated several types of analyses including hydrometric, geophysical, geochemical, isotopic, and residence time tracers to examine the hydrologic response of ground and surface water and understand the connection between distinct groundwater stores and streamflow. We compared the timing of streamflow and groundwater response to climatic drivers, quantified temporal changes in subsurface water storage, defined distinct groundwater stores, inferred recharge processes from age tracers and oxygen (δ¹⁸O) and stable hydrogen (δ²H) isotopes of water, and examined how local flow processes relate to larger-scale patterns.

2 Study site and methods

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).

$$\begin{array}{}\text{(1)}& {\displaystyle }C=C_{\mathrm{0}}{e}^{-\mathit{\lambda }t}\text{(2)}& {\displaystyle }C={\displaystyle \frac{C_{\mathrm{0}}}{\mathrm{1}+\mathit{\lambda }t}},\end{array}$$

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).

$$\begin{array}{}\text{(3)}& {\mathit{\sigma }}_t={\left[{\left({\displaystyle \frac{\mathrm{d}t}{\mathrm{d}{c}_{\mathrm{0}}}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_{c_{\mathrm{0}}}^{\mathrm{2}}+{\left({\displaystyle \frac{\mathrm{d}t}{\mathrm{d}c}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_c^{\mathrm{2}}\right]}^{\mathrm{0.5}},\end{array}$$

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.

$$\begin{array}{}\text{(4)}& {\displaystyle }t=-\mathit{\lambda }\ln \left({\displaystyle \frac{A_{\mathrm{sample}}}{A_{\mathrm{0}}}}\right)\text{(5)}& {\displaystyle }{q}_{\partial {}^{\mathrm{13}}\mathrm{C}}={\displaystyle \frac{\partial {{}^{\mathrm{13}}\mathrm{C}}_{\mathrm{measured}}-\partial {{}^{\mathrm{13}}\mathrm{C}}_{\mathrm{carbonate}}}{\partial {{}^{\mathrm{13}}\mathrm{C}}_{\mathrm{recharge}}-\partial {{}^{\mathrm{13}}\mathrm{C}}_{\mathrm{carbonate}}}}\text{(6)}& {\displaystyle }t=-\mathit{\lambda }\ln \left({\displaystyle \frac{A_{\mathrm{sample}}}{q_{\partial {}^{\mathrm{13}}\mathrm{C}}\times A_{\mathrm{0}}}}\right)\end{array}$$

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 Results

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.

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 4 Volumetric 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.

Table 1 Estimated 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.

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.

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 6 Daily 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.

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 2 Number of samples (n) and average concentrations of major ions and pH of surface waters and groundwaters. Standard deviations are shown italicized and in parentheses.

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 7 Time 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.

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.

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 9 Oxygen (δ¹⁸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).

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 3 Tritium 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.

4 Discussion

This study seeks to understand the seasonal hydrologic response of groundwater as a function of the depth below the ground surface at sites 1 and 2 and to explore how CZ architecture influences water routing and seasonal groundwater contribution to streamflow in La Jara catchment. Distinct geochemical signatures of deep groundwater from a highly fractured aquifer system and shallow groundwater that is likely a perched aquifer in a collapsed breccia deposit used in combination with geochemical and hydraulic head time series of shallow and deep groundwater and surface water from the catchment outlet enable us to decipher the effect of the contrasting CZ architecture and lithology of sites 1 and 2 on streamflow generation. The following sections discuss the dynamics of the seasonal hydrologic response, water routing, recharge and water residence times, major ion chemistry, and CZ architecture to investigate streamflow contributions over time.

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 10 Schematic 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.

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.

5 Conclusions

There are multiple separate stores of groundwater in the subsurface of the ZOB at the JRB-CZO. Major ion chemistry show that these groundwater stores are chemically distinct, while δ¹⁸O and δ²H values indicate that snowmelt is the dominant source of recharge to all groundwater. Furthermore, seasonal differences in hydrologic response to periodic precipitation inputs, including a less-pronounced response to summer monsoons, are seen in groundwaters within the ZOB. Tritium concentrations and radiocarbon analysis indicate that groundwaters are a mix of modern recharge with older waters. The path and timing of water moving through the subsurface is influenced by subsurface structure, which is exemplified by variable water contents with depth as a function of within-catchment lateral and vertical lithologic variation. The shallow aquifer at site 2 resides in a disconnected collapse breccia deposit of a different geology from the other wells, is geochemically different from La Jara stream water, and does not respond to the pressure pulse associated with increasing streamflow and rising water table in the deep site 1 aquifer system, which suggests that shallow groundwater does not contribute significantly to streamflow.

We conclude that the fractured rock site 1 architecture is most representative of the CZ volume that dominantly contributes to La Jara streamflow across the water year. The similarity in shape and timing of well 1A and surface water hydrographs results from pressure pulse propagation through the transmissive, low-storage fractured aquifer system. Furthermore, the clear similarities in major ion chemistry confirm the connection between La Jara stream and site 1 groundwater. Surprisingly, deep groundwater from site 1 wells appear to be more chemically representative of waters that contribute to La Jara stream and more representative of the structure (geology, fractured aquifer system, and greater depth to water table) and function (hydrologic response, solute fluxes, and water routing) of the CZ in the greater La Jara catchment, suggesting that deep groundwater from the fractured aquifer system, rather than shallow groundwater, sustains stream baseflow. Further, we suggest that the deep subsurface flow paths observed in the JRB-CZO are likely a signature of snow-dominated volcanic catchments transferable to other deeply fractured extrusive bedrock systems, which highlights the need to consider deeper groundwater processes in integrated hydrologic models.

The dominant contribution of deep groundwater to surface flows and the hydraulic connection between the fractured bedrock aquifer system and streamflow may suggest that deep groundwater stores in fractured bedrock aquifers are sensitive to changes in climatic drivers of streamflow, like shifts in precipitation magnitude and timing, as predicted in the southwestern United States. We assert that this study emphasizes the utility of interdisciplinary research to discern the distribution of groundwater stores, their connection to streamflow, and the underlying impact of the CZ architecture on hydrologic response to climatic drivers. Furthermore, it highlights the need to better characterize the deep subsurface of mountain systems by transferring this approach to other complex settings that challenge and advance the current understanding of subsurface hydrologic systems around the world. This study provides a template of how to bring together multiple lines of evidence to simultaneously examine both the CZ architecture and hydrology through hydrometric, geophysical, geochemical, and residence time analyses.