Using Concentration–Discharge Relationships to Identify Influences on Surface and Subsurface Water Chemistry along a Watershed Urbanization Gradient

: Urban development within watersheds impacts the hydrology and water quality of streams, but changes to groundwater–surface water interactions in this “urban stream syndrome” are not yet well understood. This study focused on three stream systems in a northern Virginia (USA) protected area with 14.2, 31.7, and 66.1% developed land in their watersheds. Surface water was sampled weekly for nutrients, dissolved metals, sulfate, ancillary water quality parameters, and discharge over two non-consecutive years with the hyporheic zone sampled during the second year. Concentration–discharge relationships revealed largely chemostatic behavior in surface water solutes in the least urbanized stream, while in the two more urbanized streams, these relationships tended to have significant positive and negative slopes, indicating diverse delivery pathways depending on the constituent. In the least urbanized stream, linear regressions between discharge and solute concentrations in hyporheic water had exclusively negative slopes, indicating source-limited delivery, while the other two urbanized streams maintained largely chemostatic behavior. Average specific conductance and nitrate + nitrite concentrations in stream surface water reflected an urbanization gradient, while sulfate, Ca, K and Sr concentrations suggested a threshold effect: the stream with a mostly forested watershed had the lowest concentrations, while the other two were higher and similar. Specific conductance indicated salinization of both surface and groundwater at the two more urban streams, possibly threatening aquatic organisms. Metal concentrations in surface and subsurface water were often positively correlated with specific conductance and negatively correlated with pH, suggesting that they may originate from road salt and/or be mobilized by acid precipitation. These results indicate the importance of monitoring both baseflow and stormflow as pathways for pollution.


Introduction
Freshwater resources are threatened globally due to increased demand for potable water, widespread pollution, species introductions, and land use change [1].Nowhere is that effect more pronounced than in urban areas, where most people now live [2,3].The collective changes to urban streams have been termed "urban stream syndrome" (USS).Statistically significant differences from forested streams have been documented after only a 3% increase in impervious surface coverage [4].Changes associated with USS generally fall into three main categories: hydrology, channel morphology, and ecosystem function [5][6][7].Water 2021, 13, 662.https://doi.org/10.3390/w13050662 Academic Editor: Thomas Boving Urban streams have increased peak discharges following precipitation events ("flashiness") and shorter lag times between the start of a precipitation event and peak discharge [8,9].Urban stream channels are more highly eroded and straightened, and their headwaters have often been drained and filled or replaced with stormwater infrastructure, reducing drainage density [10][11][12].Urban streams also frequently have much higher levels of contaminants, such as nutrients and dissolved metals, from fertilizer, industrial runoff, road deicers, and carbonate weathering from impervious surfaces [13][14][15].These changes can collectively alter stream metabolism and disrupt organisms' life histories [16][17][18].
While USS has been well documented in the literature, urbanization's effects on groundwater and groundwater-surface water interactions are not yet well understood [7,19].Groundwater-surface water interactions primarily occur in the hyporheic zone, the near-subsurface region surrounding a stream.The hyporheic zone is saturated with a mix of groundwater and surface water.This zone is critical for denitrification, and it also sustains populations of aquatic organisms by limiting temperature variability fluctuations and providing shelter [20][21][22].Straightening stream channels can disconnect urban streams from their historic floodplains and hyporheic zones by lowering the water table, influencing pollutant transport and transformations [10,[23][24][25][26].As the hyporheic zone is highly dynamic, more studies capturing temporal variability at a variety of scales are needed to elucidate changes associated with urbanization.
This study used concentration-discharge (C-Q) plots to investigate hydrology and water quality in streams and their hyporheic zones.The concentration of a constituent in a stream depends on the relative amounts of baseflow and stormflow (or snowmelt) present in a stream at a given time and on the concentrations of the constituent in baseflow and stormflow.An "up" C-Q plot with a positive slope indicates higher concentrations associated with higher discharge, or transport-limited behavior [27].In contrast, a negative slope or "down" C-Q plot indicates source-limited behavior: concentrations are higher in baseflow than in stormflow, and dilution occurs during storm events.A slope of −1 corresponds to a relationship completely driven by dilution.A "flat" C-Q plot with a slope that is not significantly different from zero indicates chemostatic behavior, where the concentration is relatively constant regardless of discharge [27,28].Chemostatic behavior can occur if the constituent in distributed homogeneously throughout the watershed, leading to similar concentrations in baseflow and stormflow [27].
In some cases, C-Q plots show different trends at high vs. low flows, and it can be informative to segment them and look at different combinations of "up", "down", and "flat" patterns [27].C-Q plots may be segmented at the median discharge [27], or at another intermediate point determined using Bayesian inference [29].Alternatively, when hundreds or more C-Q data points are available from long-term monitoring, weighted regressions on time, discharge, and season (WRTDS) [30] can be used to deal more effectively with nonlinear relationships and change over time.While C-Q relationships are an essential tool for comparing pollutant transport relationships across streams [28,31,32], this technique has less often been utilized across an urbanization gradient.
The overarching goal of this study was to characterize changes in solute sources and transport in streams across an urbanization gradient ranging from almost completely natural/vegetated to 64% developed land cover.Specifically, we investigated the following research questions: (1) How do concentrations of nutrients, metals, and other water quality parameters differ in the surface and hyporheic water of three streams with varying degrees of watershed urbanization?(2) What patterns of temporal variability in discharge and water quality occur in each stream?(3) How do C-Q relationships vary among different water quality parameters and among streams with different levels of watershed urbanization?(4) How are water quality parameters correlated in surface and subsurface water, and how do these correlations differ among streams with different levels of watershed urbanization?

Methods
This study took place in Meadowood Special Recreation Management Area (Meadowood) in Lorton, Virginia, approximately 32 km southwest of Washington, DC.Meadowood contains three streams: Giles Run (GR), South Branch (SB), and Thompson Creek (TC) with watershed areas from 1.4 to 13.1 km 2 and a gradient of watershed urbanization (Figure 1, Table 1).GR is a third-order stream while SB and TC are second-order streams.These streams sit atop the Northern Atlantic Coastal Plain aquifer system, a semi-consolidated sand aquifer.The closest U.S. Geological Survey sites (ID numbers: 383830077135502 and 383423077245901) that measure depth to the water table were 4.7 and 21.7 km away from Meadowood and indicate an average depth between 60 and 180 cm from the land surface.Rainfall is evenly distributed throughout the year with an average of 1000 mm annually, while snowfall primarily occurs during the winter months and averages around 400 mm (U.S. Climate Data).

Giles Run (GR) South Branch (SB) Thompson Creek (TC)
Watershed Area (km  .Samples were collected at the most downstream point of each stream within Meadowood (Weekly sites in Figure 1).The hyporheic zone was sampled from the streambank or streambed starting in the second year (September 2016-July 2017; n = 32 per site) using a 60 cm long PushPoint sampler (MHE Products, East Tawas, MI, USA) and syringe.The vertical hydraulic gradient was assessed by using the tube attached to the PushPoint as a manometer and visually inspecting whether the hyporheic water level was above or below the level of surface water in the stream.Spatial surveys of surface and subsurface water were completed in late October 2016 and mid-March 2017 by sampling 4-6 additional sites on the main branch and tributaries within each stream system (Spatial sites in Figure 1).
Stream discharge was quantified using the standard velocity-area method [34].We attempted to use automated level loggers (Onset HOBO U20L) combined with rating curves to obtain higher-precision discharge data, but the loggers were frequently moved by the current and/or buried with sediment, which resulted in poor data quality.Precipitation was measured at a HOBO weather station (Onset, Bourne, MA, U.S.A) installed within the park at 15-30 min intervals that were later aggregated into daily measurements.An estimated runoff ratio (m 3 s −1 mm −1 ) for each stream was calculated as the average of discharge (m 3 s −1 ) divided by 24-hour antecedent precipitation (mm) for each sampling date when rainfall occurred.
Geographic analysis of the watershed upstream of each sampling site was conducted using the USGS StreamStats tool [33] and ArcGIS v10.5.1 (ESRI, Redlands, CA).Statistical analysis was completed in RStudio version 4.0.2(RStudio Team 2020).
Data analysis corresponding to each of the main research questions was conducted as follows.We tested for significant differences in (a) surface water in the three different streams, (b) subsurface water in the three different streams, and (c) surface and subsurface water within the same stream by first calculating the difference between the two data series at each point in time, excluding any dates missing data from one or both time series.We tested for lag-one serial correlation using the Yule-Walker estimate; the null hypothesis of no serial correlation was rejected when p < 0.05.If there was no significant lag-one serial correlation, the regular calculation of standard error (SE) was used.If there was significant lag-one serial correlation, the standard error was adjusted (SEadj) using the first serial correlation coefficient r: If the 95% confidence interval calculated using the appropriate SE did not overlap with zero, we determined that the difference in means between the time series was significant.
We explored patterns of temporal variability by graphing water quality data from all three streams against time (Figures 2 and 3).This allowed us to see whether the patterns of temporal variability were similar or different among the three streams.We also binned data from each stream by season, based on the dates of the equinoxes and solstices during the study period, and tested for significant seasonal differences in surface and subsurface water quality within each stream using Kruskal-Wallis tests; the null hypothesis of no difference between groups was rejected when p < 0.05.
Similar to other papers assessing C-Q relationships [28,31], we created log-log plots of water quality parameters vs. discharge within each stream.This was done for both surface and subsurface concentrations.We hypothesized that variation in stream discharge might affect concentrations of water pollutants not only in surface water, but also in subsurface water, because stream discharge might be correlated with changes in the direction or magnitude of flow in the hyporheic zone.The magnitude and statistical significance of C-Q relationships were determined by calculating the line of best fit for each log-log plot.We assessed whether the slope of each line of best fit was significantly different from zero using the ANOVA F-test; the null hypothesis that it was not significantly different from zero was rejected when p < 0.05.
Correlations among surface and subsurface water quality variables were assessed using the Spearman correlation coefficient as most variables were non-normally distributed.Normality was tested using the Shapiro-Wilk test; the null hypothesis that the data are normally distributed was rejected when p < 0.05.

Water Quality Variation among Stream Systems
When comparing mean surface water quality values among the three stream sites, several distinct patterns emerged.Specific conductance and Ca followed the urbanization gradient, with the highest mean value at GR, intermediate at SB, and lowest at TC. N+N and Sr showed an urbanization threshold effect, with GR and SB not significantly different from each other, but both higher than TC.SO4 2− , Fe, K, Mg and Mn had significant differences that did not bear an obvious relationship to urbanization (Figure 2, Table 2).
Concentrations of SO4 2− and N+N in subsurface water displayed an urbanization threshold effect where GR and SB were not significantly different from each other, but both were higher than TC.NH4 + , Fe and Mn were significantly higher in TC subsurface water, possibly reflecting differences in geology or redox state.Specific conductance in GR subsurface water was higher than at the other two streams, suggesting urbanizationrelated salinization.However, TC subsurface water was intermediate between GR and SB, which, again, may reflect different subsurface geology in the TC watershed (Figure 3, Table 2).Subsurface chemistry differed significantly from that of surface water across all three stream systems for several variables.Surface water N+N and SO4 2− concentrations were significantly higher than subsurface concentrations across all three stream systems.Conversely, subsurface NH4 + and PO4 3− concentrations were significantly higher than surface concentrations across all three stream systems.Specific conductance was significantly higher in SB surface water compared to subsurface water, although the magnitude of the difference was small.In contrast, subsurface specific conductance was almost five times higher than in surface water, and the difference was significant, in TC.Subsurface specific conductance was also significantly higher than that of surface water for GR, though the magnitude of the difference was not as high.Mean concentrations of dissolved metals were generally higher in subsurface than surface water (Al in SB, Ba, Fe, and Mn in all streams, Ca, Mg, and Sr in TC only).However, Ca and Mg in SB were higher in surface water than in subsurface water (Table 2).We also assessed spatial variability in water quality within each stream system by sampling surface and subsurface water at 4-6 points during fall and spring spatial surveys (Figure 1, Spatial sampling points).Generally, the range of temporal variability at each stream system's weekly sampling site was larger than the range of spatial variability within each stream system, although there were some exceptions to this pattern.This indicates that the differences in water quality among weekly sampling sites reflect differences among the three stream systems, rather than resulting from uncharacterized spatial variability within each system.In other words, water quality differences between GR, SB and TC would be apparent regardless of the specific weekly sampling point chosen.

Variation in Discharge and Water Quality at Multiple Temporal Scales
No significant differences in discharge or water quality parameters between the two sampling years (June 2014-July 2015 and May 2016-July 2017) were observed in any stream system.However, significant seasonal differences were observed (Figures 2-3, Table S1).In surface water, specific conductance was significantly higher in winter versus other seasons in all streams.Discharge was significantly higher in all streams during winter and spring.Surface water NH4 + concentrations showed significant seasonal variability within all three streams, but with variable patterns.During the fall, GR had significantly lower concentrations when compared to all other seasons, while SB and TC had significantly higher concentrations during the summer versus all other seasons (Figure 2, Table S1).Surface water N+N, PO4 3− , and SO4 2− concentrations were not significantly different by season across any of the streams (Figure 2, Table S1).Substantial seasonal variation in the concentrations of dissolved metals in surface water was observed.The strongest pattern to emerge was that the concentrations of many metals were highest either in winter (Ba, Ca, Fe, Mn, and Sr at GR and Ba at TC) or in fall and winter (Mg at GR; Ba, Ca, Mg and Sr at SB; Ca, Fe, Mg, and Sr at TC).Additionally, Al at SB was higher in the summer than in the winter or fall and Mn at SB was higher in the summer than in other seasons (Figure 2, Table S1).
Subsurface variables also differed significantly by season within each of the three stream systems.Subsurface specific conductance in GR was significantly higher in the winter months versus the other seasons.Subsurface PO4 3− concentrations were highest in the winter in GR and SB.Subsurface SO4 2− concentrations were significantly higher in fall versus spring for both GR and SB, and they were highest in fall versus all other seasons in TC.Subsurface N+N and NH4 + concentrations did not differ significantly by season in any of the streams.Subsurface dissolved metal concentrations showed some significant seasonal variability but were not consistently higher in winter as surface water metal concentrations were (Figure 3, Table S2).

Concentration-Discharge (C-Q) Relationships
Our data captured baseflow and stormflow conditions in all three streams and showed that both have the potential to deliver dissolved constituents, supporting the use of C-Q relationships to distinguish their contributions.Baseflow can be a source of nutrients, metals, and other constituents if concentrations are higher in subsurface water compared to surface water and if the hydraulic gradient favors groundwater influx into the stream, rather than groundwater recharge from the stream.Assessment of the vertical hydraulic gradient between the hyporheic zone and the stream suggested that the direction of water flow was generally from the hyporheic zone into the stream (70%, 57% and 82% of sampling events at GR, SB and TC, respectively).
The amount of stormflow present in each of the three streams during weekly sampling events was highly variable and, in some cases, substantial.Discharge measured at each stream spanned 2-3 orders of magnitude, which is comparable to other studies incorporating C-Q relationships [27,28,31].For each stream, a significant percentage (61, 30, and 24%, respectively, for GR, SB, and TC) of the 88 sampling events exceeded the highest low-flow value given by StreamStats [33], and up to 25 mm of rain occurred in the 24 hours before weekly sampling events.No discharge measurements for any of the 3 streams approached the corresponding 5-year flood level, but this is not surprising given that we collected weekly measurements for 2 years.This indicates that our results are applicable to conditions ranging from baseflow to moderate rainfall events, but not to heavy storms of the type that would occur only a few times a year or less frequently.
C-Q plots with surface water concentrations showed significantly positive slopes for Fe only at GR and for Al at all three streams, indicating higher concentrations associated with stormflow conditions.C-Q plots of Ca and K (for GR only) were significantly negative, indicating source-limited behavior.Fe and Mn were significantly negative as well but at SB only (Figure 4, Table S3).The absolute value of negative slopes for Ca was ≤0.1, indicating a large chemostatic component even though some dilution was apparent.The other negative slopes ranged from −0.10 through −0.18, while the positive slopes ranged from 0.17 to 0.42 (Table S3).Other plots were consistent with entirely chemostatic behavior (Figure 4, Table S3).None of the plots showed a visual indication of a different C-Q relationship at high and low flows (Figure S1).
C-Q plots incorporating concentrations in subsurface water were just as likely to have slopes significantly different from zero as those incorporating surface water concentrations, supporting the idea that, at least in some cases, stream discharge affects hyporheic water quality with a short or negligible time lag.Patterns observed in subsurface C-Q plots reflected the urbanization gradient.C-Q plots for specific conductance, SO4 2− , K, Mg, Mn and Sr had significantly negative slopes only at TC, indicating that recent precipitation dilutes the higher concentrations generally present in hyporheic water (Figure 5, Table S4).This was further supported for specific conductance, K, Mg, Mn, and Sr because concentrations were generally higher in subsurface than surface water (Table 2).While specific conductance and K had slopes with an absolute value ≤0.10, others ranged from 0.13 to 0.36 (Table S4).Ammonium was the only constituent to have a significantly positive slope at TC, indicating higher concentrations associated with stormflow conditions.Additionally, both Fe and Mn had significantly positive slopes in SB only, indicating higher subsurface concentrations associated with stormflow conditions.All constituents in GR, as well as constituents in the other two streams not just highlighted, had slopes not significantly different from zero, indicating no significant relationship between stream discharge and hyporheic concentration (Figure 5, Table S4).

Correlations among Water Quality Variables
Correlations involving nutrients (N+N, NH4 + , and PO4 3− ) and SO4 2− in surface and subsurface water were generally weaker and less significant than those involving metals (Table 3, Table 4).SO4 2− and specific conductance had significant positive correlations in surface water across all streams, although the correlation was weaker in the more urbanized GR.Significant correlations that existed only in GR or TC may reflect the influences of urbanization and forest cover, respectively.
Surface water correlations-In surface water, Ba, Ca, K, Mg and Sr tended to have significant positive correlations with each other across all three streams, suggesting common sources and transport mechanisms.In some cases (e.g., Ca-Mg and Ca-Sr), correlations approached 1 (Table S6).Fe and Mn had a significant positive correlation to each other in all three streams, but their correlations with other metals were generally weaker and less significant.Al had significant negative correlations with Ca, Mg and Sr in all streams and with Ba, Fe and Mn in some streams, indicating different controls on its behavior.Not surprisingly, many metals had significant, negative correlations with pH, likely due to enhanced solubility and mobility under acidic conditions (Table 4).
Table 3. Spearman correlation information for GR (blue), SB (red), and TC (green) for temperature (Temp), dissolved oxygen (DO), pH, specific conductance (SC), nutrients and sulfate.Measurements from both surface water (SW) and subsurface water, collected from the hyporheic zone (HZ), are included.All correlations shown as positive (+) or negative (−) are statistically significant (p < 0.05).Non-significant correlations are shown as "0".See Table S5 for r and p values.Subsurface water correlations-Many similar correlations were observed in subsurface or hyporheic zone water.Ba, Ca, Fe, K, Mg, Mn and Sr tended to have significant positive correlations with each other across most or all streams.In contrast to surface water, Al had a significant positive correlation to Ba in subsurface water.Its correlations with other parameters tended to be weak and insignificant.Some correlations in subsurface water appeared to reflect the urbanization gradient.Significant negative correlations between pH and the concentrations of various metals (Ca, Fe, Mg, Mn, Sr) were observed only in GR, the most urbanized stream.Fe had significant positive correlations with other metals (Ca, Mg and Mn) only in TC, the least urbanized stream.This may indicate a natural soil source of these metals in TC, and an anthropogenic source that does not also contribute Fe in more urbanized watersheds (Table 4).

Variable
Surface versus subsurface correlations-Significant correlations between surface and subsurface water quality within the same stream were strongest and most numerous in GR, the most urbanized stream.N+N, NH4 + and specific conductance in surface and subsurface water were only correlated in GR.This may be indicative of fewer biogeochemical transformations as water moves through the hyporheic zone in GR.Surface and subsurface PO4 3-concentrations had a significant positive correlation across all streams.Dissolved metal concentrations in surface and subsurface water were positively correlated in both GR and TC, although the correlations tended to be stronger and more significant in GR.Few significant correlations between surface and subsurface concentrations were observed at SB.The reason for this is not clear, but it may reflect less connectivity between surface and subsurface water at that site (Table 5).

Solute Sources and Transport across an Urbanization Gradient
C-Q relationships support the idea that urban development in the GR watershed has changed the sources and transport of some dissolved constituents.While Fe in surface water had a significant positive C-Q relationship, indicating transport-limited behavior, both Ca and K had significant negative C-Q relationships, indicating source-limited behavior (Figure 4).These significant relationships occurred only in the most urbanized stream, GR, or the two most urbanized streams, GR and SB.Al had a significantly positive C-Q relationship not only in GR but in the other two streams as well.Al is mobilized by acid rain [35], which occurs in the Middle Atlantic region of the United States where Meadowood is located.Thus, the positive C-Q relationships for Al could be explained if acidic precipitation both increases stream discharge and mobilizes Al.
Conversely, when examining the relationship between subsurface concentrations and surface water discharge, most significant trends found were in the least urbanized stream, TC, where significant slopes were largely negative (Figure 5).This indicates greater infiltration of precipitation and subsequent dilution of subsurface constituent concentrations, while the observed chemostatic behavior in the other two urbanized streams likely indicates a reduction in groundwater infiltration of precipitation [12,25].It is possible that discharge affects subsurface solute concentrations with a time lag; however, since we measured discharge weekly, we cannot investigate this possibility with our data set.Moreover, because subsurface water was sampled <60 cm below the land surface either in or immediately adjacent to the stream, lags would be expected to be short.
Our results suggest an effect of watershed urbanization on N transformations.The hyporheic zone is a critical area of pollutant attenuation and, particularly, N transformation.During periods of lower DO availability, denitrification is the dominant process, whereas when DO concentrations are high, nitrification dominates [22,36].Strong significant and opposite correlations were observed in GR between N+N and DO (positive) and between NH4 + and DO (negative), suggesting both nitrification and denitrification occur, but a significant fraction of N+N remains despite these processes (Tables 2 and 3).
Incomplete consumption of N+N is common.First, denitrification processes in streams may not always proceed to completion, leaving some nitrate remaining the water [37].While this has been primarily demonstrated in streams surrounded by agricultural land [38,39], urbanization within a stream's watershed can also influence how much net nitrate is exported [37,40].Second, there could be other factors limiting denitrification, such as dissolved organic carbon, which acts as a power source for the redox reaction [26,41,42].Finally, the hyporheic zone may be connected to the stream channel only intermittently, as has been demonstrated in other urban studies, where lower groundwater tables prevented mixing with higher-N+N surface water [24,25].A disconnection between the hyporheic zone and surface water seems less likely due to numerous significant correlations between surface and subsurface water in GR, including with N+N and NH4 + .
This study revealed that several constituents other than nitrogen species were higher in subsurface compared to surface water (Table 2).Many of these solutes (PO4 3− , Ba, Fe, and Mn) were higher in subsurface water across all three streams, indicating that they may enter groundwater through rock-water interactions and be transported via groundwater flow.Their source and transport would therefore be unaffected by urbanization.Conversely, Ca, Mg, and Sr were only higher in subsurface water at TC.
Recent research has focused on biogeochemical coupling in urban streams, where concentrations of nutrients and metals are highly correlated with one another due to land use change in the surrounding watershed [2,[43][44][45][46].This coupling occurs when elements come from the same pollutant source (e.g., wastewater) and/or hydrologic "end-member" (e.g., surface runoff or subsurface groundwater) [47,48].Most correlations between surface water variables other than dissolved metals in this study indicated an urbanization threshold effect where correlations were only significant in either the most urbanized stream, GR, or least urbanized stream, TC.The strongest relationships were among pH, specific conductance, PO4 3− , N+N, and NH4 + (Table 3).PO4 3− and N+N could be transported together from wastewater; the lack of seasonal differences in N+N and PO4 3− in GR and SB may be indicative of wastewater inputs, which tend to be constant throughout the year [18].
Correlations between water quality parameters in the hyporheic zone differed from those in surface water (Table 3 Table 4).The strongest relationship between subsurface solutes was between SO4 2− and PO4 3− in GR, suggesting a common source.The correlation between SO4 2− and PO4 3− was far lower and not statistically significant in both SB and TC (Table 3).The source of SO4 2− and PO4 3− to GR could be leaky sewage pipes, septic tanks or Rainwater Landfill, located in the GR watershed.As SO4 2− could also be dissolved from sediments, a natural source could be possible for both SO4 2− and PO4 3− , with biological uptake in surface water reducing the correlation between the two in surface water [14,49].
For dissolved metals, correlations in surface and subsurface water were observed across all three stream systems and were not stronger or more prevalent at GR (Table 4).This may be because the metals in this study have similar sources and transport mechanisms within each watershed, even though they may vary among the three watersheds.However, significant, negative correlations among pH and several metals occurred in GR surface and subsurface water but generally not at the other sites (Table 4).This suggests that urbanization can lead to the coupling of water quality parameters that, in less developed watersheds, are not correlated.

Stream Health and Urban Stream Syndrome
Flashy hydrology-Our data suggest that hydrology at GR, the most urbanized stream, is flashier.GR's estimated runoff ratio (0.19 m 3 s −1 mm −1 ) was similar to that of other urban streams in Maryland [50] and much higher than those of SB and TC, which both had ratios less than 0.05 m 3 s −1 mm −1 (Table 1).These results are consistent with the findings of Arnold and Gibbons [51], who found that increases of 10-20% in impervious surface coverage (ISC) associated with urbanization can cause runoff to double relative to natural catchments.
Specific conductance-Temporal trends in specific conductance (Figure 2, Figure 3) indicate instances of very high specific conductance occurring sporadically throughout the winter.These events, which likely correspond with the timing of road salt application and/or snowmelt, are most pronounced in GR and are also apparent in the other two streams.Elevated specific conductance, indicating higher concentrations of total dissolved solids (TDS), is common in northern U.S. streams affected by urban stream syndrome due to road salt application during the winter and spring seasons [14,15,52,53].
Specific conductance for both surface and subsurface water was highest in GR relative to the other streams, suggesting the hyporheic zone and groundwater pools are becoming salinized, as has occurred in other urban streams [54,55].Surprisingly, TC had relatively high specific conductance values in subsurface water as well, but these values may result from factors other than urbanization and road salt application.Two lines of evidence support this interpretation.First, the lack of elevated surface water specific conductance suggests the ions dissolved in Thompson Creek's hyporheic zone are bioavailable, unlike chloride (Cl − ), which is the primary component in road salts.Naturally occurring ions in Mid-Atlantic streams include HCO3 − , SO4 2− , Ca, Mg, and K [14,49].SO4 2− , Ca, and Mg were measured in this study and showed positive correlations with subsurface specific conductance in TC though K, also measured, did not (Table 4).Second, subsurface specific conductance concentrations in TC are consistent year-round with no significant seasonal differences.In comparison, subsurface specific conductance concentrations in GR and SB were higher during the winter and spring (Figure 3, Table S2) during times of peak road salt application and were significantly different than concentrations in the summer and fall.
Excess salts in GR and SB surface water may pose a threat to Meadowood's aquatic ecosystems.Specific conductance in GR exceeded the EPA recommended chronic benchmark to preserve aquatic life in the Mid-Atlantic region (300 µS cm −1 ; U.S. EPA 2011) 80% of the time during this two-year study, while SB exceeded the benchmark 16% of the time and TC never exceeded it.In addition to directly increasing mortality of aquatic organisms, salinization can cause changes in growth rate, reproductive potential, and other life history traits, depending on the type of road salt applied and the duration of exposure [56][57][58].The significant correlations between Ba, Ca, Mn and Sr and specific conductance suggest that snowmelt is an important delivery pathway for those metals as well.
Nutrient and sulfate concentrations-Dissolved inorganic nitrogen (DIN) concentrations in all three streams were above the EPA chronic benchmark of 0.2 mg/L for total dissolved nitrogen [59] during most weekly sampling events (87, 83 and 67 out of 88 in GR, SB and TC, respectively).Considering that we only measured inorganic nitrogen and the criterion is for the sum of dissolved organic and inorganic nitrogen, the streams likely exceed the benchmark almost all the time.The 30 µg/L chronic benchmark for PO4 3− [59] was exceeded during some weekly sampling events (14, 28 and 17 out of 88 in GR, SB, and TC respectively), while the 51.6 mg/L chronic benchmark for SO4 2− [59] was never exceeded.These results indicate that excess loading of nutrients, especially N, could harm aquatic organisms and ecosystems in Meadowood.
In general, nutrient and SO4 2− concentrations observed in this study fell within the ranges observed in other urban streams in the eastern US [13,15,48,60,61].Flow-weighted average concentrations of N+N and SO4 2− followed the urbanization gradient, while PO4 3− and NH4 + had overlapping 95% CI for all streams.
Dissolved metals-Of the metals measured in this study, Al, Ba, Fe, Mn and Sr can be toxic to aquatic life [59,[62][63][64][65][66] while Ca, Mg and K are generally not considered toxic.For Al, no samples exceeded the EPA acute aquatic life limit, but surface and subsurface water samples exceeded the chronic aquatic life limit 7% and 13% of the time, respectively.These exceedances occurred across all three streams.There is no EPA aquatic life limit for Ba, but all samples fell well below the range that can be toxic to duckweed [62].The EPA only sets a chronic aquatic life limit for Fe, which had been 1 mgL −1 but was recently revised to be 0.25-0.5 mgL −1 for ferric iron [63].We did not determine the oxidation state of iron in these samples, but even using the 1 mgL −1 criterion, 20% of all surface water samples exceeded it, indicating that Fe concentrations in these streams could have harmful effects on aquatic life.Exceedances occurred mostly in TC, occasionally in GR and never in SB.The EPA does not set aquatic life limits for Sr, but a recent review [66] recommended 75 mgL −1 as the acute aquatic life limit and 11 mgL −1 as the chronic aquatic life limit.The highest Sr concentrations measured in this study were ~3 mgL −1 , indicating that Sr toxicity is unlikely to occur.
The term "urban stream syndrome" masks a complex array of effects experienced by streams as their watersheds urbanize.In this study, some data (e.g., specific conductance and Ca in surface water) showed a clear gradient reflecting watershed urbanization.Other data pointed to a threshold effect, where SB, the stream with an intermediate level of watershed urbanization, grouped either with the urbanized GR or the forested TC.For example, average concentrations of N+N and some metals for GR and SB were not significantly different from each other, but both were higher (or lower) than in TC (Table 2).The observation that streams do not always respond to urbanization in a clear and predictable way is not surprising, since they are influenced by numerous site-specific factors including climatic patterns, underlying geology, historical land cover, urban infrastructure material, sewage treatment, and socioeconomic factors [67].Further research at these study sites and other urban streams is needed to identify unique sources and retention of contaminants and their in-stream dynamics and transformations.

Limitations of this Study and Future Research
This study had several important limitations.One was the lack of data corresponding to intense storm events.Due to the failure of automated water level loggers, we were only able to capture the range of discharge that occurred during weekly sampling events, which included flows above the low-flow level corresponding to moderate rainfall, but not any major storms.Large storms can contribute a substantial fraction of annual nutrient and solute loads [68][69][70].C-Q plots from this study (Figures S1-S2) did not show different behavior at high and low flows, although segmented C-Q plots have been observed in other studies [27,29].It is possible that if we had more data points, especially at higher values of Q, segmented C-Q relationships would be apparent and would enable a more nuanced understanding of solute transport in these watersheds.Thus, future research should focus on quantifying concentration and discharge under a full range of flow conditions.
Additionally, the method of subsurface water sampling employed in this study resulted in a mixture of groundwater and porewater from the hyporheic zone being collected.To gain a better understanding of groundwater quality in the area, installation of semi-permanent piezometers or monitoring wells, which was beyond the scope of the present study, should be performed.Finally, it is important to recognize that the results of this study, which included three streams along an urbanization gradient, likely reflect a mixture of the general effects of urbanization and specific characteristics of the streams.The results should thus be considered within the larger context of other studies on USS.

Broader Relevance
While the field work for this study focused on three stream systems in Northern Virginia, the findings are relevant to broader environmental challenges.As the urban population increases both in the U.S. and globally, land is being converted from forest to increasing levels of urban development [71].Hydrologic connectivity between groundwater and surface water is very common, so understanding the role of groundwater in mediating the effect of urbanization on stream water quality is essential.In particular, this study showed that: (1) both overland flow and baseflow were important pathways for delivering nutrients and metals to urban streams, and (2) solute concentrations had stronger and more numerous correlations in the more urbanized stream, suggesting a common source and similar transport.Understanding whether specific pollutants enter the stream via baseflow or stormflow is additionally important because it can enable researchers and managers to better predict and contextualize temporal variation in water quality and to enact more effective pollution prevention or reduction strategies.
This study is relevant to Chesapeake Bay restoration.These streams are tributaries of the Potomac River, and ultimately Chesapeake Bay, both of which have immense cultural, recreational, and economic value to residents and tourists in the Washington, DC, area [72][73][74][75].Restoration efforts in the Chesapeake Bay have been ongoing due to severe nutrient-especially nitrogen-pollution in recent few decades [76][77][78].Reinforcing the findings of previous studies in other parts of the Chesapeake basin [79][80][81], this study suggested that high concentrations of DIN in both surface and subsurface waters are not being attenuated by transformations in the hyporheic zone, allowing for downstream runoff into the Chesapeake.Thus, efforts to prevent excess DIN from entering urban streams through surface or groundwater and to restore connections with the hyporheic zone may be necessary to address the problem.

Supplementary Materials:
The following are available online at www.mdpi.com/2073-4441/13/5/662/s1.Funding: This research was supported by a National Conservation Lands System grant to K.L.K. and Dave Culver (Award #L14AC00102) and internal American University funding to J.A.B.

Acknowledgments:
We thank John Reffit, Theresa Jefferson, Christopher Brown and other current and former Bureau of Land Management staff for support and logistical assistance.Dave Culver helped write the proposal that funded the project.Colin Casey, Raquel Freitas, Melanie Friedel, Keeli Howard, Harrison Hyde, Jude Gyeongbae Jung, Jenna Keany, Melissa Knapp, Madeeha Froze Mughal, Nikita Mukherjee, Gabriel Santos, Jesse Weinstein, and William Wong assisted with field sampling and lab work.We also thank four anonymous reviewers whose thoughtful contributions helped improve this manuscript.

Data Availability Statement:
The following raw data can be requested as a downloadable Excel book.Sheets include: Metadata, Appendix 1: Spatial Data and Appendix 2: Weekly Sampling Data.

Citation:
Balerna, J.A.; Melone, J.C.; Knee, K.L.Using Concentration-Discharge Relationships to Identify Influences on Surface and Subsurface Water Chemistry along a Watershed Urbanization Gradient.

Figure 1 .
Figure 1.Map of land use in watersheds upstream of each sampling site.Watershed boundaries were generated with the USGS StreamStats tool [33].

Figure 2 .
Figure 2. Time series plots of selected surface water variables for Giles Run (blue), South Branch (red), and Thompson Creek (green) over the two sampling periods (June 2014-July 2015 and May 2016-July 2017).Daily precipitation values are shown as black bars.The gap within each panel indicates the break between sampling periods.Alternating grey and white bars indicate the seasons: summer (Su), fall (F), winter (W), and spring (Sp).High outliers with values greater than the upper limit of the y axis are shown with an up arrow and the numerical value on the graph. 1-0.4)

Figure 3 .
Figure 3.Time series plots of measured subsurface water variables for Giles Run (blue), South Branch (red), and Thompson Creek (green) from September 2016 to July 2017.Alternating grey and white bars indicate the seasons: summer (Su), fall (F), winter (W), and spring (Sp).Black bars indicate daily precipitation.

Figure 4 .
Figure 4. Log [concentration] -log [discharge] plots of selected surface water variables for Giles Run (blue), South Branch (red), and Thompson Creek (green).Plots are only shown when at least one stream had a C-Q relationship with a slope significantly different from 0. Regression lines of best fit included when the relationship was significant at p < 0.05.See figure S1 for C-Q plots of all variables.

Figure 5 .
Figure 5. Log [concentration] -log [discharge] plots of selected subsurface water variables for Giles Run (blue), South Branch (red), and Thompson Creek (green).Regression lines of best fit included when the relationship was significant at p<0.05.Plots are only shown when at least one stream had a C-Q relationship with a slope significantly different from 0. See figure S2 for C-Q plots of all variables.

Table 1 .
Watershed Characteristics for Each Stream.SEM = Standard error of the mean.
Surface water was sampled approximately once per week during the periods June 2014-July 2015 and May 2016-July 2017 (n = 88 per site

Table 2 .
Average (min-max) of variables measured during weekly sampling within each stream system.Both surface water (n = 88) and subsurface water (n = 32) measurements are included.N/A (not available) indicates the parameters were not measured.

Table 4 .
Spearman correlation information for GR (blue), SB (red), and TC (green) for pH, sulfate, and dissolved metal concentrations.Measurements from both surface water (SW) and subsurface water, collected from the hyporheic zone (HZ), are included.Symbols are identical to Table3.See TableS6for r values.

Table 5 .
Spearman correlation coefficients for surface water-subsurface water in each stream system.