A Case Study Evaluating Water Quality and Reach-, Buffer-, and Watershed-Scale Explanatory Variables of an Urban Coastal Watershed

: Land use land cover within a watershed inﬂuences stream water quality, habitat quality, and biological community structure. As development and associated impervious surface increases in a watershed, changes in storm water and nutrient inputs generally cause declines in habitat conditions and biodiversity. The ﬁrst goal of our study was to evaluate the water quality in the Charles River watershed, in which our objective (G1O1) was to establish ten 100-meter reach-scale sampling stations and conduct physical, chemical, and biological assessments. The second goal of this study was to better understand the direct and indirect effects of hierarchical variables on water quality in the Charles River watershed. Our ﬁrst objective of our second goal (G2O1) was to calculate land use land cover percentages at the pour-point subwatershed and local 100-meter buffer scale for each of our ten 100-meter reach sampling stations. Our second objective of our second goal (G2O2) was to use path analysis to determine the direct and indirect effects of land use land cover and impervious surface on water quality in the Charles River watershed. The results of G1O1 were that habitat quality assessments ranged from “marginal” to “optimal” and biological quality assessments ranged from “fair” to “good“, indicating overall “fair” or better water quality conditions in the watershed. The results of G2O2 were that our path analysis resulted in differences in effects of development between the buffer and sub-watershed scale. At the buffer scale, water quality was inﬂuenced more negatively by the percentage of developed land area versus the percentage of impervious cover. While both buffer development and habitat quality had a direct effect on Streamside Biosurvey Macroinvertebrates, buffer development also directly hindered habitat quality, thus having an indirect effect on Streamside Biosurvey Macroinvertebrates through habitat. Streamside Biosurvey Macroinvertebrate scores were shown to be more sensitive to development within the buffer versus at the sub-watershed scale, where impervious cover was a more important indicator of stream water quality. Through this small case study of 10 stations within the Charles River watershed, we illustrated how citizen-science level water quality assessments can be combined with water chemistry and hierarchical LULC data to provide insights into potential direct and indirect effects on water quality. As the ﬁelds of landscape ecology and conservation continue to grow, so does our ability to determine changes in land development and devise management strategies aimed at improving water quality.


Introduction
Streams are organized hierarchically over spatial-temporal scales ranging from largerscale catchments (stream systems, the area in which water drains through a specific downstream pour-point), which encompass the spatially smaller and temporally shorter segment scale (segment systems, 1000s meters), reach scale (reach systems, 100s meters), meso Because of its multivariate capabilities, path analyses are often conducted in ecological research to initiate decomposition of ecological complexities.
This study was part of a collaboration with the Charles River Watershed Association (CRWA) to develop a citizen-science water quality monitoring program in the Charles River watershed (CRW). Our research group developed a citizen-science water quality program that included developing training protocols and establishing and carrying out the monitoring at 10 sampling stations in October 2012 [20]. Furthermore, we took this citizen-science water quality sampling opportunity to explore if there were any direct or indirect effects of hierarchical watershed variables on water quality in the headwater stream of the CRW that may be useful for further exploring management of the CRW. Thus, there were two goals associated with this study. The first goal was to evaluate the water quality in the CRW. Our objective associated with this goal (G1O1) was to establish ten 100-meter reach-scale sampling stations and conduct physical, chemical, and biological assessments. The second goal of this study was to better understand the direct and indirect effects of hierarchical variables on water quality in the CRW so that effective conservation and management practices can be evaluated and implemented. To achieve this second goal, we had two objectives. Our first objective (G2O1) was to calculate LULC cover percentages at the pour-point subwatershed and local 100-meter buffer scale for each of our ten 100-meter reach sampling stations. Our second objective (G2O2) was to use path analysis to determine the direct and indirect effects of LULC and impervious surface on water quality at the reach, buffer, and sub-watershed scales in the CRW. The results of this study provide information on water quality in the watershed and reach, buffer, and sub-watershed scale and the LULC influences on water quality in the watershed.

Field-Area Description
The CRW is the most densely populated watershed in New England [21][22][23]. The Charles River begins in Hopkinton, Massachusetts and meanders northeastward 134 kilometers to the Boston Harbor ( Figure 1). The area of the CRW is~798 km 2 with a total of 970 km of mapped streams [24]. The watershed is covered by 39% forest, 11% wetlands, 2.96% lakes and ponds,~47% course-grained stratified drift, 22% impervious area based on National Land Cover Database (NLCD) 2011 data, and~52% developed land based on NLCD 2011 classes 21-24 [24]. Other basin characteristics include an average mean annual precipitation of 121.7 cm, a mean annual maximum air temperature over the basin of 13.3 • C, and a mean basin slope computed from a 10-meter digital elevation model (DEM) of 5.515% [24].

Study Design
Ten sampling stations within the mainstem and 1st-3rd order tributaries of the Charles River were chosen based on the CRWA recommendations for the citizen-science monitoring program and chosen based on trade-offs between watershed coverage and accessibility (Table 1; Figure 1) [20]. We then conducted pour-point sub-watershed scale and 100-meter local buffer scale geographic information systems (GIS) analyses with the goal of spatially analyzing and representing the relationship between water quality, LULC, and impervious cover at the 10 sampling stations.

Study Design
Ten sampling stations within the mainstem and 1st-3rd order tributaries of t Charles River were chosen based on the CRWA recommendations for the citizen-scien monitoring program and chosen based on trade-offs between watershed coverage a  The above stations were visually evaluated for physical habitat along a 100-meter reach and subsequently classified using the United States Environmental Protection Agency (USEPA) Rapid Bioassessment Protocol Habitat Assessment (RBP) during the October 2012 sampling effort [9]. Variables were visually observed and assessed with the Low Gradient Streams Habitat Assessment RBP, which included the 10 variables of epifaunal substrate/available cover, pool substrate characterization, pool variability, sediment deposition, channel flow status, channel alteration, channel sinuosity, bank stability (right and left), riparian zone (right and left), and bank vegetation (right and left). Based on visual observations, a score was assigned for each variable and summed to a final score of "poor" (0-59), "marginal" (60-109), "sub-optimal" (110-159), and "optimal" (160-200) [9]. Photos were taken at each site documenting habitat conditions as well as the surrounding land use. These photos were used in verifying LULC classifications presented within GIS raster data sets.

Physical-Chemical Analyses
In situ physical-chemical measurements [i.e., temperature, dissolved oxygen (DO), specific conductance (Spc), conductivity, salinity, total dissolved solids (TDS), and pH] were taken using a YSI meter (YSI Professional Plus Multiparameter Sonde, Yellow Springs, Ohio) during the October 2012 sampling effort. There were no rain events 48 h prior to sampling. At the same time, two water samples were collected at each station and were filtered using standard protocols. Chlorophyll-a analysis was performed within the lab following the Yentsch and Menzel protocol [25]. For total suspended solids (TSS)/non-volatile suspended solids (NVSS) and particulate phosphorus (PP), analyses were performed following the protocol reported by Greenburg et al. [26]. A filtrate was used for the testing of phosphates, ammonia, nitrates, nitrites, and total dissolved nitrogen (TDN) using the following protocols. A spectrophotometer (Hach DR/2400 Portable Spectrophotometer, Hach, Loveland, Colorado) was used to analyze 10 mL of sample water following the Salicylate Method for measuring Nitrogen and Ammonia, the PhosVer®3 Method for measuring phosphorus, and the Cadmium Reduction Method for measuring nitrate. TDN was Urban Sci. 2022, 6, 17 6 of 17 analyzed using a Varian DMS 100S UV visible spectrophotometer following calculations of Crumpton et al. [27].

Macroinvertebrate Sampling and Analysis
Macroinvertebrate sampling was conducted during the October 2012 sampling effort following a 20-jab multi-habitat approach [9]. In the laboratory, macroinvertebrates were identified to class and order [11,28,29]. Calculation of the Streamside Biosurvey: Macroinvertebrates (SBM) monitoring index was performed by counting the number of observed macroinvertebrates listed within three groups: sensitive, somewhat sensitive, and tolerant [11]. Each taxon was classified as rare, common, or dominant based on their abundances, and each were weighted and summed to a water quality score. This score classified a stream using a modified scoring condition of "excellent" (59-78), "good" (40-58), "fair" (20-39), or "poor" condition (0-19) [11].

Sub-Watershed and Local Buffer GIS Data Collection
Global positioning system (GPS) locations of the mid-point of each 100-meter sampling reach station was used to establish the sub-watershed scale pour-point and the center of the local 100-meter circular buffer for the following GIS analyses. Data layers used for GIS analyses included rivers, major watersheds, and a digital elevation model and were acquired from the Massachusetts Office of Geographic Information (MassGIS). A 100-meter circular buffer was delineated and included all area within 100 meters of the mid-point of each macroinvertebrate sampling station reach. A sub-watershed also was delineated that included all areas upstream of the mid-point of each macroinvertebrate sampling station reach flowing through the sampling point. Buffer scale and sub-watershed scale LULC and impervious cover data was acquired for each site from the NLCD [24].
The LULC data was in raster format, which was produced in 2011 with a 30-meter resolution. These data were grouped into 16 classes, and each was given an ID with a "W" for watershed or "B" for buffer to indicate the scale at which each was being measured. Development was broken down into various levels including developed-open (BDevOpen, WDevOpen), developed-low intensity (BDevLow, WDevLow), developed-medium intensity (BDevMed, WDevMed) and developed-high intensity (BDevHigh, WDevHigh). Other LULC classes included barren land (BBarr, WBarr), shrub/scrub (BShrub, WShrub), grassland/herbaceous (BHerb, WHerb), and pasture/hay/cultivated crops (BPlantCult, WPlantCult). Deciduous, evergreen, and mixed forests were combined to create one forests category (BForest, WForest), and woody, emergent, and herbaceous wetlands were combined to create one wetlands category (BWetlands, WWetlands). Proportions of LULC were taken by calculating the percentage of pixels representative of each different land use classification within the sub-watershed or buffer.
The impervious cover raster dataset also was produced in 2011 with a 30-meter resolution. Data were labeled as BImpCov and WImpCov, each representing the percentage of impervious cover within the total area of each 100-meter buffer or sub-watershed.
We utilized the above data files to perform the following functions using ESRI®ArcMap®10.2.2 (ESRI, Redlands, California, USA): (1) plot sample stations on CRW map using "Display XY data" function; (2) delineate sub-watersheds using the following tools: "Fill DEM", "Sink", "Fill", "Flow Direction", and "Flow Accumulation."; (3) create 100 meter buffers around each watershed using the "Buffer" tool; (4) create new layers with LULC data clipped to the buffer region using the "Clip" tool; (5) create new layers with LULC data clipped to the sub-watershed; (6) create new layers with impervious cover data clipped to the buffer region; (7) create new layers with impervious cover data clipped to the sub-watershed; (8) calculate percentages of various land uses within each buffer and sub-watershed; and (9) calculate percentages of impervious cover by area within each buffer and sub-watershed using the Raster Properties tool.

Data Reduction
Due to the small number of stream sites (n = 10), a reduction in the number of variables to be included within the path analysis was needed. The subsequent chronological data reduction steps were followed to reduce our 49 variables discussed in the previous sections down to 11 variables. A principal component analysis (PCA) was performed on the non-dimensionalized, z-transformed physical, chemical, LULC, impervious cover, and macroinvertebrate variables using PCOrd™ v6 (MJM™ Software Designs, Gleneden Beach, OR, USA). Axis 1 represented 35.79% of the variance, axis 2 represented 15.76% of the variance, and axis 3 represented 11.94% of the variance for a total of 63.49% of the variance. Based on the PCA analysis, 23 variables significantly (tau ≥ ±0.210) representing the variance in the data were retained for additional data reduction. These variables included: DO, TDS, Spc, salinity, NVSS, PP, TSS, Chl_a, WDevMed, WDevHigh, WImpCov, BDevOpen, BDevLow, BDevMed, BDevHigh, BImpCov, Habitat, WForest, WWetland, BForest, BWetland, WHerb, and WShrub.
Further data reduction was performed by aggregating all classes of developed land (developed-open; developed-low intensity; developed-medium intensity; developed-high intensity) into one "Developed" category and an "Undeveloped" category that was created out of a sum of the forest, wetland, shrubland, and herbaceous classes. We then further reduced the dataset using Pearson Correlation Coefficients (PCC) using IBM®SPSS®Statistics Grad Pack 22.0 (IBM, Armonk, New York, NY, USA) with a P-value threshold of ≤0.05. A professional opinion was used in the removal of related variables and if any uncertainty in the relationship remained between the two variables, the variable was maintained in the dataset.

Path Analysis
Eleven variables were retained for use within the path analysis and included TSS, DO, SpC, habitat, SBM, BImpCov, WImpCov, WDev, BDev, WUndev, and BUndev. A path analysis was performed on numerous preliminary models (>4) using IBM®SPSS®AMOS 20 (IBM, Armonk, New York, NY, USA) to determine which variables were the most significant influences on water quality directly and indirectly. Standardized estimates were given within the model as regression weights and are representative of the portion of the variance of the dependent variable that can be explained by its indicator variables. Direct effects were considered significant if the p-value was <0.05. The indirect effects were given within the AMOS output and were the sum of all products of two or three direct effects. These indirect path effects were considered significant if both direct path effects used in the calculation were also statistically significant. Model fit was determined by looking at the Chi-square, probability level, the goodness of fit index (GFI), normed fit index (NFI), relative fix index (RFI), confirmatory fit index (CFI), and root mean square error of approximation (RMSEA) for each proposed model.

Water Quality
Reach-scale habitat quality scores ranged from 68 to 167 with a narrative range of "marginal" to "optimal" habitat quality ( Figure 2). Reach-scale macroinvertebrate SBM scores ranged from 20 to 40, with nine out of 10 stations having scores of 20-39 ( Figure 2). This resulted in only one station having a "good" water quality ranking, with the remaining stations falling in the "fair" category for water quality classification (Figure 2).
Reach-scale in situ measurements of temperature, DO, specific conductance, conductivity, salinity, and pH varied little across stations ( Table 2). The mean pH among all stations within the watershed was 7.22, with the highest pH at GB (9.14) and the lowest at ELDS (5.88). Temperature averaged 14.1 • C with the highest temperature at 35CS (18.9 • C) and the lowest at ELDS (10.5 • C). The mean dissolved oxygen saturation across all sta-tions was 62.5%, with CHBR having the highest saturation of 95.8% and ELDS having the lowest saturation of 8.3%. Salinity averaged 0.2 ppt with the highest salinity at CHBR (0.4 ppt) and the lowest at GB (0.1 ppt). The mean specific conductance was 429.3 µS/cm, with CHBR having the highest specific conductance (817.0 µS/cm) and GB having the lowest (120.4 µS/cm). TDS averaged 307.9 mg/L, with CHBR having the highest TDS (533.0 mg/L) and GB having lowest TDS (76.1 mg/L). For six out of 10 stations, NVSS was reported to be zero or likely below the detection limit (BDL), while it was the highest for MINE at 16.6 mg/L. TSS averaged 10.4 mg/L with STOP having the highest recorded TSS (34.6 mg/L) and ELDS having a TSS of 0.0 mg/L or likely BDL. This resulted in only one station having a "good" water quality ranking, with the remaining stations falling in the "fair" category for water quality classification ( Figure 2). Reach-scale in situ measurements of temperature, DO, specific conductance, conductivity, salinity, and pH varied little across stations ( Table 2). The mean pH among all stations within the watershed was 7.22, with the highest pH at GB (9.14) and the lowest at ELDS (5.88). Temperature averaged 14.1°C with the highest temperature at 35CS (18.9°C) and the lowest at ELDS (10.5°C). The mean dissolved oxygen saturation across all stations was 62.5 %, with CHBR having the highest saturation of 95.8% and ELDS having the lowest saturation of 8.3%. Salinity averaged 0.2 ppt with the highest salinity at CHBR (0.4 ppt) and the lowest at GB (0.1 ppt). The mean specific conductance was 429.3 µS/cm, with CHBR having the highest specific conductance (817.0 µS/cm) and GB having the lowest (120.4 µS/cm). TDS averaged 307.9 mg/L, with CHBR having the highest TDS (533.0 mg/L) and GB having lowest TDS (76.1 mg/L). For six out of 10 stations, NVSS was reported to be zero or likely below the detection limit (BDL), while it was the highest for MINE at 16.6 mg/L. TSS averaged 10.4 mg/L with STOP having the highest recorded TSS (34.6 mg/L) and ELDS having a TSS of 0.0 mg/L or likely BDL. Note: pH, temperature (Temp), dissolved oxygen (DO), specific conductance (SpC), salinity (Sal.),  Note: pH, temperature (Temp), dissolved oxygen (DO), specific conductance (SpC), salinity (Sal.), chlorophyll-a (Chl-a), total dissolved solids (TDS), non-volatile suspended solids (NVSS), total suspended solids (TSS), below detection limit (BDL). * = below detection limit; ** = likely instrument error.

GIS Analysis
The sub-watershed areas for the 10 sampling stations ranged from 21,375 m 2 at BB to 4,290,975 m 2 at 90CS. Each sub-watershed showed some level of development, with 35CS and BB being 100% developed and STOP being 10% developed ( Figure 3, Table 4). Medium intensity developed had the highest values for all sub-watersheds and buffers followed by low intensity and then high intensity. Meanwhile, at the buffer level, 35CS, CHBR, and GB were 100% developed, and ELDS was the only station that had no development within the buffer.

GIS Analysis
The sub-watershed areas for the 10 sampling stations ranged from 21,375 m 2 at BB to 4,290,975 m 2 at 90CS. Each sub-watershed showed some level of development, with 35CS and BB being 100% developed and STOP being 10% developed ( Figure 3, Table 4). Medium intensity developed had the highest values for all sub-watersheds and buffers followed by low intensity and then high intensity. Meanwhile, at the buffer level, 35CS, CHBR, and GB were 100% developed, and ELDS was the only station that had no development within the buffer.

Path Analysis
Throughout the model building process, various models were developed and analyzed for model fit. The final model that most accurately represented the path of influence of variables on reach-scale SBM water quality within the CRW included the following observed, exogenous variables: watershed-scale developed (WDev) and watershed-scale impervious cover (WImpCov) (Figure 4). The observed, endogenous variables included buffer-scale developed (BDev), buffer-scale impervious cover (BImpCov), and reach-scale habitat, SBM and DO. Unobserved exogenous variables within the model included errors for all other observed exogenous variables. Model fit was determined by looking at the Chisquare, probability level, the goodness of fit index (GFI), normed fit index (NFI), relative fix index (RFI), confirmatory fit index (CFI), and root mean square error (RMSEA) for each proposed model. The Chi-square of 2.853 was representative of the fit of the proposed model versus the saturated, or perfect, model. A chi-square closest to our degrees of freedom (3) was optimal, and the probability level needed to be ≥0.05 to demonstrate that the hypothesized model was not significantly different from the saturated model at 95% statistical confidence. The probability level of our model was 0.415, indicating that the model was not significantly different than the saturated or "perfect" model. For good model fit, the GFI, NFI, and RFI values given in the model output should be >0. 9. The values for this model were 0.913, 0.973, and 0.810, respectively. The CFI value should be close to 1 and our CFI was exactly equal to 1. The RMSEA should be <0.05 and our RMSEA was 0.000. Aside from the RFI value, all other indices indicate a good model fit.
buffer-scale developed (BDev), buffer-scale impervious cover (BImpCov), and reach-scale habitat, SBM and DO. Unobserved exogenous variables within the model included errors for all other observed exogenous variables. Model fit was determined by looking at the Chi-square, probability level, the goodness of fit index (GFI), normed fit index (NFI), relative fix index (RFI), confirmatory fit index (CFI), and root mean square error (RMSEA) for each proposed model. The Chi-square of 2.853 was representative of the fit of the proposed model versus the saturated, or perfect, model. A chi-square closest to our degrees of freedom (3) was optimal, and the probability level needed to be ≥0.05 to demonstrate that the hypothesized model was not significantly different from the saturated model at 95% statistical confidence. The probability level of our model was 0.415, indicating that the model was not significantly different than the saturated or "perfect" model. For good model fit, the GFI, NFI, and RFI values given in the model output should be >0. 9. The values for this model were 0.913, 0.973, and 0.810, respectively. The CFI value should be close to 1 and our CFI was exactly equal to 1. The RMSEA should be <0.05 and our RMSEA was 0.000. Aside from the RFI value, all other indices indicate a good model fit. At the watershed-scale, WDev, and WImpCov were found to be highly correlated with each other with a correlation coefficient of 0.95. While both datasets were compiled by the NLCD measuring level of development, they were different in that WImpCov is strictly measuring the percentage of impervious cover. For WDev, the LULC dataset provides a categorical classification based on both impervious cover as well as an unknown "X" factor related to ghd percentage of constructed materials, types of constructed materials, and vegetation.
Ten out of 17 direct, standardized effects were shown to be statistically significant ( Figure 4). The direct effect of WImpCov on SBM was -2.02, indicating that for every At the watershed-scale, WDev, and WImpCov were found to be highly correlated with each other with a correlation coefficient of 0.95. While both datasets were compiled by the NLCD measuring level of development, they were different in that WImpCov is strictly measuring the percentage of impervious cover. For WDev, the LULC dataset provides a categorical classification based on both impervious cover as well as an unknown "X" factor related to ghd percentage of constructed materials, types of constructed materials, and vegetation.
Ten out of 17 direct, standardized effects were shown to be statistically significant ( Figure 4). The direct effect of WImpCov on SBM was -2.02, indicating that for every increase of one standard deviation in WImpCov, SBM declines 2.02 standard deviations (Table 5). Similarly, the direct effect of BImpCov on SBM was 1.65, DO on SBM was 1.41, and habitat on SBM was 0.70. DO was significantly, directly affected by BDev (1.51), but not WDev, WImpCov, or BImpCov. However, the habitat was directly affected by WDev (1.20), BDev (−1.16), and WImpCov (−1.50). Significant indirect effects also were found (Tables 5 and 6). For example, WDev positively, indirectly effected SBM through its influence on habitat (0.84). In contrast, BDev indirectly effected SBM through both habitat and DO. The indirect effect between BDev and SBM through habitat was -0.812, while it was 2.213 through DO. WImpCov also significantly, indirectly affected SBM through habitat (−1.05), but BImpCov did not. combined assessments. The second goal of this study was to better understand the direct and indirect effects of hierarchical variables on water quality in the Charles River watershed and was addressed through path analyses. Four major findings emerged from our path analysis regarding the effect of LULC and impervious cover on reach-scale water quality.

Water Quality
Two of the symptoms of the "Urban Stream Syndrome" concept is that urban streams physical habitat and aquatic assemblages should be degraded compared to more pristine stream systems [5,30]. Our habitat and macroinvertebrate assemblage water quality assessments were consistent with "Urban Stream Syndrome" concept expectations. We did observe degradation-"fair" or better water quality in the watershed-with our habitat quality assessments ranging from "marginal" (i.e., "fair") to "optimal" (i.e., "good") and biological quality assessments ranging from "fair" (i.e., "marginal") to "good" (i.e., "optimal"/ "sub-optimal"). Furthermore, our CRW habitat and SBM findings were similar, but slightly higher, to findings of a similar study on the adjacent coastal urban Boston Harbor Massachusetts watershed, where the Neponset River watershed habitat scores ranged from "poor" to "optimal" and SBM scores ranged from "poor" to "fair" [31].

Impervious Cover
Our first path analysis main finding was that impervious cover had the strongest influence on water quality of all LULC variables when this analysis was conducted at the sub-watershed scale. The percentage of impervious cover had a strong negative effect on SBM, implying that higher percentages of impervious cover within the sub-watershed create a sub-optimal environment for macroinvertebrates. Schiff and Benoit studied urban streams in New Haven, Connecticut and found that as impervious cover increased between 5 and 10%, instream conditions worsened and leveled off at 10%, demonstrating that instream habitat was impaired, macroinvertebrate diversity declined, and water quality was reduced [32]. Schiff and Benoit reported that the potential mechanism behind the effect of impervious cover was mainly due to a short-circuiting of the stream drainage system at the sub-watershed scale, disrupting the natural flow throughout the hydrologic and sediment cycles [32]. At the buffer scale, they concluded that impervious cover is negatively correlated with the area of riparian corridors that limits flooding and pollution attenuation and influences habitat development by reducing shading and organic matter and by increasing sediment inputs [32].
Contrary to Shiff and Benoit's buffer impervious cover influence [32], our results indicated that BImpCov had a positive, direct effect of 1.76 on SBM, meaning that for every 1 standard deviation increase in BImpCov, there is an increase in SBM by 1.76 standard deviations. Because BImpCov is strongly influenced by WImpCov, which had a strong, negative effect on SBM, this was the opposite of our anticipated result. However, based on Schiff and Benoit's hypothesis of hydrologic alteration [32], the unknown mechanism driving the negative effect of impervious cover may break down at the buffer scale.
Another explanation arises from the idea of "effective impervious surface" presented by Zarriello and Barlow, which involves measuring the percentage of impervious surface directly connected to a stream [33]. For example, if a buffer has a high percentage of impervious surface, it is likely that there will be a drain nearby. However, this drain may not necessarily connect directly to the river within the 100-meter buffer zone and therefore the runoff may be transported and released downstream outside of the buffer or the reach in which our macroinvertebrates were sampled. If a buffer has a low percentage of impervious surface, it is likely that there is some sort of vegetation or a pervious area between the road and streambank. This impervious surface would not be considered "effective" because the land intercepts runoff and serves as an additional filtration step. Perhaps our measure of impervious cover is not portraying its true relationship with runoff and water quality in the same way that a measure of "effective impervious surface" would.

Habitat and DO Positive Direct Effects on SBM
Our second major finding associated with our path analysis supported our general understanding that habitat and DO have strong, positive direct effects on SBM. When determining habitat's causal influence on SBM, we determined that as habitat scores increased by one standard deviation, SBM scores increased by 0.70 standard deviations ( Figure 4 and Table 5). Habitat has been shown in numerous studies to positively influence macroinvertebrates. In a study by O'Connor on Australian lowland streams, specieshabitat complexity curves demonstrated that more complex habitats boast a more complex macroinvertebrate community with different body sizes, sources of nutrition, and microhabitats that they occupy [34].
While DO had a positive effect of 1.41 on SBM, there have been mixed conclusions regarding DO's effect on macroinvertebrates. The USEPA concludes that dissolved oxygen is a critical variable influencing the presence and abundance of macroinvertebrates, but many other studies, such as Kaller and Kelso's study on Louisiana streams, found that the prediction of more diverse and abundant macroinvertebrates in reaches with higher DO habitats did not hold true [35]. They concluded that macroinvertebrates are highly tolerant to low DO conditions and are often actually more abundant in low DO conditions. One potential explanation of these mixed results on the influence of DO is the seasonal and diel variability of DO. Therefore, this positive effect should be analyzed further with more DO data.

Development Not Significantly Directly Effecting SBM
Our third major finding associated with our path analysis was that development did not significantly, directly affect SBM, but it did significantly, indirectly affect SBM. BDev and WDev both had strong, direct effects on habitat and therefore also indirectly affected SBM. BDev had an indirect effect of −0.812 on SBM through habitat, meaning that as development within the buffer increased, habitat declined; therefore, there is a negative indirect effect on SBM. This relationship is understandable considering that quality of habitat was determined by the buffer-scale variables and that macroinvertebrates are more likely to be found in habitats that provide shelter, food, and DO [9,11].
In a study conducted by King et al. on coastal Maryland streams, urban streams with high percentages of proximal developed land were found to experience a loss of riparian buffers that served as sources of woody debris, provide shade, and created desirable habitats [12,17]. Therefore, a lower percentage of development within the buffer would considerably influence the quality of habitat and the macroinvertebrate communities present. In contrast, in our study, WDev had a positive, direct effect on habitat, indicating a positive, indirect effect of 0.84 on SBM through habitat. This reversal and lack of a significant direct effect between WDev or BDev to SBM could be attributed to our aggregated development data rather than breaking development out into its constituent open, low, medium, and high intensities, which may have different effects on habitat and SBM individually. Moreover, because habitat is measured by reach/buffer-scale variables, the positive effect of WDev on habitat may not actually be relevant due to the difference in scale.
Numerous studies have shown why development at the watershed scale may not be as influential on macroinvertebrate communities as development at the buffer scale. In King et al.'s coastal streams study, they concluded that macroinvertebrate composition was found to change dramatically between 21 and 32% developed area within a watershed, with a 100% likelihood of change above 32% developed [17]. However, within a 250-meter buffer around each station, there was a 100% probability that a developed area greater than 22% would negatively affect macroinvertebrates, indicating the greater sensitivity of macroinvertebrates to buffer versus watershed scale development. Like our methods and results, they used a path analysis and found that many stations with higher macroinvertebrate scores had a high percentage of development at the watershed scale (21-32%). They concluded that development's effect on macroinvertebrates stems from its contribution to riparian degradation and the reduction of woody-debris within a local habitat [17].

BDev Significant Indirect Effect on SBM
Our fourth main finding associated with our path analysis was BDev's significant, indirect effect on SBM through DO was 2.13. This result indicates that as development increases within the buffer there is also a rise in DO and an increase in macroinvertebrates. However, based on King et al.'s DO results, many poorly shaded urban streams had higher DO and lower macroinvertebrate scores [17]. As mentioned before, DO can vary throughout the day and season, and therefore its influence on SBM is inconclusive and should be studied further.

Path Analysis Summary
In summary, SBM was significantly, directly affected by habitat and DO, but impervious cover at both the buffer and watershed scales had the strongest direct effects. There also were indirect effects on SBM that proved to be significant, such as WDev's positive, indirect effect, and BDev's negative, indirect effect through habitat. This is interesting because, although WDev and WImpCov are highly correlated with one another, their effects on SBM both directly and indirectly are the opposite. Similarly, WDev positively influences BDev, but WDev's effects on SBM are opposite in sign. The positive effects of WDev and BImpCov on SBM could be an example of Simpson's paradox in which aggregation of data or the inclusion of a third "lurking" variable could result in an unexpected or reversed effect [36].

Conclusions
As the amount of developed land continues to grow, so does the need to link terrestrial area to stream water quality and to consider best management practices. Through this small case study of 10 stations within the CRW, we illustrated how citizen-science level water quality assessments can be combined with water chemistry and hierarchical LULC data to provide insights into potential direct or indirect effects on water quality. As a result of our case study, we report mostly ≤ "good" (or "optimal"/ "sub-optimal") water quality in the watershed. We also report that, at the buffer scale only, the percentage of developed land area was shown to be a better indicator of water quality, as measured by SBM, than impervious cover. Buffer development was shown to hinder quality of instream habitat, indicating that macroinvertebrates are more sensitive to changes in development at the buffer scale. Additionally, quality of habitat was illustrated as both a direct, positive influence on water quality as well as an indirect connector facilitating the effect of development on water quality. At the sub-watershed scale, the primary factor controlling overall water quality was impervious cover, which generally increases with the level of development. We hypothesize that the mechanism driving this influence was the effect of impervious cover on the hydrologic regime within the sub-watershed. As impervious cover increases, so does runoff resulting in additional sediment and nutrient inputs into a water body, ultimately degrading stream integrity. As the fields of GIS and landscape ecology and conservation continue to grow, so does our ability to determine changes in land development and devise management strategies aimed at improving water quality. We illustrate, and thus suggest, that management strategies may benefit by incorporating holistic multi-ecological scale practices focusing on reducing developed land (e.g., adding green spaces) and treating increased runoff resulting from development (e.g., adding riparian buffer strips, bank stabilization) [37].  All other data are available in tables found within this paper. Anyone interested in further details about the data may contact both authors regarding their request.