Does Stream Water Composition at Sleepers River in Vermont Reflect Dynamic Changes in Soils During Recovery From Acidification?

Stream water pH and composition are widely used to monitor ongoing recovery from the deposition of strong anthropogenic acids in many forested headwater catchments in the northeastern US. However, stream water composition is a function of highly complex and coupled processes, flowpaths, and variations in soil and bedrock composition. Spatial heterogeneity is especially pronounced in headwater catchments with steep topography, potentially limiting stream water composition as an indicator of changes in critical zone (CZ) dynamics during system recovery. To investigate the link between catchment characteristics, landscape position, and stream water composition we used long-term data (1991-2015) from the Sleepers River Research Watershed (SRRW) in northeastern Vermont. We investigated trends with time in stream water and trends with time, depth, and landscape position (upslope, midslope, and riparian zone) in groundwater (GW) and soil solution. We further determined soil elemental composition and mineralogy on archived (1996) and modern (2017) soil samples to assess changes in composition with time. SRRW is inherently well-buffered by calcite in bedrock and till, but soils had become acidified and are now recovering from acidification. Although base cations, especially Ca, decrease progressively with time in GW, riparian soils have become more enriched in Ca, due to a mixture of lateral and vertical transfers. At the same time stream water Ca fluxes increased over the past two decades, likely due to the leaching of (transient) legacy Ca from riparian zones and increased water fluxes. The stream water response therefore reflects the dynamic changes in soil chemistry, flow routing and water inputs.

Acid deposition forms when emission-derived sulfur dioxide and nitrogen oxides interact with precipitation, which impacts the entire CZ through a multitude of complex and coupled processes, including declining tree and soil quality due to low pH and base cation leaching (Matzner and Murach, 1995;Driscoll et al., 2001;Raddum et al., 2007). The ensuing stream water signals are typically sensitive to these effects and showed low pH, increased concentrations of acid anions (sulfate and nitrate), and increased effluxes of base cations (e.g., Ca and Na) from many acid-impacted systems (Newell and Skjelkvåle, 1997;Driscoll et al., 2001;Garmo et al., 2014). In the United States, legislation passed in 1970 and 1990 led to reductions in emissions that caused a slow increase in precipitation pH as well as a decrease in the acid anion content across affected regions, and stream water composition has since been monitored closely for signs of recovery (Newell and Skjelkvåle, 1997;Driscoll et al., 2001;Skjelkvåle, 2003;Rice and Herman, 2012;Rogora et al., 2013;Garmo et al., 2014;Kopáček et al., 2016).
However, because streams merge water and solutes from various sources, the attribution of changes in stream water chemistry to specific process locations without detailed investigations of the CZ is difficult. For example, one CZ process strongly affected by the presence of strong acids (mostly H 2 SO 4 and HNO 3 ) is weathering, i.e., the breakdown of bedrock and regolith through physical and chemical processes (Jin et al., 2010). Chemical weathering includes simple dissolution reactions (typical for carbonates) and hydrolysis (typical for silicates). Both processes are driven by proton availability and in the case of most silicate minerals weathering reactions are incongruent, i.e., they generate dissolved weathering products and form secondary minerals (Strawn et al., 2015). As such, weathering and soil development continuously change the composition of the solid CZ and the composition of percolating waters [e.g., soil solution and groundwater (GW)]. Another important process is the leaching of base cations from acid-impacted soils that lead to declines in soil health. These changes can be clear in headwater catchments with pronounced topography, even within just a few decades of changes in precipitation chemistry due to the transfer of soluble materials to deeper soil layers, GW (vertical transfer), and/or via lateral transfer to low-lying landscape positions (Johnson et al., 2000;Nezat et al., 2004;Brantley et al., 2007;Bailey et al., 2014;Lawrence et al., 2015;Lybrand and Rasmussen, 2015).
Stream water composition is therefore strongly affected by flowpaths that intersect different layers of a heterogeneous CZ. Accordingly, backtracking to identify specific processes requires information on soils as well as subsurface waters. We attempt this backtracking to identify signals of recovery in streams and soils using the well-studied Sleepers River Research Watershed (SRRW) in northeastern Vermont as a testbed. SRRW received strong acids through wet and dry deposition but, unlike many other watersheds in the NE, contains carbonate-bearing parent material that buffers the pH of percolating water. The impact of acid deposition on the CZ is therefore less dramatic but potentially visible due to the sensitivity of carbonates to only small changes in acid inputs.
In order to investigate how the CZ at SRRW responds to changes in precipitation chemistry in the stream signal and how compositional differences by landscape position influence this signal, we analyzed long-term stream water (1991( -present), soil solution (2004( -2013( ), and GW composition (2004( -2013 at various landscape positions. Data analyzed include: pH, acid anions (nitrate and sulfate), base cations (Na and Ca) and Si. To test for changes in soil composition with time we investigated soil elemental composition (Si, Al, Na, and Ca) and mineralogy data from archived (1996) and modern (2017) soil samples of varying landscape position.
Because silicate weathering at the watershed scale is a slow process, we hypothesize that changes in acid deposition are not recorded as significant trends of silicate weathering products in the stream over our 24-year study period. However, we hypothesize that carbonate weathering and base cation leaching slowed, leading to decreases in Ca flux in the stream since 1991. We further hypothesize that topographically induced transfer of base cations led to depletions in shallow horizons at upslope and midslope locations while near-stream soils led to enrichments in base cations, especially Ca.

Field Site
The SRRW comprises a series of nested catchments with varying land use in northeastern VT (Figure 1). The forested headwater watershed, W-9, is the focus of this study and is underlain by the Waits River Formation, which includes a quartz-mica phyllite member interbedded with calcareous granulite (Ferry, 1991;Ratcliffe et al., 2011). This formation is covered by up to 4 m of dense silty basal till from the last glaciation, (Shanley, 2000) which itself is sourced from the underlying Waits River Formation (53%) and the nearby Gile Mountain Formation (25%; Hornbeck et al., 1997). The presence of carbonates in bedrock and till contributes to the high buffering capacity and high pH in stream waters, despite historic acid deposition (Hornbeck et al., 1997;Shanley, 2000).
Annual precipitation averages 1320 mm, of which up to 30% falls as snow. The average maximum snow water equivalent approaches 280 mm, and the annual spring snowmelt is a dominant hydrologic event (Shanley et al., 2004). Decadeslong research on the hydrology and biogeochemistry of W-9 has provided insights into water routing and stream flow contribution (Hornbeck et al., 1997;Shanley et al., 2015). Saturation-excess overland flow is an important contribution to streamflow generation, stored catchment water typically controls stream response, and convergent landscape positions especially contribute to event stream waters (McGlynn et al., 1999;Shanley et al., 2015). Because this area was a till deposition zone for the SE-moving Wisconsinan glaciation, the shallow aquifer is dominated by till with fine silt texture and a low saturated hydraulic conductivity that increases from near 0.1 mm/hr within  Kendall et al., 1999;Shanley, 2000, with permission).
the till to as much as 100 mm/hr in some surficial soils . Streamflow has a relatively flashy response to precipitation, but water stored in the till sustains base flow such that W-9 has never ceased to flow, while comparable size watersheds in the region, including Hubbard Brook New Hampshire, routinely dry up in summer. This stored water helps contribute to an "old-water" dominated response in streamflow during events (Shanley et al., 2002). Elevation at W-9 ranges from 519 to 678 m, and the terrain presents a mix of slopes of varying steepness and gently dipping riparian areas and wetlands. Because of this setting, landscape position strongly influences soils in W-9: Inceptisols (most prevalent) and Spodosols dominate on well-drained upland areas such as upslope and midslope areas (Shanley et al., 2004). Histosols are mostly found in poorly drained riparian areas and are often located where GW upwelling occurs (Shanley et al., 2004). The catchment is completely forested with second-growth northern hardwoods such as sugar maple, yellow birch, white ash and American beech, but softwood species (balsam fir and red spruce) are also present (Shanley et al., 2004). Like other catchments in the Northeastern United States, SRRW has been impacted by acid deposition that caused forest decline (branch dieback, some declining crown vigor) and low soil pH (Shanley et al., 2004). However, due to emissions regulations, SRRW is experiencing reduced acid inputs (Stoddard et al., 1999;Burns et al., 2005). In particular, sulfate concentrations in precipitation have been decreasing at SRRW (Shanley et al., 2004) and the NE more generally (Miles et al., 2012). Stream water dissolved organic matter has been measured weekly at SRRW since 1991 and increased steadily by 0.027 mg L −1 yr −1 (Cincotta et al., this special issues).

Sampling Methods
In order to assess trends in soil-, ground-and stream water composition, we used existing data on discharge, pH, sulfate, nitrate, Si, Ca, and Na that have been collected weekly and during events at SRRW over 24 years (chemistry database doi: 10.5066/P9380HQG, flow database doi: 10.5066/ P929KMVK). For soil solution, we used data from 3 sets of nested zero-tension lysimeters that occupy three main landscape positions (upslope, midslope, and base of the hill) and sample at 15 cm (shallow), and 50 cm (deep). Groundwater wells screened at various depths (Table 1) occupy the same landscape positions as lysimeters and were used to monitor changes in GW composition. Both soil solution and GW compositions were determined on approximately biweekly samples between 2004 and 2013 but due to intermittently dry conditions, soil solution samples were obtained less frequently. Stream water composition for the same solutes was available from approximately weekly grab sample collection between 1991 and 2015 at the stream gauge at the base of the watershed (Figure 1).
In addition to existing data, we collected samples from lysimeters, wells and the stream in the fall of 2017 for dissolved inorganic carbon (DIC) isotope analysis to test for carbonate weathering signals. Lysimeters and wells were pumped and allowed to recharge prior to sampling and all samples were capped immediately with no headspace, placed in a cooler, and shipped to the UC Davis Stable Isotope Facility within 24 h of collection.
To investigate landscape position controls on soil composition, we took soil samples every 1-3 m in three transects that spanned from upslope to midslope to near stream areas (Figure 1). Two transects cut across the two main tributaries, and one transect cut across a few meters below the confluence of the tributaries (Figure 1). These transect samples were taken using a bucket auger in 15-cm depth increments until 90 cm or depth of refusal.
To increase the statistical power of results and to allow for an in-depth analysis of compositional changes in soil with depth and landscape position, we additionally sampled each of the three soil orders from soil pits by taxonomic horizon in various locations throughout W-9 in 2017 (Figure 1). In order to allow for comparison with soils prior to recovery, we also analyzed archived samples (stored at room temperature in the dark in sealed glass vials) collected as part of a soil survey in 1996 (n = 156). We used average values of soil composition in archived (1996) and modern (2017) samples, some of which were re-sampled in the same approximate location (Figure 1).

Sample Processing and Analyses
Soil samples prepared for XRF and powder XRD were air dried (4-5 days), sieved (<2 mm), ground and homogenized in a ball mill, and stored dry prior to analysis. Soil elemental composition was determined on both archived samples (n = 156) and modern transect samples (n = 121) as well as samples from representative landscape positions of modern (n = 43) and archived soils (n = 78) using a Thermo Scientific Niton XL3t X-Ray Fluorescence Analyzer and a Thermo Scientific ARL QUANT'X X-Ray Fluorescence Spectrometer in the Geology Department at the University of Vermont and Middlebury College, respectively. From these data we calculated mass transfer coefficients (MTC) to determine the pedogenetic behavior of selected elements (Si, Al, Na, Ca; Brimhall and Dietrich, 1987;Brantley et al., 2007). For parent material we did not use bedrock or till, but the deepest sampled horizon (i.e., partially weathered till) because parent material composition in this area varies greatly across the watershed. MTC values therefore only show depletion vs. enrichment relative to the deepest accessible horizon. We used Ti as an immobile element because of its normal distribution in the parent material. All data from each soil order of the 1996 survey as well as the representative data from 2017 were averaged by horizon to assess the variability of each horizon of each soil order. Relative MTC was calculated as: where MTC (dimensionless) is the ratio of the concentration (C) of an element of interest (subscript M) relative to an immobile element (subscript I) in the weathered soil (subscript W) and the parent material (subscript P). Soil mineralogy was identified using X-ray diffraction (XRD) on random powder mounts from modern soil samples from each resampled pit. In order to identify secondary clay minerals, samples from the lowest horizon (till) were also analyzed using oriented texture preparations of the <2 µm fraction that was separated by decantation (Poppe et al., 2001). Suspensions were pipetted on glass slide preparations and allowed to air-dry. This approach allows for an orientation of crystallites that restricts diffraction in the crystallographic c-direction and hence only the (00l) intervals are visible (Lagaly, 1993;Moore and Reynolds, 1997). The clay fraction was analyzed air-dried (AD), treated with ethylene glycol (EG) to test for presence of swelling clay, and heat-treated (HT, 400 and 550 • C) to collapse the chlorite peaks to test for overlap with the kaolinite reflections (Moore and Reynolds, 1997). XRD analysis of all samples was conducted on a Rigaku Miniflex II Powder Diffractometer with CuKα radiation, operated at 30 kV and 15 mA. Scanning parameters were set at 0.02 step width and a count time of 10 s per step between 3 and 5-65 • 2θ.

Statistical Methods
All statistical analyses were completed in R (R Core Development Team, 2017). Stream flux was calculated to evaluate trends in recovery (sulfate, nitrate) and weathering (Si, Na, and Ca). All stream flux calculations were completed using the USGS LoadEst functions (Runkel et al., 2004) implemented within the LoadFlex package and trend analyses were conducted using the Kendall and trend packages (McLeod, 2011;Appling et al., 2015;R Core Development Team, 2017;Pohlert, 2018). LoadEst requires daily time steps, therefore concentration and discharge of highly sampled storm events were entered as daily averages. The best of the 9 predefined LoadEst models was then selected based on how accurately it predicted measured values and was applied to the discharge record for annual flux calculation. Flow adjusted residuals of the concentration-discharge relationship were then evaluated and, if they were not normally distributed, a Mann-Kendall trend analysis was conducted and the Sen slope was calculated. Baseflow composition was determined by calculating the mean stream composition during discharge conditions lower than the 20th percentile and storm flow composition was determined by calculating the mean stream composition during the of discharge conditions greater than the 80th percentile. Baseflow and stormflow composition were also computed at less than the 5th percentile and greater than 95th percentile.
The significance of trends in soil solution and GW composition with time and space were evaluated using several methods. For temporal trends in GW, linear regressions on raw data were used first as this method is powerful and if significance is established no more statistical analyses are needed to draw conclusions (Meals et al., 2011). Temporal trends in GW concentration were evaluated using data from the well with the most complete record (the well at the upslope). In the case of highly variable seasonal compositions, regressions on annual averages were implemented. When large parts of the datasets (>1/3 of the record) were missing or other tests failed to establish a trend, Welch's 2 sample t tests were used to evaluate any differences in composition (Meals et al., 2011;Derrick and White, 2016). Trends in space were evaluated by calculating the average and standard error at each depth and landscape position. If averages were within error of one another, Welch's 2 sample t test was used to determine whether the data sets were significantly different (Derrick and White, 2016).

Trends in Stream Water Composition and Flux
Discharge showed large variations (ranging from 7.4 * 10 4 m 3 yr −1 in 2001 to 1.5 * 10 5 m 3 yr −1 in 2011) and increased slightly but significantly (Figure 2A). Stream water pH values showed large, mostly seasonal, variations between 1991 and 2004 with values ranging from 6.7 up to 8.6 with less variability thereafter ( Figure 2B and Supplementary Materials). Annual stream water sulfate flux decreased more than 30% since 1991 (from 870 to 580 kg yr −1 , Figure 2C) and Mann Kendall trend analysis on flow adjusted model residuals indicated a significant negative trend in sulfate concentration (τ = −0.442, Sens Slope = −0.007, p < 0.05). The annual nitrate flux showed large interannual variability (between 60 and 105 kg NO 3 − per year) but no significant increase or decrease over time ( Figure 2D).
Both Si and Na stream water flux showed large interannual variability (ranging from 750 to 1400 kg yr −1 for Si and from 180 to 330 kg yr −1 for Na, Figures 2E,F). Mann Kendall trend analyses showed no significant trends in Si and Na concentration. The flux of Ca showed similarly large variations (ranging from 3450 and 6450 kg yr −1 ) and a slight increase that was significant (τ = 0.078, Sens Slope = 0.002, p < 0.05). The great interannual variability in fluxes is in agreement with the interannual variability in discharge (Figure 2A

Temporal and Spatial Trends in Groundwater and Soil Water Composition
We also investigated groundwater and soil water composition for trends with time, depth, and landscape position (upslope vs. midslope vs. base of the hill). Temporal trends are exemplified for selected years with typical annual precipitation (i.e., years 2004, 2008, and 2012). Significance was established using linear regression on the entire time series (see Supplementary Material Figure S1 for complete time series of GW composition).
Soil solution pH increased over time especially in deep and some shallow soil solution and was generally greater at depth and with proximity to the stream ( Figure 3A). GW pH in the upslope well increased from an annual average of 7.9 ± 0.3 in 2004 to an annual average of 8.1 ± 0.2 in 2012, which is higher than the average stream water pH at base flow (7.8 ± 0.2, Table 2).
Sulfate concentrations generally decreased with time but were greater at depth and with proximity to the stream ( Figure 3B). For example, the average sulfate concentration in all shallow soil water decreased from 3.2 ± 1.5 mg L −1 in 2004 to 1.6 ± 1.1 mg L −1 in 2008 and is in agreement with the general trend of decreasing sulfate over the years. GW sulfate decreased from an average of 5.3 ± 0.7 mg L −1 in 2004 to 4.4 ± 0.3 mg L −1 in 2012. The trend of increasing sulfate concentrations with depth was observed at the upslope position where the sulfate concentration in GW was up to 2.5 times higher than in soil water. The average stream water sulfate concentration at baseflow (<20%, 8.1 ± 2.3 mg L −1 ) was the highest mean sulfate concentration of any water source ( Table 2).
Nitrate concentrations did not change significantly over time and were lower at depth and closer to the stream ( Figure 3C). Annual average nitrate concentrations in shallow soil solution ranged from as low as 0.02 ± 0.01 mg L −1 at the base of the hill in 2004 to as high as 0.26 ± 0.12 mg L −1 at the hilltop position in   Figure 3C). Compared to upslope, nitrate concentrations at the base of the hill were slightly lower and similar to base flow stream water concentrations ( Table 2). Average annual soil and groundwater Si concentrations did not show significant changes with time but, like pH and sulfate, concentrations were generally greater with depth and proximity to the stream ( Figure 3D). For example, Si in GW at the upslope position remained relatively constant since 2004 at an average value of 12 ± 0.2 mg L −1 (Figure 3D). Soil water Si concentrations were generally greater (6.7 ± 0.8 mg L −1 ) at the base of the hill than upslope. In contrast, near stream GW Si concentrations (8.4 ± 1.2 mg L −1 ) were almost 30% lower than upslope GW concentrations, (12.0 ± 0.8 mg L −1 ) which is still almost 3 times higher than stream water base flow concentrations ( Table 2).
Annual average Na concentrations decreased with time in deep soil solution at the upslope and midslope positions and were generally lowest in upslope soil waters (as low as 0.6 ± 0.1 mg L −1 ) and largest in GW, which was similar to stream water composition at base flow ( Table 2).
Ca concentrations generally decreased through the study period but were greater at depth and proximal to the stream ( Figure 3F). Upslope GW Ca concentrations decreased significantly from 18.8 ± 1.3 mg L −1 in 2004 to 16.0 ± 0.8 mg L −1 in 2012 ( Figure 3F) and were generally greater at depth (1.0 ± 0.7 mg L −1 in shallow soil water vs. 16.4 ± 2.4 mg L −1 in GW). Stream water Ca concentrations at base flow were similar to GW concentrations ( Table 2).
The isotopic composition of DIC was determined for soil solution and GW. The mean δ 13 C-DIC of soil solution (−20.8 ± 1.2) was significantly lower than that of GW (−14.5 ± 1.5), and δ 13 C-DIC values also generally increased with proximity to the stream ( Table 3).

Spatial and Temporal Trends of Selected Elements in Soil
Spatial trends on transect samples generally showed depletions at upslope and midslope locations and accumulations in riparian areas (Figure 4 and Supplementary Figures S2, S3). The exception was Si that was, relative to the parent material, variably  Figure 4). Temporal trends in soil elemental composition were investigated by comparing MTC values from archived and modern samples (Figures 5-7). The variability in typical upslope soils (mostly Spodosols) was great, as indicated by the broad shaded areas (Figure 5), but generally showed a combined accumulation-depletion pattern (Figure 5). For example, MTC values for Si indicated mild enrichments (e.g., archived mean = 0.12 ± 0.36 in the O horizon) in the upper horizons (O, A, E) and depletions in lower horizons (e.g., archived mean = −0.28 ± 0.11 for the BC horizon). Aluminum also showed an addition-depletion profile where MTC values were negative at the surface and positive in the spodic horizon. Base cations were generally enriched and highly variable in the O horizon but depleted in underlying mineral soil (e.g., archived mean up to −0.58 ± 0.19 for Na). The differences between archived and modern samples were not significant for any of the investigated elements except Ca, which was generally more enriched in the top horizons of archived soils.
Midslope Inceptisol composition in the study area was less variable compared to Spodosols (shaded areas are less broad, Figure 6). Overall, Si was depleted throughout the profile (archived means between −0.04 ± 0.25 and −0.16 ± 0.17, Figure 6), while Al in archived samples was somewhat depleted in upper horizons and was depleted in lower horizons in modern samples. MTC values for base cations, especially Ca, were highly variable in the O horizon (e.g., archived mean = 3.65 ± 4.18 for Ca). As for Spodosols, differences between archived and modern samples was only significant for Ca in the top horizons, where Ca MTC values were significantly higher in archived soils.
The variability of riparian soil composition in the study area was high for most elements, especially in archived samples. For example, MTC values were generally positive for Si, negative for Al (especially in modern samples; Figure 7). Base cation MTC values were generally positive, especially for Ca, in the two uppermost organic horizons (32.1 ± 33.6 and 19.9 ± 29.7). Again, the differences between archived and modern samples was significant for Ca in the top horizons, where Ca MTC values were significantly higher in archived soils.

Mineralogy of Till and Soils
The deepest horizons (weathered till) contained primary minerals of the Waits River and Gile Mountain Formations including micas (with a prominent peak at 10.0 Å), chlorite (14.2/7.1/3.5 Å), pyroxenes (e.g., at 3.0 Å), feldspars (e.g., at 3.1 Å), quartz (at 3.34/4.26 Å), amphibole (at 8.5 Å), oxides (rutile or spinel at  2.5 Å), and apatite (2.8 Å, see Supplementary Materials for diffractograms). Representative soils for all landscape positions (Spodosol for upslopes, Inceptisols for most midslopes and Histosols for riparian zones) showed most peaks for the primary minerals associated with contributing bedrock formations, but several of these peaks changed or disappeared toward upper horizons (Supplementary Figures S4-S6). For example, quartz peaks decreased in intensity toward shallower horizons, whereas feldspar, muscovite, amphibole and pyroxene peaks varied but became difficult to resolve in the organic horizon of all soils. Primary chlorite is present in all soils (intense peak at 14.2 Å in powder samples). The secondary swelling clay mineral smectite (confirmed by EG treatment, Supplementary Figures S4-S6) and illite were found in deeper horizons of most soils. Smectite was most prominent in Spodosols and not well resolved in the riparian Histosol. In contrast, calcite was only identified in the lowest horizons of the riparian Histosol.

Streams as Indicators of Recovery From Acidification
Changes in stream water composition can indicate changes in CZ dynamics (Meybeck, 1993;Maher, 2011;McIntosh et al., 2017). However, a multitude of complex and coupled processes in the watershed and the stream itself can make the identification of specific processes or source locations difficult. Using the acid-impacted SRRW as a testbed, we investigated the potential and limitations of stream water composition as an indicator of watershed recovery from acidification and explored whether changes in soil composition contributed to a recovery signal. Generally, decreasing stream water nitrate and sulfate fluxes and overall increases in stream water pH are typical indicators of chemical recovery at the watershed scale (Newell and Skjelkvåle, 1997;Driscoll et al., 2001;Skjelkvåle et al., 2001Skjelkvåle et al., , 2005Rogora et al., 2013;Garmo et al., 2014). At SRRW, sulfate fluxes in the stream indeed decreased significantly ( Figure 2C) despite the fact that GW sulfate concentrations are high due to pyrite weathering at greater depth (Mayer et al., 2010). Over the past decades precipitation increased an average of 14 mm/yr (Shanley et al., 2016), which is reflected in the overall increase in discharge (Figure 2A). If sulfate still had a significant source in precipitation or release from soil stores, increases in annual precipitation may have led to stream sulfate flux increases, not decreases. The actual decreases in sulfate flux indicate that the stream is responding to the substantial reduction in sulfate deposition in the area (Shanley et al., 2004;Miles et al., 2012). In contrast, stream water nitrate flux did not change significantly, a finding that is consistent with the strong control of vegetation in a system that is not N saturated and the slower decrease in nitrate deposition (Burns et al., 2005). Stream water pH never dropped below 6.7, due to the high buffering capacity of the bicarbonate-rich waters (Shanley, 2000;Shanley et al., 2004)  but the magnitude of pH variation decreased, especially after 2004 ( Figure 2B).
Because the reduction of acid deposition reduces proton availability, we had hypothesized that stream water fluxes of Ca from base cation leaching would decrease significantly. Indeed, many recovering watersheds in the region show a substantial reduction in base cation fluxes including Ca and Na (Newell and Skjelkvåle, 1997;Stoddard et al., 1999;Burns et al., 2005;Garmo et al., 2014). However, SRRW Na (and Si) fluxes remained unchanged, whereas Ca fluxes increased (Figures 2E-G).
Because GW can indicate subsurface processes more directly than stream water (Peters and Aulenbach, 2009), we investigated temporal trends in GW Ca, Na and Si concentrations and found that the concentrations of Si remained stable, Na decreased slightly and that Ca decreased by over 30% since 2004 (Figures 3D-F). The relatively high δ 13 C-DIC values measured in GW are typical for a mixture of pedogenic DIC and carbonate weathering, confirming that carbonate is still actively weathering (Cerling, 1984;Amiotte-Suchet et al., 1999). Coupled with the decreases in GW Ca concentrations on the upslope and midslope, the isotopic patterns are consistent with either decreased carbonate weathering or increased plant Ca uptake during recovery, but they do not explain the increased Ca flux in stream water. We will discuss the potential role of Ca accumulation in the riparian zone as a mechanism to explain this pattern (section "The Riparian Zone: Integrator of Flow Paths and Soil Processes").

From Subsurface to Stream: Changes in Soil and Water Composition by Landscape Position
We had hypothesized that long-term lateral transfer of leached base cations would lead to depleted shallow horizons, especially at upslope or midslope locations, while the accumulation of base cations, especially Ca, would be largest in near stream locations. Indeed, base cation and Si concentrations in soil water generally increased with depth and also increased toward lowlying landscape positions, which might reflect a combination of vertical and lateral transfer through shallow flowpaths and supply from GW when water tables rise (Sommer, 2006). Highly variable Si MTC values at all landscape positions and in all soil horizons (Figure 4) might reflect a combination of unweathered quartz high in the profile as well as transfer of Si from the weathering of less stable aluminosilicates (Figure 4). Aluminum, at upslope locations, showed a typical Spodosol depletionenrichment pattern due to pH driven changes in Al solubility (Cronan and Schofield, 1979) but had significant enrichment in riparian soil horizons. XRD analyses indicate the weathering of aluminosilicate phases (especially amphiboles and pyroxenes) in most horizons and many of these minerals were still present in mineral soil horizons (Supplementary Figures S4-S6). Together with the high soil solution concentrations of Si (and partially for Na as well) these results indicate active silicate weathering in these soils.
Solid phase analyses of Na and Ca revealed negative MTC values for these elements in upslope and midslope soils but showed, again, substantial enrichment in riparian soils (Figures 5-7), a pattern that was also present in soils taken in transects associated with these landscape positions (Figure 4). E.g., Ca and Na were generally depleted in upslope-Spodosols (except for top horizons), relatively stable in the midslope-Inceptisols and enriched in riparian-Histosols, especially in top layers (Figure 4). Organic rich horizons tend to have less Ti, and hence the normalization to this otherwise immobile element to calculate MTC values can overestimate enrichments. However, riparian Histosol Ca content reaches up to 72000 mg kg −1 , compared to a mean of 26000 ± 13000 mg kg −1 for midslope soils. XRD analyses of these soils explain most of these patterns through weathering and disappearance of amphiboles, chlorite, apatite, and pyroxenes. The lack of calcite in the till at upslope and midslope locations indicates that the carbonate weathering front is much lower than the silicate weathering front (i.e., in the deep till). The lack of significant amounts of secondary Ca-bearing minerals in near-stream locations further suggests that the large amounts of Ca are not accommodated in crystalline phases.

The Riparian Zone: Integrator of Flow Paths and Soil Processes
At SRRW, the riparian zones are areas of confluence where shallow and deep flow paths converge, GW discharges and base cations accumulate, hence the riparian zone has a high potential for impacting stream water composition for these solutes.
Si and Na concentrations were similarly high in GW at upslope and midslope locations, however, concentrations were significantly lower in riparian zone waters and even lower in stream water [e.g., base flow stream water Si concentrations were almost 6 mg L −1 lower than riparian GW (Figure 3 and Table 2)], an observation consistent with previous research . The missing Si and Na might be stored in secondary minerals, such as smectites that were identified in weathered till (Supplementary Figures S4-S6). However, as mentioned above, we cannot explain the high riparian soil Ca content through only the accumulation of secondary minerals. Cincotta et al. (this special issues) found riparian soils at SRRW were rich in both Ca and organic carbon. Ca is common in biofilms, and is the dominant exchangeable cation in organic soils with high pH and cation exchange capacity. As a divalent cation, Ca can help bridge charges between negative soil particles (Perdrial et al., 2009). Organic matter readily forms complexes with free Ca 2+ ions (Ong and Bisque, 1968;Temminghoff et al., 1995), leading to DOM coagulation that, when Ca concentrations are high, is not pH dependent (Tipping and Ohnstad, 1984;Temminghoff et al., 1995). Therefore, Ca in riparian soils is likely associated with and stabilized by organic materials in soil aggregates. Cincotta et al. (this special issues) also showed that recovery from acidification can lead to the breakup of soil aggregates and the release of associated organic matter into soil solution. Their study did not investigate Ca release, but if Ca is associated with organic matter it would be liberated in the same process. Indeed, annual discharge has slightly increased since 1992 (Figure 2A) and so has stream water DOC flux (Schuster et al., 2007;Sebestyen et al., 2008;Cincotta et al., this special issues). Thus, the slight but significant increase in stream water Ca flux ( Figure 2G) could have the same origin as the stream DOC increase, i.e., release from soil aggregates that served as transient store. Both indicate riparian processes driven by changes in aggregate stability and solubility (not increased carbonate weathering), with the result that Ca is effectively flushed into the stream. This effect is exacerbated by increased discharge (Figure 2A) and increased abundance of high intensity events (Peters et al., 2006), advectively transferring soil-derived material into the stream.
Furthermore, our results from modern vs. archived soil samples are in agreement with a gradual, recovery-induced, leaching of legacy Ca from riparian soils. The Ca concentration of the riparian soils decreased measurably between archived (1996) and modern samples (2017) throughout the study area as well (Figure 7). We illustrate the strong impact of riparian soils and near stream waters vs. midslopes on stream water in a mixing diagram of Na-normalized mole ratios (Figure 8). Here, composition of upslope and midslope soils have overall lower Ca/Na ratios than the riparian soils, which again have lower ratios than GW or stream water. These results confirm that riparian zones are highly reactive parts of the CZ and exert a strong control on stream composition (Peters and Aulenbach, 2009;Vidon et al., 2010;Lidman et al., 2017). In the case of SRRW, the complex interactions among base cation leaching, converging flow paths, GW upwelling (providing changes in pH and further accumulation of more reactants), temporal solute storage (i.e., transient legacy), and interaction with organics, produce a highly reactive zone that responds to the same driver (decrease of acid inputs) in fundamentally different ways than the midslopes.

AUTHOR CONTRIBUTIONS
JA was responsible for sample collection, preparation, and analyses as well as figure production, table production, and primary authorship. JP as MS thesis advisor aided in interpretation, authorship, field work, lab work, lab training, as well as allowed lab use. AG and JE assisted with field sampling and lab work. NP contributed through interpretations, authorship, and lab training. MC assisted with field sampling, lab work, and interpretations through research published in this issue. DR assisted in interpretations and editing. JS provided soil solution, groundwater, and stream composition data as well as archived soil samples and also introduced authors to watershed installations, aided in interpretations and provided pre-research consulting. KU aided in statistical processing and interpretations as well as provided example R scripts. PR allowed the use of and provided training on the XRF at Middlebury College.

FUNDING
Funding supporting this research was made available through the Geology Department at the University of Vermont.