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). Decades-long 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 the till to as much as 100 mm/hr in some surficial soils (Kendall et al., 1999). 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).

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

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

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

Discharge showed large variations (ranging from 7.4^∗10⁴ m³ yr^-1 in 2001 to 1.5^∗10⁵ m³ 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₃^- 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, Figure 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), showing very wet years in 1996, 2008 and 2011 and much drier years in 1995, 2001, and 2012.

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 2008 (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 δ¹³C-DIC of soil solution (-20.8 ‰ ± 1.2) was significantly lower than that of GW (-14.5 ‰ ± 1.5), and δ¹³C-DIC values also generally increased with proximity to the stream (Table 3).

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 enriched or depleted in upslope and midslope positions (MTC values ranging from -0.51 to 0.64). Al was mildly depleted from upslope soils (-0.10 to -0.39) and accumulated in the top 45 cm of the profile closest to the stream. The base cations Na and Ca were typically depleted from the upslope soils (e.g., MTC values ranging from -0.86 to -0.08 for Ca) and accumulated in near stream soils, especially in the top 30 cm (up to 31.2 for Ca, Figure 4).

Temporal trends in soil elemental composition were investigated by comparing MTC values from archived and modern samples (Figure 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.

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.

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., 2001, 2005; Rogora 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 (Figure 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 (Figure 3D–F). The relatively high δ¹³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”).

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 low-lying 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 depletion-enrichment 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 (Figure 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.

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 (Kendall et al., 1999). 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²⁺ 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.