Purpose and Scope

This report summarizes patterns in streamflow, water chemistry, and total suspended solids (TSS), P, and N loading rates across a 12-station monitoring network in the Hampton Roads region of Virginia, from water years (WY) 2016–20. This network is operated through a partnership between USGS and the HSRD, in cooperation with the HRPDC. Analyses included the assessments of:

Description of Study Area and Monitoring Network

Hampton Roads comprises three geographic subdivisions, all of which lie within the Coastal Plain PP. South Hampton Roads contains the cities of Chesapeake, Norfolk, Portsmouth, and Virginia Beach; the subdivision known as the Peninsula contains Hampton and Newport News; the Rural Southern Virginia region contains cities and counties not included in this study. The study watersheds are all in the northeastern portion of South Hampton Roads or the southern portion of the Peninsula, which are nearly exclusively urban (98 percent; U.S. Census Bureau, 2020) and have an elevation less than 33 ft above sea level. Twelve watersheds that range in size from 37 to 273 acres and drain a representative gradient of Hampton Roads urban land-use types have been included in this study (table 1).

The 6 jurisdictions included in this study all rank in the top 10 in population and numbers of housing units in the Commonwealth of Virginia, account for 76 percent of the population in the entire Hampton Roads region, and 16 percent of the total Virginia population (U.S. Census Bureau, 2020). The region is the 34th largest metropolitan statistical area in the United States (U.S. Census Bureau, 2020). The population grew by approximately 40 percent from 1970 through 2020, with most of that growth occurring in Virginia Beach and Chesapeake (U.S. Census Bureau, 2020). The increase in housing units (90 percent) was more than double the growth in population, contributing to a substantial increase of impervious land cover. The remainder of the Hampton Roads region is predominantly agricultural or forested land and not represented by this study; however, reference streamgages, which represent non-urban areas of the Coastal Plain PP, are discussed. Throughout the report, the term “Hampton Roads” is used to represent these six cities rather than the entire region. Major rivers include the James River and York River, as well as tributaries such as the Elizabeth River, Back River, and North Landing River; the Atlantic Ocean and Chesapeake Bay border the easternmost extent of the region.

The monitored watersheds are underlain by two geologic formations: the Shirley Formation in Newport News, and the Tabb Formation (Lynnhaven and Sedgefield Members) across the remainder of the network. The Shirley Formation is composed of interbedded gravel, sand, silt, clay, and peat, at altitudes up to 35 to 45 ft (McFarland and Bruce, 2006; Johnson and Hobbs, 1990). The Lynnhaven Member of the Tabb Formation is composed of pebbly and cobbly sand grading upward into muddy, fine sand and silt, at altitudes up to 15–18 ft (McFarland and Bruce, 2006; Johnson and Hobbs, 1990). The Sedgefield Member of the Tabb Formation is composed of pebbly to bouldery, clayey sand and shelly sand, at altitudes up to 30 ft (McFarland and Bruce, 2006; Johnson and Hobbs, 1990).

Study Design

A long-term monitoring network consisting of 12 data-collection stations distributed across 12 small urban watersheds representative of the Hampton Roads region was established in 2015 (fig. 1). Two stations were located in each jurisdiction: Chesapeake, Hampton, Newport News, Norfolk, Portsmouth, and Virginia Beach. The selection of these watersheds was determined using a statistically based approach to provide a representative range of urban land-use types, land-cover attributes, and watershed scales (fig. 2). Each jurisdiction provided a list of candidate monitoring locations; watersheds for the candidate locations were delineated and attributes determined, a cluster analysis of watershed characteristics was performed to categorize each watershed, and final selections were made on the basis of observations during field visits to each location. To properly represent the variety of land-use types in the region, the selected watersheds had a primary land use of either commercial (COM), high-density residential (HDR), or single-family residential (SFR). To achieve the objectives of the study, and to the extent possible, it was necessary to avoid monitoring in tidally influenced stormwater conveyances. Because most of the region is tidally influenced, the only monitoring stations that were both non-tidal and met all other selection criteria were zero- and first-order reaches within the MS4. One station (Ramsgate) was located in an open, concrete-lined trapezoidal conveyance channel, two stations (Daisy and Ludlow) were in road culverts, and the remaining nine stations were sited in buried concrete storm-sewer pipes accessible by either manhole or outfall.

The study was designed to establish background loading rates for TSS, total phosphorus (TP), and total nitrogen (TN); therefore, only watersheds with minimal stormwater best management practice (BMP) implementation were considered; the efficacy of BMPs was not evaluated. Watershed boundaries were delineated primarily on the basis of infrastructure maps provided by each jurisdiction. These maps included alignments of buried stormwater pipes and open conveyance channels, which did not always follow topographic features. Each stormwater network was physically walked to ensure the pipe network matched as-built diagrams. Data collection began at 10 monitoring stations on various dates throughout WY 2015 and at the other 2 stations in WY 2016. Real-time water-quality and streamflow data were recorded continuously at 5-minute increments throughout the study. Approximately 40 discrete water samples were collected at each station annually; these included a mix of event-targeted, high-flow “storm” samples and routinely scheduled “monthly” (typically low flow) samples. A complete annual (WY) timeseries was required for many of the analyses presented herein; therefore, data collected in WY 2015 are not included in this report.

Data Collection, Sampling, and Laboratory Analyses

The 12 monitoring stations were instrumented (fig. 3) to collect and record 5-minute measurements of water stage, WT, SC, and TB, collectively referred to as continuous data, which were transmitted to the USGS National Water Information System website (NWIS; U.S. Geological Survey, 2022) in near real time. Continuously collected 15-minute precipitation data were obtained from HRSD rain-gaging stations located near each of the monitoring stations. In both Hampton and Portsmouth, the two monitoring stations were associated with a single rain gage because of its proximity.

Standard USGS methods for continuously measuring and verifying water stage, measuring streamflow, and computing continuous timeseries of streamflow using either stage-discharge ratings or index-velocity ratings were followed (Rantz, 1982; Sauer and Turnipseed, 2010; U.S. Geological Survey, 2014a). A stage-discharge rating assumes free flow; however, this assumption was violated at many stations owing to backwater conditions. Although streamflow was computed on the basis of stage-discharge ratings at some stations, all stations were equipped with an in situ velocity meter, which was used to develop an index-velocity rating (Levesque and Oberg, 2012). Streamflow was calculated from the stage-discharge ratings at the Craneybrook, Daisy, Lucas Creek, Ramsgate, and Rivers Ridge stations, whereas index-velocity ratings were used at the Lindsley, Ludlow, Professional, and USAA stations. A combination of these methods to compute streamflow was used at the Coliseum, Garrett, and Sheppard stations, where backwater conditions occurred only during high flows, and in situ velocity meters performed poorly during periods of low streamflow. Despite attempts to avoid siting monitoring stations in tidally influenced locations, the effects of tide were observed occasionally at Coliseum and often at USAA. Reverse flow was not observed at either station, but during high tide a rise in water level downstream of the monitoring location created a backwater effect within the stormwater conveyance system, creating a tidal pattern in the stage data. This effect typically was not evident in the computed record of streamflow.

The operation and maintenance of continuous water-quality monitors was conducted by HRSD and USGS staff in accordance with published procedures (Wagner and others, 2006). During routine servicing visits to the stations, the instruments were cleaned and recalibrated when required on the basis of published thresholds. Data collected during these visits were used to apply corrections to WT, SC, and TB timeseries records for periods of suspected sensor fouling or calibration drift.

Base flow was sustained year-round at all stations by groundwater interchange through the concrete stormwater pipes. The chemistry of the water was characterized on the basis of analyses of monthly samples collected at each station in accordance with USGS methods (U.S. Geological Survey, 2006). All samples were analyzed for TSS and a suite of nutrient constituents (table 2). Monthly samples were collected by HRSD staff using a grab (dip) approach at the centroid of streamflow. All monitoring stations were sampled on the same day so as to limit variability in conditions across the stations; additionally, the order in which stations were visited was rotated each month. During each monthly sampling trip, three stations were randomly selected for the collection of an additional quality-assurance sample (two replicates, one blank). Additionally, from 28 to 40 event-targeted stormflow samples were collected at each monitoring station annually using an automated refrigerated sampler; these samples are referred to as storm samples throughout the report. The automated sampler was triggered to begin collection when a specified set of criteria were met. Criteria included a station-specific stage (water level) or streamflow threshold and a minimum time elapsed since the previous sample. Samples represented discrete points in time and were not composited to calculate an event mean concentration. For each storm event, a subset of the collected samples was manually selected by USGS personnel and retrieved by HRSD staff for delivery to the laboratory for analyses. Samples representing a range of conditions, including the rising, peak, and falling limbs of the storm hydrograph typically were collected for each event. All water-chemistry analyses were conducted by the HRSD Laboratory, as approved through the USGS Branch of Quality Systems Laboratory Evaluation Project (U.S. Geological Survey, 2014b). Samples were analyzed for concentrations of TSS, TP, orthophosphate (PO₄^3-), total Kjeldahl nitrogen (TKN), nitrate plus nitrite (NO₃^-), and ammonia plus ammonium (NH₃); concentrations of TN and total organic N (TON) were determined by calculation (TN is the sum of TKN and NO₃^-, TON is the difference between TKN and NH₃).

Statistical Analysis of Streamflow and Water Chemistry

Nonparametric analyses were used to describe statistical relations in water chemistry and streamflow data following methods published by Helsel and others (2020). Annual streamflow yield—the volume of streamflow generated per unit of drainage area—was calculated for each of the 12 monitoring stations as well as for 6 non-urban reference streamgages to inform hydrologic comparisons of watersheds of varying size and to investigate differences between land-use types. Non-urban reference streamgages in the Maryland, North Carolina, and Virginia Coastal Plain PP (see inset map on fig. 1) were selected to represent watersheds with low urban development (app. 1, table 1.1). Reference watersheds are orders of magnitude larger than those within the network because no comparable gaged watersheds were available at a sub-square mile scale.

Base-flow separation models and stream-flashiness metrics were used to investigate variability in streamflow characteristics across stations, land-use types, and seasons. Event-based analyses included those for stormflow duration, lag to peak, time to peak, stormflow volume, peak flow, runoff ratio, and rise rate, which were computed for extracted storm hydrographs. Spatial and seasonal patterns in water-chemistry data were evaluated using a nonparametric Kruskal-Wallace (Kruskal and Wallace, 1952) test and subsequent Steel-Dwass (Morley, 1982) post-hoc test for multiple comparisons (Critchlow and Flinger, 1991). The Steel-Dwass test is preferable to the more commonly used Wilcoxon method because Steel-Dwass controls for the overall experiment-wide error rate, thus reducing the probability of Type-1 errors (rejection of a true null hypothesis; Critchlow and Flinger, 1991). Monthly samples were collected on a pre-determined date and, by random chance, were occasionally storm-influenced. Summary statistics computed on monthly samples are intended to characterize base-flow conditions; therefore, these samples were instead included in the storm sample subset. The strength and statistical significance of correlation between metrics and response variables was evaluated using the nonparametric Spearman Rank correlation coefficient as defined by Spearman (1904). Where applicable, seasonal variability was tested by splitting the data by warm (April–September) and cool (October–March) seasons; in other instances, data were subset into four seasons based on solstice and equinox dates. Unless otherwise specified, the term “significant” is used throughout the text to denote a statistically significant difference based on a p-value less than or equal to (≤) 0.05.

Precipitation

Mean annual rainfall totals from the 10 precipitation stations were at or below the 5-year mean in 4 out of 5 years, with the lowest mean annual rainfall occurring in WY 2018 and the highest in WY 2016 (fig. 4; Porter, 2022). Total annual rainfall varied spatially across the 10-station precipitation network (mean range is 16.4 inches [in.]); the greatest variance occurred in the wettest year (WY 2016; minimum is 52.6 in.; maximum is 72.7 in.; fig. 5). Variance between stations reflects the large geographical extent of the study region and highlights the importance of utilizing local rainfall data when interpreting patterns in streamflow and water quality. In all years, approximately half of the total annual precipitation fell on days with high-intensity or high-volume (greater than [>]1.0 in.) events. In WY 2016, Tropical Storm Julia produced an average of 9.3 in. precipitation across the network over a 4-day period, and Hurricane Joaquin delivered 4.5 in. over a 3-day period. In WY 2017, an average of 7.7 in. of rain fell during Hurricane Michael, ranging from 6.8 in. to 9.4 in. across the network and exceeding the 100-year, 6- and 12-hour peak rainfall recurrence interval at most stations (Hampton Road Sanitation District, 2016; National Oceanic and Atmospheric Administration, 2022a). In the fall of WY 2019, Hurricane Sally produced an average of 4.7 in. of rainfall over 2 days across the network. An additional 69 unnamed storm events produced daily rainfall totals greater than 1.0 in., two-thirds of which occurred in the warm season (April–September). High-precipitation periods can have pronounced effects on hydrologic and water-quality conditions with the potential for conveying substantial loads of sediment and nutrients.

Streamflow

Median annual streamflow yields were highest at the 3 COM watersheds—Professional, USAA, and Coliseum: 242,000, 237,000, and 140,000 ft³/acre, respectively (fig. 6). High-density residential watersheds yielded more streamflow than SFR, a difference that was driven by the yields at Craneybrook and Rivers Ridge. Single family residential watersheds were the lowest flow-yielding land use, and exhibited only a weak statistical difference (p = 0.059) from the reference watersheds, despite substantially higher levels of development and impervious cover in the SFR watersheds. Seasonal variability in streamflow yields reflect meteorological patterns. Yields were greatest in the summer, a consequence of high-intensity storms that generated high volumes of runoff. Yields were lower in winter and spring, when low-pressure systems produced longer duration but lower intensity storms that resulted in moderate volumes of runoff, greater infiltration, and higher base flows. Streamflow yields were lowest in the fall, the season when rainfall typically is low.

A reduction in base flow and an increase in flashiness are the most common effects of land development on streamflow (Hirsch and others, 2010; Poff and others, 1997). Base flows are sustained in concrete stormwater conveyance systems as a result of hydraulic interchange with the surrounding saturated-unsaturated zone through the cracks, joints, and pore space in the concrete. When the water table is above (at higher elevation than) a portion of the pipe, groundwater inflow occurs because the pressure outside of the pipe exceeds that within it; alternatively, when the water table drops below the elevation of the pipe, outflow from the pipe can occur (Peche and others, 2019). Inflow of groundwater into the stormwater conveyance system can reduce recharge and result in lowering of the groundwater table by several inches (Thorndahl and others, 2016). Base flows may also contain leakage from potable-water pipes, which typically have leakage rates of 20 to 30 percent, and from sanitary sewer lines, which are commonly constructed in close proximity to stormwater conveyances (Kaushal and Belt, 2012; Garcia-Fresca, 2007). Monitored stormwater conveyances in the Hampton Roads study area had significantly lower annual BFIs than non-urban reference streams (fig. 7A). Single-family residential watersheds had significantly higher annual BFI than the COM and HDR stations, but only two of the SFR stations (Sheppard and Ramsgate) were baseflow dominant (> 0.5) in most years. Streams in the Coastal Plain PP naturally have high base flows (greater than or equal to [≥] 75 percent of total flow) because of low surface runoff, the result of high infiltration rates in sandy soils and low topographic slope (Sanford and others, 2012). In urbanized areas, disturbance factors such as soil compaction, loss of vegetation, construction of impervious surfaces, and the direct channelization of runoff to stormwater conveyance systems reduce groundwater recharge, decrease infiltration, and reduce storage and release of base-flow sustaining groundwater to streams (Rose and Peters, 2001).

Annual RBI scores were computed for each of the 12 monitoring stations and the 6 reference streamgages. Flashiness scores ranged widely across the 12-network stations from 0.83 at Lindsley to 0.25 at Ludlow. Flows at stations in HDR and COM watersheds were flashiest, with median RBI scores of 0.73 and 0.62, respectively, followed by flows at stations in SFR watersheds (0.36), and in the reference watersheds (0.02; fig. 7B). Streams in natural Coastal Plain PP watersheds are typically more stable (less flashy) than streams in other PPs because of high soil permeability, low topographic gradient, wide stream channels, and broad alluvial floodplains. High flashiness scores computed across the monitoring network are likely related to the engineered stormwater systems that serve these watersheds, which disconnect surface water hydrological processes from the natural land attributes that typically shape the hydrologic regime in the Coastal Plain PP. Flashiness is also likely related to specific design characteristics of the constructed stormwater system, watershed size, and the degree of impervious cover. These systems are engineered to export water quickly and efficiently; however, to achieve this aim channels are built on a higher gradient and with a more direct flow path than natural Coastal Plain PP streams. Differences in degree of impervious cover and in watershed size may explain much of the spatial variability observed across different land-use types. Single family residential watersheds were larger, had less impervious cover, and greater composition of land cover attributes (such as turf grass and tree canopy) that promote uptake, retention, and interception of runoff than did HDR or COM watersheds. A significant seasonal effect was present at all stations, with greater flashiness observed in warm months than in cool months, reflecting meteorological patterns—higher intensity convective systems such as thunderstorms in warmer months and lower intensity low-pressure systems in cool months.

An evaluation of flashiness across 41 USGS streamgaging stations, which included the 12 network stations in this study, CB-NTN, Coastal Plain PP reference stations, and multiple other urban monitoring stations in Fairfax County, Virginia and the City of Roanoke, Virginia, was conducted to better understand factors affecting stream flashiness across a wide range of watershed types (small versus large, urban versus non-urban). This analysis revealed strong correlations between RBI scores and both watershed size and impervious cover (fig. 8A, B; Porter, 2022). Watershed area and impervious cover are themselves highly correlated, so a stepwise linear regression analysis was used to identify the primary driver. Although watershed area likely influences flashiness, impervious cover was the strongest driver across these 41 monitored watersheds, explaining 84 percent of spatial variability. Flashiness increased by approximately 10 percent for every 10 percent increase in impervious cover, a relation that was consistent across different drainage network types (for example, earthen streams versus engineered conveyances).

Streamflows in most of the study watersheds were stormflow dominated. To gain greater insight into this component of the hydrologic regime, storm hydrographs were extracted from continuous records of streamflow and used to characterize spatial and temporal variability in runoff response. In total, 4,720 individual storm hydrographs were extracted from the data for the 12 monitoring stations (200 to 535 storms per station). For each storm, the total streamflow yield, peak flow, runoff ratio, rise rate, lag to peak, time to peak, and event duration were computed (table 3). Higher values for metrics that describe volume, magnitude, and rate of change in flow, such as event yield, peak flow, RBI, runoff ratio, and rise rate typically are related to increasing urbanization (Paul and Meyer, 2001). Generally, as magnitudes of these metrics increase, so does the potential for constituent transport, depending on the pollutant and its source (O’Driscoll and others, 2010). Conversely, longer duration, lag to peak, and time to peak commonly result in greater constituent retention and are indicative of a more stable hydrologic system.

Event yield was positively related to precipitation volume at all stations, but the slopes of these relations varied by land use. Slope was steepest in COM watersheds followed by lesser slopes in HDR and SFR (fig. 9; app. 1, table 1.2). This demonstrates that differences in runoff, and consequently differences in streamflow, are generated by an equal volume of rain falling across these three land-use types, and likely reflect differences in impervious cover, and in SFR watersheds in particular, landscape features like lawns and trees that promote infiltration and increase retention. Median runoff ratios varied widely across stations; from 10 to 72 percent of precipitation ran directly to the conveyance system as direct runoff and was positively related to watershed imperviousness (fig. 10). On average, approximately 10 percent of precipitation becomes direct runoff in forested watersheds versus 55 percent in watersheds with 75–100 percent imperviousness (Arnold and Gibbons, 1996). Runoff ratio values less than 10 percent are common in the Coastal Plain PP (Sanford and others, 2012). These differences reflect lower infiltration and evapotranspiration with increasing imperviousness.

Principal component analysis was used to explore patterns in event-based hydrologic metrics across stations, land-use types, and seasons (fig. 11). Principal components (PC) 1 and 2 explained a combined 83.3 percent of the variability in the data. The strongest loadings to PC1 included metrics describing the magnitude (peak) and rate of change (rise rate and lag to peak) in streamflow (app. 1, table 1.3). The strongest loadings to PC2 included event yield, runoff ratio, and time-based metrics like duration and time to peak.

Commercial watersheds (Coliseum, USAA, and Professional) had the highest event yields, runoff ratios, and peak flows; rise rate was also highest at Coliseum and Professional. Variability between SFR and HDR watersheds was less pronounced, though storms at SFR stations typically were longer in duration and time to peak, and had a longer lag between peak precipitation and peak streamflow. Average storm duration, time to peak, and lag to peak were longest at the two largest watersheds, Ramsgate and Ludlow, which indicate slower export of runoff. This pattern was likely affected by two factors: (1) larger watershed area, which inherently results in a longer flow path, and (2) a high proportion of pervious land cover, such as lawn turf and tree canopy (table 4). These land-cover features can serve to disconnect non-road impervious surfaces like rooftops from stormwater inlets. Building on previous analyses, these results suggest a greater volume and rate of runoff in watersheds with high impervious cover that consequently lack areas for retention and infiltration. In particular, characteristically short lag to peak and high-rise rate at COM stations suggest a short and direct path of travel from the land surface to stormwater conveyances. It is important to note that the correlations do not imply causation, as additional factors may be driving differences in streamflow metrics.

Seasonal variability was evident in the spread of warm and cool season variables on PC1, suggesting greater magnitude and rate of change in warm months, and longer duration, time to peak, and lag to peak in cool months. Similar to stream flashiness, temporal differences in event-based metrics likely are related to seasonal meteorological patterns, such as the duration, intensity, and recurrence of precipitation events. Unlike metrics describing event magnitude, timing, and rate of change, runoff ratios and event yields did not differ across seasons. This suggests that although the hydrograph was less flashy in the cool months, the percentage of precipitation that reaches the storm conveyance and consequently the volume of runoff transported out of the watershed during a storm event is consistent year-round. Consistency in yield can be explained by the seasonal nature of precipitation events in the Mid-Atlantic region of the United States—short duration-high intensity events in the warm season and long duration-low intensity events in the cool season. Both types of storms can deliver the same volume of rain across a watershed, and in doing so produce a similar volume of streamflow, but more critically, the rate at which runoff is produced differs, and may affect the magnitude of event-derived sediment and nutrient loading. These patterns highlight the importance of considering seasonal differences in storm-event hydrology as part of stormwater-management planning and evaluation.

Water Temperature

Water temperature (WT) is a fundamental property of water quality given its role in regulating chemical and biological reactions and governing the structure of aquatic communities (Demars and others, 2011; Hillebrand and others, 2010; Mulholland and others, 2001). Water temperatures in the stormflow and monthly samples ranged from 2.0 to 31.2 degrees Celsius (°C) and from -1.3 to 35.2 °C in the continuously collected data. Water temperatures were seasonally variable, with a median networkwide cool season (October–March) value of 13.9 °C and warm season value of 21.8 °C. A positive correlation (r = 0.59, p = 0.044) was observed between WT and degree of impervious land cover in a watershed. This relation has previously been linked to increased runoff, decreased base flows, and the low albedo, or reflectivity, of roads, sidewalks, and buildings that gain and hold more heat, and consequently increase the temperature of overland runoff (Galli, 1991; O’Driscoll and others, 2010). Unlike the findings in studies of earthen streams in the urban environment (Hyer and others, 2016; Jastram, 2014; LeBlanc and others, 1997; Porter and others, 2020), diel fluctuations in daily water temperature were minimal at most stations, and this pattern was consistent throughout the year (fig. 12). Diel Fluctuations were larger at the three stations most exposed to ambient air: Ramsgate (open channel), and Ludlow and Daisy (storm drains connected to road culverts). Daily fluctuations in WT at all stations were attributed to similar oscillations in air temperature and were greatest in late winter and early spring because of heating by solar radiation prior to leaf out of the tree canopy (Porter and others, 2020; Johnson, 2004). The other nine stations, which are located in buried concrete storm pipes, are more buffered against changes in ambient air temperatures. As a result, WT also was warmer at these stations during cool season months. This difference highlights an important physicochemical difference in water flowing through earthen streams versus engineered stormwater conveyances that may affect water quality. Warmer temperatures stimulate microbial respiration, which in turn promotes the uncoupling of nutrients from organic matter and can increase the export of bioavailable P and N (Demars and others, 2011). This effect may be most evident during periods of base flow when residence time is longer and therefore biogeochemical transformations of organic substrates can occur within the storm conveyance network (Kaushal and Belt, 2012).

Specific Conductance

Specific conductance (SC) is a measure of the ability of water to conduct electric current, corrected for temperature (typically 25 ^oC) and provides an indication of ion concentration in water (Hem, 1985). Streams in less developed areas of the Coastal Plain PP typically have mean SC values less than 100 microsiemens per centimeter (µS/cm) because of the relative abundance of quartz sand and lack of soluble minerals in the soils and underlying geologic formations (Sanford and others, 2012). In areas where temperature, geology, and precipitation are relatively constant, SC serves as an indicator of urban impacts on stream chemistry (Paul and Meyer, 2001).

The ionic strength of freshwater streams has increased across the United States because of anthropogenic influences in recent decades (Griffith, 2014). Freshwater salinization is a growing concern in urban watersheds, and it has been estimated that a significant change in the degree of salinization has occurred in 37 percent of streams in the contiguous United States (Kaushal and others, 2018). Salinization can degrade stream and drinking-water quality through the mobilization of N and P from soils (Haq and others, 2018; Steele and Aitkenhead-Peterson, 2013), the leaching of metals from water-conveying infrastructure (Kaushal, 2016; St. Clair and others, 2016), and numerous negative impacts on aquatic life (Corsi and others, 2010). Factors affecting ion concentrations, and thus SC, include temperature and geologic weathering (Hem, 1985), saltwater intrusion (Kaushal and others, 2021), and the discharge of wastewater to streams or groundwater (Paul and Meyer, 2001; Potter and others, 2014). Road anti-icer and de-icer applications, in particular, have been linked to elevated ion concentrations, such as sodium and chloride, in streams (Shanley, 1994). These applications are considered a major contributor to freshwater salinization in the urban environment and have been documented throughout the northeastern and southeastern United States (Kaushal and others, 2005; Moore and others, 2020).

Weathering of concrete building materials also has been linked to salinization and other changes in stream chemistry and is particularly evident in regions underlain by a dense network of concrete stormwater infrastructure (Kaushal and others, 2017). Weathering of these concrete pipes—“urban karst”—serves as an artificial geologic source of mineral leaching that alters SC and other chemical measures (Kaushal and Belt, 2012). The effect of this “geologic” source is similar to that observed in karst geologic terranes where the dissolution of carbonate rock, primarily limestone, alters groundwater and stream chemistry. Numerous small-scale controlled studies have linked increases in SC, pH, bicarbonate, and calcium to the dissolution of concrete materials after passing rainwater through concrete pipes (Davies and others, 2010c; Davies and others, 2010b; Davies and others, 2010a; Setunge and others, 2009). Likewise, Wright and others (2011) reported a doubling of electrical conductance of the water in urban streams as compared to that in non-urban reference streams in field studies (433 µS/cm and 194 µS/cm, respectively). In those studies, dissolution of the cement products used to manufacture the stormwater pipes was linked to a tenfold increase in bicarbonate concentration and sixfold increase in calcium concentration in streams.

Specific conductance typically was inversely related to streamflow networkwide (r = 0.55, p < 0.001), indicating a dilution of higher SC groundwater by lower SC surface runoff. Median SC values ranged from 191 μS/cm at Lucas Creek to 675 μS/cm at Coliseum (fig. 13). Notably, median SC at Coliseum, Ludlow and USAA was greater than 500 μS/cm, whereas median values at the remaining nine stations were less than 270 μS/cm. Specific conductance was significantly higher in COM watersheds than HDR and SFR watersheds. The soils and geology are relatively similar across the study area; therefore, differences in SC may reflect urban impacts or groundwater salinity in these watersheds. In the USAA watershed, an aesthetic pond and a wet detention pond are connected to the stormwater conveyance network. Such landscape features can trap and accumulate road salts, eventually delivering plumes of sodium and chloride to groundwater and surface waters (Snodgrass and others, 2017; Casey and others, 2012); however, data are not currently available to investigate this link in the USAA watershed. Investigating the effect of land use on SC, and by proxy salinization, was complicated by the proximity of USAA and Coliseum (COM watersheds) and Ludlow (SFR) to tidal waters. This likely contributes to higher salinity, and therefore higher SC in the surficial aquifer that interchanges with surface flow in these systems. As a natural factor, this driver of elevated SC is unrelated to differences in urban development and human activities in COM and residential watersheds. After removing these stations from the analysis, an insufficient number of COM stations remain to effectively evaluate differences related to land use. Spatial variability in SC may be related to the age of concrete pipe infrastructure, though that information was not available.

Specific conductance of streamwater was significantly higher in cool months at 8 of the 12 stations and higher in warm months at Ramsgate. Elevated SC in cool months may suggest a seasonal effect resulting from the application and subsequent wash off of de-icer and anti-icer salts. At all stations, spikes in SC values often coincided with observed or forecast winter-weather (snow, sleet, ice; National Oceanic and Atmospheric Administration, 2022b), and typically were more than an order of magnitude higher than baseline SC values. Winter-weather storms are less common in Hampton Roads than in the northern and western regions of Virginia (National Oceanic and Atmospheric Administration, 2022b); however, even infrequent applications of anti-icer and de-icer salts can alter surface and groundwater chemistry (Kaushal and others, 2021). The ionic properties of these salts favor retention in soils and shallow groundwater and slow release over decades (Robinson and others, 2017). Furthermore, Cl^- retention capacity is highest in sandy soils with high organic matter content, which are common characteristics of Coastal Plain PP soils (Robinson and others, 2017). Higher SC in cool months may be related to seasonally elevated base flows, which act to resuspend salts applied in previous winters, and can elevate SC even in years when de-icing salts are scantly applied.

Total Suspended Solids and Turbidity

Excess suspended sediment is a leading cause of water-quality and ecosystem impairment in urban streams. Although sediment transport occurs naturally and is controlled by fluctuations in streamflow, anthropogenic alterations like land development and stream channelization can increase erosional rates, decrease sediment transport time to streams, and increase the magnitude and recurrence of peak flows (Paul and Meyer, 2001; Soulsby and others, 2006). Suspended fine sediments, which account for 90 percent of the annual sediment load to Chesapeake Bay (Zhang and Blomquist, 2018), can serve as vectors by which contaminants are transported downstream (U.S. Environmental Protection Agency, 2000a). For example, sediment transported to Chesapeake Bay carries 77 and 18 percent of annual loads of TP and TN, respectively, and also conveys metals, pesticides, poly-chlorinated biphenyls, and other organic contaminants (Zhang and others, 2015). The average concentration of P and N on suspended-sediment particles in the Chesapeake Bay watershed is 1.0 milligram of phosphorus per gram and 3.6 milligrams of nitrogen per gram, respectively (Zhang and others, 2018). Additionally, fine sediments reduce light availability, causing reductions in photosynthesis and primary production (Wood, 1997), and, when deposited, can harm aquatic organisms (Luedtke and Brusven, 1976; Poulton and others, 2007; Rasmussen and others, 2009b).

Previous studies have reported higher concentrations of suspended sediment, suspended solids, and turbidity in urban watersheds than in agricultural or forested watersheds (Coulter and others, 2004; Ghane and others, 2016). Agricultural lands are more widespread than urban lands in the Chesapeake Bay watershed and account for 69 percent of the suspended sediment load to Chesapeake Bay (Brakebill and others, 2010; Langland and Cronin, 2003). The average sediment yield from urban lands, however, is approximately 6-times greater than forested yields and 3-times greater than yields from agricultural lands (Russell and others, 2017). Sediment loads typically are greater in rapidly urbanizing areas with active construction involving the removal of vegetation and excavation of soils than in older, previously developed urban regions (Carpenter and others, 1998; Langland and Cronin, 2003; Lee and Ziegler, 2010). In established, more stable impervious areas, common sources of suspended sediment and other suspended solids are dry and wet atmospheric deposition, erosion of the built environment (buildings and public infrastructure), and deterioration of vehicle parts and by-products. Soil erosion from pervious areas, especially poorly maintained lots, also contributes to the sediment load. Most of these solids are mobilized by overland runoff. These source pools are inherently limited in highly impervious landscapes; therefore, heavy rainfall events can temporarily deplete the sediment “reservoirs,” resulting in lower yields from successive storms (Lee and Ziegler, 2010).

Turbidity (TB) is a measure of the relative clarity of a body of water and typically has a strong relation with the concentration of total suspended solids (TSS; Swenson and Baldwin, 1965). Unlike TSS concentration, however, TB can be measured continuously and often can be used to estimate TSS concentration and load more accurately than models that include only a streamflow term (Jastram and others, 2009; Rasmussen and others, 2009a). In this study, the concentration of TSS was used in lieu of suspended-sediment concentration because TSS is the analytical metric used by the U.S. Environmental Protection Agency (EPA) for the Chesapeake Bay TMDL; however, it has been well established that the TSS method can result in the underestimation of suspended-sediment concentration (Galloway and others, 2005; Gray and others, 2000; Rasmussen and others, 2009a). This bias is particularly evident when sand constitutes more than 20 percent of a sample, which may be the case in Coastal Plain PP watersheds. Sand accounts for an average of 15 to 52 percent of the samples collected at four CB-NTN stations in the Virginia Coastal Plain PP (U.S. Geological Survey, 2022).

Concentrations of TSS in monthly samples, which most commonly represent base-flow conditions, generally were low (network median = 4.7 mg/L) but were substantially higher in samples collected during stormflow conditions (network median = 41.2 mg/L; app. 1, table 1.4). Total suspended solids concentration typically increased with streamflow at all stations (r = 0.631, p < 0.001; app. 1, table 1.5). At each station, a representative range of concentrations was sampled to allow for more accurate calculation of TSS loads. Median TSS concentrations were highest at Daisy (72.5 mg/L) and lowest at Ludlow (14.8 mg/L), both of which drain SFR watersheds (fig. 14). TSS concentrations were significantly higher in SFR watersheds than in HDR watersheds; concentrations were lower in COM watersheds than in both residential land-use types.

The median TSS concentration across all samples (34.0 mg/L) was similar to NSQD data collected in Hampton Roads (41.0 mg/L) and lower than the nationwide median (56.0 mg/L; app. 1, table 1.6). Concentrations reported in the NSQD were subset to reflect only residential and COM land-use types to provide a direct comparison to data collected for this study. Total suspended solids concentrations did not differ significantly between residential and COM land-use types in the NSQD; however, more generally, land use was reported as a significant factor affecting TSS (Pitt and others, 2018).

Runoff on the rising limb of a storm hydrograph often is more contaminated than the water on the recession limb. Differences result from factors such as mobilization of anthropogenic contaminants, such as excess sediment and nutrients accumulated on the land surface during the antecedent dry weather period, and the rapid transport of precipitation runoff to the storm conveyance. Whether TB peaks before or after peak streamflow can inform the source and supply of TSS or other particulate-dominant constituents in the watershed and may provide insight for how best to implement BMPs to achieve reductions in load. The timing of peak TB on a hydrograph typically is identified through analyses of TB-streamflow hysteresis patterns; however, this method is best used to investigate single storm events rather than to characterize conditions across tens to hundreds of storms. Cross-correlation analyses, which considered data from approximately 200 to 500 storm events at each station, revealed TB peaks within 1 hour of peak streamflow in 87 percent of storms, indicating that the transport of particulate material is strongly correlated to increasing streamflow (fig. 15). At all stations, TB typically peaked before streamflow, which indicates disproportionate loading during the rising limb of the storm—a first flush. The concentration-based first-flush criterion has been previously defined by high initial concentration in the early part of a rainfall-runoff event, with a subsequent rapid concentration decline (Gupta and Saul, 1996; Sansalone and Cristina, 2004). Although lags were less than one hour for most storms, some events had strong positive lags of more than one hour (12 percent), indicating TB peaked well before peak streamflow and declined as the event progressed. This pattern suggests TSS were delivered quickly to the storm conveyances from storages within the pipe network, likely those sediments deposited by previous events, and accretion on impervious surfaces (Landers and Sturm, 2013).

Cross-correlation analyses provide visual evidence of a first flush during some storm events, and mass-volume (M[V]) plots were used to quantify the occurrence of this phenomenon based on published criteria. The average maximum difference (ΔMax) in the cumulative TSS load and cumulative stormflow volume ranged from 0.20 to 0.31 at 10 of the 12 stations, exceeding Geiger’s (1987) 0.20 first-flush threshold; only the stations Daisy (0.18) and Ludlow (0.14) were below this level (table 5). Incorporating Gupta and Saul’s (1996) time term revealed that across the network, an average ΔMax of 0.23 occurred in the first 14 percent of the storm duration (which represents the first 1.21 hours of the average storm), conveying 68 percent of the TSS load in just 45 percent of the stormflow volume. The range of values across stations was minimal. The strongest average first flush occurred at Ramsgate, where ΔMax equaled 0.31 and occurred in the first 8.8 percent of the storm event, conveying 73 percent of the TSS load in just 42 percent of the stormflow volume. The weakest first-flush statistics were observed at Ludlow, where ΔMax occurred in the first 16 percent of the storm event, conveying 63 percent of the TSS load in the first 50 percent of the stormflow volume. Median ΔMax exceeded the 0.2 threshold at most stations and typically occurred in the first 20 percent of the storm duration, but there was considerable inter-storm variability (fig. 16). Variability between storm events was likely because of differences in rainfall volume and intensity, the antecedent dry weather period, and the accumulated mass of solids across the watershed and within the storm conveyance system (Bertrand-Krajewski and others, 1998; Vaze and Chiew, 2002). Given such variance, the prevalence and importance of the first-flush phenomenon is better investigated by computing the frequency of storms that exceed the 0.2 threshold. This threshold was exceeded in more than 50 percent of storms at 6 of the stations and ranged from 77 percent of the time at USAA to only 22 percent of the time at Ludlow (table 5, fig. 17). Generally, the strength of the first flush is proportional to degree of impervious cover, and inversely related to watershed size (Li and others, 2007; Ma and others, 2002). Delta-max scores were highest in COM watersheds, but scores for the two residential type watersheds were not significantly different from one another. Across the three land-use types there exists an increasing impervious gradient from COM to SFR and inverse relation with respect to watershed size; however, correlations between ΔMax and both impervious cover (r = 0.07; p = 0.829) and watershed area (r = -0.15; p = 0.649) were weak and non-significant.

Previous studies also have used M(V) curves to suggest a first flush based on the percentage of the load conveyed by a given initial percentage of the stormflow volume; however, the criteria used to define this effect varies widely. For example, Wanielista and Yousel (1993) defined a first flush when 50 percent of the pollutant mass was transported in the first 25 percent of the volume. Bertrand-Krajewski and others (1998) employed a stricter criterion of at least 80 percent of the pollutant mass in the first 30 percent of stormflow volume, whereas Geiger (1987) defined the phenomenon as occurring when the mass ratio exceeds the stormflow volume ratio at all instances through an event. Regardless of the criterion used to define this phenomenon, this analysis describes the volume of stormflow that must be treated to reduce a given percentage of the load. In this study, station averages ranged from 60 to 75 percent TSS mass transport in the first 50 percent of stormflow volume, which when evaluated, along with ΔMax values and cross-correlation, results suggest the first flush is an important factor affecting TSS loadings during runoff events in these watersheds.

At all 12 stations, TB and streamflow were strong predictors of TSS concentration and load and explained, on average, greater than 90 percent of variability in load models and 80 percent of variability in concentration models (app. 1, table 1.7). Models typically were positively biased; however, the magnitude of bias was low—6 percent on average across the 12 station-specific models—and had high efficiency scores (average Nash-Sutcliffe efficiency = 0.72), suggesting accurate predictions of TSS loads at all stations. Half of the models included a base-flow term to adjust the y-intercept for differences in low-flow versus high-flow sample populations, and similarly a base-flow-by-TB interactive term to adjust the slope of the concentration-TB relation. At stations USAA, Daisy, and Ludlow, either water temperature or a seasonal term was fit, suggesting seasonal differences in the source of TSS inputs that alter the TSS-TB relation. Model coefficients and unit-value estimates of load and yield have been published (Porter, 2022).

Over the study period, median annual yields of TSS were highest at Professional (764 pounds per acre [lbs/acre]) and Coliseum (701 lbs/acre), two COM watersheds; the lowest yields occurred at Garrett (135 lbs/acre) and Ludlow (138 lbs/acre), two SFR watersheds (fig. 18; table 6). The average annual yield across all stations was 347 lbs/acre. In contrast to TSS concentrations, COM watersheds had on average higher TSS yields than the two residential land-use types. The seemingly contradictory relation between concentration and yield results from elevated streamflow yields generated in the COM watersheds, which transport a greater mass of TSS, albeit at a lower concentration. Across all stations, an average 88 percent of TSS loads was conveyed by stormflows (fig. 19). In past studies, both runoff volume and constituent transport have been well correlated with urban factors, such as total impervious area, effective impervious area (impervious areas connected directly to stormwater inlets), and road density (Hatt and others, 2004; Janke and others, 2017; Walsh and others, 2005b). Annual TSS yield was moderately correlated with annual streamflow yield and may reflect limited pools of available sediment in these watersheds that are typically exhausted in the early phase of a storm, regardless of magnitude or duration (fig. 2.1A).

Comparison of sediment yields (in terms of either TSS or suspended sediment concentration) from Hampton Roads, Fairfax County (Porter and others, 2020), the CB-NTN (Moyer and Langland, 2020), and Gwinnett County, Georgia (Aulenbach and others, 2017) reflect key factors affecting sediment transport such as watershed size, land development, physiographic province, and differences in earthen streams versus engineered concrete stormwater conveyances. Headwater streams typically transport more sediment than higher-order streams owing to a lack of opportunity for in-channel deposition and greater energy availability for entrainment (Smith and Wilcock, 2015), and sediment yields typically are higher in developed watersheds than in forested or agricultural watersheds (Brakebill and others 2010). The Hampton Roads watersheds are an order of magnitude smaller and substantially more developed (high percent of impervious cover) than those draining the monitored Fairfax County streams, yet sediment yields were 2 to 3.5 times lower at Hampton Roads (fig. 20).

Low TSS yields in Hampton Roads watersheds may result from a lack of bank erosion, a sediment source which is absent in concrete pipes and conveyance channels. Yields also may differ between comparison networks because of unique characteristics of each physiographic province. Higher yields in the Piedmont PP have been linked to available pools of legacy sediment and high rates of soil erosion, whereas very little fine-grained sediment is transported from the characteristically sandy streams of the Coastal Plain PP (Gellis and others, 2009). In fact, from 1952 through 2001, the average annual sediment yield of Coastal Plain PP streams (106 pounds per acre per year [lbs/acre/yr]) was nearly 10 times lower than that of Piedmont PP streams (925 lbs/acre/yr), and nearly 6 times lower than streams in the Appalachian Plateaus, Blue Ridge, and Valley and Ridge PPs (Gellis and others, 2009). Coastal Plain PP streams have been described as “bottle-necks” (Phillips and Lutz, 2008) that are created by a transition from narrow, high-energy bedrock streams in the Piedmont PP to low-energy systems with broad alluvial floodplains. These latter streams have high rates of inundation and sediment attenuation and limit delivery of sediment to Chesapeake Bay (Noe and Hupp, 2009; Ensign and others, 2015). Although TSS yields from streams in the study watersheds in the Coastal Plain PP were low relative to yields from monitored urban streams in other physiographic provinces, previous research has reported a tenfold decrease in rates of sediment trapping in urbanized Coastal Plain PP streams that, like those in this study, have been artificially channelized and consequently disconnected from surrounding floodplains (Noe and Hupp, 2009). Hampton Roads watersheds are several orders of magnitude smaller and far more developed than most CB-NTN streams, but yields did not differ significantly. That yields did not differ, despite watershed size, channelization, and high levels of imperviousness, the similarity in yields suggests a lack of available substrate and mechanism for TSS entrainment in urban Coastal Plain PP watersheds.

Phosphorus

Phosphorus is an essential nutrient for the growth of plants and animals but excess P in aquatic ecosystems can lead to eutrophication and, in turn, to ecosystem degradation. Phosphorus readily binds to soils and fluvial sediment, and therefore generally enters surface waters via overland flow rather than as inflow of groundwater. Phosphorus commonly is associated with erosional processes and transported downstream via suspended solids (Correll, 1998; Denver and others, 2010); however, these matrices have a finite capacity to sorb P (Richardson and Craft, 1993). When this capacity is met, referred to as P saturation, the unsorbed P fractions may be exported to streams as dissolved P (primarily in the form of orthophosphate), and can result in degraded aquatic ecosystems (Conley and others, 2009), as this is the preferred form for uptake by aquatic organisms. Further, the release of P from particulates can be catalyzed by warmer temperatures (Duan and Kaushal, 2013), anoxia (House and Denison, 2002; Lehtoranta and Pitkänen, 2003), anthropogenic inputs (Sharpley, 1995), increases in pH (James and Barko, 2004), and changing salinity levels (Haq and others, 2018). For example, warmer temperatures stimulate microbial degradation of the organic material to which the P was previously sorbed (Duan and Kaushal, 2013). In oxic environments, P bonds to organic matter and to clays containing iron (Fe) and aluminum (Al) oxides and may remain stored in streambank and streambed sediments (Fink and others, 2016; Lehtoranta and Pitkänen, 2003; Manning and Wang, 1994). Under reduced oxygen conditions, however, those bonds can weaken, resulting in diffusion of dissolved P to the water column (Correll, 1998). Additionally, as pH increases, P is liberated from Fe and Al oxides because of competition with hydroxide ions (Fink and others, 2016). Freshwater salinization also may increase P concentrations in streams through competition for binding sites on soil particles with negatively charged ions in the saline water (Bruland and DeMent, 2009).

In undeveloped watersheds, phosphorus bound to soils is mostly immobile; in urban systems, however, overland flow across impervious surfaces can mobilize and transport sediments and plant material to storm drains and the downstream receiving body. As a result, transport of P typically is higher in urban than in undeveloped watersheds (Alvarez-Cobelas and others, 2009; Duan and others, 2012; Meybeck, 1998). Urbanization has been associated with a shift in P transport from low to high streamflow conditions that result from greater runoff volume and increased sediment mobilization (Duan and others, 2012). Common sources of phosphorus in urban watersheds include fertilizer (applied to lawns), geologically derived P, animal or pet waste, decomposition of plant and animal matter, road salt, and point-source discharges of wastewater (LaValle, 1975; Paul and Meyer, 2001). There are no permitted point-source discharges of wastewater into the monitored Hampton Roads watersheds, highlighting the importance of understanding of non-point source inputs. Turf grass and roadways typically account for most of the P in urban runoff. Grass clippings, leaf litter, and fertilizer are the predominant sources in residential areas, whereas runoff from roadways is a bigger contributor in COM and industrial watersheds (Waschbusch and others, 1999). Orthophosphate stored in sediment pore water also can be remobilized to streams and storm conveyance systems during high-runoff events (Duan and others, 2012). A phosphorus criterion specific to Virginia streams currently does not exist; however, the EPA reference condition for streams in Hampton Roads (ecoregion 63) is 0.05 mg/L (U.S. Environmental Protection Agency, 2000b).

Concentrations of total phosphorus and PO₄^3- in monthly samples, which typically represent base-flow conditions, were generally low (network median TP = 0.06 mg/L; PO₄^3- = 0.01 mg/L; app. 1, table 1.4) but substantially higher in samples collected during stormflow conditions (network median TP = 0.24 mg/L; PO₄^3- = 0.08 mg/L). In stormflow samples, median TP concentrations were highest at Daisy (0.39 mg/L) and lowest at Professional (0.07 mg/L). Concentrations of TP were higher in SFR watersheds than in HDR watersheds; concentrations at COM were lower than in both residential types (fig. 21). This pattern may reflect differences in the type or quantity of P inputs in these watersheds or may reflect greater dilution across the impervious gradient—a function of differences in stormflow yields. A likely source of P in these watersheds is the flushing of solids accumulated on impervious surfaces during the period between storm events. Inputs of P sources such as grass clippings, leaf litter, and pet waste are likely more abundant in residential than COM areas because of differences in land-cover attributes and human activities. In a previous study, TP concentration in street accumulations was 29 percent greater in residential watersheds than in COM watersheds and accumulated at five times the rate in the former type of watersheds (Sorenson, 2013). Concentration increased with streamflow at most stations (r = 0.424, p < 0.001); however, this relation was slightly negative at Ludlow and weak at the other stations that experience backwater during periods of high flow or high tide (Professional, Coliseum, USAA, and Sheppard; app. 1, table 1.5). Backwater conditions may affect concentration of P by dilution or enhanced settling of particulates because of low water velocity; however, particulates that become disentrained are likely remobilized on the rising limb of the succeeding storm.

The median TP concentration across all network samples (0.22 mg/L) was lower than NSQD stations located in residential and COM land uses (0.27 mg/L) and the subset of those stations in Hampton Roads (0.27 mg/L; app. 1, table 1.6). Likewise, the median PO₄^3- concentration in samples collected in this study (0.07 mg/L) was less than half the nationwide average reported in the NSQD (0.16 mg/L). The results for orthophosphate analyses were not available for the subset of NSQD stations in Hampton Roads, though the median dissolved phosphorus concentration (0.14 mg/L) was double that in the study watersheds. Similar to this study, the NSQD reports a significantly higher concentration of TP and dissolved P in residential than COM land uses.

The concentration and composition of TP varied both spatially and seasonally. Concentrations of PO₄^3- and the remaining portion of TP, defined as total particulate P plus dissolved organic P, were inversely related to the development gradient in the three land-use types: COM < HDR < SFR (fig. 22). Orthophosphate accounted for a significantly smaller proportion of TP in COM watersheds than in the two residential types, suggesting spatial variability in sources of P across land-use types. Seasonal shifts in composition indicate variability in sources throughout the year. Orthophosphate concentrations increased in summer and fall, and declined sharply in winter and spring, whereas the particulate plus organic fraction was highest in winter and spring (app. 1, table 1.8). Elevated PO₄^3- concentrations in summer months may be related to high rates of microbial degradation of organic matter and diffusion of P from soils and metal oxides. Additionally, annual P loadings from unbagged grass clippings left on nearby impervious surfaces such as sidewalks and curbs can be substantial, ranging from 0.35 to 6.69 lbs/acre (Soldat and Petrovic, 2008). In fall months, leaf senescence can contribute substantial inputs of organic detritus to storm conveyance systems and has been linked to high concentrations of in-stream dissolved P (Janke and others, 2017; Kluesener and Lee, 1974; Selbig, 2016). Most of the phosphorus contained in leaf litter can dissimilate in the first 24 to 48 hours of soaking (Hobbie and others, 2014; Wallace and others, 2008). In a study of urban stormwater in Madison, Wisconsin, 56 percent of the annual TP yield was related to leaf litter in the fall, compared to only 16 percent in study areas where leaf removal was implemented, and most of the P was in the dissolved form – 85 percent in fall versus less than 50 percent in spring (Selbig, 2016). Likewise, Janke and others (2017) report a strong correlation between urban tree canopy and P loads (r = 0.84). Lawn fertilizer containing P most commonly is applied in fall months during establishment or overseeding of cool-season grasses, as allowed under the current restrictions on P in fertilizers in Virginia (VA Code 3.2-3600 and 10.1-104.5). Simulations of runoff events immediately following fertilizer applications have shown substantial losses (<1 to 18 percent of the fertilizer applied), and in worst-case simulations contributed a year’s worth of P from a single event (Soldat and Petrovic, 2008). Higher amounts and concentrations of particulate P in winter and spring months may result from elevated streamflows capable of entraining P bound to sediments and organic matter as well as seasonally lower rates of P dissimilation.

Total phosphorus concentrations peaked within one hour of peak streamflow in 83 percent of storms, indicating concentration is well correlated with streamflow (fig. 23). Positive lags greater than one hour were rare, though TP peaked before streamflow in 61 percent of storms. This pattern suggests a disproportionate mass of P typically was transported on the rising limb of the storm hydrograph. Median station ΔMax values were less than the 0.2 first flush threshold at all stations except Ramsgate (0.21), where this threshold was exceeded during 54 percent of storms (table 7; fig. 24). The threshold was exceeded in less than 50 percent of events at the other 11 stations. Although this strict criterion for a first flush was rarely met, M(V) curves confirm results of cross-correlation analyses that indicate biased loading towards the beginning phase of the runoff response. On average, approximately 60 percent of the TP load was transported by the first 46 percent of stormflow volume. Delta-max scores were highest in COM watersheds, but not statistically different between the two residential types of land uses. Phosphorus loadings typically are well correlated with TSS given an affinity to bind to sediments and organic matter. Differences in the frequency of the first flush phenomenon between these two constituents likely reflects contributions of dissolved P loadings and indicates differences in the source and timing of inputs. Storage of water-soluble P in the thin zone of surface soils and organic matter can desorb and enter waterways as bioavailable PO₄^3- (Sharpley and others, 1981). Unlike with sediment wash off from impervious surfaces or resuspension of previously deposited particulates within the stormwater system, this process requires soil saturation, and consequently is less likely to occur during the initial flush of runoff.

Total phosphorous concentrations and loads generally were well predicted by models that contained turbidity and streamflow terms. Inclusion of these terms and the strength of their predictive power reflects the strong relation between sediment and phosphorous transport; on average, models explained 96 and 68 percent of variability in load and concentration, respectively (app. 1, table 1.9). On average, models were efficient (Nash-Sutcliffe efficiency = 0.82) and had low bias (2 percent), which indicates accurate predictions of TP loads. Half of the models included a base-flow term to adjust the y-intercept for differences in low-flow versus high-flow sample populations, which suggests differences in the source of P depending on hydrologic condition. Models for most stations included either water temperature or a seasonal term, indicating seasonal variance in P inputs that alter the TP-TB relation. Orthophosphate models had lower predictive power than TP models, explaining, on average, 92 and 49 percent of the variability in load and concentration, respectively (app. 1, table 1.10). Similar to findings in previous studies (Porter and others, 2020; Robertson and others, 2018), model bias in this study was higher in PO₄^3- models (average 12 percent) and efficiency of predictions was lower (average Nash-Sutcliffe efficiency = 0.67); however, these diagnostics are within the range of acceptable results for use in load computation (Lorenz and others, 2015; Nash and Sutcliffe, 1970). Model coefficients and unit-value estimates of load and yield are published in Porter (2022).

Total phosphorus and PO₄^3- yields varied across stations, though no significant effect of land use was observed for TP (fig. 25; app. 2, fig. 2.2). Over the study period, median annual TP yield was highest at Ludlow and USAA, and networkwide ranged from 0.62 to 2.39 lbs/acre. The average annual yield across all stations was 1.37 lbs/acre (table 8). The median annual PO₄^3- yield at Lucas Creek and Ludlow were approximately double and triple the network median, respectively, and yields ranged from 0.12 to 1.57 lbs/acre networkwide (app. 1, table 1.11). Annual streamflow yield was moderately correlated with TP yield but poorly with PO₄^3- (app. 2, fig. 2.1B and 2.1C). Nonetheless, TP and PO₄^3- loads were transported predominantly during periods of stormflow (> 80 percent; fig. 19). This apparent contradiction suggests the recurrence of storms was a bigger driver of annual yield than the magnitude of those events or the volume of water exported from the watershed each year. The two most likely transport pathways of P during runoff events are erosion of sediment-bound P and dissolution of soluble P stored in soils and organic matter, both on the land surface and within the stormwater conveyance system. The portion of P transported during base flows is likely dominated by PO₄^3- and dissolved organic P leached from inputs on the land surface to the underlying water table (Yang and Toor, 2018). When buried storm conveyance systems intercept the water table, these systems provide preferential flow paths (urban karst) as groundwater moves through leaks, cracks, and joints in the concrete infrastructure, and along the highly permeable gravel-filled trenches that commonly surround these systems (Bonneau and others, 2017; Kaushal and Belt, 2012; Sharp and others, 2003). The effect of urban karst is reduced attenuation and removal of phosphorus and other soluble constituents through biogeochemical cycling and sorption, leading to their increased transport.

Mean annual yields of TP and PO₄^3- from the Hampton Roads study watersheds were higher than those documented in the monitoring networks in comparison studies (fig. 26) as well as estimations for other developed watersheds in the North Carolina-Virginia coastal region (Dodd and others, 1992; Lin, 2004). Previous studies of urban watersheds in Virginia linked high rates of TP transport to similarly high yields of suspended sediment (Jastram, 2014; Porter and others, 2020). The relatively low TSS yields reported in these 12 Hampton Roads watersheds indicate that non-sediment related sources of P, such as organic matter and PO₄^3-, may dominate in this region. The primary driver of differences in TP yields between the Hampton Roads and Fairfax County networks was elevated levels of PO₄^3-, rather than particulate P, in Hampton Roads. Yields in Hampton Roads may be intrinsically higher than those found in comparison studies because of factors related to scale—stream order and watershed size. Larger watersheds and higher-order streams typically have smaller sediment and nutrient yields than small headwater streams, as a result of fluvial trapping of sediments and those nutrients (Smith and Wilcock, 2015). Furthermore, P retention may be limited by the high drainage density of the Hampton Roads watersheds, which have been artificially expanded by the engineered stormwater network of ditches, concrete channels, and buried pipes.

Phosphorus yields also may be elevated as a consequence of soil properties and geological features of the Coastal Plain PP. The degree of phosphorus retention in Coastal Plain PP soils ranges from low to moderate (USDA-NRCS, 1998) and is influenced by the P saturation level and soil texture. The Spatially Referenced Regression on Watershed attributes (SPARROW) model has similarly identified high rates of P transport in the Coastal Plain, and researchers have attributed this to soil saturation resulting from historical P applications (Ator and others, 2011). When soils become P saturated, retention is limited, and applied P may be readily lost through runoff and leaching. Retention also is influenced by soil texture; soils with high clay content have a greater P sorption capacity than do sandy soils, because of grain size and the net negative charge of the soils, which affects surface-area availability and the ability to form bonds with positively charged nutrients (Ballard and Fiskell, 1974; Olsen and Watanabe, 1957; Xu and others, 2006). Soils in the Coastal Plain PP are characteristically coarse-textured sands and silts, particularly in the upper soil horizons, as compared to the fine clays common throughout the Piedmont PP (Sanford and others, 2012; Thomas and Harper, 2009). Furthermore, Coastal Plain PP soils are more dominated by Al oxides than by Fe oxides, whereas the inverse is true of soils in the Ridge and Valley and the Piedmont PP (Penn and others, 2006). Previous studies have reported a negative correlation between dissolved P in runoff and the clay and Fe content of soils (Cox and Hendricks, 2000; Penn and others, 2006). This relation reflects the formation of stronger bonds between P ions and Fe oxides compared to those with Al oxides; therefore, Fe-associated P is less likely to desorb in runoff (Aura, 1978). Given these factors, Coastal Plain PP soils may have a greater potential to release dissolved phosphorus in runoff compared to the upland PPs of Virginia, given identical inputs of soluble P to the landscape (Penn and others, 2006).

In the Coastal Plain PP, the proximity of the surficial aquifer to the land surface, and water table that lies within it, makes it susceptible to contamination by anthropogenic inputs (McFarland and Bruce, 2006) and may promote lateral transport of dissolved P in runoff and groundwater discharging to streams and stormwater conveyances (Ballard and Fiskell, 1974; Harris and others, 1996). Elevated P yields also may be related to high organic matter content that is characteristic of the Coastal Plain PP, a product of the low topographic gradient and resulting low-velocity streams. Organic matter can affect available P in two ways: (1) microbial mineralization of organic materials such as plant debris and animal waste and the ensuing release of previously adsorbed P (Upreti and others, 2015); and (2) coating of sorption sites on metal oxides. Negatively charged substances in organic matter can interact with positively charged oxides, such as Al and Fe, creating competition for binding sites and consequently reducing the potential for P adsorption (Liu and others, 1999) and increasing desorption (Hinsinger and others, 2011).

Nitrogen

Nitrogen (N) is a naturally occurring element that is essential to the growth of plants and animals, but inputs of N exceeding ecosystem processing and storage capacities may be exported to downgradient waters (Hem, 1985; Paul and Meyer, 2001). Excess N can cause eutrophication, excessive algal growth, hypoxia, reduced water clarity, and harmful algal blooms, which degrade aquatic ecosystems and may pose risk to human health (Galloway and others, 2003; Hem, 1985; Hobbie and others, 2017). Common anthropogenic sources of N in the urban landscape include atmospheric deposition, lawn fertilizer, leaking wastewater infrastructure, septic systems, and pet waste (Bettez and Groffman, 2013; Hyer and others, 2016; Paul and Meyer, 2001); natural sources include plant decomposition and biological fixation (Carpenter and others, 1998; Wanielista and others, 1977). The export of N associated with these inputs is positively correlated with urban development (Shields and others, 2008; Tasdighi and others, 2017), though N retention varies with watershed size (Groffman and others, 2004).

Nitrogen loading in urban watersheds is composed of both inorganic (nitrate [NO₃^-], nitrite [NO₂^-], ammonia [NH₃], and ammonium [NH₄]) and organic (both dissolved and particulate) forms, though many studies report organic N as the primary component in urban runoff (Jani and others, 2020; Lusk and Toor, 2016; Seitzinger and others, 2002; Seitzinger and Sanders, 1997; Taylor and others, 2005). Leaf litter, grass clippings, and other organic detritus and particulate matter are the primary sources of organic N in urban watersheds (Bratt and others, 2017; Janke and others, 2017; Selbig, 2016). Organic N, previously considered largely biologically unavailable, has been reported in several recent studies to be a significant source of reactive N for the phytoplankton and cyanobacteria linked to water-quality impairment (Bratt and others, 2017; Jani and others, 2020; Petrone, 2010; Wiegner and others, 2006). Coastal Plain PP streams typically have elevated concentrations of dissolved organic N, a consequence of characteristically high rates of organic matter accumulation, but relatively low concentrations of dissolved inorganic N (Gold and others, 2019). In urban streams, inputs of dissolved inorganic N may be elevated by anthropogenic activities, and changes in watershed hydrology following urban development have been shown to increase loadings of all forms of N (Gold and others, 2019; Petrone, 2010). A nitrogen criterion specific to Virginia streams has not been established; however, the EPA reference condition for streams in Hampton Roads (ecoregion 63) is 0.87 mg/L for TN and 0.04 mg/L for NO₃^- (U.S. Environmental Protection Agency, 2000b).

The network median TN concentration in base-flow samples (1.28 mg/L) was similar to the median stormflow concentration (1.43 mg/L) but the concentration and composition of TN varied between stations (app. 1, table 1.4). Concentrations were highest at Lindsley (3.43 mg/L) and Lucas Creek (3.93 mg/L), but typically were low (approximately 0.80 mg/L) at the three COM stations. Spatial variance in TN was moderated during stormflows (range = 1.07 mg/L), reflecting low inter-station variability of NO₃^- and TON in runoff, and highlighting the occurrence of NO₃^- in groundwater discharges as an important driver of overall differences between stations. The network median NO₃^- concentration of 0.46 mg/L in base-flow samples was similar to the median concentration in the surficial aquifer of the Coastal Plain PP (0.2 mg/L), reported by McFarland (2010). Some wells sampled by McFarland (2010) that were located on agricultural lands exceeded 10 mg/L of NO₃^-; however, most of the 482 wells measured were less than 1 mg/L, similar to the monitored Hampton Roads stormwater conveyances.

Concentrations of total N varied significantly across the three land-use types—highest in SFR watersheds and lowest in COM watersheds (fig. 27). This pattern may reflect differences in the type or quantity of N inputs across these land-use types or greater dilution across the impervious gradient, a function of differences in stormflow yields. The median TN concentration for the 12 study watersheds (1.41 mg/L) was lower than that for NSQD residential and COM stations (2.63 mg/L); TN concentrations were not available for the subset of NSQD stations within the Hampton Roads study area. Likewise, mean TKN and NO₃^- concentrations in samples collected in this study were lower than those in NSQD data collected in the Hampton Roads area and nationwide (app. 1, table 1.6). Similar to values from this study, TN, TKN and NO₃^- concentrations in the NSQD were higher in watersheds draining residential than COM land uses.

Total nitrogen was composed of approximately 70 percent TON, 20 percent NO₃^-, and 10 percent NH₃ across all three land-use types. Concentrations of NO₃^- and NH₃ were lowest in COM watersheds, and TON concentrations varied across the full land-use gradient (fig. 28). Composition of TN also varied with hydrologic conditions; during base flow, TON constituted only 29 percent of TN but became 69 percent of TN during stormflow. In contrast, NO₃^- accounted for 57 percent of TN in base flow and only 15 percent in stormflow. This compositional shift signals a change in the source of N from NO₃^--enriched groundwater discharges during dry weather periods to TON dominant runoff during storm events. On average, TN was not well correlated with streamflow (r = -0.002, p < 0.921), a consequence of the contrasting positive relation to TON and negative relation to NO₃^- that cancel out when aggregated as TN. In most watersheds, which were more strongly dominated by TON, weak but significant positive correlations exist (app. 1, table 1.5). Conversely, at Lucas Creek and Lindsley, where NO₃^- is the dominant form of nitrogen, the negative relation suggests dilution with increasing streamflow.

The concentrations of the different forms of nitrogen varied seasonally, with the highest concentrations of TN in the spring, driven primarily by seasonal fluctuations in TON (app. 1, table 1.12). Organic N typically is bound to sediments and organic matter; therefore, elevated concentrations in springtime may be related to similarly high concentrations of TSS. Seasonal differences in NO₃^- were minimal and suggest year-round contributions from one or more sources, or low denitrification potential in local soils. Terrestrial denitrification (reduction of NO₃^– to atmospheric N²) often is the dominant nitrogen loss term in N budgets of mid-Atlantic watersheds (Van Breemen and others, 2002), and typically results in greater NO₃^– loss in warmer months because of increased heterotrophic respiration (Deek and others, 2010; Meng and others, 2018; Porter and others, 2020; Tufford and others, 1998). Rates of NO₃^- loss vary in the Coastal Plain PP as a result of heterogenous lithology but are low in watersheds with predominantly sandy sediments where water remains oxic (Ator and others, 2011; Ator and Denver, 2015). Further, concrete stormwater infrastructure largely disconnects streamflow from the hyporheic zone, where microbial mediated processes such as denitrification occur (Merill and Tonjes, 2014). Denitrification potential was not measured as part of this study but may be a relevant factor driving the lack of a seasonal variation in NO₃^- concentrations across the network. The lack of seasonality in NO₃^- may also be related to the absence of periphyton and other algae in these stormwater systems. In natural streams, these organisms play an important role in improving water quality by temporarily sequestering N and P during periods of high growth, but owing to a lack of sunlight, are unlikely to thrive in subterranean stormwater pipes (Santhana Kumar and others, 2017).

Concentrations of TN typically peaked within one hour of peak streamflow (77 percent of storms; fig. 29). The timing of peak concentration reflects the dominant particulate composition of TN in stormflows and suggests that N transport may occur via similar pathways as do TSS and TP during periods of runoff. Total N peaked before streamflow in 72 percent of storms, suggesting a disproportionate mass of TN typically is transported on the rising limb of the storm hydrograph. This occurs from either the remobilization of deposits within the storm conveyance system or from sources nearby or directly connected to sewer inlets. Longer positive lags occurred in COM and HDR than in SFR watersheds, which may result from differences in the areal coverage and connectivity of impervious surfaces to stormwater infrastructure. These landscape features typically promote a more rapid first flush of constituent loadings to the stormwater conveyance system because of the directness of the flow path (Li and others, 2007). Exceptions were observed at Lucas Creek and Lindsley, where a negative lag occurred (TN peaked after streamflow) in more than 20 percent of storm events. Base-flow NO₃^- was substantially higher in these two watersheds, and a pattern of negative lags suggests (1) delayed contributions of NO₃^- from groundwater or interflow; or (2) transport of N sources farther from the monitoring station.

Although M(V) curves confirm that loadings of N were greater in the beginning phase of the runoff response, strict criterion (ΔMax ≥ 0.2) of a first flush was rarely met (table 9; fig. 30). On average, the ΔMax criterion was exceeded only about 20 percent of the time across all 12 stations, mobilizing about 50 percent of the TN load in the initial 39 percent of stormflow volume. Delta-max scores were highest in COM watersheds, but not statistically different from scores for the two residential types of land use. Relatively low occurrences of a first flush for TN suggest that sources are less directly connected to the stormwater system. For example, NO₃^- and dissolved organic N may be flushed out of soils in residential lands and onto impervious surfaces connected to stormwater infrastructure; however, this process first requires soil saturation. Further, loadings of the various forms of N—dissolved versus particulate and organic versus inorganic—can be staggered in time and therefore reduce the first-flush effect. For example, previous studies have reported a rapid first flush of particulate organic N followed by decreasing concentrations at peak streamflow because of dilution. Concentrations of dissolved N (primarily NO₃^-) then peak on the receding limb of the storm hydrograph and remain elevated for several subsequent days after the event, responding to higher inputs from interflow (Fletcher and others, 2004; Taylor and others, 2005).

The predictive power of TN load and concentration models developed in this study was weaker than that of models for TP or TSS but stronger than the predictive power of similar models for urban stream monitoring stations in Virginia (Porter and others, 2020; Jastram, 2014). The relative strength of these models likely is related to the high proportion of particulate organic nitrogen, which is well predicted by streamflow and turbidity terms. Most models also contained an SC term, a proxy for dissolved solids (such as NO₃^-), and a seasonal term to account for temporal fluctuations of inputs and speciation. Total nitrogen models accounted for on average 95 percent of the variability in loads and 53 percent in concentrations (app. 1, table 1.13). Total Kjeldahl nitrogen, TON, and NO₃^- models explained a similar degree of variance in loads and slightly greater variance in concentrations (> 60 percent; app. 1, tables 1.14, 1.15, and 1.16). On average, TN, TKN, TON, and NO₃^- models were efficient (Nash-Sutcliffe > 0.7) and had low bias (< 3 percent), suggesting accurate predictions of loads. Model coefficients and unit-value estimates of load and yield are published (Porter, 2022).

Yields of TN, TKN, and TON varied across stations and with land-use types, though no effect of land use was observed for NO₃^- yields (fig. 31; figs. 2.3, 2.4, and 2.5). The average annual yield across all stations for TN, TKN, TON, and NO₃^- was 10.2 lbs/acre, 7.3 lbs/acre, 6.3 lbs/acre, and 3.0 lbs/acre, respectively (table 10; tables 1.17, 1.18, and 1.19). Over the study period, median annual TN yields were highest at Professional (14.8 lbs/acre) and USAA (14.6 lbs/acre), despite having two of the three lowest median annual concentrations; the lowest median yields occurred at Garrett and Lindsley (5.9 lbs/acre). In COM watersheds, low concentrations are offset by high rates of runoff, producing elevated constituent yields. Total nitrogen yields determined from samples collected in this study were composed primarily of organic N; therefore, the distribution of annual TKN and TON yields across the 12 stations were similar to those of TN. Total nitrogen, TKN, and TON yields were strongly correlated to annual streamflow yield (figs. 2.1D, E, and F), highlighting the importance of particulate wash off in the transport of nitrogen in all its forms. These patterns are consistent with those found in other studies of urban streams in which substantial decreases in N retention were observed between dry and wet years (Duncan and others, 2017; Jastram, 2014; O’Driscoll and others, 2010; Porter and others, 2020). Nitrate was poorly related to annual streamflow yield (app. 2, fig. 2.1G). Annual streamflow yields were predominantly affected by the volume and recurrence of stormflows, whereas 53 percent of the annual NO₃^- load, on average, was conveyed during base flows.

The median annual TN yield for the 12 monitoring stations in the Hampton Roads study (10.0 lbs/acre) was similar to that for monitored streams in Fairfax County, Virginia (8.5 lbs/acre), and significantly higher than median values for monitored streams in Gwinnett County, Georgia (5.4 lbs/acre), and the CB-NTN stations (5.4 lbs/acre; fig. 32). Conversely, NO₃^- yields were significantly lower in Hampton Roads (2.7 lbs/acre) than in the monitored Fairfax County streams (4.4 lbs/acre) and the CB-NTN (3.3 lbs/acre), but similar to those in Gwinnett County (2.5 lbs/acre). Yields of neither TON nor TKN were available from the comparison networks; however, average annual TN and NO₃^- yields from those networks can be used to infer regional differences in the composition of N loads. As with other constituents, care must be taken when comparing yields across these networks because of dissimilarities in stream order and watershed size. Elevated TON yields in Hampton Roads may also highlight an important difference in carbon processing between buried stormwater conveyances and natural earthen streams. Natural streams typically consist of a series of erosional and depositional zones. Depositional areas commonly serve as “hotspots” for the microbial decomposition of allochthonous materials and transformation of nutrients (Kaushal and Belt, 2012). Traditional engineered stormwater systems are designed to discourage the deposition and retention of materials to prevent the buildup of materials that reduce flow capacity and efficiency. Further, the lack of light and natural habitat is not conducive to microbial communities or benthic macroinvertebrates, which can process allochthonous nutrients and carbon inputs (Santhana Kumar and others, 2017). These factors may result in greater export of organic matter and organic N from lands drained by engineered stormwater conveyances than by natural streams.