# CRITICAL ZONE (CZ) EXPORT TO STREAMS AS INDICATOR FOR CZ STRUCTURE AND FUNCTION

EDITED BY : Julia Perdrial, Ashlee Dere, Nicole West and Pamela L. Sullivan PUBLISHED IN : Frontiers in Earth Science and Frontiers in Environmental Science

#### Frontiers eBook Copyright Statement

The copyright in the text of individual articles in this eBook is the property of their respective authors or their respective institutions or funders. The copyright in graphics and images within each article may be subject to copyright of other parties. In both cases this is subject to a license granted to Frontiers. The compilation of articles constituting this eBook is the property of Frontiers.

Each article within this eBook, and the eBook itself, are published under the most recent version of the Creative Commons CC-BY licence. The version current at the date of publication of this eBook is CC-BY 4.0. If the CC-BY licence is updated, the licence granted by Frontiers is automatically updated to the new version.

When exercising any right under the CC-BY licence, Frontiers must be attributed as the original publisher of the article or eBook, as applicable.

Authors have the responsibility of ensuring that any graphics or other materials which are the property of others may be included in the CC-BY licence, but this should be checked before relying on the CC-BY licence to reproduce those materials. Any copyright notices relating to those materials must be complied with.

Copyright and source acknowledgement notices may not be removed and must be displayed in any copy, derivative work or partial copy which includes the elements in question.

All copyright, and all rights therein, are protected by national and international copyright laws. The above represents a summary only. For further information please read Frontiers' Conditions for Website Use and Copyright Statement, and the applicable CC-BY licence.

ISSN 1664-8714 ISBN 978-2-88963-623-5 DOI 10.3389/978-2-88963-623-5

#### About Frontiers

Frontiers is more than just an open-access publisher of scholarly articles: it is a pioneering approach to the world of academia, radically improving the way scholarly research is managed. The grand vision of Frontiers is a world where all people have an equal opportunity to seek, share and generate knowledge. Frontiers provides immediate and permanent online open access to all its publications, but this alone is not enough to realize our grand goals.

#### Frontiers Journal Series

The Frontiers Journal Series is a multi-tier and interdisciplinary set of open-access, online journals, promising a paradigm shift from the current review, selection and dissemination processes in academic publishing. All Frontiers journals are driven by researchers for researchers; therefore, they constitute a service to the scholarly community. At the same time, the Frontiers Journal Series operates on a revolutionary invention, the tiered publishing system, initially addressing specific communities of scholars, and gradually climbing up to broader public understanding, thus serving the interests of the lay society, too.

#### Dedication to Quality

Each Frontiers article is a landmark of the highest quality, thanks to genuinely collaborative interactions between authors and review editors, who include some of the world's best academicians. Research must be certified by peers before entering a stream of knowledge that may eventually reach the public - and shape society; therefore, Frontiers only applies the most rigorous and unbiased reviews.

Frontiers revolutionizes research publishing by freely delivering the most outstanding research, evaluated with no bias from both the academic and social point of view. By applying the most advanced information technologies, Frontiers is catapulting scholarly publishing into a new generation.

#### What are Frontiers Research Topics?

Frontiers Research Topics are very popular trademarks of the Frontiers Journals Series: they are collections of at least ten articles, all centered on a particular subject. With their unique mix of varied contributions from Original Research to Review Articles, Frontiers Research Topics unify the most influential researchers, the latest key findings and historical advances in a hot research area! Find out more on how to host your own Frontiers Research Topic or contribute to one as an author by contacting the Frontiers Editorial Office<EMAIL_ADDRESS>
# CRITICAL ZONE (CZ) EXPORT TO STREAMS AS INDICATOR FOR CZ STRUCTURE AND FUNCTION

Topic Editors:

Julia Perdrial, University of Vermont, United States Ashlee Dere, University of Nebraska Omaha, United States Nicole West, Central Michigan University, United States Pamela L. Sullivan, University of Kansas, United States

Credit: (*modified after Chorover et al. 2007, with permission*)

Citation: Perdrial, J., Dere, A., West, N., Sullivan, P. L., eds. (2020). Critical Zone (CZ) Export to Streams as Indicator for CZ Structure and Function. Lausanne: Frontiers Media SA. doi: 10.3389/978-2-88963-623-5

# Table of Contents

*04 Editorial: Critical Zone (CZ) Export to Streams as Indicator for CZ Structure and Function*

Julia Perdrial, Pamela L. Sullivan, Ashlee Dere and Nicole West

Justin B. Richardson and Elizabeth K. King

*40 Solute Fluxes Through Restored Prairie and Intensively Managed Critical Zones in Nebraska and Iowa*

Ashlee L. Dere, Andrew W. Miller, Amy M. Hemje, Sara K. Parcher, Courtney A. Capalli and E. Arthur Bettis

Scott W. Bailey, Kevin J. McGuire, Donald S. Ross, Mark B. Green and Olivia L. Fraser

*103 Hysteretic Response of Solutes and Turbidity at the Event Scale Across Forested Tropical Montane Watersheds*

Adam S. Wymore, Miguel C. Leon, James B. Shanley and William H. McDowell

*116 Soil Aggregates as a Source of Dissolved Organic Carbon to Streams: An Experimental Study on the Effect of Solution Chemistry on Water Extractable Carbon*

Malayika M. Cincotta, Julia N. Perdrial, Aaron Shavitz, Arianna Libenson, Maxwell Landsman-Gerjoi, Nicolas Perdrial, Jesse Armfield, Thomas Adler and James B. Shanley

## Editorial: Critical Zone (CZ) Export to Streams as Indicator for CZ Structure and Function

#### Julia Perdrial 1 \*, Pamela L. Sullivan2 , Ashlee Dere3 and Nicole West 4

*1 Department of Geology, University of Vermont, Burlington, VT, United States, 2 College of Earth, Ocean, and Atmospheric Science, Oregon State University, Corvallis, OR, United States, 3 Department of Geography/Geology, University of Nebraska at Omaha, Omaha, NE, United States, 4 Department of Earth and Atmospheric Sciences, Central Michigan University, Mt. Pleasant, MI, United States*

Keywords: critical zone, stream, system, interdisciplinary, catchment

#### **Editorial on the Research Topic**

#### **Critical Zone (CZ) Export to Streams as Indicator for CZ Structure and Function**

The critical zone (CZ) is commonly defined as the zone of life, spanning from the tops of the tree canopy down to the actively cycled ground water (Brantley et al., 2007; Chorover et al., 2007; Richter and Mobley, 2009). Considering the CZ as an open thermodynamic system, where energy and matter can exchange freely, precipitation is one of many important CZ inputs, while stream exports dissipate energy and matter (Chorover et al., 2011).

Traditionally, investigations of the Earth's surface are limited to separate parts of the CZ (e.g., the biosphere, lithosphere or hydrosphere) that are analyzed outside of their context. However, streams are integrators of the physical structure and connected nature of subsurface pore space in the CZ as well as the biogeochemical reactions that play out along the multitude of flow paths. In order to decipher such pathways and related CZ processes, we need to investigate CZ processes within the catchment context.

For example, we know that relief and soil thickness play a large role in determining hydrologic flowpaths, fluid residence time, and soil drainage class (Maher, 2011). These processes can result in changes in soil dynamics, including biogeochemistry, that cascade through the system and ultimately control stream water composition (Perdrial et al., 2015).

To probe stream water dynamics, high-frequency datasets are now increasingly generated and long-term datasets are used to identify changes in CZ processes (Crowther et al., 2016; Robinson et al., 2016; Hirmas et al., 2018). For example, many streams show long-term changes in fluxes of dissolved organic carbon (Raymond and Saiers, 2010; Perdrial et al., 2014; Zarnetske et al., 2018), weathering solutes (McIntosh et al., 2017), and nutrients (Björnerås et al., 2017; Jarvie et al., 2018), which have been related to changes in climate drivers, land use, and overall soil conditions. In particular, investigations of the relationships between solutes and their discharge regimes have shown strong potential to understand the impacts of extreme perturbations (Chorover et al., 2017; McIntosh et al., 2017) and the sensitivity of CZ processes to these changes (Abbott et al., 2018; Sullivan et al., 2018). The goal of this special topic issue is to highlight opportunities as well as limitations of streams as indicator of CZ structure and processes.

Several contributions highlight how relief, topography, and soil thickness affect CZ evolution and are reflected in stream composition. Cincotta et al. and Armfield et al. found that CZ

#### Edited and reviewed by:

*Alexandra V. Turchyn, University of Cambridge, United Kingdom*

\*Correspondence: *Julia Perdrial<EMAIL_ADDRESS>
#### Specialty section:

*This article was submitted to Biogeoscience, a section of the journal Frontiers in Earth Science*

Received: *03 December 2019* Accepted: *04 February 2020* Published: *25 February 2020*

#### Citation:

*Perdrial J, Sullivan PL, Dere A and West N (2020) Editorial: Critical Zone (CZ) Export to Streams as Indicator for CZ Structure and Function. Front. Earth Sci. 8:37. doi: 10.3389/feart.2020.00037*

**4**

composition varies drastically with landscape position in a northeastern United States forested headwater catchment, and that near stream environments disproportionally affect stream water response. Bailey et al. found that, along stream segments, the connectivity to shallow and deep upslope soils strongly impacted stream water geochemistry. Resulting elluvial soil processes were found to control stream chemistry in shallow soil environments while illuvial soil processes persisted in deep soils. Using a refined model and two contrasting forested sites, Billings et al. explored linkages between upland and depositional environments and the impact of soil organic carbon (SOC) erosion and burial on carbon budgets within watersheds on decadal to centennial time scales. They noted that the reactivity and distribution of SOC across watersheds is necessary to predict carbon storage across the landscape.

Several Critical Zone Observatory (CZO) studies reveal the importance of biotic processes and storm event magnitude in controlling elemental behavior. Olshansky et al. found that concentration-discharge patterns during snowmelt in the Jemez River Basin CZO in New Mexico indicate biologically mediated silicate weathering that propagates from the shallow CZ down through groundwater and strongly impacts stream water composition. A logical extension of this work is that processes impacting root/microbial respiration and biologically generated acids will alter both silicate weathering and stream water composition. Using sensor data from the Luquillo CZO in

#### REFERENCES

Puerto Rico, Wymore et al. found transport vs. supply limitation using concentration-discharge patterns of specific conductance and turbidity. The authors note that these relationships are especially influenced by the magnitude of events (including the category 5 Hurricane Maria) and CZ lithology.

Variations in bedrock lithology are also a theme in a cross site comparison of oxyanion dynamics at several CZOs (including Boulder Creek, Calhoun, Luquillo, and Southern Sierra). The study by Richardson et al. emphasizes the importance of deep bedrock and transport processes within the CZ as controls for stream water dynamics. The role of land use (agricultural and restored prairie) is investigated in the contribution by Dere et al., where solute fluxes through soils and to streams in Nebraska and Iowa are quantified. Their results indicate that changes in land use propagate to streams and drive weathering product export less than 50 years after land use change.

#### AUTHOR CONTRIBUTIONS

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

#### FUNDING

This research was supported by the National Science Foundation, grant no. EAR-1724171.

may intensify water cycle. Nature 561, 100–103. doi: 10.1038/s41586-018- 0463-x

alternative stable states of soil moisture. Sci. Rep. 6:20018. doi: 10.1038/srep 20018

**Conflict of Interest:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2020 Perdrial, Sullivan, Dere and West. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

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

#### Yaniv Olshansky1 , Alissa M. White2 , Bryan G. Moravec1 , Jennifer McIntosh2 and Jon Chorover1 \*

1 Department of Soil, Water and Environmental Science, The University of Arizona, Tucson, AZ, United States, 2 Department of Hydrology and Atmospheric Sciences, The University of Arizona, Tucson, AZ, United States

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, United States). High frequency, concurrent sampling of soil solution and gasses, groundwater, and surface water enabled identification and sourcing of five groups of solutes with distinct C-Q behavior. Non-hydrolyzing 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 biologically induced 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 gasses provided evidence for biologically induced silicate weathering reactions that initiated in the soil subsurface and propagated down through groundwater to streams.

Keywords: critical zone, mineral weathering, groundwater, cluster analysis, PCA

#### Edited by:

Pamela L. Sullivan, The University of Kansas, United States

#### Reviewed by:

Daniel Enrique Ibarra, Stanford University, United States Matthew Winnick, University of Massachusetts Amherst, United States Elizabeth Herndon, Kent State University, United States

#### \*Correspondence:

Jon Chorover<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Biogeoscience, a section of the journal Frontiers in Earth Science

Received: 02 June 2018 Accepted: 12 October 2018 Published: 09 November 2018

#### Citation:

Olshansky Y, White AM, Moravec BG, McIntosh J and Chorover J (2018) Subsurface Pore Water Contributions to Stream Concentration-Discharge Relations Across a Snowmelt Hydrograph. Front. Earth Sci. 6:181. doi: 10.3389/feart.2018.00181

### KEY POINTS

feart-06-00181 November 9, 2018 Time: 15:13 # 2

#### 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., 2009, 2015; Vaughan 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, 2017; Herndon 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.

#### MATERIALS AND METHODS

#### 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, 2018; Vázquez-Ortega et al., 2015, 2016; Zapata-Rios et al., 2016; McIntosh et al., 2017). La Jara is a 3.3 km2 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 km2 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

groundwater monitoring wells are indicated.

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 Stations1 (**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 PCO2, 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 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+, Mg2+, Al, Si, K +, Ca2+, Mn, Fe, Sr2+, 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−, NO3 −, SO4 2−, and PO4) 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

1https://wrcc.dri.edu/vallescaldera

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

$$\mathbf{C} = \mathbf{a} \times \mathbf{Q}^b \tag{1}$$

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 (CVC/CVQ) was calculated as

$$\frac{CV\_C}{CV\_Q} = \frac{\mu\_Q \sigma\_C}{\mu\_C \sigma\_Q} \tag{2}$$

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 CVC/CVQ (Godsey et al., 2009; Thompson et al., 2011; Musolff et al., 2015). Chemostatic behavior was defined as b = 0 and CVC/CVQ ≤ 0.2, whereas chemodynamic behavior is defined as |b| ≥ 0 and/or CVC/CVQ > 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 (Qi,norm) and concentration (Ci,norm) values were calculated as

$$Q\_{i,norm} = \frac{Q\_i - Q\_{\min}}{Q\_{\max} - Q\_{\min}} \tag{3}$$

$$C\_{i,norm} = \frac{C\_i - C\_{\min}}{C\_{\max} - C\_{\min}} \tag{4}$$

where Qi and Ci are the discharge and concentration value at time i, and Qmin, Cmin, Qmax, and Cmax are the minimum and maximum discharge and concentration values. The Cnorm values for each limb were interpolated at increments of 0.01 of Qnorm 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, CVC/CVQ 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 PCO2 and alkalinity (estimated by the diffrence between non-hydrolyzing cations and strong acid anions, Stumm and Morgan, 1981).

#### RESULTS

#### 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 m3 h −1 and 450 m3 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 m3 h −1 in ZOB and LJ streams, respectively. As snowmelt resumed, a second hydrograph rise was observed (up to 20 and 170 m3 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+, SO4 2−, Mg, Ca2+, and Sr2+ were significantly higher in Well 2D than in the LJ stream. For example, DIC, Mn, Ca2+, SO4 2− and Sr2+ mean concentrations were 18, 16, 15, 7 and 7 times higher in Well 2D groundwater than in the LJ stream, respectively. Highest NO3 − 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 Mg2+ was higher in Well 2C. Wells in site 1 (Well 1A, and Well 1B) had up to 10-fold lower Ca2+, SO4 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 CO2(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 PCO2 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 PCO2 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, Ca2+, 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) Ca2+ 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

Mean values in µM except for Zr, Mn, and REE, which are reported in nM. The 95% confidence intervals are shown in parentheses. †Letters display differences between means examined using Tukey's-HSD.

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 PCO2 measured in the ZOB instrumented pedons (**Supplementary Figure S4**). Linear correlation with soil PCO2 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 PCO2 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 CVC/CVQ (<0.2). Besides Na+, non-hydrolyzing cations and strong acid anions (K+, Ca2+, Mg2+, Sr2+, SO4 2−, Cl−, and NO3 −) exhibited slight positive C-Q relations (0.05 < b < 0.2) with low CVC/CVQ (<0.2) during R1. D1 showed lower positive or chemostatic trends for the anions (SO4 2−, NO3, and Cl−). Other cations showed similar positive trends in R1 and D1 (Mg2+, Ca2+, Sr2+, 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 CVC/CVQ). 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. Mg2+, Ca2+, Sr2+, DIC, and Cl− had similar hysteresis values for both hydrograph peaks, whereas NO3 −, and SO4 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 (CVC/CVQ), 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+, Ca2+, Mg2+, Sr2+, SO4 2−, Cl−, and NO3 −), (3) DOC and REEs (4) Al, Mn, Zr and Ti, and, (5) Fe. The first three principle components (PC) explained 78% of

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 CVC/CVQ, were separated based on the higher value of CVC/CVQ 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

#### TABLE 2 | Summary of power law fit (Eq. 1) and coefficient of variance (Eq. 2) for C-Q relations in the LJ stream.

†Asterisks represent level of significance with ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001. ‡Poor fit to Eq. 1 (p > 0.05) are marked in italics. Power exponents (b) are presented with 95% confidence interval in parentheses.

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

each plot are the mean and standard deviation of the colloidal fraction over the course of the experiment.

In contrast to Si, concentrations of Ca2+, Na+, Mg2+, Mn, and SO4 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 SO4 2−, Ca2+, 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 Ca2+ 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

are shown in brown arrows. Ellipsoid probability (95%) encapsulates the groups.

calcite. Hence, the Na+ peak in Well 2D during that same period may derive either from the weathering of primary silicate minerals, or from Ca2+-induced displacement of Na+ from clay surfaces. In the current study, Ca2+ 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 Ca2+ 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, 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, SO4 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, 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.

TABLE 3 | Summary of mean hysteresis indices (HI) calculated for pairs of hydrograph limbs.

95% confidence intervals are presented in parentheses.

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, SO4 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 Ca2+ to Well 2D groundwater samples. Since calcium carbonate can contain significant amount of SO4 2− (ca. 10 mg kg−1 , Gill et al., 2008; Balan et al., 2014), the presence of calcite can explain the elevated SO4 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 PCO2 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 PCO2 are known to accumulate under snowpack (Brooks et al., 1996), consistent with an active heterotrophic microbial community during the cold season whose respiratory products (CO2 and DOC with high HIX values) are discharged to surface water during snowmelt, as demonstrated by the close similarity in PCO2 between streams water and soil gas. Demonstrating important mechanism of CO2 transport from soil to streams. It is important to note, however, that higher PCO2 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 PCO2), 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 PCO2 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 (Ca2+, Mg2+, K+, NO3 −, SO4 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 CVC/CVQ), 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 PCO2 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

#### REFERENCES

of such effects would be useful for the ongoing development of reactive transport models of critical zone evolution (Li et al., 2017b).

#### AUTHOR CONTRIBUTIONS

YO performed the data analysis, developed the conceptual framework, and wrote the manuscript with support from JC, JM, AW, and BM. JC and JM conceived the original idea. JC supervised the project.

#### FUNDING

This research was funded by the National Science Foundation, Grant No. EAR-1331408, which supports the Catalina-Jemez Critical Zone Observatory, and through the Binational Agricultural Research and Development (BARD) program, which awarded a postdoctoral fellowship to YO (Grant No. FI-534-2015) and the National Science Foundation, Grant No. EAR-1331408, which supports the Catalina-Jemez Critical Zone Observatory.

#### ACKNOWLEDGMENTS

All data utilized in this paper are available through the criticalzone.org website at http://criticalzone.org/catalina-jemez/ data/. We are grateful to Mary Kay Amistadi, Rachel Nadine Burnett, Mark Losleben, and Adam Killebrew for assistance with laboratory chemical analyses and field work.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart. 2018.00181/full#supplementary-material

during reactive transport. Vadose Zone J. 13:7. doi: 10.2136/vzj2013.10. 0179

and solute generation in forested tropical montane watersheds, Water Resour. Res. 53, 6279–6295. doi: 10.1002/2016WR020016

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Olshansky, White, Moravec, McIntosh and Chorover. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Regolith Weathering and Sorption Influences Molybdenum, Vanadium, and Chromium Export via Stream Water at Four Granitoid Critical Zone Observatories

#### Justin B. Richardson1,2 \* and Elizabeth K. King3

1 Department of Geosciences, University of Massachusetts Amherst, Amherst, MA, United States, 2 Earth and Atmospheric Sciences, College of Agriculture and Life Sciences, Cornell University, Ithaca, NY, United States, 3 Department of Earth, Ocean, and Atmospheric Sciences, University of British Columbia, Vancouver, BC, Canada

#### Edited by:

Julia Perdrial, University of Vermont, United States

#### Reviewed by:

Angelica Vazquez-Ortega, Bowling Green State University, United States Dragos George Zaharescu, Georgia Institute of Technology, United States

\*Correspondence:

Justin B. Richardson<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Biogeoscience, a section of the journal Frontiers in Earth Science

Received: 20 June 2018 Accepted: 19 October 2018 Published: 09 November 2018

#### Citation:

Richardson JB and King EK (2018) Regolith Weathering and Sorption Influences Molybdenum, Vanadium, and Chromium Export via Stream Water at Four Granitoid Critical Zone Observatories. Front. Earth Sci. 6:193. doi: 10.3389/feart.2018.00193 Understanding the fate of oxyanions in the Critical Zone is important because of their biological significance and the potential for their use as geochemical tracers in terrestrial environments and subsurface systems. This study assessed the partitioning and transport of a suite of oxyanion metals (Mo, V, and Cr) in regolith profiles and stream waters from four granitoid Critical Zone Observatories (CZOs) (Boulder Creek, Calhoun, Luquillo, and Southern Sierra). For regolith profiles, we compared Mo, V, and Cr in total digestions and two extractions targeting oxyanions adsorbed to organic matter and amorphous oxides (H2O2 + 0.1 M acetic acid) and secondary Fe oxides (citrate– bicarbonate–dithionite). Total Mo, V, and Cr ranged from 0.4 to 2.5 mg kg−1 , 16 to 208 mg kg−1 , and 0.2 to 55 mg kg−1 , respectively. The greatest concentrations of the oxyanions did not occur in surface soil samples, nor deepest regolith samples (7–10 m in depth), but instead in subsurface peaks that corresponded with secondary Fe oxides and total organic carbon. The average organic and amorphous oxide bound phase was 0.1–3.5% while the secondary Fe oxide fraction was 4–27% of the respective total concentrations for oxyanions, suggesting that secondary Fe oxides were an important phase across the regolith profile. Stream water Mo, V, and Cr concentrations ranged from 0.02 to 0.25 µg L−1 , 0.2 to 1.8 µg L−1 , and 0.08 to 0.44 µg L−1 , respectively. Our results demonstrate that the deep regolith (2–7 m in depth) play an active role in both sourcing and retention of oxyanions. In addition, we observed that increased weathering intensity at warmer, wetter climates does not always lead to increased depletion in regolith or stream water export, which implies the importance of transport processes within regolith. Further quantification of oxyanion export from regolith can aid in developing their use as geochemical tools for global weathering.

Keywords: deep regolith, oxyanions, stream export, Critical Zone, weathering profile, soil biogeochemistry

### INTRODUCTION

feart-06-00193 November 8, 2018 Time: 16:14 # 2

The Critical Zone (CZ) is the thin region on Earth's surface where the atmosphere, hydrosphere, lithosphere, and biosphere interact through chemical, biological, and physical processes (e.g., Brantley et al., 2007; Chorover et al., 2007; Giardino and Houser, 2015). Exploring the natural processes that govern the CZ is important as they ultimately sustain terrestrial life and exert controls on the biotic and abiotic cycling of energy and nutrients (Amundson et al., 2007). However, the CZ exhibits substantial variations in climate, lithology, and biogeography, both across landscapes and with depth. It is paramount to understand and quantify the mechanisms, particularly the poorly understood processes occurring at depth, that influence the CZ's development and function (Brantley and Lebedeva, 2011; Giardino and Houser, 2015).

The biological importance of oxyanion metals is twofold as they can serve as micronutrients in a number of life-sustaining processes, but may also have toxic effects on organisms at elevated concentrations. For example, plant growth can be negatively impacted by total V concentrations greater than 0.8 mg L −1 in soil solution (Larsson et al., 2013) and total Cr water concentrations greater than 5 mg L−1 is considered hazardous to plants and animals (Singh et al., 2013). In this study, we investigate how lithology and climate influence the retention of oxyanions in the regolith and mobilization to stream water. We focus on molybdenum (Mo), chromium (Cr), and vanadium (V) because of their role as essential micronutrients and potential toxicity for organisms.

In addition to toxicity, the unique geochemistry and speciation patterns of oxyanions allows them to be proxies for various environmental parameters, particularly oxidationreduction reactions and weathering. Molybdenum is a key component for a number of enzymatic processes that occur during nitrogen and sulfur metabolism and is also responsible for the biosynthesis of plant hormones (Kaiser et al., 2005). Of particular interest, Mo is a cofactor for the nitrogenase enzyme required for nitrogen fixation and can become a limiting micronutrient in forest soils (e.g., Barron et al., 2008). Low concentrations of Mo are required by plants, with typical foliar values ranging between 0.03 and 1.61 mg kg−1 (Adriano, 2001). At concentrations exceeding 100 ppm in soils Mo can become toxic for plants, and at concentrations exceeding 10 mg kg−1 in vegetation humans and animals are susceptible to molybdenosis (Gupta, 1997; Adriano, 2001).

Vanadium (V) is another trace metal oxyanion that enhances plant growth and serves as a micronutrient for a variety of biological processes. At low concentrations (<2 mg kg−1 ), V has been shown to enhance chlorophyll synthesis and nitrogen fixation in plant tissues; however, higher concentrations may cause chlorosis and stunted growth (Kabata-Pendias, 2011). Due to similarities in their ionic radii, V can substitute for Ca and Mg, inhibiting chlorophyll generation, protein production, ATP formation, cellular pH, and enzymatic activities at high abundance within plant tissues (Cantley et al., 1977, 1978; Marschner et al., 1986; Crans et al., 1995; Larsson et al., 2013). Nevertheless, V toxicity in plants are uncommon due to its poor translocation from belowground to aerial tissues in plants (Adriano, 2001). There is no significant data to conclude that V is a necessary part of nutrition for humans (Imtiaz et al., 2015).

The biological importance of Cr is highly dependent on its oxidation state. Chromium(VI) does not serve any known roles for plants or animals, is a highly toxic carcinogen for humans (Kota´s and Stasicka, 2000; Dayan and Paine, 2001; Levina et al., 2003). Chromium (III), however, is essential for humans as it has been shown to be a cofactor for biologically active molecules that enhance the effects of insulin on target tissues; Cr(III) deficiencies include abnormal glucose utilization and increased insulin requirements (see Kota´s and Stasicka, 2000; Dayan and Paine, 2001). Chromium toxicity in plants is possible in very acidic soils with <5 pH, but is uncommon in soils with >5 pH due to low bioavailable concentrations (Adriano, 2001).

The retention and mobility of oxyanions (Mo, V, and Cr) are largely controlled by climate, weathering, and sorption reactions, all of which are variable in the Critical Zone. Molybdenum concentrations across all rock types range from 0.3 to 300 mg kg−1 (Adriano, 2001), and all soil types range from 0.01 to 17 mg kg−1 (Kabata-Pendias, 2011). Molybdenum exists in oxic soils primarily as the molybdate oxyanion (MoO4 2−, +VI oxidation state) with molybdic acid (H2MoO4) and its oxyanion (HMoO4 −) becoming increasingly present below pH 4, (Adriano, 2001; Alloway, 2013). Vanadium concentrations in rocks range from 20 to 250 mg kg−1 while soils have a broader range of 3–500 mg kg−1 depending on the initial concentration of their parent material (Adriano, 2001; Kabata-Pendias, 2011). Vanadium can be present as V(IV), mainly as VO2 (vanadyl), whereas V(V) exists as the oxyanions H2VO4 −1 or as HVO4 −2 (collectively referred to as vanadate), depending on the pH and oxidation status. Vanadate is considered the most toxic vanadium species (Martin and Kaplan, 1998; Larsson et al., 2013). Chromium concentrations in soils range from 2 to 150 mg kg−1 and is found mainly in two oxidation states, the trivalent Cr(III) and the hexavalent Cr(VI) (Adriano, 2001).

The terrestrial geochemistry of Mo, V, and Cr is a function of sorption processes and chemical speciation. Under oxic conditions, Mo, V, and Cr oxyanion sorption onto soil solids is determined by Fe and Al oxides, clay minerals and organic matter, of which most are pH dependent (e.g., Gäbler et al., 2009; Pérez et al., 2014; Wisawapipat and Kretzschmar, 2017). As an example, adsorption of MoO4 2− ions from soil solution onto positively charged metal oxides occurs between pH 4 and 5 (Riley et al., 1987; Gupta, 1997; Xu et al., 2013). In acidic soils, molybdate, vanadate, and chromate anions are sorbed onto positively charged Fe, Mn, and Al oxides as well as clay minerals and organic matter (Grove and Ellis, 1980; Pérez et al., 2014; King et al., 2016; Larsson et al., 2017; Wisawapipat and Kretzschmar, 2017). In addition, sorption to primary minerals such as inherited clay minerals like illite can also occur. In three Swedish soils, V(IV) was predominantly bound to primary minerals (Gäbler et al., 2009; Larsson et al., 2015; Wisawapipat and Kretzschmar, 2017). Oxyanion sorption to organic matter, Al and Fe oxides, and clay minerals in soils is important as it directly controls the retention and export of Mo, Cr, and V to stream waters (Kerr et al., 2008).

The influence of deep regolith—the unconsolidated layer of physically and chemically weathered rock that overlies unweathered bedrock—on oxyanion biogeochemistry in the Critical Zone is largely unconstrained below 2 m depth. Previous studies of Mo, Cr, and V sourcing to stream waters have primarily focused on surface soils <0.3 m depth (e.g., Salvador-Blanes et al., 2006; Gäbler et al., 2009; Wichard et al., 2009; Pérez et al., 2014; Morrison et al., 2015; Larsson et al., 2017; Wisawapipat and Kretzschmar, 2017) or only aqueous systems (e.g., Kerr et al., 2008; Gardner et al., 2017). However, regolith underlying soils at depths >1 m may be an important source of metals to surface waters and may also influence their immobilization. Here, we investigated the importance of regolith on oxyanion geochemistry at four Critical Zone Observatories (CZOs) developed on monolithic granitoid parent material, spanning a diverse climatic gradient (wet, hot tropical forest to cold, dry mountain-steppe). We hypothesized that CZOs with more secondary Fe oxide phases will have greater oxyanion sorption deeper within their profiles (>1 m depth) beyond the typical depth of surface soils (<0.3 m depth). Illuminating the role of the deep CZ is essential to capture its control on soil and surface water geochemistry.

#### MATERIALS AND METHODS

#### Critical Zone Observatories (CZOs) Boulder Creek CZO

Boulder Creek CZO is located in Arapahoe National Forest at the front range of the Rocky Mountains in north central Colorado, (**Figure 1**). Regolith samples were collected by borehole drilling in June 2014 from Gordon Gulch, a valley with low topographic relief, within the Boulder Creek CZO. Eight regolith samples were collected over a 7.0 m depth profile with samples collected every 0.5 m for the top 2.5 m and one deep sample at 7.0 m. Gordon Gulch is characterized by a montane-steppe climate with a mean annual temperature (MAT) of 4◦C and a mean annual precipitation (MAP) of 735 mm year−1 , with the majority of precipitation occurring between March and May (**Table 1**) (Dethier and Bove, 2011; Aguirre et al., 2017). The lithology of Gordon Gulch is dominated by biotite-gneiss (Dethier and Bove, 2011) and is forested with evergreen vegetation, primarily lodgepole pine (Pinus contorta)

and ponderosa pine (Pinus ponderosa) (Dethier and Bove, 2011; Barnes et al., 2012). Despite the low MAP, Boulder Creek CZO regolith is moderately weathered, on the basis of Ca, K depletion and presence of neoformed kaolinite (Condie et al., 1995; Dethier and Bove, 2011). Stream water samples were collected under baseflow conditions from the permanent flowing river in Gordon Gulch in March and April of 2017, using acid washed HDPE bottles following USGS water sampling techniques (United States Geological Survey [USGS], 2006).

#### Calhoun CZO

The Calhoun CZO is located in the Southern Piedmont in South Carolina (**Figure 1**), within the Calhoun Experimental Forest, and has been extensively studied to inform researchers about land and water degradation caused by agricultural uses (e.g., Richter and Markewitz, 2001). The sampling location at Calhoun CZO is on a low relief hill, derived from granitic-gneiss. The climate is subtropical with a MAT of 16◦C and a MAP of 1,250 mm year−1 (**Table 1**). Vegetation is dominated by loblolly pine (Pinus taeda), mixed hardwoods (primarily oak, Quercus spp.), and hickory (Carya spp.) (Richter and Markewitz, 2001). This study focused on samples collected every 0.6 m within a 7.6 m weathering profile collected from Research Area 1. Details of the collection and processing of drill core samples can be found in Bacon et al. (2012). In brief, soil samples (0–6.0 m) were collected from three cores while regolith was sampled with a three-wing auger bit (6.1–18.3 m). Samples were crushed, sieved to <2 mm, and air-dried. Bacon et al. (2012) has noted regolith at Calhoun CZO is intensively weathered, on the basis of Ca, K depletion in the upper 10 m. Stream water samples were collected under baseflow conditions from the nearest flowing surface water source, Watershed 4 stream, in March and April, 2017 using acid washed HDPE bottles following USGS water sampling techniques (United States Geological Survey [USGS], 2006).

#### Luquillo CZO

The Luquillo CZO is located in the El Yunque National Forest in northeastern Puerto Rico (**Figure 1**), and has also been studied as a Long-Term Ecological Research site and a USGS Water, Energy, and Biogeochemical Budgets (WEBB) site. The Luquillo CZO is a tropical montane forest and has a MAT of 21◦C (Murphy et al., 2012). Precipitation increases with elevation from 2,500 to 4,500 mm year−1 (Garcia-Martino et al., 1996). This paper focuses on the Guaba study site, a montane ridge that receives a MAP of 4,200 mm year−1 and is located 680 m.a.s.l. in the Rio Icacos watershed (**Table 1**) (Murphy et al., 2017). The most abundant tropical vegetation at the Guaba study site were Dacryodes excelsa, Prestoea montana, Sloanea berteroana, Cordia borinquensis, Manilkara bidentata, Schefflera morototoni, Cecropia schreberiana, Micropholis garciniifolia, Henriettea squamulosa, and Quararibea turbinate (Weaver and Gould, 2013). Surface soil (top 0.3 m of regolith) and upper regolith samples were collected in 2017 down to 1.8 m depth. Deeper regolith samples were subsampled in 2010 from archived material from site LG1 (Buss et al., 2005, 2017) and borehole LGW1 (Buss and White, 2012; Orlando et al., 2016). In total, eleven regolith samples were analyzed. Samples were spaced

TABLE 1 | Climate, ecosystem, and geologic properties of sampling sites for the four Critical Zone Observatories (CZO).

every 0.4 m down to 1.8 m depth, with three deeper samples collected at 3.1 m, 4.2 m, and 8.6 m depth. An intact regolith sample was collected from 1.8 m depth, at the surface soil-regolith interface. Regolith at Luquillo CZO is intensively weathered, on the basis of Fe and base cation depletion throughout the upper 6 m (Buss et al., 2005, 2017). Stream water samples were collected under baseflow conditions (White et al., 1998; Pett-Ridge et al., 2009) from the Rio Icacos in February and March of 2017 using acid washed HDPE bottles following USGS water sampling techniques (United States Geological Survey [USGS], 2006).

#### Southern Sierra CZO

The Southern Sierra CZO is located in Fresno County, California (**Figure 1**), within the Kings River Experimental Watershed (KREW), a long-term research site established by the Pacific Southwest Research Station of the United States Forest Service (Hunsaker and Eagan, 2003). The Southern Sierra CZO lies outside the limits of recent glaciation and is characterized by a sequence of parallel ridges and valleys, with alternating steep and gentle terrain of granitic bedrock. We focused on site P301 at the head of Providence Creek, which has a montane-steppe climate with a MAT of 9◦C and a MAP of 1,100 mm year−1 (**Table 1**). The topography is a 0.75 km wide valley with moderate to steep relief. Vegetation is dominated by a mixed-coniferous forest consisting of Abies concolor, Pinus ponderosa, Pinus jeffreyi, Quercus kelloggii, Pinus lambertiana, and Calocedrus decurrens, with intermixed chaparral. Eleven samples were collected using a Geoprobe R every 1.0 m down to 11.0 m in 2011 (Holbrook et al., 2014). Based upon the shallow depth of soil and high variability in regolith depth consistent with chemical weatheringlimited regolith formation, we consider regolith at Southern Sierra CZO to be weak to moderately weathered (Bales et al., 2011). Stream water samples were collected above baseflow conditions from Providence Creek 301 in March and April of 2017 (see Godsey and Kirchner, 2014) using acid washed HDPE bottles following USGS water sampling techniques (United States Geological Survey [USGS], 2006).

#### Elemental Analyses

#### Extractions and Total Digestions of Regolith Samples

A two-step sequential extraction was performed to quantify organic bound and secondary oxide fractions. Subsamples of regolith from the four CZOs were dried at 105◦C and then pulverized in a boron carbide mortar and pestle to <100 µm. Two grams of regolith was weighed into 50 mL centrifuge tubes. First, organic matter and amorphous oxide bound elements were collected by adding 5 mL of 30% H2O2 to oxidize organic matter and then extracting the soil solution with 10 mL of 0.01 M acetic acid to target amorphous and short-range-ordered Al-, Mn-, and Fe- (oxyhydr)oxides. The solution was centrifuged at 7,000 rpm for 45 min and decanted. This was repeated with two washes of 10 mL of deionized water, which were added to the total extract. For the secondary Fe oxide fraction of Mo, Cr, and V, the remaining soil was treated with citrate–bicarbonate– dithionite (CBD) following the methods described in Mehra and Jackson (1960), which focuses on the dissolution of crystalline Fe oxides such as goethite and hematite without dissolution of Febearing silicates. In brief, 40 mL of analytical grade 0.3 M sodium citrate and 5 mL of analytical grade 1 M sodium bicarbonate was added to the dry samples. The tubes were heated to 60◦C and 1 g of analytical grade sodium dithionite was added. After 45 min of reaction, samples were centrifuged at 3,500 rpm for 30 min and the supernatant was collected, filtered to <0.4 µm, and diluted with deionized water. For every batch of samples, a blank, duplicate, and a standard reference material (SRM) was included. Montana Soil 2711 from the National Institute of Standards and Technology was used as a standard reference material (SRM).

For total sample digestion, 100 mg of dried, pulverized regolith were dissolved in 5 mL of 25 M HF acid and 5 mL of distilled 16 M HNO3 (Trace metal grade, VWR Analytical, Radnor, PA, United States) at 130◦C in sealed 30 mL Teflon PFA vials for 48 h or longer as needed with quartz-rich samples. After complete digestion of visible materials, samples were dried to a powder and redissolved in 0.5 M HNO3. For every batch of samples, a blank, duplicate, and two SRMs [W-2 diabase and G-2 granite from the United States Geological Survey (USGS)] were included with every 20 samples to monitor reproducibility.

#### ICP-MS Analyses

Regolith digests and extracts were diluted with deionized water and analyzed for trace metals with an Element 2 ICP-MS (Thermo-Scientific, Waltham, MA, United States). Recoveries for total digests of Mo, V, Cr, Fe, and Ti were 81–103% of certified values. The variation between intra-sample duplicates was <6% and Ti

concentrations in the preparation blank samples were <0.01 µg g−1 .

Stream water samples were collected from each CZO under baseflow conditions between February and April 2017. Stream water samples were collected using acid-washed LDPE plastic bottles, filtered to <0.45 µm and acidified to pH 1 with trace metal grade HNO3. Stream water samples were analyzed for trace metals using an Agilent 8800 triple quadrupole (QQQ) ICP-MS at the Dartmouth College Trace Element Analysis Laboratory (Agilent Technologies, Santa Clara, CA, United States). The QQQ-ICP-MS was calibrated using NIST traceable single and multi-element standards containing Sc, Ge, Y, In, Rh, Te, and Bi as internal standards. Ten point calibration curves were constructed for each analyte with correlation coefficient criteria >0.99. The calibration was followed by an Initial Calibration Blank (IBC) and an Initial Calibration Verification (ICV). Acceptance criteria for the ICV and CCV are ±10%. The CCV was run after every 10 samples. Analytical duplicates and spikes were performed every 20 samples. Spike recoveries were 93–101% for Mo, 91–99% for Cr, and 95–102% for V. Coefficient of variations (CVs) for duplicates were <5% for Mo, Cr V, Fe, and Ti.

#### Micro X-Ray Fluorescence Mapping

The goal of the synchrotron-based Micro X-ray Fluorescence analysis (µ-XRF) was to elucidate whether V and Cr were located in primary minerals or secondary Fe oxide minerals. Molybdenum was not imaged due to its low concentration. Intact regolith samples were only collected from the Guaba Ridge site of the Luquillo CZO in February 2017 at the soil-regolith interface (1.8 m depth). The regolith sample was oven dried at 70◦C and cemented in place using a room-temperature-curing epoxy-resin (EPO-TEK 301-2FL). The hardened soil columns were cut into 15 mm thick sections and wet polished using silicon carbide starting at 240 grit. Minimizing the amount of deionized water and decreasing the total amount of time samples spent wet-polished were taken to limit dissolution or disturbance of minerals, particularly the weathered rind, and water soluble elements. The weathering zone of the regolith was analyzed by µ-XRF imaging at the Cornell High Energy Synchrotron Source (CHESS) on the F3 bending-magnet beamline. At CHESS, a 20 µm beam (incident energy of 17 KeV) with an incident flux of ∼3 × 109 photons s−1 at the sample was used. The detector was mounted at 90◦ with respect to the incident X-ray beam and samples were mounted at 45◦ to the incident beam and scanned through the X-ray beam. The Vortex silicon drift detector with a Digital X-ray Processor (DXP) for digital signal processing to detect fluorescent X-rays. X-ray spectra were fitted for elements and quantified using the PyMCA software package (Solé et al., 2007).

#### Total Carbon Quantification

Soil carbon data was obtained for each CZO: Calhoun (Fimmen et al., 2008; Bacon et al., 2012), Luquillo (Buss et al., 2005), and Southern Sierra (Tian, 2018). Total organic C (TOC) concentrations in soil and regolith samples were measured for Boulder Creek CZO samples using a Carlo-Erba Elemental Analyzer. In brief, 6 ± 1 mg subsamples <0.5 mm were combusted. Every 20 samples included one blank, one Atropine SRM, and a duplicate. TOC concentrations in Atropine SRMs were within 3% of its certified value and <10% relative percent difference. Soils were treated with 0.1 mM HCl acid to remove C from carbonates. Values for Boulder Creek CZO were comparable with surface soils from Eilers et al. (2012).

#### Data Analyses

We tested the associations of TOC and secondary Fe oxide with total Mo, V, and Cr concentrations and concentrations within secondary Fe oxide phase using linear regressions for each CZO. All regressions were performed using Matlab (Matlab, 2016, Mathworks, Natick, MA, United States). For elemental concentrations, error bars were propagated based upon the coefficient of variation for each sample measurement and the standard error measured through duplicate and triplicate samples.

Mass transfer coefficients (τ) of Mo, V, and Cr were calculated to assess the degree of chemical weathering for specific elements using Ti as an index element (Brantley and Lebedeva, 2011). Mass transfer coefficients are defined for the element of interest (Mo, V, and Cr) (j) in the soil or regolith sample normalized to an immobile index element Ti (i) in the parent material:

$$\pi\_{\mathbf{j},\mathbf{w}} = \frac{C\_{\mathbf{j},\mathbf{w}}C\_{\mathbf{i},\mathbf{p}}}{C\_{\mathbf{j},\mathbf{p}}C\_{\mathbf{i},\mathbf{w}}} - 1 \tag{1}$$

where Cj,w and Ci,w are total concentrations of the element of interest (Mo, V, and Cr) and Ti in the soil and regolith samples, respectively, and Cj,p and Ci,p are the element of interest and index element concentrations in the least weathered regolith material (e.g., Brimhall and Dietrich, 1987; Kraepiel et al., 2015). Negative τ values indicate a net loss of an element from a specific soil depth due to leaching. Positive τ values indicate a net enrichment of an element. Although Ti can be mobilized at strongly acidic, heavily weathered tropical soils (Kurtz et al., 2000), application of Ti-based mass transfer coefficients is effective for estimating relative losses and gains in the CZOs of our study.

Stream water concentrations were grouped across collection dates and compared across the four CZOs using the Kruskal–Wallis non-parametric tests. Using measured watershed area and mean annual stream discharge, we calculated the export of Mo, V, and Cr at each of the four CZOs. To accomplish this, we obtained the total watershed area, MAP, and mean annual discharge for each of the four CZOs from existing online data1 (accessed May 2018). Stream water Mo, V, and Cr concentrations were multiplied by total discharge to estimate export rates for each element.

1 criticalzone.org

### RESULTS

#### Regolith Profiles

feart-06-00193 November 8, 2018 Time: 16:14 # 6

Regolith profiles exhibited a range of TOC and secondary oxide Fe with depth across the four CZOs. Overall, %TOC decreased with depth from approximately 1% down to 0.1% with subsurface peaks present at all four CZOs (**Figure 2**). Calhoun CZO had <0.1 %TOC below 2 m depth, which was lower than the other three CZOs (**Figure 2**). Total Fe was relatively constant with depth at Boulder Creek and Southern Sierra CZOs. Subtropical Calhoun CZO had low total Fe concentrations in the surface soils

FIGURE 2 | Regolith profiles of Mo, V, and Cr concentrations in total digests and two extractions of secondary Fe oxide bound (citrate–bicarbonate–dithionate) and organic + amorphous oxide bound at the four CZOs. Error bars are coefficients of variation for each measurement. Please note each scale is different for each phase fraction.

while tropical Luquillo CZO had higher total Fe concentrations above 2 m depth (**Figure 2**). Secondary oxide Fe concentrations were highest throughout the regolith profile at Luquillo CZO by nearly two orders of magnitude compared to Boulder Creek and Southern Sierra CZOs (**Figure 2**). Calhoun CZO had a subsurface peak of secondary oxide Fe concentration between 1 and 1.5 m depth (**Figure 2**).

Total Mo profile concentrations ranged from 0.4 to 2.5 mg kg−1 across the four CZOs (**Figure 3** and **Supplementary Table 1**). The highest total Mo concentrations at the Calhoun and Luquillo CZOs were in the surface soils, at 0.6 m and 0.4 m in depth, respectively. The colder, higher elevation Boulder Creek and Southern Sierra CZOs had their highest total Mo concentrations deeper in the regolith, at 2.4 m and 4 m in depth, respectively. Total V profiles exhibited two orders of magnitude variation in concentrations across regolith profiles. Total V concentrations were highest in the top meter for Boulder Creek, Luquillo, and Calhoun CZOs, and were generally consistent with depth at the Southern Sierra CZO. Total Cr concentrations at the Calhoun CZO had subsurface peaks of total Cr between 2.0 and 4.0 m. Boulder Creek exhibited a significant total Cr peak between 1.0 and 2.0 m, which decreased with depth. Southern Sierra and Luquillo CZO exhibited nearly homogeneous (±3 mg kg−1 ) total Cr profiles with depth.

In the regolith, organic matter and amorphous oxide bound Mo, V, and Cr were substantially lower than their respective secondary oxide and total concentrations. Across the four CZOs, organic matter and amorphous oxide bound Mo ranged between 0.01 and 0.13 mg kg−1 (**Supplementary Table 2**), which on average was 0.5–0.9% of the total Mo concentration (**Supplementary Table 4**). The organic matter and amorphous oxide bound fraction of Mo were greatest near the surface (<1 m depth) for Calhoun and Southern Sierra CZOs while Boulder Creek and Luquillo CZOs had subsurface peaks (>2 m depth). Vanadium varied between 0.11 and 3.1 mg kg−1 (**Supplementary Table 2**), which on average was 0.3–2.1% of the total V concentration (**Supplementary Table 4**). The organic matter and amorphous oxide bound fraction of V were greatest near the surface (<1 m depth) for Luquillo and Southern Sierra CZOs while Boulder Creek and Calhoun CZOs exhibited subsurface peaks (**Figure 3**). Organic bound Cr concentrations ranged between 0.01 and 1.07 mg kg−1 (**Supplementary Table 2**), which on average was 0.9–3.5% of the total Cr concentration (**Supplementary Table 4**). Organic matter and amorphous oxide bound Cr was greatest in the top 2 m of regolith at Luquillo and Southern Sierra CZOs but were greater below 2 m in depth at Boulder Creek and Calhoun CZOs (**Figure 3**).

Secondary Fe oxide Mo concentrations ranged between 0.02 and 0.68 mg kg−1 (**Supplementary Table 3**), which on average was 6.6–40.1% of the total Mo concentration (**Supplementary Table 5**). Vanadium concentrations in the secondary Fe oxide ranged between 1 and 52 mg kg−1 (**Supplementary Table 3**), which on average was 6.0–39.7% of the total V concentration (**Supplementary Table 5**). Secondary Fe oxide Cr concentrations varied between 0.2 and 8.7 mg kg−1 , which on average was 13.0–40.5% of the total Cr concentrations (**Supplementary Table 5**).

Using Pearson correlations, we observed that total Mo, V, and Cr were correlated with secondary oxide Fe concentrations [a proxy for long-range ordered to crystalline Fe oxy(hydr)oxides] (**Table 2**). At Calhoun CZO, secondary oxide and total Mo, V, and Cr concentrations were positively correlated with secondary oxide Fe concentrations. For Southern Sierra CZO, secondary oxide V and Cr were positively correlated with secondary oxide Fe concentrations. At Boulder Creek CZO, secondary oxide and total concentrations of V and Cr were negatively correlated with both secondary oxide Fe and TOC. At Luquillo CZO, secondary oxide Mo and V were positively correlated with %TOC.

#### Mass Transfer Coefficients and µ-XRF Mapping

Calculated tau values for Mo, V, and Cr show both net enrichment and depletion across the CZOs (**Figure 4**). At the montane coniferous Boulder Creek CZO, regolith down to 2.0 m depth were depleted in Mo and V, but Cr was enriched below 1.0 m

TABLE 2 | Pearson correlation coefficients for linear regressions between secondary Fe oxide and %TOC with secondary Fe oxide and total Mo, V, and Cr at each of the four CZOs.

∗∗ indicates P < 0.01 and ∗ indicates P < 0.05.

depth. At the subtropical Calhoun CZO, V was depleted while Mo and Cr were enriched throughout the regolith profile. The tropical Luquillo CZO showed a wide variation but overall Cr was depleted throughout the entire profile while Mo and V were enriched at certain depths. The Southern Sierra CZO regolith profile was depleted in V and Cr but enriched in the surface and subsurface between 2.0 and 5.0 m depth.

Analysis of the µ-XRF spectra for the weathered regolith from Luquillo CZO showed Fe, V, and Cr coinciding together in hornblende and biotite primary minerals (Buss et al., 2005, 2017). This is particularly observable in the upper left-hand corner of the images and the outlines of the hornblende minerals are distinguishable in all four images. Moreover, Fe, V, and Cr also correspond with each other in the Fe oxide weathering rind on the right side of the image (**Figure 5**). More information on the weathering rinds of quartz diorite from Luquillo CZO is available in Buss et al. (2008).

#### Stream Water Concentrations

Overall, stream water concentrations of oxyanions were generally comparable across CZOs (**Table 3** and **Figure 6**), but there were a few significant differences. Stream water V concentrations were similar across CZOs and ranged between 0.04 and 0.70 µg L−1 (**Figure 6**). Stream water Cr concentrations at Calhoun and Boulder Creek CZOs were significantly than at the Southern

Sierra and Luquillo CZOs (**Table 3** and **Figure 6**) using the nonparametric Kruskal–Wallis Test (p < 0.05). Stream water Mo concentrations were significantly higher at Boulder Creek and Calhoun CZOs compared to Southern Sierra CZO (p < 0.05, **Table 3** and **Figure 6**).

We compared the concentration of oxyanions in the stream water to the deepest, least weathered regolith to determine if variations in lithology explained the differences in stream water concentrations across CZOs and among oxyanions. In **Figure 7**, Mo stream water to regolith ratios were significantly greater at subtropical Calhoun and tropical Luquillo CZOs than montane Southern Sierra and Boulder Creek CZOs. In addition, Mo stream water to regolith ratios at Calhoun, Luquillo, and Southern Sierra CZOs were significantly greater than V and Cr ratios. Stream water to regolith ratios for V were not significantly different for Boulder Creek, Calhoun, and Luquillo CZOs but were greater than at Southern Sierra CZO. In addition, Stream water to regolith ratios for V were generally lower than Mo and Cr across CZOs (**Figure 7**). Stream water to regolith ratios for Cr were significantly higher at Calhoun CZO than the other three CZOs. Similarly, stream water to regolith ratios for Cr were significantly higher at Boulder Creek CZO than Luquillo and Southern Sierra CZOs.

#### DISCUSSION

#### Weathering and Partitioning of Mo, V, and Cr

Our results support the hypothesis that deep regolith play a substantial role in sorption of oxyanions. Overall, both surface soils (<0.3 m depth) and deeper regolith served as important sources and retainers of oxyanions. We observed many peaks in Mo, V, and Cr in both the organic matter and amorphous oxide and the secondary Fe oxide extractions at depths are far deeper than typical 0.3–1 m depth investigated for oxyanion sorption and weathering in soils (e.g., Shiller and Mao, 2000; Wichard et al., 2009; Wisawapipat and Kretzschmar, 2017). For example, at the montane regolith of Southern Sierra CZO had >39% of Mo and Cr were in the secondary Fe oxide fraction from 3 to 4 m in depth and the tropical forest Luquillo CZO, 52% of Cr and 23% of V were in the secondary Fe oxide fraction at 3 m in depth.

In addition, enrichment and depletion of oxyanions >2 m in depth indicated that regolith played both a significant role in the sorption and retention of oxyanions (**Figure 4**). Except for Cr, Boulder Creek CZO demonstrated the most depletion of oxyanions, extending down 2 m in depth (**Figure 4**). This is particularly important as the mountain-steppe Boulder Creek CZO has the least intensive weathering environment with the lowest MAP and MAT, but due to its long exposure, has been moderate to extensively weathered. The more tropical regolith profiles at Calhoun and Luquillo CZOs show much more complicated patterns of accumulation and depletion (**Figure 4**), likely due to variations in retention as described later in this section.

These results support other weathering observations by Dethier and Bove (2011) at Boulder Creek CZO, Bacon et al. (2012) at Calhoun CZO, and Buss et al. (2005) at Luquillo CZO. Subtropical Calhoun CZO and tropical Luquillo CZO were the most weathered regolith profiles and had the highest accumulations of secondary Fe oxides (**Figure 2**), consisting of hematite and goethite as observed by Buss et al. (2005); Fimmen et al. (2008), and Dethier and Bove (2011). In addition, neoformed kaolinite was present at the previously mentioned CZOs at >4% by mass (Fimmen et al., 2008 at Calhoun CZO; Dethier and Bove, 2011 at Boulder Creek CZO; Buss et al., 2017 at Luquillo CZO) demonstrating their extensive weathering history. The extent of weathering corresponded with the sorption of oxyanions throughout the regolith profiles. Southern Sierra had

Averages are given ± one standard error. †The flux estimates are based upon baseflow oxyanion concentrations. Here, we assumed oxyanions were isostatic, which may not necessarily be applicable for each element at all of the CZOs.

lower secondary Fe oxide fractions of Mo, Cr, and V compared to the other three CZOs. Moreover, secondary oxide Mo, Cr, and V were greater in regolith with higher secondary Fe oxide concentrations at much warmer and wetter Calhoun and Luquillo CZOs when compared to the colder and drier Boulder Creek CZO. Thus, we conclude that weathering has influenced the formation and retention of oxyanions with secondary Fe oxides.

A generalized trends across the three oxyanions were present, demonstrating some linkages between weathering intensity and climate with the sorption of oxyanions to secondary Fe oxides. Total Mo, secondary Fe-bound Mo and organic and amorphous oxide bound Mo concentrations were comparable across the four CZOs, suggesting weathering intensity nor lithology did not strongly influence Mo in the soil and regolith profiles. Secondary oxides did not appear to play a dominant role for Mo retention at Luquillo CZO as total Mo was not correlated with secondary Fe oxides. Instead secondary oxide Mo was correlated with %TOC despite the fact that organic matter and amorphous oxide bound Mo was not a substantial fraction of the total Mo concentrations at any of the CZOs, generally ranging below 5%. This largely agrees with previous studies that the sorption and retention of Mo is dominated by both, secondary Fe oxides and organic compounds (e.g., Wichard et al., 2009; King et al., 2016, 2018). Due to the relatively small range in %TOC, there may not have been enough variability to capture the relationship between Mo and TOC.

Total V concentrations varied significantly across the CZOs (**Supplementary Table 1**), suggesting lithologic controls on the sourcing of V to the soil and regolith profiles. Secondary Fe oxide fraction of V, however, coincided with weathering intensity. Average secondary Fe oxide concentrations were highest for the more tropical Calhoun (14 mg kg−1 ) and Luquillo CZO (26 mg kg−1 ) compared to mountain-steppe Southern Sierra (9 mg kg−1 ) and Boulder Creek CZOs (3 mg kg−1 ). Vanadium exhibited significant linear relationships with secondary Fe oxide at Calhoun, Boulder Creek, and Southern Sierra CZOs. The substitution of V(III) for Fe(III) in minerals has been widely observed (Schwertmann and Pfab, 1996; Shiller and Mao, 2000; Wisawapipat and Kretzschmar, 2017) and V(II) and V(III) for Fe(II) and Fe(III) in primary Fe-silicate minerals during magmatic crystallization (e.g., Goldschmidt, 1954). Even though this relationship was not clearly observed at Luquillo CZO as with the other three CZOs using linear regressions with secondary Fe oxide, we did observe this relationship in µ-XRF maps of regolith samples. Our µ-XRF imaging of V in an intact weathering rind sample from Luquillo demonstrates its spatial distribution with Fe, both in primary minerals like hornblende and biotite and in secondary Fe oxides within regolith at Luquillo CZO.

Chromium concentrations exhibited a pattern similar to V. Total Cr concentrations were varied strongly with lithology but Cr concentrations in secondary Fe oxides appeared to be climate influenced. The more tropical Calhoun and Luquillo CZOs had higher % Cr in the secondary Fe oxide fraction than the montane-steppe Southern Sierra and Boulder Creek CZOs, suggesting intensive weathering promotes Cr sorption by secondary Fe oxides. The sorption of Cr by secondary Fe oxides and organic matter has been well observed by previous studies (e.g., Gustafsson et al., 2014; Pérez et al., 2014; Saputro et al., 2014).

Some potential variables that may explain the variation of oxyanions not explored here are the effects of erosion, atmospheric inputs, hydrological properties, and other source of sorption such as clay minerals and other oxides. Luquillo CZO differs from the other three CZOs due to its substantially greater annual precipitation (**Table 1**), which could not only flush and dissolve solutes at rates that may prevent sorption (see Dosseto et al., 2012) but also lead to increased atmospheric inputs of oxyanions in the surface soils (Marks et al., 2015; King et al., 2016). This may also influence τ values by enriching oxyanions at the surface and a net depletion due to leaching losses with depth. Nevertheless, the effect of atmospheric inputs on the deep regolith is likely minimal. In addition, Luquillo CZO undergoes more erosion in the form of landslides and surface soil erosion than the other three CZOs (Larsen, 2012). Lastly, Luquillo CZO also had much greater clay content in its soils and upper regolith than Calhoun, Boulder Creek, and Southern Sierra CZOs. Clays can provide an unaccounted source of oxyanion sorption capacity (Goldberg et al., 1996).

Our study observed low fractions of oxyanions in the organic bound and amorphous oxide phase and poor correlations

between %TOC and oxyanion concentrations in the deep regolith. However, this conclusion should not be extrapolated to surface soils. Many previous studies have observed a significant relationship in surface soils and oxyanion concentrations (see Grove and Ellis, 1980; Wichard et al., 2009; King et al., 2014, 2016, 2018; Pérez et al., 2014; Marks et al., 2015; Siebert et al., 2015; King, 2017; Larsson et al., 2017; Wisawapipat and Kretzschmar, 2017). Therefore, sorption of oxyanions with organic matter may be most important spatially but not with depth into deep regolith (>2 m).

#### Stream Water Export of Mo, V, and Cr

Oxyanion concentrations and export from the Critical Zone via stream water varied substantially across CZOs. The most intensively weathered regolith had higher average Mo

concentrations in stream water (**Table 3** and **Figure 6**). However, Mo concentrations were generally comparable across CZOs when removing the outlier samples at Calhoun and Luquillo CZOs. Calhoun and Luquillo CZOs had higher Mo stream water to regolith ratios compared to V and Cr (**Figure 7**). This implies that Mo had lower relative retention in regolith, increasing its export by stream water compared to V and Cr at Calhoun and Luquillo CZOs. Alternatively, the stream water to regolith ratios increase with increasing MAP across the CZO (**Table 1** and **Figure 7**). This suggests that atmospheric inputs of Mo via precipitation (Marks et al., 2015; Siebert et al., 2015; King et al., 2016) may lead to increased export of Mo from the regolith. There may be additional transport processes influencing Mo transport such as colloidal Al, Fe, and Mn oxides and dissolved organic matter (Pearce et al., 2010; Trostle et al., 2016; Prouty et al., 2017), which can enhance its transport through soil. Furthermore, fine colloids <0.45 µm in solution may also be responsible for the two elevated Mo concentrations at Calhoun and Luquillo CZOs. When these high concentrations are not considered, Mo concentrations were similar. Thus, there is not a strong trend in Mo concentrations with weathering intensity but potential variability due to colloids was greatest at the most intensively weathered CZOs.

Few studies have examined the controls on V concentrations in stream water (Shiller and Boyle, 1987; Shiller, 1997; Shiller and Mao, 2000; Strady et al., 2009). Stream water V concentrations were not significantly different across the four CZOs (**Table 3**). The similarity in stream water to regolith ratio for V suggests that regolith at Boulder Creek, Calhoun, and Luquillo CZOs are weathering and exporting V at consistent rates. Further, the stream water to regolith V ratios were lower than Mo and Cr, suggesting that V is the least mobile of the three oxyanions during weathering (**Figure 7**). The elevated stream water V concentration measured in one sample at Luquillo CZO may be due to colloids, particularly DOC and Al and Fe particulates. Stream water to regolith ratios for V was not significantly different than Fe for each CZO, suggesting V export is directly

controlled with Fe movement. This may be as either dissolved or fine colloids <0.45 µm, the size of the filter pores used in this study (e.g., Pokrovsky and Schott, 2002; Gardner et al., 2017). Stream water Cr concentrations exhibited the greatest variation among the four CZOs. Boulder Creek CZO had the highest stream water Cr concentrations (**Table 3**). Moreover, Boulder Creek CZO had significantly higher stream water to regolith ratio, despite having the lowest MAT and lowest MAP compared with the other four CZOs in this study (**Figure 7**).

Using measured watershed area and mean annual stream discharge, we calculated the export of Mo, V, and Cr at each of the four CZOs. Luquillo CZO had Mo, V, and Cr fluxes (105 to 106 mg year−1 ) up to four orders of magnitude greater than the other CZOs (102 to 104 mg per year) due to 102 to 103 times greater stream discharge than the other three CZOs. We observed that Boulder Creek CZO had greater Mo, V, and Cr fluxes than Southern Sierra CZOs, despite having lower MAP. Even when scaled for different watersheds size, Boulder Creek had >10× greater export of oxyanions per km2 . This highlights the importance that more weathered systems (based upon interpretation of mass transfer coefficients) export more oxyanions at greater rates than fresher regolith. It must be noted these approximations are based upon oxyanion concentrations in baseflow conditions. While this assumption might not hold for each of the oxyanions at the four CZOs, previous studies such as Shiller and Mao (2000), Aguirre et al. (2017) have found oxyanions can be isostatic with respect to concentrationdischarge (CQ) dynamics. However, this provides an estimate in the variation of Mo, V, and Cr export via stream waters across different watershed with contrasting climates.

#### CONCLUSION

Our study has shown that retention and sourcing of Mo, V, and Cr from regolith extends far below the depth that most typical soil studies examine. Depletion of Mo, V, and Cr extends several meters in depth at all four CZOs, from the montanesteppe Boulder Creek CZO to the tropical rainforest Luquillo CZO. Moreover, secondary Fe oxides containing Mo, V, and Cr are present beyond surface soils, extending down >8 m in depth and were important fraction of the total percentage. In particular, >25% of total Cr and V were in the secondary Fe oxide fraction at Luquillo and Calhoun CZO between 0 and 2 m depth.

Each of the three oxyanions exhibited a unique pattern in stream water samples. Total concentrations and regolith to

#### REFERENCES

stream water ratios of Mo and Cr exhibited significant variation among the four CZOs while V concentrations were similar across CZOs. Our regolith and stream water results highlight that each oxyanion has different geochemical pathways but share unifying principles across weathering and climate gradients. Identifying the upland causes for different oxyanion retention is key for developing Mo, V, and Cr as geochemical tools for understanding contemporary and paleo-terrestrial and marine global biogeochemical cycles.

#### AUTHOR CONTRIBUTIONS

JR conducted laboratory experiments, analyzed, and interpreted data and shaped discussion. EK analyzed and interpreted data and shaped discussion.

#### FUNDING

This work was made possible by the National Science Foundation Grant (NSF-1360760) to the Critical Zone Observatory Network National Office. Funding from the National Science Foundation for each CZO was also instrumental for each CZO studied in this project: Boulder Creek (NSF-1331828), Calhoun (NSF-1331846), Luquillo (NSF-1331841), and Southern Sierra (NSF-1331939).

#### ACKNOWLEDGMENTS

We are thankful for help with obtaining samples from Erin Stacy (UC Merced), Stephen Hart (UC Merced), Scott Hynek (Penn State), Sue Brantley (Penn State), and Andy Kurtz (Boston University). We are thankful for the UC Angelo Coast Range Reserve and the analytical help from Brian Jackson (Dartmouth College). CHESS is supported by the NSF and NIH/NIGMS via NSF award DMR-1332208, and the MacCHESS resource is supported by NIH/NIGMS award GM-103485.

#### SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart. 2018.00193/full#supplementary-material

Amundson, R., Richter, D. D., Humphreys, G. S., Jobbágy, E. G., and Gaillardet, J. (2007). Coupling between biota and earth materials in the critical zone. Elements 3, 327–332. doi: 10.2113/gselements.3.5.327

Bacon, A. R., Richter, D. D., Bierman, P. R., and Rood, D. H. (2012). Coupling meteoric 10Be with pedogenic losses of 9Be to improve soil residence time estimates on an ancient North American interfluve. Geology 40, 847–850. doi: 10.1130/G33449.1

Bales, R. C., Conklin, M. H., Kerkez, B., Glaser, S., Hopmans, J. W., Hunsaker, C. T., et al. (2011). "Sampling strategies in forest hydrology and biogeochemistry," In Forest Hydrology and Biogeochemistry: Synthesis of Past Research and Future

Directions, eds D. F. Levia, D. Carlyle-Moses, and T. Tanaka (Dordrecht: Springer), 29–44. doi: 10.1007/978-94-007-1363-5\_2

of bacterial and archaeal communities in soil. Soil Biol. Biochem. 50, 58–65. doi: 10.1016/j.soilbio.2012.03.011

Goldschmidt, V. M. (1954). Geochemistry. Oxford: Oxford University Press, 730.

**Conflict of Interest Statement:** The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2018 Richardson and King. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Solute Fluxes Through Restored Prairie and Intensively Managed Critical Zones in Nebraska and Iowa

Ashlee L. Dere1 \*, Andrew W. Miller2 , Amy M. Hemje2 , Sara K. Parcher1 , Courtney A. Capalli3 and E. Arthur Bettis III3

1 Department of Geography/Geology, University of Nebraska Omaha, Omaha, NE, United States, 2 Department of Chemistry, University of Nebraska Omaha, Omaha, NE, United States, 3 Department of Earth and Environmental Sciences, University of Iowa, Iowa City, IA, United States

#### Edited by:

Alexandra V. Turchyn, University of Cambridge, United Kingdom

#### Reviewed by:

Justin Byron Richardson, University of Massachusetts Amherst, United States Joel Moore, Towson University, United States

\*Correspondence:

Ashlee L. Dere<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Biogeoscience, a section of the journal Frontiers in Earth Science

Received: 04 September 2018 Accepted: 04 February 2019 Published: 21 February 2019

#### Citation: