Subsurface Pore Water Contributions to Stream Concentration-Discharge Relations Across a Snowmelt Hydrograph

This study investigated the concentration discharge (C-Q) patterns of selected elements transported to streams during spring snowmelt through an instrumented mixed-conifer forested catchment in rhyolitic terrain in the Jemez Mountains (NM, USA). High frequency, concurrent sampling of soil solution and gases, groundwater, and surface water enabled identification and sourcing of five groups of solutes with distinct C-Q behavior. Non-hydrolysing cations and strong acid anions, had mostly positive C-Q relations and a clockwise hysteresis pattern related to flushing of a limited reservoir of solutes accumulated in soils under snowpack. Rare earth elements (REEs) and dissolved organic carbon (DOC), demonstrated large positive C-Q relations and a clockwise hysteresis pattern, consistent with their co-transport as metal-ligand complexes, and signaling bio-weathering reactions in the soil. Silicon and dissolved inorganic carbon (DIC) exhibited chemostatic C-Q trends and an anti-clockwise hysteresis pattern consistent with sourcing from deep groundwater. Hydrolyzing metals (Mn, Al, Ti and Zr) with high coefficients of variance (CV) for concentration relative to CV values for stream discharge and with no significant C-Q pattern, were found to be transported mainly as filterable colloids. Fe C-Q behavior was similar to this hydrolyzing metals group, but complexation with DOC was also important for Fe during the initial stage of the snowmelt hydrograph. Investigation of time-series of solutes and gases provided evidence for biologically-induced silicate weathering reactions that initiated in the soil subsurface and propagated down through groundwater to streams.


INTRODUCTION
Solute concentration-discharge (C-Q) relations in streams contain useful information on flow paths and biogeochemical processes that occur during water transport through the critical zone (Chorover et al., 2017). During transport, water interacts with multiple CZ compartments, including soil pores, groundwater, and hyporheic zones, accumulating and redistributing solutes and suspended matter. The chemical signature of the water changes as a function of contact times and composition in these compartments as a result of thermodynamic and kinetically limited reactions (e.g., ion exchange, mineral transformation, organic matter decomposition), and mixing of different sources. Streamwater C-Q relations integrate these multiple upgradient CZ processes (Evans and Davies, 1998;Pokrovsky et al., 2006;Godsey et al., 2009;Bouchez et al., 2017;Winnick et al., 2017).
Concentration-discharge relations are defined as chemostatic or chemodynamic, depending on whether variation in concentration is negligible or significant, respectively, relative to variation in discharge (Godsey et al., 2009;Herndon et al., 2015). Large solute reservoirs and high chemical reaction rates relative to advective velocity are hypothesized to generate chemostatic C-Q relations (Godsey et al., 2009;Thompson et al., 2011;Vaughan et al., 2017), whereas supply limited conditions (i.e., small mass stores and low reaction rates) result in solute dilution (Maher, 2011). Solute accretion (i.e., increased concentration with discharge) has been attributed to flushing of a concentrated source (e.g., soil pore water, hyporheic zones) during high flow events (Hoagland et al., 2017;Vaughan et al., 2017). Positive C-Q relations may also result from colloid mobilization during high through-flux events (Trostle et al., 2016).
Hysteresis in C-Q relations between rising and falling hydrograph limbs provides additional information on processes controlling solute export to the stream. Higher solute concentrations on the rising versus falling limb (i.e., clockwise hysteresis), suggests limited solute supply from a proximal source (Evans and Davies, 1998;Chanat et al., 2002;Bowes et al., 2009Bowes et al., , 2015Vaughan et al., 2017), whereas the opposite trend suggests longer transport time and a distal solute source, or a change in rate of groundwater relative to near surface water discharge to the stream.
Alternatively, non-stationary transit time distributions of waters through the subsurface may also result in hysteresis in the C-Q relation (Benettin et al., 2015;Winnick et al., 2017).
Most studies on C-Q relations have focused solely on measurements of stream water itself and general catchment characteristics (i.e., size, relief, and geology). To improve mechanistic interpretations, it should be beneficial to concurrently assess time-dependent hydrochemical behavior of upgradient sources (i.e., soil pore water and groundwater). Few prior studies have simultaneously collected time-series data on the chemistry and hydrologic response of soil water, groundwater and streams (Kim et al., 2014(Kim et al., , 2017Herndon et al., 2015;Hoagland et al., 2017). Furthermore, while stream hydrochemistry can vary significantly over the course of individual hydrologic events, sampling across a wide range of discharge is frequently limited to grab samples taken at random discharge levels often separated in time by days to weeks (Evans and Davies, 1998;Godsey et al., 2009;McIntosh et al., 2017). Several recent studies have indicated the need for higher (sub-daily) sampling frequency to better resolve system response to hydrologic events (Bowes et al., 2009;Kirchner and Neal, 2013;McIntosh et al., 2017). Prior work at the Jemez Critical Zone Observatory (NM, United States) has indicated that most of the groundwater recharge and the largest variation in stream discharge occurs during the annual snowmelt period, which is characterized by a single large stream hydrograph that extends over several weeks (Zapata-Rios et al., 2016). In the present study, we collected time series data across a snowmelt hydrograph at daily minimum and maximum streamflow, to test the hypotheses that (i) event-based changes in soil and groundwater chemistry would be reflected in stream C-Q relations and (ii) this would enable us to better resolve mechanistic controls over stream hydrochemistry. Therefore, the aim of this study was to link stream C-Q relations to hydrochemical dynamics in up gradient sources, by simultaneously probing soil, groundwater and stream hydrologic and hydrochemical response.

Site Description
The study site is part of the Jemez River Basin-Critical Zone Observatory (JRB-CZO) (Chorover et al., 2011), and the La Jara catchment studied here has been described in several prior publications (Perdrial et al., 2014(Perdrial et al., , 2018Vázquez-Ortega et al., 2015Zapata-Rios et al., 2016;McIntosh et al., 2017). La Jara is a 3.3 km 2 catchment extending in elevation from 3432 to 2728 m, draining the south facing slope of Redondo dome under montane mixed conifer forest within the Valles Caldera National Preserve (VCNP) in the Jemez Mountains of northern New Mexico (Figure 1). The upper headwaters include a highly instrumented 0.15 km 2 zero order basin (ZOB). The La Jara catchment (LJ) and its nested ZOB are the focus of the current study. The site was impacted by the Thompson Ridge wildfire in May-June 2013 with medium to high severity burned areas distributed in a patch-wise fashion (Reale et al., 2015).
The instrumented ZOB and larger LJ catchment contain a mixture of the upper Tshirege member of Bandelier Tuff and caldera collapse breccias. While the LJ catchment is dominated by Bandelier Tuff, in the ZOB, these two geological substrates are bisected by a series of normal faults that lie parallel to the catchment axis. The Bandelier Tuff is rhyolitic welded to non-welded vitrophyre primarily made up of plagioclase, sanidine, quartz, pyroxene and biotite in a glass matrix intra-stratified by ash-fall deposits. Asymmetrical jointing in tuffs was observed from the surface to a depth of at least 50 m in drill cores excavated from the ZOB in June, 2016. Evidence of discontinuous hydrothermal alteration exists at FIGURE 1 | Digital elevation model (DEM) of the La Jara catchment (LJ, outline in white) based on 2010 LiDAR flight data. LJ drains the southeastern slope of Redondo Peak within the Valles Caldera National Preserve. The mixed conifer zero order basin (ZOB, also outlined in white and shown in expanded form in upper inset), is an instrumented sub-catchment nested within the LJ catchment. Sampling locations of meteorological stations, flumes, instrumented pedons and groundwater monitoring wells are indicated.
Frontiers in Earth Science | www.frontiersin.org the surface on some portions of the Bandelier Tuff, while upwelling of hydrothermal fluids has altered the volcaniclastics. Caldera collapse breccias consist primarily of reworked rhyolitic volcaniclastics and pre-caldera sedimentary rocks (Hulen and Nielson, 1986;Chipera et al., 2008).

Field Data and Sample Collection
Water dynamics in the LJ catchment during spring 2017 were measured using two meteorological stations, six instrumented pedons and five monitoring wells located in the ZOB, and two flumes located at the outlet of the ZOB and at the base of LJ catchment (Figure 1).
Meteorological data (i.e., air temperature, precipitation and snow depth) were obtained from the Valles Caldera National Preserve Headquarters Redondo Climate Stations 1 (Figure 1). The Redondo station is located in the upper elevation of the LJ catchment close to the ZOB (hereafter upper LJ) whereas the Headquarters station is located closer to the LJ flume at the lower elevation of the catchment (hereafter lower LJ).
To probe the soil compartment of the CZ, six instrumented pedons were installed in the ZOB. Two pedons were located on the south-west facing slope, two were situated on the south-east facing slope, and two were placed in the convergent hollow (Figure 1). Soil volumetric water content (VWC) was measured continuously at three depths (10, 30, and 60 cm) in each of the six pedons, using 5TE sensors (Decagon, Pullman, WA, United States). Soil P CO2 , and temperature were measured continuously at depths of 2, 10, 30, and 60 cm, with GMP222 (Vaisala Inc., Houston, TX, United States). Biweekly soil water samples were collected from each pedon at four soil horizon depths (10-130 cm), using Prenart SuperQuartz soil-water samplers (2 µm pore size, Prenart Equipment ApS, Frederiksberg, Denmark). A vacuum pump applied constant suction of 40 kPa to collect the samples into acid-washed high density polyethylene (HDPE) bottles.
Groundwater was monitored with five wells located on the south-west (site 1) and south-east (site 2) facing slopes of the ZOB (Figure 1). Well depths were 43, 39, 47, 32, and 7 m below ground surface (mbgs, corresponding to elevations of 2989, 2992, 2983, 2990, and 3014 m above mean sea level) for Well 1A, Well 1B, Well 2A, Well 2C, and Well 2D, respectively. Well 2D probes a shallow perched aquifer on the south-east facing slope. The groundwater table was monitored using vibrating wire piezometers (Model 4500, Geokon Inc., Lebanon, NH, United States). Water samples were collected bimonthly from March to November 2017 from wells 1A, 1B, 2A, and 2C into acid-washed HDPE bottles using a vacuum pump (as described for the soil water samples above). Daily samples were collected from well 2D into acid-washed 1 L polypropylene bottles using an automatic water sampler (Teledyne ISCO, NE, United States). The outlets of the ZOB and LJ catchment were installed with Parshall flumes fitted with pressure transducers (InSitu Level TROLL 500 15 psig Model: 89010) that measured stage heights continuously at 15 min intervals over the study period. Discharge was estimated using standard stage-discharge relationships for 12 or 18 inch Parshall flumes. Water samples from these 1 https://wrcc.dri.edu/vallescaldera flumes, were collected in 1 L polypropylene bottles, during the highest and lowest daily discharge, using an automatic water sampler.

Aqueous Sample Processing Prior to Analysis
All water samples were packed on ice and sent overnight to UA to be filtered within 48 h of collection using 0.45 µm nylon filters. Sample aliquots for dissolved inorganic and organic carbon [dissolved inorganic carbon (DIC) and dissolved organic carbon (DOC), respectively] were filtered through combusted 0.7 µm glass-fiber filters. A subset of samples was used to evaluate colloidal contribution to element concentrations. These samples were cascade filtered through a 1.2 µm polycarbonate filter, and then through a 0.025 µm cellulose acetate/nitrate filter. Element concentrations in the "colloidal fraction" were calculated as the difference in concentrations measured for the 1.2 and 0.025 µm filtrates.

Sample Analyses
Sample chemical analysis was performed at the University of Arizona Laboratory for Emerging Contaminants (ALEC). Sample aliquots for major and minor solutes (Na + , Mg 2+ , Al, Si, K + , Ca 2+ , Mn, Fe, Sr 2+ , and Zr) as well as lanthanide series elements (rare earth elements, REEs, hereafter designated as light, medium and heavy, LREE: La to Nd, MREE: Sm to Tb, HREE: Dy to Lu) were preserved in acid-washed HDPE bottles by acidification to ca. pH 2 using trace metal grade nitric acid. All water samples were stored at 4 • C prior to analysis by inductively coupled plasma mass spectrometry (ICP-MS, Perkin Elmer ELAN DRC II ICP-MS, Shelton, CT prior to 2013; and thereafter Agilent 7700x) equipped with either a dynamic reaction (Perkin Elmer) or collision (Agilent) cell to eliminate molecular isobaric interferences. Rhodium (Rh) was included as an internal standard for Fe, Al, Mn, and P, and indium (In) was used as an internal standard for REE based on similarities in first ionization potential. Calibration curves included at least five points, and correlation coefficients were >0.995. Creation of a satisfactory calibration curve was followed by the analysis of an initial calibration blank (ICB) and calibration verification (ICV) solution, with a concentration in the low to midrange of the calibration curve. A continuing calibration blank (CCB) and continuing calibration verification (CCV) solution were measured after every 20 samples and at the end of each batch. Each batch also included measurement of at least one quality control (QC) solution from a second source, such as the NIST 1643e. Acceptable QC responses were between 90 and 110% of the certified value. DIC and DOC were measured using a Shimadzu TOC-VCSH analyzer (Shimadzu Scientific Instruments, Columbia, MD, United States). Major anions (Cl − , NO 3 − , SO 4 2− , and PO 4 ) were measured by ion chromatography (Dionex ICS-1000 equipped with an AS22 Analytical Column, Sunnyvale, CA, United States). Organic carbon molecular characterization was conducted by collecting UV-Vis absorbance spectra (190-655 nm) (Shimadzu Scientific Instruments UV-2501PC, Columbia, MD, United States), and fluorescence excitation-emission matrices (EEM) (FluoroMax-4 equipped with a 150 W Xe-arc lamp source, Horiba Jobin Yvon, Irvine, CA, United States). The EEMs were acquired by scanning excitation (Ex) from 200 to 450 nm and emission (Em) from 250 to 650 nm in 5 nm increments. Spectra were collected with Ex and Em slits at 5 and 2 nm band widths, respectively, with an integration time of 100 ms. Ultrapure water blank EEMs were subtracted and fluorescence intensities were normalized to the area under the water Raman peak, collected at an excitation of 350 nm. Inner-filter correction was performed based on the corresponding UV-Vis scans (Murphy et al., 2013). Details on the described methods, precision, and reproducibility for analysis performed in ALEC are reported elsewhere (Perdrial et al., 2014;Vázquez-Ortega et al., 2015;Trostle et al., 2016).

Data Processing and Statistical Analysis
All data processing and statistical analyses were performed using R (R Core Team, 2017) and the following packages: CAR, lsmens, nlme, and ggplot2 (Wickham, 2009;Fox and Weisberg, 2011;Lenth, 2015;Pinheiro et al., 2017). C-Q relations for each limb of the hydrograph were fitted to a power law (Eq. 1) using non-linear least-squares estimates (Gauss-Newton).
Where C is the solute concentration, Q is the discharge, b is the fitted power exponent and a is a fitting parameter. The ratio of coefficients of variance for concentration and discharge (CV C /CV Q ) was calculated as where µ and σ are the mean and standard deviation of the solute concentration (C) and discharge (Q), respectively. Solute export behavior was defined on the basis of the power law exponent, b (Eq. 1), and the ratio CV C /CV Q (Godsey et al., 2009;Thompson et al., 2011;Musolff et al., 2015). Chemostatic behavior was defined as b = 0 and CV C /CV Q ≤ 0.2, whereas chemodynamic behavior is defined as |b| ≥ 0 and/or CV C /CV Q > 0.2 (Hoagland et al., 2017). The hysteresis index (HI), which quantifies the extent of offset between C-Q trends for each limb in the hydrograph was calculated following Vaughan et al. (2017). Normalized discharge (Q i,norm ) and concentration (C i,norm ) values were calculated as where Q i and C i are the discharge and concentration value at time i, and Q min , C min , Q max , and C max are the minimum and maximum discharge and concentration values. The C norm values for each limb were interpolated at increments of 0.01 of Q norm by fitting local polynomial regression (LOESS function, stat package R Core Team, 2017). The hysteresis indices were then calculated by the difference between normalized solute concentration pairs at equivalent discharge. C-Q analysis were performed only for LJ stream and not for the ZOB stream due to low data resolution on volumetric discharge for the latter.
To further study the relations between solute C-Q behaviors, a hierarchical cluster analysis approach was applied. First, the fitted values of the power law exponents, CV C /CV Q values, and mean HI indices for each solute and limb were combined into a single matrix. In cases of a poor fit to Eq. 1 (i.e., p > 0.05), b was assigned a value of zero. Columns were standardized by subtracting the mean and dividing by the standard deviation of each column. The Euclidian distance matrix was calculated using the "dist" function in R. Hierarchical clustering was than calculated from the distance matrix using the "hclust" function with the "ward.D2" method. Principle component analysis (PCA) was performed on the standardized matrix, using the "prcomp" function to evaluate the different behaviors between clusters.
Differences between solute concentration means in water sources (i.e., soil, groundwater and streams) were assessed using linear mixed effect model for repeated measures (lme function, nlme package Pinheiro et al., 2017) one-way ANOVA and Tukey's honest significant difference (HSD) test (ANOVA and lsmean function car, and lsmeans packages, respectively, Fox and Weisberg, 2011;Pinheiro et al., 2017). Humification indices (HIX) were calculated as the ratio between the sum of fluorescence intensity in the 435-480 nm region to the region of 300-345 nm for a specific excitation wavelength of 255 nm (Ohno, 2002). Visual Minteq 3.1 was used for geochemical speciation modeling and saturation indices calculation. Carbonate speciation in streams and groundwater was calculated from pH and DIC data, whereas for soil pore water, carbonate speciation was calculated using P CO2 and alkalinity (estimated by the diffrence between non-hydrolyzing cations and strong acid anions, Stumm and Morgan, 1981).

Water and Snow Dynamics
Snowpack accumulation began in late November 2016, as (daily mean) air temperature dropped to −4 • C, below freezing, and it reached a maximum depth of ca. 1000 mm in upper LJ catchment in mid-January 2017 (data not shown). The onset of snowmelt was measured in the upper LJ climate station on March 1, 2017 as average air temperature increased above 0 • C (Figures 2A,B). During the first 3 weeks (March-1 to March-22), snow depth in upper LJ decreased from 930 to 440 mm. A subsequent drop in air temperature decreased the rate of snowmelt and this was accompanied by further snow accumulation (ca. 200 mm) as a result of new precipitation events from March 22 to April 5. Air temperature increased again and snowmelt resumed until completion on April 21. Therefore, this snowmelt hydrograph can be considered as a large melting event superimposed by a sub-event with small snow accumulation and melting. It is important to note that air temperature measured in the lower LJ climate station, was higher by 2.5 ± 0.42 degree, and most of snowpack at this location had melted by mid-March.
Following snowmelt trends, mean VWC in the ZOB, averaged across the six instrumented pedons ( Figure 2C), increased, signaling infiltration into the near surface soil layers. The duration of elevated VWC increased with soil depth, with surficial (10 cm depth) soils showing a relatively short-lived period of elevated water content between mid-March and early April, whereas the deeper horizons (60 cm depth) showed sustained saturation into early May, and the intermediate depth (30 cm) showed an intermediate behavior ( Figure 2C). Previous work indicated minimal overland flow over still frozen ground (Liu et al., 2008;McIntosh et al., 2017).
The perched groundwater table elevation, represented here by Well 2D, showed a steady response to snowmelt ( Figure 2D). Overall the perched water table increased by 1320 mm over the snowmelt period. Variation in the rate of groundwater elevation increase corresponded to the timing of snowmelt. During the initial melt period (March 1 -March 22) this rate was 1.6 × 10 −3 ± 4 × 10 −6 m h −1 . During this period snow depth was negatively correlated with groundwater table elevation (slope = −1.25 m snow m water table −1 , R 2 = 0.97). Then as the rate of snowpack ablation decreased with decreasing temperature (and new snow accumulated between March 22 and April 5), the groundwater table continued to rise, but at a greatly reduced rate of 6.5 × 10 −4 ± 2 × 10 −6 m h −1 . For this period there was a positive correlation between snow depth and groundwater table elevation (slope = 0.61 m snow m water table −1 , R 2 = 0.59). The rate of the groundwater table rise increased again to 8 × 10 −4 ± 1 × 10 −6 mh −1 during the second period of snowmelt (April 5-24). Snow depth was negatively correlated with the groundwater table elevation (slope = −0.45 m snow m water table −1 , R 2 = 0.98). Maximum water table elevation was measured on April 24, 6 days after the second peak in the LJ hydrograph.
In contrast to the Well 2D, the deep groundwater table elevation measured in the south-west facing slope (Well 1A), exhibited two pronounced peaks. During the first peak (March 1-22) the groundwater table elevated by 2 m, corresponding to the initial period of snowmelt. Following the first peak, the water table decreased by 1.3 m. The minimum between the two water table peaks (April 8) was observed 3 days after the second snow accumulation peak. The second max in groundwater table occurred on April 10. This peak precedes the maximum water table elevation measured in Well 2D by 2 weeks, and coincides with the decrease in mean ZOB VWC at 10 and 30 m depth.
The hydrograph of the ZOB stream exhibited two pronounced peaks during spring snowmelt, March 20 to April 16, 2017 ( Figure 2E). A similar pattern was observed in the hydrograph of the LJ stream ( Figure 2F), however, maximum discharge consistently occurred 2 days later in LJ (March 22 and April 18), despite the fact that the ZOB is a higher elevation stream and subject to colder air temperatures (Figure 2A). During the first snowmelt period, ZOB and LJ stream discharge increased from ca. 0 to 30 m 3 h −1 and 450 m 3 h −1 , these maximum values corresponding to specific discharges of 2 × 10 −4 and 1.36 × 10 −4 m h −1 , respectively (Figures 2E,F). The accumulation of the snowpack (March 22 to April 5) resulted in decreased discharge to 0 and 170 m 3 h −1 in ZOB and LJ streams, respectively. As snowmelt resumed, a second hydrograph rise was observed (up to 20 and 170 m 3 h −1 , for ZOB and LJ streams, respectively). Both ZOB and LJ hydrographs contain daily peaks and troughs associated with diurnal fluctuations in temperature, and hence melt rates.

Bulk and Time Series of Water Chemistry
Mean solute concentrations in soil solution, groundwater and streams (LJ and ZOB) are presented in Table 1. The pH in Well 2D groundwater was higher by 0.5 units than those of soil and stream waters. Soil solutions were significantly enriched in Zr, Al, DOC, and REEs. Mean concentrations of DIC, Si, Mn, Na + , SO 4 2− , Mg, Ca 2+ , and Sr 2+ were significantly higher in Well 2D than in the LJ stream. For example, DIC, Mn, Ca 2+ , SO 4 2− and Sr 2+ mean concentrations were 18, 16, 15, 7 and 7 times higher in Well 2D groundwater than in the LJ stream, respectively. Highest NO 3 − concentrations, which were measured in the LJ stream, were 2-3 times higher than concentrations in the ZOB stream, ground water, and soil water.
Deeper ground water at site 2 (Well 2C) showed similar chemistry to Well 2D, except that Mg 2+ was higher in Well 2C. Wells in site 1 (Well 1A, and Well 1B) had up to 10-fold lower Ca 2+ , SO 4 2− , and DIC concentrations, whereas Mn was higher in this site compared to site 2 (Supplementary Figure S2). Groundwater was slightly oversaturated with respect to calcite (mean saturation index, SI = 0.22 ± 0.009) while streams and soil pore water showed consistent under-saturation (mean SI = −2.77 ± 0.08, −2.42 ± 0.08, and −3.26 ± 0.10 for LJ, ZOB streams, and soil pore water, respectively, Supplementary Figure S1). Since some degassing of CO 2(aq) from water samples may have occurred before sample analysis, calculated SIs for calcite may overestimate the actual values in the streams and Well 2D waters. To address this issue calcite SI values were calculated for Well 2D water in equilibrium with three P CO2 concentrations, 0.4, 4, and 40 mL L −1 . The calculated SIs were, 1.35 ± 0.01, 0.44 ± 0.01, and −0.55 ± 0.01. Therefore, even for high P CO2 values, Well 2D water remained close to equilibrium with calcite. It is important to note that some calcite precipitation may have occurred between sample collection and analysis, and this may result in underestimation of calcite SI in Well 2D.
Time series of Si, DIC, Ca 2+ , and Na + in soil water, groundwater and ZOB and LJ streams are shown in Figure 3. Si concentrations in both ZOB and LJ surface waters exhibited similar trends, and concentration ranges (ca. 200 to 1000 µM). Two pronounced peaks in surface water Si concentrations occurred during the first and second falling limbs of the snowmelt hydrographs. The first set of peaks in both surface water systems occurred March 24 to April 4 when Si concentration increased from ca. 370 µM to ca. 500 µM. The second set of peaks occurred May 8-18 and was substantially larger than the first; Si concentrations increased to 1030 µM. It is noteworthy that both surface water Si peaks follow the observation of elevated Si concentrations in groundwater (Well 2D). This is especially evident for the second peak in streamwater Si concentration that follows a rapid increase in Well 2D concentration that started April 26. The increased concentration of Si in groundwater occurred in two steps the first from ca. 500 to 800 µM (April 26) and then to ca. 1300 µM (May 5). This increase in Si concentration was also observed in the deeper groundwater wells located on both slopes of the ZOB (wells 1A and 2C, Supplementary Figure S2). Groundwater Si concentration decreased to 720 µM by June 28 (Supplementary Figure S2), however, due to a gap in data between May 17 and June 28, it is unclear precisely when the Si concentrations declined in Well 2D. Interestingly, the peaks in Si concentration in groundwater coincided with increases in concentrations of DIC, Na + and (to a lesser extent) Ca 2+ in the groundwater (Figure 3). Furthermore, these increases in groundwater solute concentrations preceded by a few days the increases in surface water Si concentrations.
The ZOB stream showed a decrease in DIC during the first hydrograph rise followed by an increase during the second descending limb of the hydrograph (Figure 3). In contrast to the ZOB stream, DIC concentrations in LJ stream exhibited smaller variability. Notably, the humification index (HIX) measured for dissolved organic matter (DOM) collected from the ZOB stream, was correlated with DIC (R 2 = 0.53, Supplementary Figure S3).
In addition, DIC concentrations in the ZOB stream were correlated in time with P CO2 measured in the ZOB instrumented pedons (Supplementary Figure S4). Linear correlation with soil P CO2 explained up to 76% of the variance in ZOB stream DIC concentrations (Supplementary Figures S5, S6). It is important to note that stream water and soil gas P CO2 values were similar for LJ and ZOB streams, and for soil gas (1.28 ± 0.33, 2.18 ± 0.33, and 2.48 ± 36, mL L −1 , respectively).

C-Q Relations
To increase the information obtained from the collected data, each limb (Rise 1, Descent 1, Rise 2, and Descent 2 are denoted as R1, D1, R2, D2, respectively) of the snowmelt hydrograph was analyzed separately (Figure 4). This approach allowed us to elucidate sequential trends in snowmelt hydrochemical response that would otherwise be missed in bulk analysis of C-Q data. C-Q trends differed among solutes and hydrograph sections (Figure 4 and Table 2). For both R1 and R2, Si, DIC, and Na exhibited chemostatic or slight dilution trends (b ≤ 0, Eq. 1) with relatively low CV C /CV Q (<0.2). Besides Na + , non-hydrolyzing cations and strong acid anions (K + , Ca 2+ , Mg 2+ , Sr 2+ , SO 4 2− , Cl − , and NO 3 − ) exhibited slight positive C-Q relations (0.05 < b < 0.2) with low CV C /CV Q (<0.2) during R1. D1 showed lower positive or chemostatic trends for the anions (SO 4 2− , NO 3 , and Cl − ). Other cations showed similar positive trends in R1 and D1 (Mg 2+ , Ca 2+ , Sr 2+ , K + ). During R2, a shift from positive or chemostatic to dilution (b < −0.1) was observed for all non-hydrolyzing cations and strong acid anions except K + , which retained a positive C-Q trend. During D2, all ions except Na + showed a positive power law slope (b > 0.2). The largest b values were calculated for DOC and REEs in all of the four hydrograph limbs ( Table 2). REEs exhibited higher positive C-Q trends during R2 and D2 (b > 0.7), while DOC showed highest b values on the two descending limbs (0.49 and 0.65 for D1 and D2, respectively). The metals Mn, Al, Zr, and Ti exhibited scattered chemodynamic behavior (i.e., low b values and poor fit to power law function and high CV C /CV Q ). The power law slopes for Fe were similar to those for these other metals for the last three hydrograph limbs, but in the first rising limb, Fe concentrations were much higher than the subsequent hydrograph limbs and showed a tighter positive correlation with discharge (Figure 4 and Table 2). All of these polyvalent metals (Mn, Al, Zr, Ti, and Fe) also exhibited a large colloidal fraction in both ZOB and LJ streams (Figure 5), whereas other elements did not show a significant colloidal contribution.

C-Q Hysteresis
The La Jara stream snowmelt hydrograph exhibited different hysteresis patterns for the tested solutes. Since the snowmelt hydrograph in this work can best be described as one large snowmelt event superimposed with a smaller sub-event of snow accumulation and melting, hysteresis indices were calculated for the R1-D2 limbs as well as each of the sequential sub limbs (i.e., R1-D1, D1-R2, and R2-D2). Hysteresis indices (HIs) calculated as described in Section "Data Processing and Statistical Analysis, " are presented in Table 3 and Supplementary Figure  S7. Pronounced anticlockwise hysteresis was observed between the R1 and D2 for Si, DIC, Zr, and Al. In contrast, major ions, DOC, and REEs exhibited overall clockwise hysteresis between R1 and D2. Strong clockwise hysteresis was observed for Fe for the first hydrograph peak (R1-D1) with no hysteresis for the second peak (R2-D2). The elements Na + , K + , Fe, DOC, and REEs, exhibited larger hysteresis on the first hydrograph peak. Mg 2+ , Ca 2+ , Sr 2+ , DIC, and Cl − had similar hysteresis values for both hydrograph peaks, whereas NO 3 − , and SO 4 2− had larger hysteresis on the second hydrograph peak. Si showed anticlockwise hysteresis for the first hydrograph peak.

Cluster and Principle Component Analyses
Cluster analysis and PCA reduction methods were applied to separate solutes into groups with similar C-Q behaviors. As described in Section "Data Processing and Statistical Analysis, " values of the power law exponent (b in equation 1), coefficients of variance ratio (CV C /CV Q ), and mean HIs, were combined into a single matrix. Results of the cluster analysis and PCA performed on this matrix are presented in Figure 6. The solutes clustered into five groups: (1) Si and DIC, (2) non-hydrolyzing cations and strong acid anions (Na + , K + , Ca 2+ , Mg 2+ , Sr 2+ , SO 4 2− , Cl − , and NO 3 − ), (3) DOC and REEs (4) Al, Mn, Zr and Ti, and, (5) Fe. The first three principle components (PC) explained 78% of FIGURE 4 | La Jara flume C-Q relations for selected elements. Solids lines are best fits to power law (Eq. 1) for each of the four hydrographs limbs: first rising (blue), first falling (green), second rising (yellow), and second falling (orange). See Table 2 for extracted parameters and associated statistics. Group numbers in each plot are linked to Figure 6. the data variability, and the eigenvectors of these first three PCs are presented in Supplementary Table S1. Biplots of the first and second PCs (Figure 6), demonstrate that the hydrolyzing metals (groups 4 and 5) that exhibited higher CV C /CV Q , were separated based on the higher value of CV C /CV Q loadings in PC1 and PC2. DOC and REEs (group 3) were separated by the high power law exponent value (negative correlation with b values was dominant in PC2). Scores of group 1 plotted in the opposite direction of the power law exponent and HI loadings. K + and Na + were situated on the extremes of the major ions ellipsoid; K + scores plotted closer to group 3 (DOC and REEs) whereas Na + scores plotted closer to group 1 (Si and DIC). There was no strong separation observed between groups 3 and 4 in the biplot of PC1 and PC3, whereas Fe was strongly separated from the other elements based on its high hysteresis in the first hydrograph peak. Strong correlations were observed between DOC and REEs (R 2 > 0.5, p < 0.001) for all hydrograph limbs except D1.

DISCUSSION
The aim of this study was to link hydrogeochemical processes in the critical zone pore waters to their manifestation in stream C-Q relations on the time scale of the snowmelt hydrograph. Spring snowmelt 2017 in the LJ Jara catchment, characterized by two distinct hydrograph peaks resulting from snow melt and accumulation periods, provided a unique opportunity to test the hypothesis that stream hydrochemical response can be resolved in terms of soil solution and groundwater contributions. Results presented in Section "Water and Snow Dynamics, " suggest that soil solution, groundwater and surface waters all responded rapidly to snow melt with a short lag between the upper ZOB stream and the LJ catchment outlet. Variation in groundwater table rise correlates with snow melt and stream discharge trends (Figure 2), consistent with connectivity between soil water, groundwater and streamwater on the hydrograph time scale. The following sections discuss the hydrogeochemical processes observed in the La Jara catchment during spring snowmelt of 2017. The first section interprets relationships among time series of the different CZ compartments, and associated analysis performed on the bulk data. The second section discusses the analysis of C-Q and its link to short term CZ processes.

Time Series Data
Elevated Si concentrations in groundwater have been postulated as the source for contributing the bulk of Si to LJ surface waters (McIntosh et al., 2017). Data on Si concentrations in soils, groundwater and surface waters (Table 1 and Figure 3) support this hypothesis. During D1, Si values in LJ and ZOB streams reached similar values as found in groundwater (Well 2D), suggesting that at this period streamflow was associated with a progressively decreasing contribution of laterally transmitted soil water that diluted groundwater-derived Si during the high flow periods (R1 and R2). This suggests that pressure wave propagation from snow melt through the soil profile and into groundwater exhibits a lag where stream waters are heavily influenced by discharging groundwater later in the hydrograph pulse, even after soils ( Figure 2C) have largely transferred drainage waters to the saturated subsurface ( Figure 2D).
In contrast to Si, concentrations of Ca 2+ , Na + , Mg 2+ , Mn, and SO 4 2− did not elevate in the streams during D1, even though their concentrations were larger (by factors of 7-16, Table 1) in Well 2D groundwater relative to surface water. This suggests that stream water Si may have an additional or even different source than that represented solely by Well 2D. Indeed, concentrations of SO 4 2− , Ca 2+ , and DIC in Wells 1A and 1B were lower than observed in Well 2D (Figures 3  and Supplementary Figure S2), and are more consistent with concentrations measured in LJ and ZOB surface waters, which suggests that groundwater permeating the welded tuff on the eastern side of the ZOB is likely an important source to the ZOB and LJ streams. In any case, the relative magnitude of deeper groundwater contributions appears to increase during the descending limbs of the hydrograph, which helps to explain both C-Q and hysteresis behaviors, and the separation of solutes by source, as discussed further in Section "C-Q Analysis." The increase in Ca 2+ and DIC concentrations in Well 2D during the first peak of the snowmelt hydrograph (Figure 3) is consistent with calcite dissolution. This is with agreement with saturation indices calculated from geochemical modeling; these indicate groundwater solutions are near equilibrium with calcite. Hence, the Na + peak in Well 2D during that same period may derive either from the weathering of primary silicate minerals, or from Ca 2+ -induced displacement of Na + from clay surfaces. In the current study, Ca 2+ concentrations in Well 2D remained high (Table 1 and Figure 3) regardless of changes in water table depth. This observation is different from the rapid decrease in Ca 2+ concentration in groundwater due to dilution with fresh rain water during the wet periods as reported for an argillite bedrock in the Eel River Critical Zone Observatory (Kim et al., 2014(Kim et al., , 2017. This suggests that in the Jemez CZO study site, water percolating thorough the unsaturated zone to groundwater accumulates sufficient amounts of solutes (from relatively rapid mineral dissolution and cation exchange reactions) to maintain the high concentration in groundwater. A second silicate weathering pulse was observed as the water table in Well 2D began to recede. A similar Si pulse was captured at the same time in Wells 1A and 2C (Supplementary Figure S2). Moreover, in Well 1A this pulse was coincident with pulses of Na + , Mn, SO 4 2− , Cl − , DOC, and HIX. The increase in DOC and HIX, suggests that reactants (e.g., carbonic acid, DOM) giving rise to the weathering pulse propagate from the biologically active soil layers down through the vadose zone and into groundwater. The resulting fast Si pulse likely reflects dissolution of volcanic glass, which is prevalent in the ZOB subsurface. This rapid silicate dissolution during high water flow through the CZ is in agreement with previous observations for an argillite bedrock site (Kim et al., 2014(Kim et al., , 2017. Furthermore this rapid pulse of Si in groundwater and streams offers an additional mechanism for rapid Si recharge of groundwater, contrasting with the slow silicate dissolution mechanism postulated previously for this study site (McIntosh et al., 2017), and is in agreement with the assumption of the dispersive mixing model. On-going investigations into the mineralogy of cores extracted during the construction of wells 1 and 2 will help to elucidate the solid phases contributing to weathering pulses, and to explain the different trends in Na, Mn, SO 4 2− , Cl − , and DOC concentrations between these wells on the eastern (Wells 1, Bandelier Tuff) and western (Wells 2, caldera breccia) slopes of the ZOB. For example, relative to the eastern slope, the western slope contains threefold higher mass concentration of Fe (oxy)hydroxide (ca. 30 g kg −1 in the 1-15 m depth increment), predominantly as ferrihydrite, which we postulate acts as a reactive adsorbent for DOC (Gu et al., 1994;Eusterhues et al., 2005;Vazquez-Ortega et al., 2014). The western slope also contains calcite inclusions, not observed in the eastern slope boreholes, that likely contribute elevated Ca 2+ to Well 2D groundwater samples. Since calcium carbonate can contain significant amount of SO 4 2− (ca. 10 mg kg −1 , Gill et al., 2008;Balan et al., 2014), the presence of calcite can explain the elevated SO 4 2− concentrations in Well 2D water samples. In addition, we postulate that groundwater recharge on the eastern slope is dominated by fracture flow, whereas on the western slope it is dominated by matrix flow.
Analyzing multi-year time series on LJ surface water carbon dynamics, Perdrial et al. (2014) concluded that DOC with a terrestrial fluorescence signature and high humification index (HIX) originated from lateral water flow through near-surface soils, whereas DIC, as well as DOC with a stronger microbial fluorescence signature, were sourced from groundwater. Results of the present study indicate DIC dynamics in the ZOB stream are strongly influenced by soil biological activity. Specifically, DIC concentration and HIX (a reflection of partial microbial degradation of organic matter) in ZOB surface waters were highly correlated with each other (Supplementary Figure S3), as well as with P CO2 in the soil profile (Supplementary Figures S5,  S6) implying that soils serve as the dominant source of both DIC and DOC to groundwater and stream. Higher P CO2 are known to accumulate under snowpack (Brooks et al., 1996), consistent with an active heterotrophic microbial community during the cold season whose respiratory products (CO 2 and DOC with high HIX values) are discharged to surface water during snowmelt, as demonstrated by the close similarity in P CO2 between streams water and soil gas. Demonstrating important mechanism of CO 2 transport from soil to streams. It is important to note, however, that higher P CO2 value occurs during the warm season (Brooks et al., 1996;Perdrial et al., 2014). This activity enables the processing of organic matter and accumulation of the degradation products in the soil profile. As water starts to move through the soil profile it transports this processed organic material, carbonic acid, and DOM to the subsurface and ZOB stream. Increased soil saturation during the snowmelt (Figure 2C), limits microbial activity (i.e., decrease in soil P CO2 ), probably due to oxygen limitation. Thus, degradation of soil organic matter and transport of processed DOM to the stream decreases (i.e., lower HIX). As soils begin to drain, microbial and plant root activity increases (indicated by increased P CO2 and HIX values). Hence, the soil hydrologic response from snow melt affects heterotrophic microbial activity. This, in turn, drives carbonic acid and complexing ligand (DOC with high HIX) production. These reactants propagate downward into the deeper subsurface to affect the second weathering pulse observed in the wells on both sides of the ZOB. These results provide support for the analytical model proposed by Calabrese et al. (2017), which suggests that transport of carbonic acid from bioactive soil exerts a principal control over subsurface bedrock chemical weathering rates. Results of the current study suggest that organic ligand production also controls these weathering reactions and needs to be included in future models.

C-Q Analysis
Principle component analysis is commonly used to develop mixing models (Nils Christophersen, 1992;Hooper, 2003;Liu et al., 2008). Wymore et al. (2017) used PCA, to reduce dimensionality of catchment variables, enabled them to compare C-Q parameters between multiple sites. In the current study, cluster analysis and PCA enabled us to sort the large set of quantified major and trace solutes into five distinct groups on the basis of their hydrochemical behavior (Section "Cluster and Principle Component Analyses, " Figure 6). These groups represent different coupled processes that occur in the CZ during spring snowmelt, which in this system represents the largest through-flux of fresh water. The first group, Si and DIC, represents mobile products of silicate weathering reaction. Elements of this group demonstrated chemostatic or dilution C-Q trends and anti-clockwise hysteresis patterns. Their higher concentrations in groundwater (Table 1), and the Si concentration peak during the first hydrograph recession (Figure 3), suggest that Si and DIC derive from a relatively large groundwater reservoir that is displaced by new meteoric waters entering through surficial soils. Such an assessment is also consistent with the observed groundwater table rise ( Figure 2D) that coincides with high downward fluxes through saturated soils ( Figure 2C). Groundwater Si is, therefore, potentially recharged by pulsed upgradient weathering events, as observed in this study.
Non-hydrolyzing cations and strong acid anions formed the second group. With the exception of Na + , ions in this group exhibited slightly positive power law exponents in the R1, D1, and D2 hydrograph limbs (Table 2 and Figure 4). Most ions of this group exhibited dilution trends during the second hydrograph rise (R2), indicating lower generation rates of these ions during the second part of the snowmelt hydrograph, consistent with flushing and removal during the first hydrograph peak. The high frequency sampling method provide increased resolution for these trends and shifts in C-Q trend, compare to previous study that demonstrated chemostatic behavior of these elements (McIntosh et al., 2017). Potassium was unique in retaining a positive trend during the second rise of the hydrograph. All solutes in this group also showed clockwise hysteresis for both peaks in the hydrograph (Table 3 and Supplementary Figure S7). These trends are consistent with a flush of a limited pool of nutrient ions (Ca 2+ , Mg 2+ , K + , NO 3 − , SO 4 2− ) that accumulate under snowpack as a result of both wet deposition followed by downward ion migration through the snowpack during freeze-thaw cycles as well as organic matter decomposition under the snowpack (Brooks et al., 1996;Schmidt and Lipson, 2004). Trend of positive or chemostatic C-Q relations for these elements also been observed in the Susquehanna Shale Hills CZO, and explain by the large contribution of soil later flow and mobilization accumulated elements to stream (Bao et al., 2017;Li et al., 2017a). Na + , which is not a nutrient cation, diverged from the rest of the group in its display of slight dilution to chemostatic behavior. Likewise, on the PC biplots (Figure 6), Na + scored closest to the Si/DIC group. This difference is likely the result of lower biocycling of Na + (and hence less accumulation in decomposing organic matter of surface soils) relative to the other ions in this group. Moreover, groundwater obtained from wells on both sides of the ZOB was more enriched in Na + than was soil or stream water. Therefore, Na + discharge may be dominated by release from the deeper groundwater weathering zone (similar to the case for Si and DIC), as has been suggested in previous work from this site (Zapata-Rios et al., 2015). The other borderline member of this group was K + , which did not change the power exponent sign from positive to negative during the second hydrograph rise, R2 (Table 2). Interestingly, on the PC biplot K + scored closer to group 3 solutes (DOC and REEs). Since soil is the primary source of DOC to streams, our analysis suggests that K + discharge is dominated by this source as well, and indeed it is also known to be mobilized readily by leaching of biomass (Tripler et al., 2006).
Prior observations of correlations between REE and DOC concentrations has led to the suggestion that fractionation of the lanthanide series resulting from variation in metal-ligand stability is a useful probe of ligand-promoted biological weathering in the CZ (Vázquez-Ortega et al., 2015). A prior study in a granitic catchment showed that these elements are transported as colloidal material in streams, and that this process explained positive power law exponents in C-Q relations (Trostle et al., 2016). Results of the present study indicate that, in the volcanic terrain of La Jara catchment, REEs in streamwater occur as aqueous organic complexes that are tightly correlated with and potential controlled by DOC transport.
The last two groups (4 and 5) exhibited strong coefficients of variance in concentrations compared to discharge (i.e., high values of CV C /CV Q ), and also showed large colloidal fractions throughout the snowmelt hydrograph (Table 2 and Figure 5). In contrast to group 4 (Ti, Zr, Mn, and Al), Fe showed a high concentration flush and a positive C-Q trend in the first hydrograph limb (R1, Table 2 and Figure 4), and a strong clockwise hysteresis in the first hydrograph peak, which led to it being distinguished as a separate solute group. The positive correlation of Fe with DOC during R1 (R = 0.56, p < 0.001), and large positive hysteresis pattern for the same period, suggests a rapid flush of small Fe-DOM colloids that accumulated over winter in the soil or riparian zones. Overall, the C-Q relations of elements in groups 4 and 5, are consistent with a dominant colloidal (<1.2 µm) control on aqueous filtrate concentrations, as was confirmed by direct cascade filtration (Figure 5). Previous studies have suggested increased fluid shear stress during lateral soil water flow is the dominant mechanism for colloid detachment and transport from soils and hyporheic zones to the stream (Bradford et al., 2002;Wörman et al., 2002;Trostle et al., 2016). Therefore, elements that transport to streams as colloids may demonstrate pronounced positive C-Q relations (Trostle et al., 2016). The fact that components of groups 4 (Ti, Zr, Mn, and Al) and 5 (Fe) did not show consistent and clear positive power law exponents in the present study (Table 2 and Figure 4) despite being dominantly present in colloidal form ( Figure 5) is consistent with the fact that the Jemez CZ is deeper and more complex than that of the Catalina site studied by Trostle et al. (2016). In the LJ catchment, colloids mobilized from soils may be subsequently immobilized in pores and fractures of the deep vadose and saturated zones.

CONCLUSION
The La Jara catchment is representative of the high elevation, mountain headwaters of the Rio Grande, which serves as the principal water source for Albuquerque, El Paso and other cities. Monitoring of pore waters and surface waters in the LJ catchment during spring snowmelt identified five distinct groups of C-Q behaviors that are distinguished on the basis of their sourcing and transport behavior in the CZ. Probing multiple CZ compartments, including soil water, groundwater and surface water, enabled the establishment of linkages between C-Q behavior of these groups and CZ processes. Si and DIC (group 1) derive from a groundwater source, as their contribution to surface waters become more dominant when shallow subsurface flows (soil) decreased during the descending limb of the hydrograph. Lanthanide series elements and DOC (group 3), had stronger positive C-Q trends and clockwise hysteresis patterns that were more pronounced compared with other tested elements. The C-Q behavior of these solutes suggest that they are useful signals of soil biologically influenced weathering reactions (e.g., strong positive chemodynamic behavior). Non-hydrolyzing cations and strong acid anions (group 2) demonstrated mostly positive C-Q trends and a clockwise hysteresis pattern. These trends imply flushing of a limited reservoir of solutes that accumulated under the snowpack. Endmembers of this group included Na + , which showed some similarity to DIC and Si, suggesting contribution from deep subsurface weathering, and K + , which exhibited C-Q behavior transitional to REEs and DOC, consistent with its sourcing from decomposing organic materials and active bio-cycling role. Groups 4 and 5 comprise metals that exhibit higher coefficients of variance (CV) for concentration relative to discharge, with generally poor fits to the power law function. These groups also had high colloidal fractions in surface water, demonstrating that adhesion and shearing effects that control colloid transport in complex geological settings are not strongly correlated with discharge. Although Fe had much in common with these other hydrolyzing metals and also comprised a significant colloidal component, it was represented by a separate group (5) due to its apparent complexation with DOM and transport to the stream as a metal-organic complex during the first rise of the hydrograph.
Frequent sampling of surface and ground water further revealed a large pulse of Si in streams and wells of the LJ catchment that occurred during the second descending limb of the hydrograph. This pulse was linked to enhanced bioactivity (increased soil P CO2 and HIX values) and vertical downward flow of DIC and DOC that promoted fast weathering of silicate minerals, likely including volcanic glass, which is an abundant and highly weatherable mineral present in the CZ at this volcanic site. Given that soil biological activities likely drive the DIC and DOC impacts on silicate weathering observed here, consideration of such effects would be useful for the ongoing development of reactive transport models of critical zone evolution (Li et al., 2017b).