Wastewater reuse and predicted ecological risk posed by contaminant mixtures in Potomac River watershed streams

A wastewater model was applied to the Potomac River watershed to provide (i) a means to identify streams with a high likelihood of carrying elevated effluent‐derived contaminants and (ii) risk assessments to aquatic life and drinking water. The model linked effluent discharges along stream networks, accumulated wastewater, and predicted contaminant loads of municipal wastewater constituents while accounting for instream dilution and attenuation. Simulations using 2016 data suggested that nearly 30% (8281 km) of streams were wastewater impacted. Low‐ to medium‐order streams had the largest range of accumulated wastewater (ACCWW%) values. ACCWW% exceeded a 1% threshold at >39% of drinking‐water intakes (varied by temporal condition). Risk assessments of municipal wastewater‐contaminant mixtures indicated that 22% (1479 km) of streams impacted by municipal wastewater (5.5% of all reaches modeled) may pose high risk to aquatic organisms under mean‐annual conditions, with fish more susceptible to chronic‐exposure effects relative to other taxa. Risk varied temporally and by stream order, with the greatest risk occurring in the summer in small streams. These findings suggest that wastewater may be an important factor contributing to environmental degradation in the Potomac River watershed.

Many biologically active chemicals introduced through domestic and industrial activities are not completely removed or transformed during wastewater treatment processes and are discharged to surface waters in the effluent. Effluent-derived contaminants have been detected in rivers across the United States and span a wide range of organic and inorganic compounds, including hormones, analgesics, hygiene and personal-care products, agrochemicals, and trace metals, among others Bradley et al., 2017;Focazio et al., 2008;Kolpin et al., 2002). Although often occurring at trace concentrations due to stream dilution, many synthetic organic compounds (e.g., pesticides and pharmaceuticals) are designed for their potency on biological pathways, resulting in potential adverse effects following exposure to environmental levels. Unlike pollutants from other point sources or the landscape, these contaminants are continuously discharged to receiving waters. Complex mixtures of effluent-derived contaminants in streams have been shown to increase the risk to aquatic organisms through interactive or cumulative effects (Barber et al., 2022;Kortenkamp, 2007;Vajda et al., 2008). Furthermore, low-flow conditions may increase stream vulnerability to these contaminants when instream dilution is reduced Rice & Westerhoff, 2015.
Assessing the risk posed by wastewater to aquatic life and drinking-water quality is a global environmental health priority (Malaj et al., 2014;Ramírez-Malule et al., 2020). Sampling and laboratory evaluations of contaminants, however, can be a costly and time-consuming endeavor.
Mechanisms to identify streams with a high likelihood of carrying elevated amounts of effluent and contaminants are needed to target sampling, assessments, and possible mitigations. Although hydrologic models like the U.S. Geological Survey (USGS) SPARROW model (Schwarz et al., 2016) typically consider WWTP discharges as one of several inputs, model outcomes often focus on water availability or contributing effects from discharged nutrient loads. Only a few models exist to assess wastewater impacts on streams across large geographic scales and broad categories of organic and inorganic contaminants. National studies assessing de facto wastewater reuse (the presence of treated wastewater in the water-supply source) resulted in some of the first geospatial models capable of predicting wastewater-impacted streams (U.S. Environmental Protection Agency, 1980; Rice et al., 2013). The models calculate accumulated wastewater (percentage of streamflow from total upstream WWTP discharges) at the intake locations of public water systems (PWSs) supplied by stream waters. Screening-level exposure models have been developed that use wastewater discharges to calculate predicted environmental concentrations (PECs) of contaminants in the effluent and receiving streams. These approaches incorporate chemical-specific information (e.g., per capita use and environmental fate and transport behavior) to estimate effluent contaminant loads discharged to surface waters, and account for cumulative upstream contributions, dilution, and instream decay during transport (GREAT-ER model, Koormann et al., 2006;PhATE model, Anderson et al., 2004;HydroROUT model, Lehner & Grill, 2013;iSTREEM® model, Kapo et al., 2016; Shenandoah Accumulated Ratio Mapper, Barber, Rapp, et al., 2019).
Wastewater models provide insight into the potential occurrence of wastewater-derived contaminants in streams. Accumulated wastewater estimates allow for examination of the potential impacts of discharges on downstream water quality and drinking-water sources, such as effects on disinfection byproducts (Weisman et al., 2019(Weisman et al., , 2021 and occurrence of pathogens (Soller et al., 2019). Estimates of PECs from wastewater models can be leveraged to provide screening-level assessment of risk within a stream reach or portion of the watershed when they are compared to known toxicity thresholds (e.g., predicted no-effect concentrations, PNECs) and used alongside risk-ranking approaches.
One common risk-ranking approach is evaluation of risk quotients (RQs), the ratio of the constituent concentration (measured or predicted) to a toxic-effect concentration (Barber et al., 2022;Donnachie et al., 2014;Johnson et al., 2017), where contaminants with larger RQs are assumed to present more risk. This approach links constituent concentrations to potential biological effects and allows comparison of individual contaminants, both within and across contaminant classes. This provides an improved understanding of relative risk so that chemicals can be prioritized for focus or regulation. RQs can also be used to examine the cumulative instream risk presented by multiple pollutants to which aquatic taxa may be simultaneously exposed. This is done by aggregating the individual contaminant RQs into a combined risk index (RI) to assess synergistic (de Zwart & Posthuma, 2005) or additive effects (Barber et al., 2022;Bliss, 1939).
Wastewater models vary in extent, scale, accessibility, flexibility, and purpose. Often, models are limited to select watersheds, chemical constituents, or model outcomes (e.g., hydrologic accumulation of wastewater but no chemical estimates). This article presents a flexible approach for wastewater reuse modeling, PEC modeling, and hydrological and chemical characterization to assess potential risk to aquatic health at the basin scale. The approach integrates wastewater and hydrologic datasets that can be applied to basins throughout the United States. Similarly, the approach to estimate PECs and perform risk-ranking assessments can be applied to any constituent for which reliable use, chemical, biological, and ecotoxicological benchmark parameters are available. This study focuses on the Potomac River watershed and is an expansion of studies on the Shenandoah River watershed (a Potomac River watershed subbasin) to examine de facto reuse and wastewater

Research Impact Statement
Modeling approaches used for screening-level assessment of wastewater reuse and risks to aquatic life provide a flexible riskmanagement tool applicable to multiple geographies and contaminants.
The objectives of this study were to (1) provide a comprehensive description of the modeling approach, inputs, and equations; (2) apply the approach to the Potomac River watershed to develop simulations of accumulated wastewater for municipal-plus-industrial wastewater contributions and PECs for select municipal wastewater-derived contaminants; (3) examine de facto reuse at surface-water PWS intakes in the context of potential risk associated with disinfection byproduct formation; and (4) conduct screening-level ecological risk assessments using RQs and additive effects of contaminants. These approaches are straightforward to implement and interpret. The purpose of this model is to identify river reaches that are impacted with accumulated wastewater or have predicted ecological risk that may require further attention by resource managers, either through targeted contaminant or biological monitoring and sampling or WWTP treatment upgrades to improve contaminant removal.

| Study area
The Potomac River watershed is the second largest watershed within the Chesapeake Bay watershed. It spans four states (Maryland, Pennsylvania, Virginia, and West Virginia) and the District of Columbia and encompasses several major cities (including Frederick, Maryland, Alexandria, Virginia, and Washington, D.C.). In 2010, the watershed supported a population of over six million people (U.S. Census Bureau, 2010;U.S. Geological Survey, 2014). Three major tributaries converge to form the main stem Potomac River (a Strahler seventh-order stream, Horton, 1945;Strahler, 1952Strahler, , 1957Figure 1a Table S1) with a period of record ≥10 years as of September 30, 2017 (U.S. Geological Survey, 2021). These records were used to calculate flow-exceedance probabilities for each of the 101 gages. Flow-exceedance probabilities represent the percent of time that specific streamflows were equaled or surpassed during a given period, where the 25th percentile flow-exceedance probability represents high-flow conditions, the 50th percentile represents median-flow conditions, the 75th percentile represents low-flow conditions, and flows exceeding the 90th percentile represent drought conditions (Krstolic & Ramey, 2012;Searcy, 1969).

| Modeling approach
The process-based model used literature-derived data (Table 1; Figure 2) to calculate the following estimates for 14,885 National Hydrography Dataset Version 2.1 (U.S. Environmental Protection Agency, 2012a) stream segments in the Potomac River watershed: (1) accumulated wastewater from municipal-plus-industrial effluent flows as a percent of total streamflow (ACCWW%); and (2) PECs for select municipal WWTP effluent-derived chemical constituents. Table S2 provides a list of key model definitions and equations that are discussed in additional detail in the methods section. Variation in ACCWW% was examined throughout the watershed and at surface-water PWSs as a measure of de facto reuse (Table S3; extracted from the Safe Drinking Water Information System, U.S. Environmental Protection Agency, 2020c). The PECs were used to calculate RQs and RIs to examine the potential ecological risk posed by the wastewater-derived chemical mixtures. The conceptual diagram (Figure 2) of key wastewater model parameters integrated with the NHDPlus V2 flow network depicts how model estimates of ACCWW% and PEC values at a stream reach represent the combined total upstream wastewater discharges and direct discharges into the segment. As in previous studies, this approach assumed complete lateral and depth-integrated mixing of water bodies Rice et al., 2013). In cases of network flow divergences, wastewater flow was routed through the major divergence as defined by NHDPlus V2. Wastewater modeling was performed using open-source Python (version 2.7.16; Python Software Foundation, 2017) and R (version 4.0.4; R Core Team, 2020) scripting languages. Statistical analyses were performed using R and OriginPro (version 2020b, OriginLab Corporation, 2020). Geographic information system tasks were performed using ArcGIS Pro (version 2.7.1; Esri Inc, 2020).

| Wastewater data
Wastewater inputs consisted of average daily National Pollutant Discharge Elimination System (NPDES) permitted WWTP effluent discharges to surface waters during 2016 (U.S. Environmental Protection Agency, 2022). Effluent discharge volumes were obtained from the Chesapeake Assessment Scenario Tool (CAST) wastewater report (Chesapeake Bay Program, 2020) with the substitution of 79 Shenandoah River watershed facilities with previously published and quality-assured 2016 discharge data (Weisman et al., 2021). The CAST tool is specific to Chesapeake Bay WWTPs; however, it synthesizes data from state databases and the U.S. Environmental Protection Agency (USEPA) Enforcement and Compliance History Online database (ECHO; U.S. Environmental Protection Agency, 2020b), which provides discharge volumes and outfall locations for facilities across the nation. Permits in CAST are classified as: (1) either industrial (e.g., non-publicly owned treatment works, NON-POTW; wastewater is derived from manufacturing, mining, and refining) or municipal (e.g., publicly owned treatment works, POTW; wastewater is derived from homes, institutions, and small businesses); and (2) as either significant or non-significant, based on jurisdictional-specific criteria (Chesapeake Bay Program, 2020). Significant WWTPs are those with large flow rates that are required to meet total maximum daily load targets for nutrient regulations. Due to these extensive monitoring and reporting requirements, information for significant facilities is typically better characterized than for non-significant facilities (Chesapeake Bay Program, 2020; Supporting Information). Significant and non-significant classification is not based on a ranked importance of water-quality impacts.
Outfall locations for each permit were linked to an NHDPlus V2 common identifier (COMID) using the following approaches, in order of best available information: (1) previously published NHDPlus V2 discharge locations Kandel et al., 2019); (2) reach codes provided by the USEPA ECHO database (U.S. Environmental Protection Agency, 2020b); and (3) spatial linkage to the nearest stream segment or catchment based on outfall locations determined by the Chesapeake Bay Model Phase 6 (Chesapeake Bay Program, 2020), TA B L E 1 Key model parameters and sources used in the calculation of accumulated wastewater percentages and predicted environmental concentrations in the Potomac River watershed.

Contaminant
Use, environmental fate and transport Daily per capita use (mg/capita/day): consumption (and percent removal with human metabolism), emissions in influent, or emissions in effluent; WWTP treatment removal efficiency, in percent; instream removal during transport (s −1 ) Source(s): Barber et al. (2022) Ecotoxicological data Toxicological endpoint (

F I G U R E 2
Conceptual diagram of wastewater model parameters and integration with the NHDPlus V2 (U.S. Environmental Protection Agency, 2012a) flow network. Each WWTP outfall is assigned to a NHDPlus V2 stream segment COMID. Percent accumulated wastewater (ACCWW%) and predicted environmental concentration (PEC) values are produced for each COMID. These ACCWW% and PEC values represent accumulated flow from upstream and direct discharges and incorporate instream decay (where applicable for select constituents) and instream dilution. PECs are only produced from municipal effluent flows and may incorporate additional attenuation factors (e.g., removal during WWTP treatment and metabolic removal following human consumption).
USEPA ECHO, or quality-assured USEPA Discharge Monitoring Report data (Williams et al., 2021). Multiple adjacent outfalls belonging to a single permit were combined into a single location and their discharge amounts summarized. Discharge locations were checked for all significant facilities using independent sources (e.g., facility website and satellite imagery) and fine-tuned as necessary. Non-significant WWTPs typically had fewer independent sources to leverage for quality assurance and linkage to a NHDPlus V2 COMID primarily relied on spatial linkages of the discharge coordinates with the nearest stream segment or catchment.

| Hydrologic data
Hydrologic data obtained from NHDPlus V2 consisted of stream segment flow routing and connectivity (hydroseq, used to trace network flows; divergence, used to identify divergent network flows and major and minor paths therein; and startflag and terminalfl, used to identify headwaters and terminal flow features, respectively) and gage-adjusted Enhanced Runoff Method (EROM) streamflow and velocity estimates based on 1971-2000 flows, available for both mean-annual and mean-monthly conditions (U.S. Environmental Protection Agency, 2012a). For most segments, velocity and length attributes provided information for calculating time of travel used to estimate instream decay of wastewater contaminants where applicable. Relative to measured conditions at streamgages, EROM streamflow conditions were important for the PEC and subsequent RQ and RI calculations because the modeling of reactive constituents requires streamflow and time of travel estimates at all upstream reaches in addition to at the reach of interest. Artificial flowlines through reservoirs and lakes do not have velocity attributes in the NHDPlus V2 from which to calculate a time of travel. A method was developed to estimate monthly time of travel for 688 artificial flowlines through reservoirs and lakes in the Potomac River watershed. This method assumed that time of travel is equal to the water volume replacement time and used the mean-annual travel time in days (TOTMA) value provided in NHDPlus V2 to calculate stream segment volume in cubic meters. Monthly travel times were derived from the estimated volume in cubic meters by dividing by the monthly EROM streamflows. In all, 17 reservoir/lake segments did not have a TOTMA value provided in NHDPlus V2. Because these segments were not wastewater impacted, the modeled travel times were set to zero without affecting model results.
Tidal segments and permits discharging to tidal segments were excluded from the model. Non-initialized segments (those with an unknown or unattributed flow direction in NHDPlus V2) could not be used to track and accumulate wastewater flows and were excluded from the model.

| Accumulated wastewater ratios and PECs
Point discharges of effluent leaving facilities were summarized by COMID and routed through the NHDPlus V2 flow network, resulting in ACCWW% and PEC values for each modeled stream segment. The ACCWW% value was calculated as the proportion of accumulated wastewater to total streamflow (Equation 1; Barber, Rapp, et al., 2019;Barber et al., 2022): Modeled ACCWW% values using municipal-plus-industrial effluent allowed examination of total wastewater reuse and occurrence in watershed streams. The ACCWW% was examined at 68 PWSs with surface-water intakes in the Potomac River watershed to assess de facto reuse during mean-annual and quarter-year EROM streamflow conditions. Quarter 3 was used to represent low-flow conditions. Quarter-year evaluation aligns with PWS monitoring for disinfection byproducts and other water-quality constituents to assess compliance based on annual averages (U.S. Environmental Protection Agency, 1998).
Due to difficulties in obtaining chemical data inputs for industrial wastewater effluent emission concentrations, PECs were modeled only for municipal effluent inputs. Municipal wastewater-derived contaminant PECs were modeled using a compilation of literature-derived consumption, fate, and transport data (Table S4; Barber, Rapp, et al., 2019;Barber et al., 2022). Values for PECs were derived as follows: (1) calculation of contaminant loads in micrograms per day (μg/day) leaving each municipal WWTP (Equations 2 and 3; Barber, Rapp, et al., 2019;Barber et al., 2022) and (2) application of dilution and attenuation factors to contaminant loads as they accumulate and travel downstream (Equation 4). Two attenuation factors could be applied to adjust the WWTP load (WWTP load ) where applicable: (1) fraction of unmetabolized drug excreted in urine and feces (F ex ; many pharmaceuticals are administered orally and may be partially or fully metabolized prior to excretion); and (2) fraction of the chemical removed during WWTP treatment (F WWTP ). Loading factors derived from literature data on emissions in influent (Table S4)  (1) ACCWW % = Σ upstream and incoming WWTP discharge EROM streamflow + Σ upstream WWTP discharge × 100. (2) where P (capita) was population served and U (mg/capita/day) was a per capita chemical use loading factor (derived from data on consumption, emissions in influent, or emissions in effluent).

The 2012 Clean Watersheds Needs Survey (CWNS; U.S. Environmental Protection Agency, 2016) and Shenandoah Accumulated Ratio
Mapper Kandel et al., 2019) provided facility-specific data on population served and treatment type for select facilities. Treatment codes could be assigned to each facility as (1) primary (e.g., basic removal of solids, such as a lagoon); (2) secondary (e.g., biological removal of contaminants by activated sludge, oxidation ditch, or rotating biological contactor); or (3) tertiary (e.g., additional removal of inorganic and organic compounds by enhanced nutrient removal or filtration/adsorption). Each treatment code could be assigned a removal rate specific to a given contaminant. For facilities with multiple treatment processes (e.g., secondary and advanced), the more advanced treatment process was assigned. For facilities without treatment or population data, primary treatment was assumed, and the population served (P estimated ) was estimated by standardizing to 378.5 L/capita/day (100 gal/capita/day; Equation 3; Barber, Rapp, et al., 2019): where Q WWTP is WWTP discharge (L/day). This standardization value was obtained from the water-supply plans for the Northern and Upper Shenandoah Planning Districts (Central Shenandoah Planning District Commission, 2011; Shenandoah Valley Regional Commission and Technical Advisory Committee, 2011). This may be a conservative estimate for some urban facilities, as minimum design ranges for WWTPs in U.S. cities have been reported at 946.4-1514.2 L/capita/day (250-400 gal/capita/day; Merritt, 1983). In general, population is a driver of municipal WWTP flow, where discharge increases or decreases in proportion to changes in the community's population.
Following the load calculations, a PEC (in micrograms per liter, μg/L) was estimated (Equation 4; Barber, Rapp, et al., 2019;Barber et al., 2022) after accounting for dilution factors and instream attenuation (Rapaport, 1988): where T is the travel time (s), R stream is the removal rate during stream transport (s −1 ), and Q stream is the streamflow (L/day).
Some of the constituents modeled are naturally occurring and loading from anthropogenic sources such as WWTPs represent only a fraction of background levels (Dyer & Caprara, 1997). Boron and lithium PECs in the Potomac River watershed were corrected for background con-

| Ecological risk assessment
Ecological risk to aquatic organisms was assessed using RQs as previously described (Barber et al., 2022;Donnachie et al., 2014;Johnson et al., 2017). Using RQs allows for the direct comparison of exposure and risk by producing a standardized value across contaminants, where contaminants with a higher RQ present more risk (Equation 5; Barber et al., 2022): Two aquatic toxicity datasets were used to develop RQs for wastewater compounds (Table S4). The first approach used predicted no-effect concentration equivalents (PNEC eq ) normalized from toxicological endpoints for the most sensitive aquatic taxa (fish, amphibian, invertebrate, or plant) using toxicity assessment factors (Barber et al., 2022). Since empirical PNEC data are not always available for many important contaminants of concern, Organization of Economic Co-Operation and Development guidelines (OECD, 1995) were used to extrapolate from a variety of experimentally determined toxicological (or effects) endpoints to a PNEC eq as a common estimator of environmental concern (Equation 6; Barber et al., 2022): where the toxicological (or effects) endpoint is an aquatic concentration in μg/L, and the toxicity assessment factor is a dimensionless numerical adjustment of 1, 10, or 1000 (Table S4). Explanations of criteria for toxicity assessment factors used to normalize more than 10 different toxicological endpoints into a PNEC eq (e.g., PNEC eq = lowest observable effect concentration/10) are described in Barber et al. (2022).
(3) P estimated = Q WWTP 0.264172 × 1, 000, 000 100 , (6) PNEC eq = Toxicological endpoint Toxicity assessment factor , The second approach to develop RQs involved calculating taxa-specific RQs for fish (fathead minnow), invertebrates (daphnid), and plants (green algae) using chronic-toxicity values from the EPI Suite™ Ecological Structure Activity Relationships (ECOSAR) model, without toxicity assessment factors applied (U.S. Environmental Protection Agency, 2012b). In both approaches, individual RQ values for all constituents with toxicity values for which PECs were simulated for a given stream reach were summed to create an aggregated RI, as previously used for assessments of chemical mixtures (Barber et al., 2022;Covert et al., 2020;Nowell et al., 2018). The RI values for each stream segment were used to represent cumulative risk from contaminants across the watershed, assuming an additive effect of the chemical mixtures. Although this approach does not consider different modes-of-action or synergetic effects, Cedergreen (2014) suggests that synergistic effects in typical exposure scenarios are relatively small compared to additive effects. Additive-effect toxicities can yield predictions within a factor of three among heterogenous mixtures (Belden et al., 2007;Faust et al., 2003). Values for individual compound RQs and RI values calculated from PNEC eq (RI PNEC ) fell into one of four categories (Barber et al., 2022): (1) "low risk" (>0 ≤ 0.1); (2) "potential risk" (>0.1 ≤ 1.0); (3) "risk" (>1 ≤ 100); or (4) "high risk" (>100). RI values calculated using the ECOSAR data (RI chronic ) allowed comparison of relative risk from the contaminant mixture among each aquatic taxon but were not categorized by the risk thresholds defined above. Because no toxicity assessment factors were applied, RI chronic values were smaller than RI PNEC values.  (2) the mean-annual PEC calculated using the mean-annual EROM streamflow values (Table S5).

| Comparison of PECs with measured environmental concentrations
Duplicate samples reported in more than one database and samples collected at facility outfalls were removed from the data compilation. Because the approach presented here allows the prediction of contaminant concentrations below the method detection limit ( value was used to scale the results to a range of (0 to 1) where nNSE = 1 indicates ideal model performance and nNSE ≤0.5 suggests that the observed mean is a better predictor than the model simulation. The sensitivity of the NSE parameter to large dependent variable outlier values (PECs) has been well documented (Krause et al., 2005). The RMSE is the standard deviation of the errors between observed and predicted concentrations and nRMSE compares datasets with different scales (Golmohammadi et al., 2014), where a value of zero suggests optimal model performance and large positive values suggest poor model performance. Both streamflow-normalized and raw concentration paired PEC/MEC datasets were analyzed for stream segments with predicted municipal wastewater influence.
Compounds with PECs close to or lower than the MDL are less suitable for the PEC/MEC comparison. It is important to identify compounds where MDL may interfere with the availability of MEC data and therefore limit evaluation of model performance. A measure of the percent occurrence of PEC values above the MDL (%PEC above MDL ; Equation 7) was used to assess the relative amount of detection-limit interference in MEC data availability: where PEC above MDL is the number of COMIDs with PECs greater than the MDL. High values (70 to 100%) indicate that most of the stream reaches with municipal wastewater contained simulated concentrations that were above MDLs. Low values (0 to 10%) indicate that the majority of stream reaches containing municipal wastewater had PEC values that were below MDLs, which may limit the availability of MEC data for comparisons.
Temporally, ACCWW% values followed a bimodal seasonal hydrograph wherein the highest ACCWW% values occurred in the summer months (lowest streamflow estimates) and the lowest ACCWW% values occurred in the winter months (highest streamflow estimates) (Figure 4). Mean-August EROM conditions typically represented the lowest streamflow month (and thus the month with the highest ACCWW%) and mean-March EROM conditions typically represented the highest streamflow month (lowest ACCWW%). During mean-August conditions, 41% of wastewaterimpacted streams exceeded the mean-annual third-quartile value of 2.3% ACCWW% (Figure 3d). Below the confluence of the Monocacy and Potomac Rivers north of Washington, D.C., ACCWW% was 1.3% during mean-March conditions and increased to 5.9% during mean-August conditions. The ACCWW% also was evaluated at streamgages using measured streamflows to illustrate the range of potential combinations of streamflow and wastewater conditions that can occur ( Figure S1). Streams had a larger range of ACCWW% values as they decreased in stream order, with some experiencing >30% ACCWW% both during mean-August low-flow EROM conditions ( Figure 4) and historical low-flow conditions (represented by high-flow percentiles).
Of all modeled stream kilometers in the Potomac River watershed, first-order streams comprise 63%, second-order 17%, third-order 9.0%, fourth-order 7.0%, fifth-order 3.0%, sixth-order 1.0%, and seventh-order <1.0%. The first-order streams received the majority of direct wastewater discharges. Of the 1837 outfalls, 71% discharged wastewater to first-and second-order streams and 27% discharged to third-, fourth-, and fifth-order streams. Despite receiving only 2.0% of direct discharges, sixth-and seventh-order streams had the highest mean ACCWW% values (6.8% and 6.1% during mean-August, respectively) relative to other reaches (mean-August: first-order, 5.3%; second-order, 3.9%; third-order, 3.6%; fourth-order, 4.6%; fifth-order, 5.3%), and all sixth-and seventh-order streams considered in the model were wastewater impacted. Relative to other streams, there was a smaller range of ACCWW% in sixth-and seventh-order streams, which did not exceed 21% under low-flow EROM conditions and did not exceed 50%, even under the lowest historic conditions at streamgages, except for site number 01643000 (Monocacy River at Jug Bridge near Frederick, Maryland in HUC 0207009).

| Wastewater reuse
The 68 PWSs ( Figure 5; Table S3) considered collectively serve roughly 4.5 million people (75% of the Potomac River watershed population) and draw their source waters from a range of stream sizes (first-order, 15 intakes; second-order, 4 intakes; third-order, 10 intakes; fourth-order. 13 intakes; fifth-order, 9 intakes; sixth-order, 11 intakes; seventh-order, 6 intakes). In all, 65 PWSs are community water systems, two are non-transient non-community water systems, and one is a transient non-community water system. In all, 46 (68%) intakes were impacted by industrial or municipal wastewater. At the PWS intakes, mean-annual ACCWW% values ranged from <0.01% to 16% (Figure 5a). The majority (7) % PEC above MDL = PEC above MDL Total number of wastewater-impacted reaches × 100, (89%) of PWSs drawing from low-order streams had no de facto reuse (ACCWW% values of zero). PWSs drawing from high-and medium-order streams generally had the highest levels of de facto reuse (elevated ACCWW%). For each intake, median ACCWW% values from the monthly estimates within each quarter were used to obtain the quarter-year values ( Figure 5b; Table S3). Median ACCWW% during Quarter 3 was much larger than the other quarters, with a maximum of 31% at one intake (PWS JJ in HUC 02070010, Middle Potomac-Anacostia-Occoquan).

| Predicted and measured environmental concentrations
Over a quarter (27%) of Potomac River watershed stream segments had upstream or direct municipal effluent discharges contributing wastewater-derived contaminants. In these streams, PECs were calculated for 69 contaminants, consisting primarily of pharmaceuticals (n = 37; 54% of the total), using previously compiled per capita consumption and environmental fate data (Table S4; Barber et al., 2022).
The remaining contaminants generally were categorized as perfluorinated alkyl substances (6), consumer product chemicals (6) of their relation to municipal WWTP discharge, PECs and ACCWW% generally were positively correlated ( Figure S2A). However, PECs may be subject to attenuation factors not considered in the ACCWW% calculation that can affect relations between these two values among individual constituents. For example, estriol (naturally occurring hormone) and oxycodone (opioid analgesic pharmaceutical) have different WWTP removal and instream decay rates (estriol has high instream decay and wastewater treatment removal rates relative to oxycodone, Figure S2B).

F I G U R E 4
Boxplots showing the interquartile range, median, and outliers of municipal-plus-industrial ACCWW% across all wastewaterimpacted stream segments in the Potomac River watershed as a function of month and Strahler stream order. The total number (n) of municipal-plus-industrial wastewater-impacted segments are listed. The dashed lines represent 1% ACCWW%; higher-order streams are less sensitive to flow variation, resulting in a smaller range of ACCWW% values relative to lower stream orders. Mean ACCWW% in lower stream orders, however, is generally smaller relative to higher stream orders (Streamflow conditions represent mean monthly gage-adjusted EROM streamflow from 1971 to 2000, U.S. Environmental Protection Agency, 2012a). n = 30) and metals/metalloids (log-transformed; ρ = 0.89, r 2 = 0.74, n = 97). These classes generally followed a quasi-linear relationship between observed and mean-monthly predicted values. The dataset without streamflow normalization showed similar trends (Table S6; Figure S3).
Compounds classified as being derived from both consumer and industrial mixed-use sources (bisphenol A, EDTA, nonylphenol), plasticizers/fire retardants (tributyl phosphate) and hormones (17-alpha-ethinylestradiol, estrone) frequently had PEC values lower than detection limits (low %PEC above MDL ) and thus lower in magnitude than streamflow-normalized MEC values, for which values below the method detection limit were excluded (Table S6; Figure 6b). The MEC/PEC pairs that were not streamflow-normalized highlighted additional statistically relevant relationships with the increased sample size for compounds that did not have associated streamflow data ( Figure S3). Constituents with low availability of streamflow measurements with category 2 or 3 designations included pharmaceuticals (anhydroerythromycin, metformin, ranitidine) and hormones (17-alpha-ethinylestradiol, estriol).

| Ecological risk
Potential ecological risk of the contaminant mixtures on aquatic organism health was assessed using previously published ecotoxicity datasets (Barber et al., 2022; Table S4). The RI PNEC values produced using the literature-based approach involved converting empirically derived data with various toxicity endpoints (e.g., PNEC, predicted no-effect concentration; NOAEC, no observable adverse effect concentration; EC50, concentration having effect on 50% of population; LC50, concentration lethal to 50% of population; and others represented in  Table S4). The ECOSAR model uses log octanol-water partition coefficients and semi-quantitative structure activity relationships to estimate chronic-toxicity benchmarks for organic chemicals when empirical toxicity data are not available. The RI values in Table S7 for each stream reach COMID were calculated using both approaches.
Of the Potomac River watershed stream segments receiving municipal effluent, 1479 km (22%) had mean-annual RI PNEC values >100, indicating "high risk" of potential effects for the most sensitive aquatic species (Figure 7a,b). These streams represent 6% of the total stream kilometers in the Potomac River watershed. Over a third (598 km; 36%) of municipal wastewater-impacted first-order streams were in the "high-risk" category under mean-annual EROM streamflow conditions. The percentage of "high-risk" stream kilometers decreased with ascending stream order (Figure 7a). The reverse trend occurred in the "risk" category (RI PNEC > 1 ≤ 100), which comprised 4849 km (72%) of impacted Potomac River watershed stream segments.
Temporal variation was largest within the "high-risk" and "risk" categories in small-to medium-order streams during mean-August lowflow and mean-March high-flow EROM streamflow conditions. Although 752 km (46%) of impacted first-order streams were "high risk" during mean-August conditions, only 217 (13%) were "high risk" during mean-March conditions. During mean-March conditions, these streams shifted to the "risk" category from "high risk" under mean-August conditions. Sixth-order stream kilometers (consisting of the Cacapon, North Branch Potomac, Middle Potomac, Shenandoah, and Monocacy Rivers) in the "high-risk" category ranged from 35 (9%) to 0 under mean-August and mean-March conditions, respectively.
The only seventh-order stream in the Potomac River watershed is the Middle Potomac-Catoctin and it is the integration point of streamflow from throughout the watershed just before the Potomac River empties into the Chesapeake Bay (HUC 02070008). The entire length of the Middle Potomac-Catoctin was in the "risk" category (RI PNEC > 1 ≤ 100) throughout the year regardless of streamflow condition. Whereas a subset of smaller streams had RI PNEC ≤ 1 ("low risk" or "potential risk"), model simulations suggested that all sixth-and seventh-order streams carried wastewater contaminants year-round at concentrations sufficient to present "risk" or "high risk" to aquatic organisms. At the 8-digit HUC level under mean-annual EROM streamflow conditions, >90% of impacted stream kilometers had RI PNEC > 1 representing "risk" or "high risk" (Table 2; Figure 7b). The Monocacy River (HUC 02070009) had the largest percentage of impacted stream kilometers in the "high risk" category (46%) followed by the Shenandoah River (HUC 02070007; 40%). The RI chronic results indicated that the relative order of susceptibility to chronic-exposure effects was fishes > invertebrates > plants (Figure 8a), a pattern that was consistent regardless of stream order (Figure 8b).

| Wastewater reuse
The use of treated wastewater for potable applications is important for sustainable water reuse programs, particularly where growing populations and climate change place stress on limited water resources (National Research Council, 2012). Identifying the extent and significance of de facto reuse, however, continues to be a national priority (U.S. Environmental Protection Agency, 2020a). In contrast to planned reuse paradigms, de facto (unintended) reuse may not be officially recognized or accounted for, with unknown consequences on consumers.
Potential risks include exposure to pathogens, such as Giardia, and other effluent-derived contaminants, such as pharmaceuticals. De facto wastewater reuse also may introduce elevated levels and types of precursors of chlorinated DBPs like HAAs and THMs, which have been associated with increased bladder cancer risk in epidemiological studies (Villanueva et al., 2004(Villanueva et al., , 2007Weisman et al., 2019Weisman et al., , 2021. Relative to national studies (Rice & Westerhoff, 2015), PWSs in the Potomac River watershed are more frequently impacted by upstream wastewater discharges (70% compared to 50% nationally) and were more likely to have de facto reuse above the 1% threshold (40% compared to 25% nationally). Of note, the national model used by Rice and Westerhoff (2015)

| Predicted and measured environmental concentrations
Surface waters contain natural constituents from geologic materials as well as anthropogenic contaminants, such as pharmaceuticals and personal-care products. Major anthropogenic sources of contaminants in streams include continuous discharges from municipal and industrial WWTPs (National Research Council, 2012;Rice & Westerhoff, 2015), baseflow groundwater discharge (Schilling & Wolter, 2001), and runoff affected by agricultural practices (Stone et al., 2014). Additional contaminants may be introduced to receiving waters from other non-point sources, such as urban stormwater runoff (Masoner et al., 2019), and point sources, such as landfill leachate and combined-sewer overflows (Barnes et al., 2004;Weyrauch et al., 2010).
The PEC calculations in this study reflected contributions of contaminants from municipal WWTPs only and were controlled by modeled streamflow conditions and the input parameters assigned to calculate facility loads and instream removal rates. Because municipal WWTPs represent only one pathway of contaminants to the environment, agreement between MECs and PECs was not expected for all compounds.
Aggregated comparisons across general contaminant classes indicated that the best model performance was for pharmaceuticals and consumer product chemicals relative to other constituent types. Excretion of human pharmaceuticals via urine and feces, and subsequent discharge from municipal WWTPs in treated wastewater, is the primary loading source of these contaminants (Winker et al., 2008). Similarly, consumer product chemicals generally are released with household wastewater following everyday use, such as cleaning and showering, and have limited alternative pathways into the environment.
Predicted environmental concentrations for most compounds were lower than MECs, which may reflect contaminant-loading sources separate from municipal wastewater (e.g., industrial effluent, landscape inputs, and agricultural runoff). MDL interferences, inaccurate massloading estimates from municipal outfalls, or other contaminant sources can result in data points plotting below the 1:1 correlation line ( Figure 6a; Figure S3A). Low correlation can occur for a variety of reasons, including unsuitable first-order decay coefficients, additional contaminant sinks or sources in the watershed, errors in sample-collection protocols, laboratory analytical errors, MDL limitations, and reporting or database inaccuracies. High correlation with low predictive accuracy can occur when many or all of those conditions are satisfied, but predictions of mass entering the system are too low, such as from unsuitable estimates of per capita use or overestimates of removal occurring during treatment. The predictions herein were made using the best available chemical data compiled from multiple studies across different regions and leveraging different methods. These estimates have inherent uncertainty and can vary both temporally and spatially.
Highly correlated PEC/MEC pairs were observed for streamflow-normalized contaminants such as N,N-diethyl-meta-toluamide, trimethoprim, tris-(2-chloroethyl) phosphate, and triphenyl phosphate, indicating that model input parameters (hydrology, dilution, and first-order decay rate coefficients) adequately described many constituent sources, transport, and fate within the watershed.
Observed differences in streamflow-normalized hormone PECs and MECs may reflect inputs from livestock and agricultural fields treated with manure (Gall et al., 2011;Lange et al., 2002). Most hormone MECs included in the streamflow-normalized dataset were compiled from studies in the Shenandoah River watershed, where agricultural land use is prevalent . The hormones frequently occurred at environmental concentrations below the MDL (mean %PEC above MDL of 6.9%), which complicates data interpretation. The availability TA B L E 2 Percent of municipal wastewater-impacted Potomac River watershed 8-digit HUC stream kilometers (km) falling within each predicted no-effect concentration equivalent risk index (RI PNEC eq ) category for high-flow (mean-March), mean-annual, and low-flow (mean-August) EROM streamflow conditions (U.S. Environmental Protection Agency, 2012a).

| Wastewater condition variability and ecological risk
Aquatic organisms typically are exposed to instream contaminants in the form of complex mixtures that are potentially toxic to multiple taxa through various modes-of-action (Brain et al., 2008;Fent et al., 2006). National evaluations of surface waters across the United States by Kolpin et al. (2002) detected a median of seven organic wastewater contaminants in a given water sample, with as many as 28 detectable compounds at some sites. A follow-up study by Bradley et al. (2017) further characterized contaminant-mixture complexity in 38 streams across the nation and detected complex surface-water organic-contaminant mixtures even in undeveloped and uninhabited watersheds. In the Shenandoah River watershed, a subbasin of the Potomac River watershed, 273 inorganic and organic constituents were detected in water samples collected between 2014 and 2016 (Barber et al., 2022). The potential impact of complex chemical mixtures is not well characterized, particularly at the ecosystem level. At the organismal level, exposure to wastewater has been associated with complex effects that include endocrine disruption, reduced fecundity, mouthpart deformities, immunosuppression, increased mortality, and phytotoxicity Brain et al., 2004;Ford & LeBlanc, 2020;Iwanowicz et al., 2009;Vajda et al., 2008;Watts et al., 2003). Application of USEPA EPI Suite™ ECOSAR (U.S. Environmental Protection Agency, 2012b) chronic-toxicity data suggested more potential risk from chronic exposure to chemical mixtures simulated in this study in fishes than in invertebrates and plants. This is consistent with the evolutionary distance between each taxon to humans, for whom many of the active ingredients in pharmaceuticals and consumer products are designed (Gunnarsson et al., 2008).
The ACCWW%, PECs, RQ, and RI are affected by instream dilution, which is a function of stream order and temporal changes in streamflow. As observed in other studies, a subset of low-order streams may experience higher levels (and greater fluctuations) of wastewater relative to high-order streams due to streamflow volumes that generally are lower (Rice et al., 2013;Rice & Westerhoff, 2015). Chemical mixtures that may present "high risk" to sensitive aquatic organisms were observed more frequently in small streams impacted by municipal wastewater regardless of streamflow condition. All high-order streams were wastewater-impacted and carried contaminants at levels sufficient to present "risk" or "high risk" throughout the year, resulting from a combination of direct wastewater discharges and contributions from wastewaterimpacted upstream tributaries.
Low-order streams are critical hotspots of aquatic biodiversity and habitat variability (Biggs et al., 2017;Freeman et al., 2007;Meyer et al., 2007), while medium-order streams have been shown to support rich fish assemblages (Rapp et al., 2020). In the Potomac River watershed, a combined 5387 km of first-, second-, and third-order streams were impacted by direct or upstream discharges of municipal and/or industrial wastewater (23% of combined first-, second-, and third-order stream kilometers). A combined 2433 km of fourth-and fifth-order streams were impacted (92% of the total). As observed in other studies, lower-order streams had higher streamflow variability and lower dilution factors (Rice et al., 2013;Rice & Westerhoff, 2015), placing the diverse habitats and fauna in these reaches at the highest risk of negative impacts. Increasing wastewater percentages and instream contaminant concentrations in early to late summer may coincide with early life stages for many species, which may be more sensitive to toxicants relative to the adult organism (Hutchinson et al., 1998). Wastewater also may combine with alternative inputs of contaminants to worsen negative effects. For instance, aquatic organism exposure to agricultural runoff generally is greatest in the spring (Smalling et al., 2021;Thurman et al., 1991). In smallmouth bass, it has been hypothesized that risk of endocrine disruption may be greatest following the spring flush of agricultural chemicals, but such effects may increase in severity when individuals subsequently are exposed to additional contaminants (Kolpin et al., 2013). This may create a scenario where increased wastewater risk in the summer threatens already compromised organisms.

| Tidal streams
Simulation of transport times, mixing dynamics, and salting-out effects of contaminants in tidal reaches requires multi-dimensional consideration of variation in water movements, wind, and baroclinic forcing outside the scope of the current model (Gao et al., 2020;Shen et al., 2016).
The results presented in this article are for the non-tidal portion of the watershed and exclude the 1675 NHDPlus V2 tidal stream segments (38% of which are wastewater-impacted) beginning near Washington, D.C. The tidal segments include 64 industrial and 26 municipal WWTPs ( Figure 9) discharging a combined 16 m 3 /s of treated effluent (99% from municipal WWTPs). The total estimated mass load of the 69 contaminants entering tidal waters from the 26 municipal WWTPs was 694 kg/day (Table S8). Contaminant loading from WWTPs to the tidal Potomac likely will persist and may intensify due to expanding metropolitan populations and industry. Consequently, hydrodynamic modeling techniques to improve understanding of interfaces between freshwater and estuary systems and to predict the impacts of tidal loading and chemical transport phenomena warrant further exploration. Additional uncertainties in the model input parameters include WWTP information (e.g., effluent flows, population served, treatment processes, and discharge location), data used to calculate PECs (e.g., municipal effluent emissions and human consumption), and streamflow.

| Methodological considerations
The model assumed constant wastewater discharge rates over time, when wastewater flows can vary based on factors such as water use, industrial activities, and precipitation (Merritt, 1983). Per capita chemical use data were compiled from literature sources spanning a range of time periods and regions. The ACCWW% and PEC values were based on mean streamflow from long periods of record (1971 to 2000), but streamflow variations occur based on multiple environmental and anthropogenic factors (e.g., precipitation, water withdrawals, and vegetation; Mays, 2011). Flow-exceedance percentiles at streamgages illustrate the variability in streamflow conditions that can occur. The advantage of using mean EROM streamflows compared to measured streamflow at streamgages, is that EROM streamflows allow consideration of stream dilution and time of travel at each ungaged reach in the watershed. Refinements to this approach have been demonstrated in the Shenandoah River watershed using local streamgage information to estimate ACCWW% for measured historic flows at ungaged reaches (Weisman et al., 2021).
Limitations to the literature-based approach to calculate RI PNEC values include range of taxa and species assessed, data availability (where some constituents are not well studied), and the range of biological endpoints and reporting units. Limitations to using the empirical ECOSAR toxicological endpoints to calculate RI chronic values are that it applies only to organic compounds and there are uncertainties associated with the estimated values (Card et al., 2017). Computational approaches, however, are common in environmental toxicity studies and provide a method of risk screening where empirical data are sparse (Card et al., 2017;Sanderson et al., 2004). In general, risk assessments require multiple assumptions and have inherent uncertainty. Thus, results should be considered relative rather than absolute.

| CON CLUS ION
This study examined variations in simulated industrial-plus-municipal wastewater conditions, evaluated de facto reuse at PWS drinking-water intake locations, and generated PECs for select municipal wastewater-derived contaminants in Potomac River watershed streams. The PECs were evaluated against ecotoxicity benchmarks and aggregated to estimate cumulative-exposure risk conditions for aquatic life across the watershed. This study expands on previous wastewater reuse and risk-ranking approaches (Barber et al., 2022;Barber, Rapp, et al., 2019) and can be adapted for different geographic areas or contaminants to (1) evaluate potential impacts of constituents with concentrations below detection limits in surface waters; (2) identify areas with high percentages of accumulated wastewater or elevated predicted ecological risk; and (3) spatially integrate widely available hydrologic, wastewater discharge, and contaminant data.
In the Potomac River watershed, reports of declines in fish health (e.g., intersex fish and fish mortalities; Blazer et al., 2010) have raised concerns about aquatic populations and drinking-water quality. The screening-level approach for assessing risk to aquatic organisms (stream water) and humans (drinking water) presented in this article suggests that wastewater may be an important factor contributing to environmental degradation of Potomac River watershed streams. Modeled de facto reuse at over a quarter of the PWS intakes examined exceeded levels previously associated with the formation of halogenated DBPs in some systems (Weisman et al., 2019(Weisman et al., , 2021. Assessment of potential toxicity from a chemical mixture consisting of 51 contaminants indicated that 22% of the reaches (1479 km) receiving municipal wastewater (5.5% of all reaches modeled) exceeded "high-risk" thresholds of potential ecological risk for aquatic organisms. Municipal wastewater-based chemical mixtures may result in more negative effects from chronic exposure in fishes relative to invertebrates and plants.
The results presented for the Potomac River watershed may be beneficial as a screening-level risk assessment for resource managers for developing comprehensive water-sampling plans targeting stream reaches with a high likelihood of containing biologically active, toxic compounds contributing to increased risk to aquatic organisms. These results also could inform water-supply permitting decisions in medium-to high-order streams, which could consider the potential increased risk levels resulting from increased total water withdrawals. In addition, results could be used to target WWTP upgrades to reduce the discharge of certain anthropogenic contaminants to improve downstream water quality (Barber et al., 2012. Actual risk will vary temporally and spatially with changes in streamflow and wastewater flow. Future risk modeling efforts could be improved through (1) targeted sampling validation efforts to create a consistent MEC dataset for comparison with predictions and (2) the integration of other sources contributing to riverine contaminant loading, which is impeded by data gaps such as accurate and current chemical input parameters and treatment efficiencies for industrial effluent-derived chemicals.

This research was supported by the U.S. Geological Survey Environmental Health Program (Toxic Substances Hydrology and Contaminant
Biology). We thank Karen Rice of the U.S. Geological Survey and Jacelyn Rice-Boayue of the University of North Carolina, Charlotte for their reviews, which strengthened this report. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. government.

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors declare no conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available in ScienceBase at https://doi.org/10.5066/P9CB2YM7.