Bull Run flows through four regional parks: Manassas National Battlefield Park (20 km²), Johnny Moore Stream Valley Park (less than 0.005 km²), Bull Run Regional Park (6 km²), and Hemlock Overlook Regional Park (1.6 km²) (Manassas, 2021; Hemlock Overlook Regional Park, 2022; Johnny Moore Stream Valley Park, 2022; Northern Virginia Regional Park Authority, 2022). Synoptic sampling began at Bull Run Regional Park and Hemlock Overlook Regional Park. The confluence of Bull Run and a rapidly urbanizing tributary, Cub Run, occurs within Bull Run Regional Park.

The Upper Occoquan Service Authority (UOSA) operates a water reclamation facility which discharges 54 million gallons of treated wastewater effluent per day into Bull Run within Bull Run Regional Park (Onyango et al., 2014). The discharge flows to the Occoquan Reservoir and is part of an indirect potable reuse (IPR) system (Cubas et al., 2014). The UOSA treatment process includes a conventional treatment, chemical and physical advanced treatment, disinfection, and digestion and sludge handling before discharging to Bull Run (Upper Occoquan Service Authority, 2022). The wastewater treatment process can remove biological nitrogen, therefore discharging reclaimed water with lower phosphorus and organic matter concentrations (Cubas et al., 2014). Nitrified effluent is discharged to Bull Run in the warmer months to combat the start of anoxic conditions in the reservoir; in the winter months when the reservoir is oxygenated, UOSA must remove nitrogen in the wastewater and follow regulations for point source nitrogen loads (Cubas et al., 2014). UOSA is the only wastewater treatment plant within the studied watersheds.

Four longitudinal synoptic surveys were conducted along the length of Bull Run. Each synoptic included portions of Bull Run above and below the location where the treated wastewater is discharged from UOSA. The first two synoptic sampling events were conducted during baseflow conditions and extended from Bull Run Regional Park to the Occoquan Reservoir. The last two synoptics were sampled during snow events and at a high-frequency along the flowpath (Table 2, Supplementary Figure S3); we were not able to sample the full length of the previous synoptics due to safety concerns.

Rock Creek flows through Rock Creek Regional Park in Maryland and directly into Rock Creek National Park in Washington, D.C. Rock Creek Regional Park, at about 7.3 km², is managed by the Maryland-National Capital Park and Planning Commission (Montgomery Parks, 2022). Downstream of the regional park, the 7.1 km² Rock Creek National Park in Washington D.C. was established in 1890 by the National Parks Service (National Parks Service, 2022). While the river flows through two geographical regions, about 80% of Rock Creek’s drainage area is in Maryland and the remaining 20% in D.C. (Cintron, 2016). Rock Creek has been the focus of numerous water quality studies, addressing the contribution of the stormwater drain runoff and the 217 combined sewer outfalls to pathogen, metal, and organic matter pollution in the stream (Miller et al., 2013; Cintron, 2016).

Five longitudinal synoptic surveys were conducted along Rock Creek, each with varying distances downstream. All synoptics began in Rock Creek Regional Park in Maryland and ended in Rock Creek National Park in Washington, D.C. (Table 2). Road salt event 1 was sampled 4 days after a snow event and the second road salt event was collected 3 days after a smaller snowstorm (Supplementary Figure S3).

Scotts Level Branch is a suburban stream with a narrow riparian buffer (Wood et al., 2022). Two floodplain reconnection projects were completed on Scotts Level Branch (Scotts Level Branch Stream Restoration Project, 2019). The first restoration project at McDonogh Road was completed in the spring of 2014. Around 2,000 linear feet of the stream was restored, starting 2.6 km downstream of the headwaters (Scotts Level Branch Stream Restoration Project, 2019). The second project began in the summer of 2019 at the headwaters, a large storm drain, near Marriottsville Road and extended 1,500 feet to Tiverton Road (Scotts Level Branch Stream Restoration Project, 2019). Unlike at Rock Creek and Bull Run, sampling at Scotts Level Branch was conducted at the same location during each synoptic. Scotts Level Branch was the shortest stream length with restricted access due to the residential areas, therefore, we could only access certain locations. Samples were collected from the headwaters of Scotts Level Branch originating at a stormwater drain to the confluence of Scotts Level Branch and the Gwynns Falls.

The relationship between specific conductance and multiple ion mobilization suggests that chemical cocktails associated with salinization are formed in urban streams (e.g., Kaushal et al., 2018a; Kaushal et al., 2019; Morel et al., 2020; Galella et al., 2021). PCA was used to diagnose the composition of the chemical cocktails which formed seasonally over an annual monitoring cycle (Figures 2, 3). Dimension 1 of the PCA (PC1) (Figure 4, see Supplementary Table S2 for the discriminant matrix) captures patterns in salt ions and trace metals (∼40% variance explained), whereas dimension 2 of the PCA (PC2) captures patterns in redox sensitive elements (∼21% of the variance explained). Water samples collected during road salt pulses are characterized by elevated concentrations of Na⁺, K⁺, Cu, and Mn (see pink circles in the top right quadrant of Figure 4A; small symbols indicate the individual samples collected, large symbols indicate their central tendency). Most of these ions, except for Cu, have a positive significant relationship with increasing specific conductance (Supplementary Table S2). Water samples collected during summer baseflow and stormwater events had lower concentrations of Na⁺, K⁺, Cu, and Mn (see blue symbols loading opposite the road salt cluster, Figure 4A).

A second chemical cocktail can be seen in the lower right quadrant of Figure 4A, characterized by elevated concentrations of Ca²⁺, Mg²⁺, Sr²⁺, and TDN. These elements had a positive significant relationship with increasing impervious surface cover (Supplementary Table S2). Stormflow and baseflow samples collected in the fall (brown symbols), spring (green symbols), and summer (blue symbols) had variable concentrations of constituents in these chemical cocktails (note the variability in the locations of samples from these groups with respect to PC1). Only baseflow and stormflow samples collected in the winter exhibited chemical cocktails that were similar to those evident during road salt events. Notably, you can trace a cycle from one chemical cocktail towards the next by following the seasons, starting with road salt in the upper right quadrant, shifting down and left (away from chemical cocktail 1 as winter transitions to spring and then to summer), relaxing towards chemical cocktail 2 (lower right quadrant) transitioning into fall. The stark differences between chemical cocktails by season illustrate the importance of factoring seasonal changes into our understanding of the mobilization of chemical cocktails by freshwater salinization. Moreover, the relationship between different chemical cocktails and variables like land cover (imperviousness) or management practices (road salt application) highlight the importance of these factors of seasonal change.

Chemical cocktails formed longitudinally following different pollution pulse events. This is evident in Figure 4B, where water samples from the longitudinal sampling events at Bull Run, Rock Creek, and Scotts Level Branch were grouped based on the particular types of pollutant pulses. Samples that contained more Na⁺, Ba²⁺, Mg²⁺, Cu, Sr²⁺, Ca²⁺, and K⁺ tended to be associated with road salt (red and blue symbols) or wastewater treatment (purple symbols; positive PC1) (sensu Amrhein et al., 1992; Acosta et al., 2011; Galella et al., 2021). Samples that contained less Na⁺, Ba²⁺, Mg²⁺, Cu, Sr²⁺, Ca²⁺, and K⁺ tended to be associated with stream baseflow conditions (pink and green symbols) or stormwater runoff (yellow symbols; negative PC1). Relative to road salt pulses, wastewater pulses were less enriched in Fe and more enriched with Ca²⁺, Sr²⁺, and K⁺ (purple symbols shifted down towards negative PC2). Summer baseflow samples were enriched in Fe relative to winter baseflow samples, perhaps because Fe was mobilized in the summer during low oxygen conditions (note absence of pink symbols from positive PC2).

In Bull Run, different road salting events along the flowpath exhibit different distinct chemical cocktails. Both are elevated in salt ions, however, road salt event 1 was more influenced by Mg²⁺ and Sr²⁺ (Mg²⁺ mobilizes Sr²⁺) and road salt two was influenced by Na⁺, Cu, and Ba²⁺ (Na⁺ mobilizes Cu) (Galella et al., 2021).

Water samples collected near and upstream to the wastewater treatment plant have a chemical cocktail characterized by high Ca²⁺, K⁺, and Sr²⁺ and low Mn and Fe (blue symbols; positive PC2, Figure 5A). While wastewater effluent can act as a source of Na⁺ (Bhide et al., 2021), individual road salt events increase Na⁺ concentrations in the stream more than wastewater does. This chemical signature aligns with the longitudinal trends in Figure 6—at the wastewater treatment plant, the in-stream concentrations of K⁺ increased in the stream by an average of 228% and the in-stream Ca²⁺ concentrations increased by an average of 44%. This makes Ca²⁺, K⁺, and Sr²⁺ the chemical signature of ions from the wastewater treatment effluent in Bull Run (Jalali et al., 2008). The baseflow samples had low concentrations of all salt ions and trace metals measured. The summer baseflow samples had concentrations that were influenced by the wastewater treatment plant.

A stronger divergence between the two road salt events is evident along the flowpath from Rock Creek Regional Park to Rock Creek National Park. The chemical cocktails from road salt event 1 had high concentrations of salt ions, particularly Na⁺, Ca²⁺, and Mg²⁺, whereas road salt event 2 had higher concentrations of Cu (Figure 5B). These differences are likely due to the conditions during the road salting event. During that snow event, the stream’s conductivity at the USGS gauge in Rock Creek National Park peaked around 11,700 μS/cm with a discharge of 1.9 ft³/s (Supplementary Figure S3) (peak specific conductance of 6,230 μS/cm, discharge of 10.5 ft³/s), which may have transmitted salt ions (K⁺, Na⁺, Ca²⁺, Mg²⁺) out of the watershed more quickly. During baseflow, Rock Creek exhibited stormwater and baseflow samples with elevated ion concentrations, in particular Fe (Figure 5B).

The composition of the chemical cocktails formed during baseflow and storms depended on the impervious surface cover in the watershed. Scotts Level Branch (Figure 5C) has the highest impervious surface cover as well as the smallest watershed area. The stormwater and baseflow events grouped on the biplot by the seasons, though all were elevated in salt ions. In contrast to the other two streams, both road salt events at Scotts Level Branch were associated with similar chemical cocktails, containing high concentrations of a highly correlated chemical cocktail comprised of Cu, Na⁺, K⁺, Ca²⁺, Mn, and Sr²⁺ (Figure 6C). The summer stormwater and summer baseflow events were correlated with chemical cocktails heavily enriched in Fe. During the storm events in the spring and the fall, the chemical cocktail was enriched in both salt ions and Fe.

Wastewater effluent transiently increased concentrations of Ca²⁺, K⁺, and Na⁺ in spring and summer baseflow at Bull Run. Concentrations of these ions subsequently declined along Bull Run’s flowpath, the majority of which is located within a riparian buffer zone that ultimately feeds into Hemlock Overlook Regional Park (green region, Figure 6A). Longitudinal patterns in Mg²⁺ concentrations were somewhat different (i.e., they were lower in wastewater effluent than the uppermost regions of the riparian buffer zone, and only began to decline approximately 3.5 km downstream after Little Rocky Run; see dashed line at ∼3.5 km, Figure 6A). The rapid decline observed across all salt ions prior to Hemlock Overlook Regional Park reflects the proximity of a large tributary, Popes Head Run (see leftmost green dashed line, Figure 6A). Approximately 86% of the Popes Head Run drainage area was classified as residential-conservation in 1982 to preserve water quality in the Occoquan Reservoir (Fairfax County, 2005).

For both road salt events, Na⁺ and Mg²⁺ peaked higher up in the watershed (i.e., at the confluence of Cub Run and Bull Run; leftmost dashed line, Figure 6A). Concentrations of these salt ions declined along the flowpath to the wastewater treatment plant, and then increased somewhat (or stayed relatively stable) thereafter. A similar pattern occurred for K⁺ during the second road salt event, but not the first. During the first road salt event for K⁺ and both road salt events for Ca²⁺, concentrations were highest at the wastewater treatment plant and subsequently stayed stable or declined (i.e., the longitudinal pattern observed was comparable to baseflow conditions).

The dominant pattern in chemical cocktails for summer and winter baseflow at Bull Run (i.e., PC1 from Figure 5A; positive PC1 indicates higher salt ion concentrations) was well explained by an MLR model including riparian buffer width (ion concentrations were lower when riparian buffer width was higher; see partial effects plot, Figure 7A) and sampling event (ion concentrations were higher in summer baseflow samples than winter baseflow samples; see partial effects plot, Figure 7B). The fitted model explained 66% of the observed variance in PC1 at Bull Run. Approximately two-thirds of this explanatory power was attributable to sampling event (45 of 66%; averaging over ordering, Supplementary Table S4), suggesting that observed chemical cocktails may differ more between summer and winter sampling events (i.e., over time) than along the Bull Run flowpath for a given sampling event. The remainder was attributable to riparian buffer width (21 of 66% variance explained, Supplementary Table S4). Attenuation of chemical cocktails with increasing buffer width was more pronounced in summer than winter baseflow samples (note elevated scatter in winter baseflow samples (celadon points) relative to summer baseflow samples (orange points) in Figure 8A).

In Rock Creek, salt ion concentrations during summer baseflow events tended to increase along the flowpath through the regional park (white region, upper portion of the synoptic) and then decrease through the national park itself (green region, lower portion of the synoptic; Figure 7B). During stormflow and road salt events (purple and blue lines, Figure 6B), ion concentrations tended to increase along the entire flowpath. This was most evident for Ca²⁺, K⁺, and Na⁺. Mg²⁺ was more variable, particularly during road salt event 1. During this event, concentrations of Mg²⁺ increased until approximately 16 km downstream of Rock Creek’s headwaters. Mg²⁺ then rapidly declined over a short ∼5 km stretch only to begin increasing again throughout the national park. Changes in discharge and specific conductance were evident throughout each road salt pulse as indicated in Supplementary Figure S4.

The dominant longitudinal pattern in road salt events at Rock Creek (i.e., PC1 from Figure 5B; positive PC1 indicates higher salt ion concentrations) was well explained by an MLR model including watershed drainage area (ion concentrations were higher when drainage area was higher; see partial effects plot, Figure 8C) and sampling event (ion concentrations were higher in the first road salt event than the second; see partial effects plot, Figure 8D). The fitted model explained 95% of the observed variance in PC1 at Rock Creek. As observed with Bull Run, the majority of this explanatory power was attributable to the sampling event (84 of 95%; averaging over ordering, Supplementary Table S5), which suggested that the observed differences in chemical cocktails between road salt events are greater than any longitudinal differences observed within road salt events. The remaining 11% of the explainable variance was captured by drainage area (Supplementary Table S5). This suggests that the increase in salt ion concentrations evident during road salt events across the entire Rock Creek flowpath (recall Figure 7B) may reflect increasing downstream salt inputs from successively larger drainage areas.

In Scotts Level Branch, salt ion concentrations during baseflow and stormwater events tended to decrease along the flowpath, leveling out past the second restoration site at approximately 3 km (Figure 6C). This pattern was most notable for divalent base cations (Ca²⁺, Mg²⁺). Different attenuation patterns were observed across all ions (except Mg²⁺) during road salt events. In road salt event 1, salt ion concentrations sharply declined through both restoration regions, but rebounded shortly thereafter, approaching pre-restoration levels (blue curve, Figure 6C). During the second road salt event, salt ion concentrations appear to increase through the first restoration area and decline only marginally through the second. Given that these road salt events occurred within a week of each other, one possible explanation for this pattern is that the second salting event simply overwhelmed restoration areas, and the restoration areas ceased to provide further water quality improvements.

Treated wastewater formed a unique, distinct chemical cocktail of Ca²⁺, K⁺, and Sr²⁺ (Figure 5A) in this study. During baseflow events, the concentrations of Na⁺, Ca²⁺, and K⁺ increased at the wastewater treatment plant while Mg²⁺ declined (Figure 7A). Some common elements and ions in wastewater are total nitrogen, K⁺, Na⁺, Cl⁻, and SO₄ ⁻² with Mn, Ni²⁺, Zn²⁺, Cr, Cu, As, Pb, or Cd²⁺ (Kaushal et al., 2020). While Na⁺ is typically a dominant ion in treated wastewater effluent, the excess of Na⁺ during the road salt events causes Ca²⁺ and K⁺ to be a unique signature from the effluent. Although not measured in this study, wastewater is a source of other nutrients, organic matter, emergent pollutants, and pathogens to streams (Castelar et al., 2022), all of which can create chemical cocktails that can be transported and transformed along a flowpath.

Bull Run receives the highly treated wastewater from the wastewater treatment facility and is a part of the indirect potable reuse (IPR) management strategy. Indirect potable reuse is a sustainable management practice that discharges highly treated wastewater to the surface water and the groundwater that can, in turn, be used as a drinking water source for a community downstream (Rodriguez et al., 2009). The wastewater treatment facility, UOSA, discharges effluent into Bull Run as surface water recharge which flows to the Occoquan Reservoir, the receiving waters for the Fairfax Country Griffith drinking water facility (Onyango et al., 2014). Therefore, the watershed management system is established upon the ability for Bull Run to attenuate the additional nutrients and ions, particularly total dissolved nitrogen.

Salts from the untreated wastewater are not removed and can be compounded from wastewater treatment plant and directly discharged into the stream during a road salt event (Bhide et al., 2021). During dry and median weather conditions, the treated wastewater effluent is the majority of the sodium mass load that enters the reservoir (Bhide et al., 2021). The sodium mass loadings from the wastewater treatment facility to Bull Run can vary from 60% to 80% in dry weather conditions (Bhide et al., 2021). This is shown in the PCA—the summer baseflow chemical cocktails found at Bull Run were most similar to wastewater effluent during the low flow conditions, indicating that the wastewater treatment facility was a large contributor to Bull Run’s discharge. The capacity to dilute, biogeochemically process, and attenuate these nutrients and ions from the wastewater treatment plant along a stream’s flowpath is dependent on the flow conditions, residence times of the ions, ion exchange capacity of sediments, and seasonality (Castelar et al., 2022; Kaushal et al., 2022). Depending on the climatic conditions and the discharge, Bull Run can act as a source, sink, or transporter of ions downstream of the wastewater treatment plant. Future work should continue to investigate the transport and biogeochemical transformation of various pollutants and contaminants before and after wastewater effluent enters streams.

Differences among the formation of salt ion-enriched chemical cocktails from road salt events could depend on the timing and intensity (i.e., the quantity of road salt application per impervious surface cover) of road salt application, the surrounding land use, and the timing of the sampling event within the watershed. We observed trace metal comobilization with base cations during road salt events both temporally and longitudinally, potentially due to cation exchange (Figure 7) (Kaushal et al., 2019; Galella et al., 2021; Kaushal et al., 2021). As the specific conductance rises and peaks from the influx of road salt, there is a strong relationship between the base cations and the trace metals (Galella et al., 2021). After the initial peak in road salt pulse, the declining specific conductance suggests a weaker relationship between the base cations and trace metals (Galella et al., 2021). However, the relationships between chemical cocktails downstream may be due to varying watershed characteristics such as the surrounding land use and the management strategies implemented into the watershed (Kaushal et al., 2023a).

There were not consistent longitudinal downstream patterns of dilution or attenuation along the flowpath during road salt events, even within conservation and restoration areas. Through all the management areas (i.e., Bull Run Regional Park, Rock Creek Regional and National Park, and the restoration sites along Scotts Level Branch) there were variable patterns of retention and release with Na⁺, Ca²⁺, Mg²⁺, and K⁺ (Figure 6). These variable patterns are potentially due to the hydrologic connectivity of the stream to surrounding land use through subsurface storm drains and exposed sewer lines.

The surrounding watershed land use around Bull Run has been managed and regulated since the early 1970s (there is newer development in northern Virginia compared with our urban Maryland sites). Fairfax County in Virginia has required stormwater ponds as a management strategy to protect the drinking water reservoir downstream (Fairfax County, 2023). The Cub Run tributary and Bull Run watersheds have over 420 stormwater ponds and 26% of the watershed area has been conserved and restored to improve the peak flows and pollution runoff (Fairfax County, 2023). Therefore, there are more measures to reduce stormwater and road runoff pollution implemented in Bull Run than Rock Creek, leading to differences amongst the patterns during road salt events.

Around Rock Creek, densely urban Washington D.C. surrounds Rock Creek National Park. Although the riparian buffer and forest cover increase within the park, the drainage area of the watershed and the sampling event (Figures 8C, D) were the parameters that best explained the chemical cocktails captured by PC1. The PCA in Figure 5B shows that the two road salt events had different chemical cocktails, with the first road salt event in the positive PC1 and negative PC2 space (higher concentration of salt ions, lower concentration of redox sensitive elements) and road salt event 2 tracking with the positive PC2 space (higher concentration of redox sensitive elements). The partial effect of higher drainage areas with more salt ions indicates the influence of the sewer lines and storm drains entering the stream, with some as combined storm and sanitary sewers that were built in the 1880s (Hartman, 2001). During rainfall or snowfall runoff events, some of the storm and wastewater flows to a wastewater treatment plant in the region but the remaining untreated combined sewage overflow is discharged into streams directly, including about 4 million gallons that empty into Rock Creek (Hartman, 2001; Groves and Wolfe, 2019).