Removal of Fecal Indicator Bacteria by River Networks

: Fecal contamination is a signiﬁcant source of water quality impairment globally. Aquatic ecosystems can provide an important ecosystem service of fecal contamination removal. Understanding the processes that regulate the removal of fecal contamination among river networks across ﬂow conditions is critical. We applied a river network model, the Framework for Aquatic Modeling in the Earth System (FrAMES- Ecoli ), to quantify removal of fecal indicator bacteria by river networks across ﬂow conditions during summers in a series of New England watersheds of different characteristics. FrAMES- Ecoli simulates sources, transport, and riverine removal of Escherichia coli ( E. coli ). Aquatic E. coli removal was simulated in both the water column and the hyporheic zone, and is a function of hydraulic conditions, ﬂow exchange rates with the hyporheic zone, and die-off in each compartment. We found that, at the river network scale during summers, removal by river networks can be high (19–99%) with variability controlled by hydrologic conditions, watershed size, and distribution of sources in the watershed. Hydrology controls much of the variability, with 68–99% of network scale inputs removed under base ﬂow conditions and 19–85% removed during storm events. Removal by the water column alone could not explain the observed pattern in E. coli , suggesting that processes such as hyporheic removal must be considered. These results suggest that river network removal of fecal indicator bacteria should be taken into consideration in managing fecal contamination at critical downstream receiving waters.


Introduction
Fecal contamination of surface waters is a worldwide concern due to its potential to contain pathogens that cause illness [1,2]. Waterborne diseases such as giardiasis and cryptosporidiosis [3] are responsible for the death of more than 2 million people globally each year [4]. Pathogens are the most common pollutant leading to water quality impairments throughout the United States, and Escherichia coli (E. coli) is the most frequent cause of impairments [5,6]. Fecal indicator bacteria such as E. coli are used as a proxy for detecting the presence of pathogens and identifying potential health risks [7] since monitoring pathogens is expensive and time-consuming [8]. Water quality impairment by fecal contamination at any specific location is the product of bacterial sources, transport, and transformations in the upstream watershed. Key factors that must be considered include spatially and temporally varying sources to surface waters, transport and net attenuation (i.e., removal or die-off) throughout river systems [9,10]. Thus, a river network perspective is essential to improve management of fecal pollution sources.
River networks provide important regulating ecosystem services of dilution and removal that can help maintain downstream water quality [11,12]. The "sanitation ecosystem service" of fecal waste dilution and filtering by rivers has been largely neglected [13]. Natural removal of fecal contamination during downstream transport in surface waters is potentially an important factor that influences fecal contamination fluxes at critical downstream locations [14]. Processes contributing to fecal contamination removal include bacterial die-off in the water column [15] and filtering by and eventual mortality in the streambed [16], which are further influenced by hydrological conditions and water temperature [17]. These processes together constitute a potentially important aquatic ecosystem service that can improve downstream water quality [18,19]. Better understanding of this ecosystem service and how it is influenced by watershed size, distribution of sources, and flow variability [16,[20][21][22] is vital, particularly as land use and climate change.
Physical and biological processes interact to control fecal contamination removal by river networks. Physical factors include temperature [23], solar radiation [24], and hydrologic factors [25]. Hydrologic factors that influence removal include hydrological conditions and routing, which influence residence times, exchange between water column and sediments [21,26,27] and exchange with floodplains [28]. The portion of the stream bed sediments that interacts with surface water through water exchange is also known as the hyporheic zone [29]. Biological processes include mortality in the water column and in sediments. These processes result from predation by viruses and protozoa [30]. The decrease in fecal indicator bacteria concentration with time can generally be described by first-order kinetics [31]. Removal rates of E. coli, a fecal contamination proxy, have been estimated at the whole-reach scale that accounts for all channel processes, including both water column die-off and filtration by sediments [16]. This study has shown that exchange with the hyporheic zone contributes more than decay in the water column to whole channel removal. While empirical estimates of water column die-off rates have been implemented in watershed scale models, including the Soil and Water Assessment Tool (SWAT), Hydrologic Simulation Program-FORTRAN (HSPF), and Storm Water Management Model (SWMM) [32][33][34][35][36][37][38], few studies have accounted for the role of non-water column removal at watershed scales [28]. The novelty of our study is represented by quantifying the role of hyporheic removal at watershed scales.
Removal of E. coli by hyporheic zone sediments is determined by exchange of surface water containing E. coli into the sediments and the subsequent fate of E. coli. Hyporheic exchange is the process of exchange of surface channel water with pore water in sediments beneath the stream. Streambed filtration occurs as suspended particles enter the streambed and are strained by pores that are too small to pass or settle due to slower water velocities [39]. The hyporheic zone has an important influence on many biogeochemical fluxes (nitrate, dissolved organic carbon) as water flows downstream [26,40]. Microorganisms and other fine particles can also be trapped in the hyporheic zone [16,41,42]. The function of the hyporheic zone is controlled by the combination of hydrologic exchange with surface water and the physical, chemical and biological processes that often differ substantially from the overlying water column [43].
Water quality impairment due to fecal contamination is commonly associated with storm events [2]. Generally, managers use a rainfall threshold as an indicator for impairment [44,45] without consideration of antecedent moisture conditions, spatial variability of potential source areas within watersheds, effectiveness of removal by river networks, and their interactions. Much of the storm response is due to mobilization of sources [46]. At the same time, the capacity of river networks to remove biogeochemical constituents (nitrogen, carbon) is also strongly influenced by flow [20,22]. Hydrologic variability controls the removal efficiencies of materials by influencing water depth, residence time, and the proportion of flowing water entering the hyporheic zone. For example, increased streamflow has a profound impact on the proportion of fecal indicator bacteria die-off in the water column due to decline in the residence time [28]. High flows also reduce the proportion of surface water exchanging with the hyporheic zone due to changes in the ratio of cross-sectional area of transient storage zone (A s ) to main channel zone (A) and α (the per time exchange rate coefficient) [47]. A more mechanistic understanding of the response of the flux of fecal indicator bacteria to storm events would help optimize management (e.g., beach postings, shellfish bed closures) and identify opportunities for mitigating microbial pollution. Such an understanding requires a river network perspective that explicitly considers the effect of sources, transport, and removal of fecal indicator bacteria concentrations [10].
We integrated field measurements and river network scale modeling of E. coli to answer the research question: What is the ability of river networks to regulate E. coli fluxes to downstream water bodies across flow conditions? This is an important issue that affects strategies for managing upstream pollution sources where downstream drinking water, irrigation, and recreational uses may be impaired [16,48]. We hypothesized that, under low flow conditions, large river networks are able to significantly reduce E. coli pollutants from reaching downstream water bodies through removal processes that include filtering by sediments and water column die-off. We further hypothesized that the effectiveness declines with increasing flow conditions during storm events because of declining hyporheic exchange with sediments and reduced water column residence times [49]. Finally, we hypothesized that non-water column removal is an important contributor of E. coli removal during downstream transport.

Materials and Methods
We applied a spatially distributed modeling approach to assess the temporal variation of E. coli removal by a number of New England river networks draining watersheds that varied in size. The model includes two major components: terrestrial loading and aquatic removal. A terrestrial loading function was determined using new empirical measurements. The model simulated both water column (WC) and hyporheic zone (HZ) E. coli removal using parameters derived from previous studies (Table S3) [50,51]. Simulated E. coli concentrations were compared to independent observations at the mouth of each watershed. The model application focused on summer (1 June-31 August), when fecal contamination is of the greatest concern in most water bodies.

Study Sites
River networks draining six watersheds across New England, USA, were studied. These networks had a range of drainage areas, source land uses, and spatial arrangements of sources, including the Oyster R., NH; Winnicut R., NH; Cocheco R., NH; Lamprey R., NH; Ipswich R. MA; Merrimack R. NH/MA; and Penobscot R., ME (Table S1 and Figure 1). Selected watersheds include small watersheds (Oyster R. and Winnicut R., drainage Area, <50 km 2 ), medium watersheds (Lamprey R., Ipswich R., and Cocheco R., 400-600 km 2 ), and large watersheds (Merrimack R. and Penobscot R., >10,000 km 2 ). The percentage of the watershed classified as urban/suburban ranges from 2% in the Penobscot to 36% in the Ipswich. Agriculture is relatively unimportant (<10%) and is typically dominated by hay fields. Developed land use is skewed towards the mouth of the watershed in five of the watersheds (skewness index < 1, see Section Skewness index analysis) and towards the headwaters in the Winnicut R. and Ipswich R. (skewness index > 1). Annual precipitation averages 1100-1200 mm for all watersheds. Mean annual air temperature is 8.4, 5.6, and 9.5 • C for Merrimack R., Penobscot R., and other coastal watersheds, respectively.

E. coli River Network Model Description
To understand the fate and transport of E. coli, we used a river network model, the Framework for Aquatic Modeling in the Earth System (FrAMES), that routes flow and constituents throughout river systems. FrAMES is a gridded, fully spatially distributed hydrology and biogeochemical river network. FrAMES was chosen for this study because: (1) the model has been well tested and was able to predict streamflow and hydraulic conditions with low errors at the annual [52,53] and seasonal (summer low-flow periods) scales [53,54], (2) the model has the required characteristics for an E. coli model, including temporal resolution (daily time step) and key environmental processes (e.g., hyporheic exchange [26]) [55]. The model has been applied to assess the relative impact of water column and hyporheic zone on pollutant removal [26], which is an essential model characteristic to answer our research question, (3) the model has been applied to investigate dilution and and in-stream biogeochemical process at the whole river network to predict nitrogen removal and export [52,53,56], chloride concentrations [53,57], dissolved organic carbon [54], and stream temperature [53], and (4) the model is a fully spatially distributed model, which can provide better spatial representation of hydrology and E. coli dynamics better than lumped or semi-distributed models. We added an E. coli component to FrAMES (FrAMES-Ecoli) that simulates the sources, transport, and fate of E. coli at the daily time step and accounts for in-stream die-off of E. coli in both the water column (WC) and the hyporheic zone (HZ) ( Figure S1).

E. coli River Network Model Description
To understand the fate and transport of E. coli, we used a river network model, the Framework for Aquatic Modeling in the Earth System (FrAMES), that routes flow and constituents throughout river systems. FrAMES is a gridded, fully spatially distributed hydrology and biogeochemical river network. FrAMES was chosen for this study because: (1) the model has been well tested and was able to predict streamflow and hydraulic conditions with low errors at the annual [52,53] and seasonal (summer low-flow periods) scales [53,54], (2) the model has the required characteristics for an E. coli model, including temporal resolution (daily time step) and key environmental processes (e.g., hyporheic exchange [26]) [55]. The model has been applied to assess the relative impact of water column and hyporheic zone on pollutant removal [26], which is an essential model characteristic to answer our research question, (3) the model has been applied to investigate dilution and and in-stream biogeochemical process at the whole river network to predict nitrogen removal and export [52,53,56], chloride concentrations [53,57], dissolved organic carbon [54], and stream temperature [53], and (4) the model is a fully spatially distributed model, which can provide better spatial representation of hydrology and E. coli dynamics better than lumped or semi-distributed models. We added an E. coli component to FrAMES (FrAMES-Ecoli) that simulates the sources, transport, and fate of E. coli at the daily time step and accounts for in-stream die-off of E. coli in both the water column (WC) and the hyporheic zone (HZ) ( Figure S1).
E. coli in small headwater streams was used to estimate terrestrial loading due to the proximity to terrestrial source and short water residence time and minimal possibility for biogeochemical transformations [58]. To parameterize terrestrial E. coli inputs across flow conditions, we used existing E. coli data from three headwater streams (<2 km 2 ) in E. coli in small headwater streams was used to estimate terrestrial loading due to the proximity to terrestrial source and short water residence time and minimal possibility for biogeochemical transformations [58]. To parameterize terrestrial E. coli inputs across flow conditions, we used existing E. coli data from three headwater streams (<2 km 2 ) in southeastern NH collected by the New Hampshire Department of Environmental Services between 2008 and 2014 and data from our fieldwork in summer of 2015. The three headwater streams differ in land use, including (1) College Brook (urbanized, impervious = 28.4%), (2) Dube Brook (forested, impervious = 4.8%), and (3) Chesley Brook (intermediate, impervious = 6.2%). A continuous storm sampling encompassing the prior baseflow and the entire hydrograph is required to examine E. coli dynamics/behaviors through storms. Therefore, during June 2015 a storm event was sampled for E. coli, including baseflow just prior to the storm and 5-6 time points on both the rising and falling limb using an autosampler (Sigma) in all three headwater streams. Samples were kept on ice and shielded from light until processed within 24 h. E. coli concentrations were quantified using the IDEXX Colilert-18, Quanti-Tray/2000 (IDEXX) method. The mouth of the Oyster R. watershed (~50 km 2 ) was also sampled throughout the June 2015 storm, but data was reserved for model validation purposes.
E. coli loading to river networks is represented as a function of flow variability and land use, thereby accounting for temporal and spatial variability of inputs [59]. The relationship of E. coli concentrations measured in headwater streams and environmental variables was established by Partial Least Squares Regression (PLSR). PLSR has been shown to perform better than multiple linear regression because it eliminates the problem of multicollinearity in the independent variables [60] and has been recommended for predicting E. coli levels [61]. The spatial and temporal distribution of E. coli inputs to the river network (colony forming unit (CFU)/100 mL d −1 ) were as follows: log(EC) = 0.87 + 0.049 Precip + 0.046 T Air + 0.014 DEV + 0.0052 FOREST (1) where EC is E. coli concentration (CFU/100 mL d −1 ), Precip is precipitation (mm d −1 ), T Air is air temperature ( • C), DEV is the percentage of urban developed area in the watershed, and FOREST is the percentage of forested area in the watershed (R 2 = 0.60, n = 90). Our loading function can explain similar or greater E. coli variation compared to previous studies (R 2 = 0.27-0.62, Table S2). The net effect of this equation is higher loading concentrations when there is more urban area, warmer temperatures, and wetter conditions. In-stream removal of E. coli is simulated for each stream reach within the river network. The river network is represented by a series of grid cells with flow direction [62]. We accounted for E. coli removal in the water column and in the hyporheic zones of each stream reach (=grid cell) as: where R i (dimensionless) is the total proportional removal of E. coli within stream reach i, R WC,i (dimensionless) is the proportional removal of E. coli in the water column (WC) and R HZ,i (dimensionless) is proportional removal of E. coli that has entered the subsurface hyporheic zone (HZ). TE HZ,i (dimensionless) is the fraction of water and E. coli in the water column entering the hyporheic zone (HZ) within stream reach i. E. coli is assumed to travel with the water. R WC estimates the fraction of E. coli removal in the water column (WC) as a function of water travel time, die-off rate, and water temperature [31]: where K 20 is die-off rate of E. coli in the water column (WC) at 20 • C (d −1 ), τ i (d) is the residence time of water in the water column (WC) of the stream reach, T i is water temperature ( • C), and θ (dimensionless) is temperature adjustment factor.
where L i is the reach length (m), V i is the mean flow velocity (m s −1 ), and 86,400 is for unit conversion to days. T i and V i are predicted by FrAMES [63], while L i is determined by the gridded river network [26]. The proportion of water flowing into a stream reach that is exchanged between the water column (WC) and hyporheic zone (HZ) is determined as in Stewart et al. (2011) [26] and Mulholland and DeAngelis [64]: where α i is the exchange coefficient (s −1 ), A cross,i is the channel cross-sectional area (m 2 ), L i is reach length (m), and Q i is discharge (m 3 s −1 ). The downstream flux of E. coli from stream reach i (Flux i ) is calculated as: where Upstream i (CFU d −1 ) is the sum of E. coli input into stream reach i from the channel(s) immediately upstream, and Local i (CFU d −1 ) is the total input from land that drains directly to stream reach i based on the empirical loading function (Equation (1)). This study did not consider floodplains because of minimal lateral connectivity during low-flow summer periods [65].

E. coli Model Parameterization and Testing
Key parameters in the model include the E. coli die-off rate in the water column (K 20 ), the temperature adjustment factor (θ), the water column to hyporheic zone exchange coefficient (α HZ ), and the proportional removal of E. coli in the hyporheic zone (R HZ ) (Table  S3). K 20 and θ were parameterized a priori (K 20 = 0.8, θ = 1.07) from previous modeling studies of E. coli removal by the water column [32,66] or from laboratory experiments [67]. Hydraulic parameters that influence residence time and exchange with sediments, and hence removal, were estimated based on empirical functions of width, depth and velocity vs. discharge (Equations (S1)-(S6)).
The exchange rate coefficient of water mass transfer (HZ exchange parameter, Time −1 ) represents the rate at which surface water enters the transient storage [43]. The hyporheic zone (HZ) exchange parameter, α HZ , is based on measurements in first-through fifth-order stream segments within the Ipswich and Parker rivers, Massachusetts [51], and previously applied at river network scales to quantify the role of hyporheic zone (HZ) on nitrogen removal [26]. The mean α HZ from that study (9.5 × 10 −6 s −1 ) was applied in all watersheds of this study and was assumed to remain constant across flow conditions [26]. Little information exists to isolate the die-off rate in the hyporheic zone at the river network scale. Thus, we initially assumed all E. coli entering the hyporheic zone (HZ) would settle and die (R HZ = 1). As a result, variability in removal by the hyporheic zone (HZ) over time is controlled by hydrologic exchange between the water column (WC) and hyporheic zone (HZ). We also explored scenarios where R HZ is not 1 through sensitivity analysis.
Sensitivity analyses were performed to evaluate how river network scale E. coli removal changes with adjustments in hydraulic characteristics and removal rates in water column (WC) and hyporheic zone (HZ) at river network scales. We applied the sensitivity analysis for the Oyster R. watershed to examine the model responses to adjusting parameters. Each parameter was changed by +25% and −25% (except for Temperature adjustment factor (θ) and R HZ ) while leaving all others constant and quantifying the change of E. coli removal at the river network scale. Temperature adjustment factor (θ) was adjusted to the maximum and minimum values from the literature [33,34] since 25% changes exceeded the range. R HZ was only tested by −25% because R HZ cannot exceed 1.

Hypothesis Testing
To test the hypothesis that network removal is a potentially important control on downstream E. coli concentration patterns, and that the hyporheic zone (HZ) is a major contributor to network scale E. coli removal, we applied the model with three different network removal scenarios: (1) E. coli behaves conservatively and is transported with dilution only (Mixing), (2) E. coli is removed by water column processes only (RWC), and (3) E. coli is removed by both water column and hyporheic zone processes (RWC + RHZ). The Mixing scenario assumes E. coli is non-reactive and all E. coli entering surface waters is transported to the basin mouth, with concentration patterns only affected by mixing and dilution depending on land use in different tributaries, along river corridors and their weather conditions. The other two scenarios (RWC and RWC + RHZ) assume that E. coli is removed by streams and rivers. Differences in removal between RWC and RWC + RHZ represent the effect of the hyporheic zone (HZ). We applied the model, parameterized a priori, to seven New England watersheds in order to simulate E. coli input and removal during the summer periods of 2011-2014. We extended the model period to 2015 for five of seven watersheds to test the model with additional observations at the basin mouth during 2015. Model predictions of E. coli concentrations were tested against measurements at the basin mouths from June to August, 2015 across flow conditions (data sources: field measurements from the New England Sus-tainability Consortium project (http://ddc-bacteria-assessment.sr.unh.edu/about/project_ description/, accessed on 12 July 2021) and published data from New Hampshire Department of Environmental Services (https://www.waterqualitydata.us/, accessed on 12 July 2021). These measurements were not used to calibrate the model. We calculated the model percentage error (MPE) of simulated E. coli concentrations at the basin mouth for each removal scenario (Mixing, RWC, and RWC + RHZ).

Skewness Index
To understand how watershed E. coli source distributions and their associated flow path distances in the river network affect the E. coli removal efficiency at river network scales, we estimated a skewness index (SI) [56] for each watershed. SI characterizes the spatial distribution of source areas within watersheds relative to the average river network distance that a water molecule travels within the watershed. An SI > 1 means that on average E. coli enters the network further upstream than the typical water molecule, whereas an SI < 1 means the average E. coli enters closer to the basin mouth. To distinguish the contributions of watershed sizes and loading distributions on E. coli removal, we evaluated the effect of skewness on E. coli removal in watersheds with similar sizes but different SIs (Oyster R. vs. Winnicut R.; Cocheco R. vs. Ipswich R.).

E. coli Exceedance
To examine the benefit (i.e., ecosystem service) that network removal has on fecal indicator bacteria levels at the watershed mouth, we examined E. coli concentrations and exceedance probabilities of water quality standards in scenarios with and without E. coli removal. There are two classes of impairment for surface waters determined by E. coli thresholds: Class A is for water potential for public water supply, and Class B is for recreational water use. These differ from state to state, so we used standards set for NH, USA, as a demonstration. The Class A surface water quality standard for E. coli is defined by a single sample with a count of 153/100 mL. The Class B surface water quality standard for E. coli is a count of 406/100 mL [68]. We compared the number of days each standard is exceeded in each watershed based on the time series of E. coli concentrations using the Mixing and RWC + RHZ scenarios. The Mixing scenario represents the water quality without the ecosystem service of E. coli removal, while RWC + RHZ scenario represents the water quality with the full ecosystem services for E. coli removal.

Model Validation
FrAMES simulated summer discharge with reasonable accuracy (simulated vs. observed discharge R 2 : 0.46-0.73; Percentage bias (PBIAS): −16.9-+25.9%). For E. coli concentrations, median model percent error (MPE) was the highest for the Mixing scenario (112-960%) and declined when adding water column (WC) removal in the RWC scenario (44-418%) (Figure 2). MPE was the lowest and relatively unbiased after adding hyporheic zone (HZ) removal for the RWC + RHZ scenario (−73-+27% for RWC + RHZ scenario). Across the entire summer period, the model provided reasonable predictions of E. coli concentrations at all the basin mouths, particularly for the scenario with both water column and hyporheic zone removal (RWC + RHZ) (Figure 2). Consistently positive MPE values of Mixing and RWC scenarios indicate overestimation bias, consistent with the need to consider the role of the hyporheic zone (HZ) in order to understand network scale E. coli removal. To evaluate the model's performance for E. coli concentration prediction during storm events, the predicted E. coli concentration in three removal scenarios (Mixing, RWC, RWC + RHZ) were compared with the 5-day sampling data at the Oyster River basin mouth. FrAMES-Ecoli also predicted the E. coli concentration response at the basin mouth during the storm period reasonably well at the Oyster R. basin mouth during the 2015 storm sampling, with the RWC + RHZ scenario coming closest to observations (Figure 3). The comparison of predictions against observations at the basin mouth represent an indepen- understand network scale E. coli removal. To evaluate the model's performance for E. coli concentration prediction during storm events, the predicted E. coli concentration in three removal scenarios (Mixing, RWC, RWC + RHZ) were compared with the 5-day sampling data at the Oyster River basin mouth. FrAMES-Ecoli also predicted the E. coli concentration response at the basin mouth during the storm period reasonably well at the Oyster R. basin mouth during the 2015 storm sampling, with the RWC + RHZ scenario coming closest to observations (Figure 3). The comparison of predictions against observations at the basin mouth represent an independent test of model prediction, since these data were not used to calibrate the model. We conclude that the model is reasonable for exploring factors controlling E. coli removal by river networks.  concentration prediction during storm events, the predicted E. coli concentration in three removal scenarios (Mixing, RWC, RWC + RHZ) were compared with the 5-day sampling data at the Oyster River basin mouth. FrAMES-Ecoli also predicted the E. coli concentration response at the basin mouth during the storm period reasonably well at the Oyster R. basin mouth during the 2015 storm sampling, with the RWC + RHZ scenario coming closest to observations (Figure 3). The comparison of predictions against observations at the basin mouth represent an independent test of model prediction, since these data were not used to calibrate the model. We conclude that the model is reasonable for exploring factors controlling E. coli removal by river networks.

Sensitivity Analysis
To investigate the uncertainty of parameters, we applied sensitivity analysis on parameters that influence aquatic removal (A cross , α HZ , R HZ , K 20 , θ) to the Oyster R. watershed. Sensitivity results indicated that hyporheic zone (HZ) E. coli removal was higher than water column (WC) E. coli removal for all parameter sets (Table S4). Both R WC and R HZ increased with increasing A cross (and associated lower water velocities) due to the increase of hydraulic residence time and hyporheic storage. Water column (WC) percentage E. coli removal was most sensitive to the temperature adjustment factor (θ), which was consistent with previous results reported by Coffey [69], and Sowah et al. [70]. R WC was 3% and 16% for θ = 1.12 and θ = 1.02, respectively. The increase in R WC occurs with θ = 1.02 because the die-off rate in the water column (WC) is higher in this scenario across the range of water temperature. K 20 had a slight effect on R WC . Reducing α HZ and R HZ from 100% to 75% caused network scale RHZ to drop from 20% to 16%. Changes in parameters for one of the compartments result in slight increases in removal in the other compartment, indicating compensatory effects (more supply and, as a result, removal in downstream reaches).

Whole River Network E. coli Removal Using Full Removal Scenario
River networks removed 26-87% of the E. coli inputs to surface waters during transport over the entire summer in the full removal scenario RWC + RHZ (Figure 4). Most removal in this scenario occurred in the hyporheic zone (HZ), which removed 19-72% of E. coli inputs, while water column (WC) removal accounted for 6-16% (Figure 4). E. coli removal increased with watershed size. The smallest proportion removed was in the Oyster R. (26%, 44 km 2 ) and Winnicut R. watersheds (28%, 40 km 2 ) and the largest proportion removed was in the Penobscot (87%, 22,691 km 2 ). There was a limited amount of E. coli removed in small watersheds due to short residence times and limited potential to exchange with the hyporheic zone (HZ).
higher than water column (WC) E. coli removal for all parameter sets (Table S4). Both RWC and RHZ increased with increasing Across (and associated lower water velocities) due to the increase of hydraulic residence time and hyporheic storage. Water column (WC) percentage E. coli removal was most sensitive to the temperature adjustment factor (θ), which was consistent with previous results reported by Coffey et al. (2010) [32], Niazi et al. (2015) [33], Kondo et al. [69], and Sowah et al. [70]. RWC was 3% and 16% for θ = 1.12 and θ = 1.02, respectively. The increase in RWC occurs with θ = 1.02 because the die-off rate in the water column (WC) is higher in this scenario across the range of water temperature. K20 had a slight effect on RWC. Reducing αHZ and RHZ from 100% to 75% caused network scale RHZ to drop from 20% to 16%. Changes in parameters for one of the compartments result in slight increases in removal in the other compartment, indicating compensatory effects (more supply and, as a result, removal in downstream reaches).

Whole River Network E. coli Removal Using Full Removal Scenario
River networks removed 26-87% of the E. coli inputs to surface waters during transport over the entire summer in the full removal scenario RWC + RHZ (Figure 4). Most removal in this scenario occurred in the hyporheic zone (HZ), which removed 19-72% of E. coli inputs, while water column (WC) removal accounted for 6-16% (Figure 4). E. coli removal increased with watershed size. The smallest proportion removed was in the Oyster R. (26%, 44 km 2 ) and Winnicut R. watersheds (28%, 40 km 2 ) and the largest proportion removed was in the Penobscot (87%, 22,691 km 2 ). There was a limited amount of E. coli removed in small watersheds due to short residence times and limited potential to exchange with the hyporheic zone (HZ).

Effect of Watershed Size
E. coli removal was higher in large than in small watersheds (Figure 4), particularly at high flows ( Figure 5). Watershed size had a greater influence on network scale removal during high flow than during low flow conditions. The fraction of E. coli removed during low flow periods was high and relatively unaffected by watershed area (ranging from 68-99% from small to large watersheds). In contrast, removal was lower during high flows, with a greater range of removal across different watershed areas, ranging from 19-85% from small to large watersheds ( Figure 5). The reduced removal efficiency at high flow was due mainly to declines in hyporheic zone (HZ) removal ( Figure 5). at high flows ( Figure 5). Watershed size had a greater influence on network scale removal during high flow than during low flow conditions. The fraction of E. coli removed during low flow periods was high and relatively unaffected by watershed area (ranging from 68-99% from small to large watersheds). In contrast, removal was lower during high flows, with a greater range of removal across different watershed areas, ranging from 19-85% from small to large watersheds ( Figure 5). The reduced removal efficiency at high flow was due mainly to declines in hyporheic zone (HZ) removal ( Figure 5).

Skewness Index Analysis
Distribution of land use within the watershed had a modest influence on network scale removal proportions. Watersheds with source areas (developed land use) skewed towards the distant headwaters (i.e., high SI) removed roughly 2-16% more E. coli than

Skewness Index Analysis
Distribution of land use within the watershed had a modest influence on network scale removal proportions. Watersheds with source areas (developed land use) skewed towards the distant headwaters (i.e., high SI) removed roughly 2-16% more E. coli than those where E. coli sources are skewed towards the mouth ( Figure 6). For example, the Oyster R. watershed (44 km 2 ) and Winnicut R. watershed (40 km 2 ) both have similar watershed sizes, yet the Winnicut R. network (SI = 1.06) removed 2% more E. coli than the Oyster R. network (SI = 0.77) at high flows and 16% more at low flows ( Figure 6). A similar response occurs in medium-sized watersheds ( Figure 6). The Ipswich R. network (SI = 1.09) removed 8% more E. coli than the Cocheco R. network (SI = 0.72) at the high flows and 13% more at the low flows.
Oyster R. watershed (44 km 2 ) and Winnicut R. watershed (40 km 2 ) both have similar watershed sizes, yet the Winnicut R. network (SI = 1.06) removed 2% more E. coli than the Oyster R. network (SI = 0.77) at high flows and 16% more at low flows ( Figure 6). A similar response occurs in medium-sized watersheds ( Figure 6). The Ipswich R. network (SI = 1.09) removed 8% more E. coli than the Cocheco R. network (SI = 0.72) at the high flows and 13% more at the low flows.

Prediction of E. coli Level Exceedance Probabilities
Removal of E. coli by river networks reduces the probability of exceeding water quality thresholds at downstream locations compared to if river network removal was not considered. The RWC + RHZ scenario had a much lower probability of exceeding water quality thresholds at the basin mouth (mean of 19% exceedance for class A, and 9% for class B) for all watersheds than when assuming conservative mixing of sources alone (mean of 83% and 36%, Figure 7). In all watersheds but the Penobscott R., without E. coli removal by the river network, E. coli levels exceed the Class A standard more than 90% of the time, including during low flows ( Figure 7A), and even the loading from the landscape is much lower at this time (e.g., Figure 3). With E. coli removal, only high flows show exceedance, consistent with storm observations (e.g., Figure 3) and the rule of thumb management of closing shellfish beds following precipitation events.
The response is similar for the Class B standard, though the mixing-only scenario results in less severe frequency and greater variability compared to the higher Class A exceedances ( Figure 7B). The Class B exceedance assuming mixing only ranged from < 10% of days in the Penobscot and >80% of days in the highly developed Ipswich watershed (Figure 7). Exceedance probabilities tended to decrease as watershed size increased for the RWC + RHZ scenario (grey bars in Figure 7). The Penobscot R. watershed showed low exceedance probability even without the ecosystem service (Mixing scenario) because the watershed is relatively undeveloped (Table S1). The Ipswich watershed is the most developed watershed in this study, but shows the greatest difference in exceedance probability between Mixing and RWC + RHZ scenarios. In this watershed, most sources are located in the headwaters (the highest SI of all the watersheds) and, as a result, a large proportion of fecal indicator bacteria are removed during transport.

Prediction of E. coli Level Exceedance Probabilities
Removal of E. coli by river networks reduces the probability of exceeding water quality thresholds at downstream locations compared to if river network removal was not considered. The RWC + RHZ scenario had a much lower probability of exceeding water quality thresholds at the basin mouth (mean of 19% exceedance for class A, and 9% for class B) for all watersheds than when assuming conservative mixing of sources alone (mean of 83% and 36%, Figure 7). In all watersheds but the Penobscott R., without E. coli removal by the river network, E. coli levels exceed the Class A standard more than 90% of the time, including during low flows ( Figure 7A), and even the loading from the landscape is much lower at this time (e.g., Figure 3). With E. coli removal, only high flows show exceedance, consistent with storm observations (e.g., Figure 3) and the rule of thumb management of closing shellfish beds following precipitation events.
The response is similar for the Class B standard, though the mixing-only scenario results in less severe frequency and greater variability compared to the higher Class A exceedances ( Figure 7B). The Class B exceedance assuming mixing only ranged from <10% of days in the Penobscot and >80% of days in the highly developed Ipswich watershed (Figure 7). Exceedance probabilities tended to decrease as watershed size increased for the RWC + RHZ scenario (grey bars in Figure 7). The Penobscot R. watershed showed low exceedance probability even without the ecosystem service (Mixing scenario) because the watershed is relatively undeveloped (Table S1). The Ipswich watershed is the most developed watershed in this study, but shows the greatest difference in exceedance probability between Mixing and RWC + RHZ scenarios. In this watershed, most sources are located in the headwaters (the highest SI of all the watersheds) and, as a result, a large proportion of fecal indicator bacteria are removed during transport.

Discussion
Hydrologic networks do not act as neutral pipes for many biogeochemical fluxes [71][72][73]. We demonstrate that this is likely to be true also for E. coli. River networks are able to remove a significant proportion of E. coli inputs from land and reduce the impact on downstream water bodies. They have much greater removal efficiency at lower flows and in large river networks. The fact that most observed contamination of water bodies (e.g., shellfish bed closures; beach postings) is associated with storm events and high flows [74,75] is likely in part to be due not only to additional sources at high flows, but also the lower capacity of river networks to remove inputs during high flows. Removal of fecal indicator bacteria by river networks is thus an important ecosystem service.
The hyporheic zone is potentially a significant influence on E. coli removal at river network scales, and it can act as the rivers' "liver" [76], processing anthropogenic pollutants and preventing them from being transported downstream [77]. If water column-sediment interaction is ignored, predicted concentrations often greatly exceed

Discussion
Hydrologic networks do not act as neutral pipes for many biogeochemical fluxes [71][72][73]. We demonstrate that this is likely to be true also for E. coli. River networks are able to remove a significant proportion of E. coli inputs from land and reduce the impact on downstream water bodies. They have much greater removal efficiency at lower flows and in large river networks. The fact that most observed contamination of water bodies (e.g., shellfish bed closures; beach postings) is associated with storm events and high flows [74,75] is likely in part to be due not only to additional sources at high flows, but also the lower capacity of river networks to remove inputs during high flows. Removal of fecal indicator bacteria by river networks is thus an important ecosystem service.
The hyporheic zone is potentially a significant influence on E. coli removal at river network scales, and it can act as the rivers' "liver" [76], processing anthropogenic pollutants and preventing them from being transported downstream [77]. If water column-sediment interaction is ignored, predicted concentrations often greatly exceed measured concentrations, particularly at low flows (e.g., Rehmann et al. [78]). Our model results account for water exchange rates between the water column and the hyporheic zone [51] and E. coli removal by the hyporheic zone. Drummond et al. [16] using field measurements reported that the hyporheic zone is the dominant control of fecal indicator bacteria removal at whole stream reach scales, suggesting that this is a reasonable assumption, although they did not specifically isolate removal in each compartment. The sensitivity analysis further suggested that, even if removal of fecal indicator bacteria once they enter the sediments is reduced (R Hz = 0.75, with the rest remobilized), the hyporheic zone would still contribute significantly to network scale removal. Thus, despite the uncertainty in hyporheic removal estimation, a significant amount of E. coli removal by river networks is likely to occur in the hyporheic zone.
Hydrologic transport, storage, and reaction rates interact to control fecal indicator bacteria [79]. Hydrological variability controls dynamics of both E. coli input and removal. A disproportionate input of E. coli occurs at higher flow conditions at the same time that E. coli removal efficiency decreases. Removal efficiency decreases more rapidly for hyporheic zone removal than water column removal during storms because the proportion of flowing water that exchanges with the hyporheic zone declines rapidly due to increasing depths and declining residence times [41]. The decline in water column removal is not as steep because water column removal is represented as a volumetric process (i.e., per time decay), and therefore is only affected by the decline in residence time [80]. The net effect is a reduction in removal efficiency by river networks due mainly to a decline in filtering ability of the hyporheic zone.
Watershed size and land use distribution are also important factors influencing E. coli removal by river networks. Removal efficiency is higher in large watersheds than small watersheds across all flow categories because of an increased proportion of E. coli retained in the hyporheic zone and a longer residence time for E. coli to be removed in the water column. The decline in removal efficiency at the river network scale with increasing flow was less in larger watersheds, indicating they have a greater buffering effect on water quality [81,82]. Previous studies also found that large rivers can play a disproportionate role in the function of entire river networks [81,83]. Land use (i.e., source) distribution was also a factor controlling E. coli removal as previously shown for nutrients [56], though had less of an impact than other factors described above. Nevertheless, the watershed with E. coli sources (i.e., suburban land use) furthest from the mouth (the Ipswich R.) showed a much greater drop in exceedance due to removal by the river network compared to a similarly sized watershed with E. coli sources closer to the mouth (the Cocheco R. watershed) ( Figure 6). Development in larger watersheds that is skewed towards the distant headwaters therefore has less impact on bacterial impairments. Conversely, a focus on fecal contamination management should target sources in small watersheds with land use skewed near critical water bodies. Refinement of rainfall threshold for shellfish bed closures depending on the type of watershed and where rain falls within a watershed is also a possibility.
The model presented in this study is a useful tool for understanding the potential role of river networks in removing E. coli. However, a number of uncertainties remain. The role of the hyporheic zone in E. coli removal is based on an assumption regarding fate once E. coli are in the hyporheic zone, since there is little empirical information to inform decay of E. coli once in this zone. We assume that hydrologic exchange is the limiting process for E. coli removal, and once in the hyporheic zone, 100% of the E. coli remains out of the water column. This proportion is likely to vary depending on sediment characteristics. Further, the hydrologic exchange parameter (α) is based on a single watershed that focused on baseflows [51]. Little is known what controls variability in this parameter across space and as a function of flow. Despite coming from animal gastrointestinal tracts, fecal-borne bacteria can persist in aquatic ecosystems [84]. To avoid intensive data requirements and model complexity, we only include the major controlling factors, such as streamflow, upstream-downstream linkages, vertical exchange between water column and hyporheic zone, and water temperature. Water temperature is typically considered the primary factor influencing E. coli removal [85], and it has more significant effects than nutrient depletion and total suspended sediment concentration [23,31,86]. Other factors are also likely to affect die-off rates in both the water column and the sediments (redox conditions, nutrients, temperature, solar radiation, predation), but were not explicitly considered here. To the best of our knowledge, there is no watershed scale model that includes predation into E. coli removal estimation (Table S5). However, predation might be indirectly reflected by temperature. Laboratory experiments have found that predation rates in environmental waters were positively correlated to temperature [87,88]. This could be explained by the presence of a larger community of protozoa when the temperature increases [88]. We assumed that sunlight exposure has limited effects on E. coli removal. Intermediate streams or floodplains might receive direct solar radiation. However, there are substantial reductions in solar radiation in small streams (due to riparian shading) [89] or large rivers (due to water depth) [85]. Although canopy cover can reduce solar radiation exposures of rivers [90], especially in forest-dominated watersheds, the water column removal in shallow streams might be underestimated without considering solar radiation [85]. Interaction with sediments, suspended particles and organic debris can result in persistence of fecal-borne bacteria [84]. Resuspension of sediments could re-introduce these bacteria to the water column and be sources for surface waters [91,92]. The significance of these factors can be taken into account by future models when measurements are available.
The loading function in this study was based on limited data from three headwater sites. Although the model performed with reasonable accuracy (Figure 2), our results are not predictive, but more heuristic, so that we could demonstrate the potential role of water column and hyporheic zone removal by entire river networks. If the loading function were biased high in the model (E. coli loading is actually lower), then water column removal alone could be sufficient to match observed E. coli concentrations. However, given that we parameterized loading and removal in the model a priori and tested against an independent data set of observations through time at the basin mouth, we have confidence in the model mechanisms. Other models, including SWAT and HSPF, tended to over-predict E. coli concentrations when observed concentrations are low (i.e., low flows) and under-predict E. coli concentrations when observed concentrations are high (i.e., high flows) [31]. These biases may have occurred because these models did not account for the role of the hyporheic zone, which has a greater effect at base flows [29]. FrAMES-Ecoli accounts for sedimentation processes through advective hyporheic zone exchange.
Mitigation of the economic impacts of E. coli on receiving waters, including closures of shellfish beds and pollution postings at recreational beaches, requires an understanding of E. coli removal between sources in the landscape and water bodies of interest. The current water quality management agencies apply a rainfall threshold for beach closure decision making to all watersheds, despite the different E. coli removal capacity [42,43]. Without the consideration of the different E. coli removal capacity among watersheds (Figure 4), E. coli concentration might be underestimated or overestimated (Figure 2). The regulation of pollutants by aquatic ecosystems is critical to help protect downstream water bodies [93]. Without the ecosystem service of E. coli removal performed by river networks, even low flows could result in high E. coli levels that exceed standards. Anthropogenic changes such as channel straightening, dredging, and erosion could contribute to declines in network removal by reducing hydraulic complexity and hyporheic flux exchange [94], and lower the threshold of flows over which river networks are able to control E. coli fluxes. Lower downstream export due to E. coli removal contributes to improved human health. The major factor controlling removal in freshwater ecosystems, flow, is sensitive to change in climate, land use and management [95]. Small coastal watersheds, especially those with urban land cover skewed toward the basin outlet, have less removal capacity than large watersheds, making them particularly sensitive to changing climate. Decision makers should prioritize management of these watersheds to mitigate the risk of fecal contamination in coastal areas.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/w14040617/s1, Figure S1: The conceptual model of the E. coli fate and transport in the FrAMES model, Table S1: Characteristics of study watersheds, Table S2: Comparison of E. coli loading function, Table S3: A priori parameter values for the FrAMES-Ecoli module, Table S4: Sensitivity analysis of parameter impacts to WC, HZ, and total E. coli removal by the Oyster R. network compared to the Base scenario as defined in Table S3, Table S5: E. coli removal parameters included in previous studies.