Upstream Basin Characteristics

The dataset of upstream basin characteristics was developed by delineating stream basins upstream from the intersections between roads and streams in southern New England and analyzing selected basin properties by using GIS software. Spaetzel and others (2020) generated this dataset of basin properties above road-stream crossings by using the intersections of roads as defined by the USGS National Transportation Dataset and streams as defined by the USGS StreamStats modified National Hydrography Dataset. Although all the detected road crossings were within southern New England, the delineated basins include basin properties and roadway characteristics from areas of New York, Vermont, and New Hampshire that are within the delineated basins. The selected basin characteristics are drainage area (in square miles), 10-85 slope (in feet of elevation change per mile), longest flow path (in miles), number of road crossings by road class, impervious cover (in percent), length of roads by road class (in miles), and length of streams (in miles; Spaetzel and others, 2020). The 10-85 slope and longest flow path in this dataset correspond to the slope and main-channel length in SELDM. The number of road crossings, length of roads by type, and impervious cover were determined to assess variations in the magnitude of development in southern New England stream basins. The length of streams was determined so that it could be used with drainage area to estimate the stream density (in miles per square mile). The stream density commonly is used to develop streamflow estimates (Bent and others, 2014), and it can be used to estimate the length of overland flow from drainage divides to tributary stream channels within the upstream basin (Horton 1945; Carlston 1963; Jeznach and Granato, 2020).

The geographic analysis by Spaetzel and others (2020) resulted in 53,131 basin delineations in southern New England with drainage areas ranging from 0.000032 to 1,938.8 square miles. Because delineation of very small basins and determining their characteristics is highly uncertain, the 48,466 basins delineated upstream from paved roads with a minimum drainage area of 0.025 square miles (16 acres) and a main-channel length, main-channel slope, and a drainage density greater than zero were selected for further analysis (fig. 2). Only 5,545 (about 11 percent) of these basins were delineated upstream from arterial roads (The National Map Functional Road Classification definitions in table 1); these basins will be described herein as the arterial-upstream basins. About 22 percent of these 48,466 basins contain one or more upstream arterial road crossings.

The road-crossing basin count may seem large for southern New England, but most delineated basins are nested within larger basins. For example, the Blackstone River Basin above Interstate 95 in Providence, Rhode Island is 475 mi², has an imperviousness of about 12 percent, and has 2,340 upstream road crossings; the Charles River Basin above Interstate 93 in Boston, Massachusetts is 313 mi², has an imperviousness of about 23 percent, and has 1,365 upstream road crossings; and the Park River Basin above Interstate 91 in Hartford, Connecticut is 77.2 mi², has an imperviousness of about 27 percent, and has 539 upstream road crossings. These delineated basins, by definition, do not represent the confluence of tributary streams, but the drainage-area pattern is similar to Giusti’s law, which indicates that the number of upstream basins of any size is about 0.3 times the ratio of the basin area to the selected subbasin area (Giusti and Schneider, 1965). For example, this relation would indicate that a 250-square-mile basin would be expected to have 30, 300, and 3,000 tributary stream subbasins with drainages areas of 2.5, 0.25, and 0.025 mi², respectively.

The characteristics of southern New England stream basins are shown in figure 2 and table 3. The values of basin characteristics for all the basins and the arterial-upstream basins vary by almost three (for basin length) to almost five (for drainage area) orders of magnitude. The coefficient of variation (COV) for drainage areas, which is the standard deviation divided by the average, also indicates the large variability of the basin properties. For example, the COV of the drainage areas is 6.65 for all basins and 5.04 for arterial-upstream basins. Although the values are wide ranging, most basins have small drainage areas. For example, the median, average, and geometric mean drainage areas for all basins are 0.455, 7.65, and 0.6 square miles, respectively. The median, average, and geometric mean drainage areas for arterial-upstream basins are almost twice the size with values of 0.721, 22.0, and 1.11 square miles, respectively.

Information about relations between basin properties is needed to guide the choice of a limited but representative set of values for simulating the potential effect of runoff on receiving waters. To this end, an analysis of correlations between basin properties was done by calculating the nonparametric rank correlation coefficient (Spearman’s rho) and the product-moment correlation coefficient (Pearson's R) for the logarithms of data (table 4). Spearman’s rho is calculated by ranking the data and calculating the correlation coefficients between the rank values rather than the data values (Haan, 1977; Helsel and Hirsch, 2002). Spearman’s rho indicates the strength of the relation regardless of the linearity of the relation between variables. Correlations among the drainage area, the main-channel (10-85) slope, the main-channel length, and the basin-lag factor (BLF) are moderate (correlation coefficients with an absolute value greater than or equal to 0.5 and less than 0.75) to strong (correlation coefficients with an absolute value greater than or equal to 0.85, table 4). Drainage density and imperviousness of the basins may be considered to be random variables with respect to the other basin properties because they have weak correlations (correlation coefficients with an absolute value less than 0.5) with all the other basin variables.

Correlations for the basin-lag factor (BLF), which is the main-channel length divided by the square root of the main-channel slope, also were calculated (table 4) because the BLF is the controlling variable used to calculate the basin lagtime that determines the timing of runoff from the upstream basin (Granato, 2012, 2013). Although correlations of the BLF to length and slope are strong because the BLF is a function of length and slope, the correlation between the BLF and drainage area indicates the potential for using drainage area as the master variable for other basin properties. The correlations between drainage area and BLF in this study are very similar to the correlations calculated by Granato (2012) using National datasets with hundreds of sites. Although the correlations between drainage area and main-channel slope are only moderately strong (correlation coefficients with an absolute value greater than or equal to 0.5 and less than 0.75), correlations between drainage area and length are strong. This indicates that potential relations between drainage area and main-channel slope are less influential than the relation between drainage area and length for simulating basin lagtimes in southern New England. Granato (2012) also determined that drainage area was almost as strong a predictor for basin lagtime than the BLF and imperviousness, which further indicates that drainage area is the master variable for simulating the timing of stormflow from the upstream basin.

Regression relations were developed to select representative values of main-channel length and slope from drainage area (fig. 3, table 5). Because the potential effects of high-leverage outliers in datasets ranging over several orders of magnitude on regression relations can be large, the Kendall-Theil robust line method (Granato, 2006) was used to develop these equations. Because use of the Kendall-Theil robust line method on the full 48,466 basin dataset would require the calculation of about 1.2 billion slopes in arithmetic and 1.2 billion slopes in logarithmic space, the full dataset is too big to process in the KTRLine software (Granato, 2006). Therefore, a subsample of 6,923 basins was used to develop the regression equations. The subsample was created by sorting the dataset by basin size, and the data from every seventh basin out of the 48,466 were selected. The basins were selected if the remainder, or modulus, of the index number divided by 7 equals zero. Experiments conducted by shifting the index number and repeating the regression analysis indicated that the regression equations in table 5 are representative of the whole dataset.

The total impervious area (TIA) is an important variable for simulating runoff because it is used to calculate runoff coefficients and basin lagtimes in SELDM (Granato, 2010, 2012, 2013; Granato and Jones, 2014; Jeznach and Granato, 2020). The TIA of the delineated basins ranges from 0 to 85.8 percent with a median of about 3.84 percent among all basins and 10.43 percent among the arterial-upstream basins (fig. 2, LC16IMP in table 3). About 45 percent of all delineated basins and 66 percent of arterial-upstream basins exceed the TIA threshold of 5 percent that is commonly used to indicate the lower limit of substantial stream ecologic degradation (Jeznach and Granato, 2020). About 18 percent of all delineated basins and 30 percent of arterial upstream basins exceed the TIA threshold of 20 percent that is commonly used to indicate complete degradation of natural stream ecology. Because correlations between TIA (LC16IMP) and other basin variables are very weak (table 4), this variable must be considered as a random variable with respect to the other basin properties.

The stream density, which is the length of all streams in the basin divided by the drainage area, has a smaller range than the other basin characteristics in this study, and the differences between stream density for all the basins and the arterial-upstream basins are relatively minor. One-half of the reciprocal of the stream density can be used to estimate the length of overland flow from drainage divides to tributary stream channels; this estimate is known as the Horton half-distance (Horton 1945; Carlston 1963; Jeznach and Granato, 2020). The average stream density in the study area is 2.52 miles per square mile, and the reciprocal of this value is about 0.4 mile. Therefore, the Horton half-distance for overland flow in southern New England would be about 0.2 mile, or 1,056 feet.

Highway Site Characteristics

SELDM is nominally a highway-runoff model, but it can be used to simulate runoff for any site of interest by using the characteristics of the site of interest and representative water quality. Because SELDM is a lumped-parameter model, the basin characteristic values chosen for the highway, urban, or other developed areas that are simulated as the site of interest can be literal or interpretive (Granato, 2013; Stonewall and others, 2019; Jeznach and Granato, 2020). A literal site would be simulated by using the particular characteristics of an individual drainage pathway. For example, a literal site may be a section of roadway draining to a stream or a developed area with a trunkline drainage system. The basin characteristics for a literal site may be derived from actual drainage plans or estimated by using online tools like StreamStats and Google Earth. Interpretive sites may have multiple drainage pathways. Interpretive sites are used to simulate the net effect of multiple outfalls on the receiving water quality downstream from a point of interest (Granato, 2013; Granato and Jones, 2017a; Stonewall and others, 2019; Jeznach and Granato, 2020). For example, an interpretive site may represent a bridge (or two highway bridges in parallel) with many individual scupper drains, a bridge with two approach sections that discharge to a stream through multiple outfall locations, a road paralleling a stream with multiple outfall locations, or an agglomeration of developed areas that drain to a point of interest along a stream. An interpretive site may be simulated by selecting representative basin properties that produce the volume and timing of runoff characteristic of the entire simulated drainage area. Because decisionmakers commonly need information about the net effect of multiple stormwater outfalls on the receiving water quality at a given location, interpretive sites commonly are used to simulate stormwater quality (Granato and Jones, 2017a; Smith and others, 2018; Stonewall and others, 2019; Weaver and others, 2019; Jeznach and Granato, 2020).

If simulations are done to develop annual total maximum daily load (TMDL) yields, then the timing of runoff during individual events is not of concern and the site may be simulated by using an area of 1 acre and a representative TIA value; the remaining basin properties may be specified by using generic values (Granato and Jones, 2017b, Stonewall and others, 2018; Granato and Friesz, 2021a). For TMDLs, the yields can be applied to the areas of different road classes and to the developed impervious or land-use areas to estimate loads from simulated yields (Granato and Jones, 2017b, Stonewall and others, 2018; Lantin and others, 2019; Granato and Friesz, 2021a).

In this study, available information about roadway geometry and drainage-system characteristics were used to simulate runoff from hypothetical, but representative sites. Runoff from roadways was simulated by using the paved area rather than the area of the entire right-of-way because roadway-runoff quality data collected in southern New England were collected from paved areas (Smith, 2002; Smith and Granato, 2010; Smith and others, 2018; Granato, 2019a), and the grassy swales and strips of the shoulders and medians alter the flows, concentrations, and loads from pavements in ways that are unique to each site (Granato, 2014; Taylor and others, 2014; Granato and others, 2021). Roadway geometry and drainage characteristics were estimated by using the AASHTO (2001) policy on geometric design of highways and streets and Federal Highway Administration Hydraulic Engineering circulars (Young and others, 1992; Federal Highway Administration, 1993, 2013). The State DOTs in southern New England follow these standards with minor modifications (Massachusetts Highway Department, 2004, 2006; Rhode Island Department of Transportation, 2008; Connecticut Department of Transportation, 2020).

The road data incorporated into StreamStats provides information about the lengths of various road classes above any given point on a stream, but information about road widths is needed to estimate the drainage areas of roads within a delineated basin. The AASHTO (2001) guidelines indicate that roadway widths are commonly specified by road class, rated speed limits, and traffic volume. Individual travel-lane widths commonly range from 9 to 12 feet. Safety is the prime consideration for travel-lane widths, but there are many design considerations such as access for pedestrians or bicyclists and parking that may come into play. Although 12-foot lane widths have become a design standard, 9-foot lanes are considered acceptable for low speed, low volume rural and residential roads. Lane widths less than 12 feet may be legacy widths or the result of right-of-way constraints in urban areas. The width of each roadway also includes the width of paved shoulders. A minimum shoulder width is 1 foot for drainage, 2 feet to help protect the integrity of the pavement edge, and up to a full 12-foot shoulder to permit emergency parking along multilane limited-access roadways. A minimum width of 4 feet is the design standard for the shoulder where vertical barriers or guard rails are present. The minimum parking-lane width for residential areas commonly is 8 feet. Parking lanes may be 10 feet wide on connecting roads and full access arterial roadways. Bicycle-lane design standards are a minimum of 4 feet wide in open areas and 6 feet in commercial areas. The desired width is 8 feet to accommodate multiple bicyclists within the lane. In southern New England where rights-of-way are constrained, bicycle lanes commonly are created by reducing motor vehicle lane widths rather than widening the paved roadway area.

The area of paved roads is needed to calculate the road-runoff flows and loads. Typical road widths may be estimated based on the AASHTO (2001) guidelines and information about the number of lanes from the road census (Federal Highway Administration, 2022a, b) and NBI (Federal Highway Administration, 2020) for southern New England listed in table 1. On average, local roads and minor collectors have 2 travel lanes. Major collectors have an average of 2.1 travel lanes, which indicates that about 94 percent of these roads have 2 travel lanes and 6 percent have 4 or more lanes. On average, minor arterial roads have 2.2 to 2.5 travel lanes indicating that about 75 to 90 percent of these roads are 2-lane roads. On average, principal arterial roads have 2.6 to 3 travel lanes indicating that about 50 to 70 percent of these roads are 2-lane roads. Road widths of principal arterials are similar to minor arterials, but there is a greater proportion of principal arterial roadways with multiple lanes in each direction. The divided limited-access highways including freeways and expressways and interstates were combined in table 1 to calculate the number of lanes in both directions. The average number of lanes for freeways and expressways are about 4.2 lanes, indicating that about 90 percent of these road types have 2 travel lanes in each direction within southern New England. Using these estimates of width and the number of lanes for each road type, the estimated pavement area may range from about 2.4 acres per mile for a 2-lane local road without roadside parking to 7.8 acres per mile for each roadway of an 8-lane limited-access arterial roadway (table 6).

Because the cost of building and maintaining drainage systems to manage runoff are large, direct drainage to the local land surface is used where possible for infiltration. Highway drainage-design guidelines specify use of grass strips and swales rather than storm sewer systems wherever practical (American Association of State Highway and Transportation Officials, 2001). Consequently, only a small part of the road network may drain directly to receiving waters. Available information (Granato and others, 2022) on the impervious areas draining to roadway stormwater-conveyances in Massachusetts and Rhode Island indicates that the distributions of designed drainage networks in these States are similar to each other (fig. 4A). The values shown for Massachusetts are the percent distribution of impervious roadway area sizes contributing to the structural stormwater best management practices, and the values shown for Rhode Island are the percent distribution of impervious roadway area sizes contributing to stormwater conveyances, which may include areas draining to a single stormwater inlet. These roadway drainage areas range from about 0.01 to 32 acres, with a median of 0.63 acres for Massachusetts and 0.57 acres for Rhode Island (fig. 4A). These distributions are similar, which may be the result of the similarities in hydrology within southern New England and use of hydrologic design guidelines based on National standards.

The NBI (Federal Highway Administration, 2020) also provides information that can be used to estimate roadway-runoff source areas; the NBI can be considered as a random sample of road characteristics across each State. Precipitation that falls on bridges may be directly discharged by using bridge scuppers or routed to stormwater treatment facilities adjacent to the receiving water body (Federal Highway Administration, 1993). Figure 4B shows the distribution of areas of 4,973 roadway bridge decks over water in southern New England that have the bridge deck width and structure length values in the NBI, which are needed to calculate the bridge-deck areas, and the functional class, which is needed to identify different road types (there are another 514 roadway bridge decks over water that do not have all of these values in the NBI). These bridge-deck areas range from 0.007 to 13.98 acres with a median of 0.045 acres. Median bridge-deck areas for local roads, collectors, full-access arterials, and limited-access arterials are 0.0275, 0.0386, 0.0686, and 0.236 acres, respectively. The minimum areas may be smaller than would be expected given the average road widths and areas shown in table 6, but the NBI includes bridges and large culverts, which both may have spans as short as 20 feet.

In SELDM, the length of the drainage flow path is used with its slope to simulate the timing of runoff from the highway or urban site to the point of interest. The site of interest may have two lengths; the physical length for calculating drainage-basin area, and the main-channel drainage length used for calculating the hydrologic basin lagtime for the site. The main-channel drainage length for the site of interest is estimated as the characteristic drainage length that controls the timing of runoff from the drainage divide to the stormwater outfall. For example, when simulating runoff from a bridge with direct-discharge scuppers, the physical length of the bridge may be used to calculate area, and the average distance from the crown of the road to the nearest scupper, which is the average length of the flow path of precipitation on the bridge, may be the hydraulic length. Similarly, on a long stretch of highway with multiple drainages of varying lengths to a parallel stream, the length of that road segment may be used to calculate area, and the average distance from the crown of the road through the drainage system to the nearest outfall may be the hydraulic length. If a highway site stretches across the entire hydrologic stream basin and there is one outfall where it crosses the stream, then the divide-to-divide distance would be used to calculate the roadway area, and the hydraulic length would be the longer distance from the stream to one of the divides. The physical length of highway conveyances can be estimated from the information in table 6 and figure 4. For example, given a median area of about 0.6 acre, the length of a local road without parking can be estimated as about 1,320 feet, and the length of a two-lane full-access arterial road can be estimated as about 930 feet. The length of a divided 4-lane limited access highway may be estimated as about 660 feet if one of the roadways of the highway drains to the conveyance, and 330 feet if both roadways of the divided highway drain to the same conveyance. Similarly to finding physical lengths of parking lot conveyances, the area of a parking lot may be calculated from a physical length, and the hydraulic length may be the route that the main drainage pipe follows from one edge of the parking lot to collect water from the pavement and discharge it to the receiving water body.

Highway drainage slopes can be estimated by using information from roadway design guidelines and hydraulic design circulars (table 7). Highway and drainage design guidelines commonly specify slopes by using percentages or dimensionless ratios, but SELDM uses the watershed slope convention of (vertical) feet of elevation change per (horizontal) mile. Pavement cross slopes may represent the first segment of the flow path from the drainage divide to the drainage-system outfall. Because SELDM uses a representative slope calculated from the elevations of points at 15 and 85 percent up the main channel from the outlet, the pavement cross slopes may not be critical except in the case of a bridge deck for calculating the timing of runoff from the crown of the road to the nearest scupper. The minimum and maximum of longitudinal drainage slopes and road grades may be used to estimate drainage slopes for drainage pathways that follow the highway to the waterway; these slopes commonly represent the longest distance of the flow path. The selected slope may depend on the components of the drainage system. The minimum storm-drain slopes, which are designed to maintain the self-cleaning flow velocity of 3 feet per second, represent minimum slopes for closed drainage systems (table 7). The maximum unlined channel slope specification of 2 percent, which is equivalent to about 106 feet per mile (ft/mi), represents the open channel velocity at which unwanted erosion of roadside swales may begin. The estimated range of roadway drainage slopes in table 7 is about 4 to 900 ft/mi; this is well within the range of 0.046 to 2,186 ft/mi for main-channel slopes of stream basins above of road crossings in southern New England (table 3). When road drainage slopes are being calculated, any vertical drops (such as from the road surface to the catch basin outflow or from a bridge deck scupper or an overhanging drainage outlet to the stream) should not be included in the slope. This is because vertical drops are almost instantaneous and so do not contribute to the basin lagtime of the highway (or urban) runoff drainage pathway.

The drainage characteristics, which include drainage area, length, slope, and imperviousness, for other developed land covers also can be estimated from StreamStats, highway design information, and other sources. SELDM can be used to simulate runoff from a particular site or the upstream drainage areas can be aggregated into a site by lumping the areas and using representative hydraulic characteristics (Stonewall and others, 2019; Jeznach and Granato, 2020). If specific sites are to be simulated, then actual basin characteristics may be derived by using local GIS data (Granato and Friesz, 2021a). The drainage area for simulating runoff from developed areas may be estimated by using the imperviousness, the percent developed area, or land-cover areas (Stonewall and others, 2018, 2019, Jeznach and Granato, 2020; Granato and Friesz, 2021a). Studies of the components of impervious surfaces in developed areas consistently indicate that, on average, off-street parking, roofs, roads, and other anthropogenic surfaces comprise about 35, 32, 25, and 8 percent of the TIA in these areas, respectively (Tilley and Slonecker, 2006; Roy and Shuster, 2009; Wang, 2013). If, as indicated in table 2, the State DOTs own about 15 percent of the road network (about 17.2, 8.1, and 18.3 percent of roadways in Connecticut, Massachusetts, and Rhode Island, respectively), then DOT owned roads would represent about 3.8 percent of the total impervious area in urban areas. However, the percentage of imperviousness composed of local roadways and State-owned roadways may be much higher outside developed areas than in developed areas with a high proportion of off-street parking, roofs, and other anthropogenic surfaces. Granato and Friesz (2021a) determined that the imperviousness of developed areas increase with increasing percentages of developed area because of urban intensification. In southern New England, the areas of State roadways can be estimated from StreamStats (Spaetzel and others, 2020) and subtracted from the total impervious area to produce a more robust estimate of DOT and non-DOT runoff areas. The drainage length may be measured if a literal site is being simulated but must be estimated for interpretive sites. Jeznach and Granato (2020) used the Horton half distance calculated from stream density as the drainage-length to simulate the timing of urban runoff because it is the average distance from the local drainage divide to the nearest stream segment. Roadway and drainage-system slopes in table 7 also may be used to simulate urban runoff because many of the same drainage design constraints influence urban drainage design. As with highway sites, a basin development factor equal to –1 can be specified to use the basin lagtime equation that is based on the imperviousness of the site.

The annual average daily traffic (AADT) volume, which is a count of the number of vehicles using the roadway per day, is commonly viewed as a basin characteristic of roadway sites that is indicative of runoff quality. AADT data is primarily collected to measure and plan roadway capacity needs, but it has been used, with mixed success, as an explanatory variable for estimating highway-runoff quality (Driscoll and others, 1990; Granato and Cazenas, 2009; Smith, and Granato, 2010; Wagner and others, 2011; Granato and Friesz, 2021a). In the National highway-runoff monitoring study by the FHWA (Driscoll and others, 1990) water-quality monitoring sites were categorized as being “rural” if they had an AADT value of less than 30,000 vehicles per day (VPD), and were categorized as “urban,” if they had an AADT greater than or equal to 30,000 VPD based on statistical differences in runoff quality.

State DOTs run small numbers of continuous traffic monitoring stations and supplement these stations spatially by using many more short-period counting locations, which are used to estimate AADT values (Krile and others, 2015). Studies show that the uncertainty in AADT estimates from short-term monitoring stations commonly is on the order of ±20 percent and as high as ±50 percent for low AADT roads (less than 1,000 vehicles per day) and that estimates are highly uncertain for all traffic volumes with measurement durations of less than a full day (Krile and others, 2015).

The NBI (Federal Highway Administration, 2020) was used as a sample of roadway locations in southern New England to assess the AADT population characteristics of roadways near stream crossings. Figure 5 shows the distribution of AADT of the population of all bridges and State-maintained bridges over water in southern New England. The State-maintained bridge population has higher AADTs (median of 11,700 VPD) than the AADTs (median of 3,995 VPD) for the population of all bridges in this area because State-maintained bridges carry higher capacity motorways. About 9.2 percent of all bridge crossings and 20.3 percent of State-maintained bridge crossings in southern New England have AADT values over 30,000 vehicles per day, which is the traditional rural to urban water-quality threshold known as the Strecker number (Driscoll and others, 1990). In comparison, the population of Highway-Runoff Database monitoring sites (Granato and Cazenas, 2009; Granato, 2019a; Granato and Friesz, 2021b) in southern New England has a median of 61,534 VPD, with 69 percent of sites having AADT values greater than 30,000 VPD (fig. 5).

The population of all southern New England bridges in the NBI (Federal Highway Administration, 2020) was used to examine AADT by road class (fig. 6). While the road-class labels in the NBI are slightly different than in table 1 and the roadways on the bridges are identified as rural or urban, the descriptors in table 1 broadly apply to the categories in figure 6. The NBI definitions of rural or urban are based on the census designation for the location of the bridge, rather than being based on the traffic volume or the traditional 30,000 VPD runoff-quality threshold (Federal Highway Administration, 2020). The AADT values are shown as VPD per lane in figure 6 to normalize the values for comparison across roads with different numbers of lanes. Even when normalized for lane count, AADT values increase from category to category and from rural to urban within categories. For example, among rural road classes, the median per-lane traffic volume for interstate arterials is about 68 times the median traffic volume for local-roads (fig. 6A). In comparison, the median per-lane traffic volume for interstate arterials is about 24 times the median traffic volume for local-roads among urban road classes (fig. 6B). Comparison of the urban to rural traffic volumes for similar road classes indicates that the urban road medians are about 1.4 times and 4.7 times the rural road medians for principal arterials and local roads, respectively. Traffic volumes per lane for urban local-roads are comparable to rural collectors; volumes for urban collectors are comparable to rural minor arterials; volumes for urban minor arterials are comparable to rural principal arterials, and volumes for urban principal arterials are comparable to rural interstate highways (fig. 6). Therefore, given the overlapping traffic volumes between urban roads and rural roads, uncertainty in individual AADT values, differences in traffic patterns (more starting and stopping on urban roads), background air quality between urban and rural areas, and other factors, road class and traffic volume may have limited use as a predictor variable for the quality of roadway runoff at the watershed scale (Granato and Friesz, 2021a, b).

Precipitation Statistics

SELDM uses precipitation statistics to stochastically simulate a large series of runoff-generating events. Storm-event precipitation statistics define the characteristics of each simulated storm event and the number of events in the simulation (Granato, 2013, Risley and Granato, 2014, Stonewall and others, 2019, Weaver and others, 2019). SELDM also uses precipitation statistics to aggregate events into annual-load accounting years, which can be used to assess long-term annual loads or yields that can be used for TMDL analyses (Granato, 2013; Granato, and Jones, 2017b; Smith and others, 2018; Stonewall and others, 2018; Lantin and others, 2019; Granato and Friesz, 2021a). SELDM uses the EPA definition of a runoff-generating event, which is based on hourly precipitation values, a minimum precipitation volume of 0.1 inch (in.), and a minimum inter-event period of 6 hours between events (Driscoll, Palhegyi and others, 1989; Granato, 2010, 2013). To simulate the events, SELDM uses the event volume (in inches), duration (in hours), and the time between event midpoints (in hours). SELDM uses the event duration and the time between event midpoints to group random collections of events into the annual-load accounting years; subsequent events are assigned to a year when the accumulated hours equal 365 or 366 days. The number of runoff-generating events per year specified from the selected precipitation statistics is used to calculate the minimum number of events to be simulated in each run (Granato, 2013).

Regional simulations were done by using precipitation statistics for three EPA Level III ecoregions that include parts of Massachusetts, Connecticut, or Rhode Island, and statistics for southern New England (table 9, fig. 7). The ecoregion statistics are the median of statistics for all NOAA hourly precipitation data stations in the ecoregions, which cover areas inside and outside the southern New England States. The statistics for southern New England were calculated as the median of values from precipitation data stations within and adjacent to Connecticut, Massachusetts, and Rhode Island (table 10, fig. 7). Precipitation stations outside but adjacent to southern New England were selected to better characterize conditions within this region even though the additional area of the bounding box reduced the station density in comparison to ecoregion 59 (the Northeastern Coastal Zone), which covers most of southern New England (table 9).

In this study, precipitation statistics for individual NOAA precipitation data stations also were used to do TMDL simulations (a level 2 analysis) to provide information about variations in long-term average yields and to do sensitivity analyses on the effect of variations in precipitation statistics on flows, concentrations, and loads of runoff-constituents of concern within the region. Statistics for the 45 precipitation stations within and adjacent to southern New England are shown in table 10 and are compared to the southern New England median (shown as a red line in figure 8) and the ecoregion medians (shown as black diamonds in figure 8). The median statistics for the Northeastern Highlands (ecoregion 58) are offset from the estimated median of selected southern New England stations in this ecoregion because only a small part of the Northeastern Highlands ecoregion lies within the Southern New England States. The median statistics for the Northeastern Coastal Zone (ecoregion 59) are similar to the estimated median of selected southern New England sites because most stations within the Northeastern Coastal Zone ecoregion are located within or near Connecticut, Massachusetts, or Rhode Island. The median statistics for the Atlantic Coastal Pine Barrens (ecoregion 84) also are similar to the estimated median of selected precipitation stations in this ecoregion within southern New England, even though most of the precipitation monitoring stations in the Atlantic Coastal Pine Barrens ecoregion are located in New Jersey and New York. Among the 45 stations within and adjacent to southern New England, the event volumes, durations, and time between midpoints in table 10 vary from the minimum to maximum values by a factor of about 1.41, 1.77, and 1.48 respectively. Among the 45 precipitation stations in southern New England, the event volume is not strongly correlated to the duration or time between midpoints (with Spearman’s rank correlation coefficient values of −0.24 and 0.37). The duration and delta statistics are moderately correlated (with a Spearman’s rank correlation coefficient value of −0.714) because half the storm-event durations before and after each time between midpoint are a component of that value.

Prestorm Streamflow Statistics

SELDM uses streamflow statistics to stochastically simulate a large series of prestorm streamflow volumes from the basin upstream from the point of interest (Granato, 2013; Risley and Granato, 2014; Stonewall and others, 2019; Weaver and others, 2019). SELDM uses the frequency factor method to generate a population of nonzero streamflows from the average, standard deviation, and skew of the logarithms of nonzero flows; and conditional probability methods to simulate the fraction of zero flows (Granato, 2013). Prestorm streamflow, which may include base flow (generally defined as groundwater discharge) and stormflow from a previous storm, is one component of the total stormflow from the upstream basin. The prestorm streamflow is simulated as the instantaneous flow at the beginning of a storm, which is added to the current storm runoff for the duration of the runoff or BMP discharge from the site of interest. This is used to simulate the total flow available for dilution in the current storm. The prestorm streamflow is simulated by using the average, standard deviation, and skew of the logarithms of streamflow data. Some proportion of prestorm streamflows may equal 0 if the stream is intermittent or ephemeral. Granato (2010) provides a detailed discussion of the methods and data used for estimating prestorm streamflows for use with SELDM. Estimates of prestorm streamflow in receiving waters are important for assessing risks of adverse effects of runoff on water quality because prestorm streamflow can be a substantial proportion of total stormflow.

Streamflow statistics for three EPA Level III ecoregions that include areas within and outside of Massachusetts, Connecticut, or Rhode Island, (the Northeastern Highlands, Northeastern Coastal Zone, and Atlantic Coastal Pine Barrens ecoregions; tables 8, 11) were considered for use in streamflow simulations. The selected ecoregion statistics were the median of statistics for all USGS streamgages included in the SELDM database within each ecoregion. As table 11 indicates, there is much more variation within each ecoregion than there is between the regional medians in this area of the country.

Because these three ecoregions also include large areas outside of southern New England, statistics from three other streamflow datasets also were evaluated as alternatives for simulating prestorm streamflows (table 11). The first “SELDM” dataset includes all streamgages within Connecticut, Massachusetts, or Rhode Island that were originally included in the SELDM database (Granato, 2013); this dataset includes 106 streamgages with at least 24 complete years of record during the period 1960–2003 with drainage areas ranging from 10.6 to 497 square miles. The second, the “1901–2015” dataset, includes all 385 streamgages in Massachusetts, Connecticut, or Rhode Island with 1 or more complete years of record during water years 1901–2015 with drainage areas ranging from 0.35 to 9,660 square miles (Granato, 2017; Granato and others, 2017). The third, the “Index” dataset, includes streamgages that are within or adjacent to southern New England and are commonly selected as index streamgages for characterizing minimally altered streamflows over a common period in this area (Granato and others, 2022). This index dataset includes 73 streamgages with drainage areas ranging from 0.49 to 404 square miles. The period of record for these streamgages was water years 1961 through 2017 (57 years of record); 32 of these streamgages had complete periods of record and 41 streamgages had records that were extended to cover the full period by using the maintenance of variance type 1 (MOVE.1) technique (Granato, 2009; Granato and others, 2022).

Although the drainage-area distributions and periods of record are different, statistics for the three datasets are similar (fig. 9). Some of the similarity in flow statistics is a result of the fact that these datasets share many streamgages. The SELDM dataset shares 92 stream gages with the 1901–2015 dataset and 31 streamgages with the Index dataset; the Index dataset shares 59 streamgages with the 1901–2015 dataset. The 1901–2015 dataset has the greatest variability in statistics because it is composed of many streamgages with short-term records that may not represent long-term conditions. The Index dataset is slightly more variable than the SELDM dataset because it covers a wider geographic area and has a much wider range in drainage areas. The drainage areas for all three streamgage datasets (fig. 9) are larger than many of the drainage-basin areas above road crossings in southern New England (fig. 2). The minimum drainage area among the streamgages in the three streamflow datasets is 0.35 square miles; about 45 percent of basins above roadways have drainage areas less than this value (fig. 2). Only 8.5 percent of basins above roadways have drainage areas greater than 10 square miles (fig. 2); in comparison, 100 percent of SELDM dataset gages, 67 percent of the1901–2015 dataset gages, and 64 percent of Index gages have drainage areas greater than 10 square miles (fig. 8).

Information about relations between the average, standard deviation, and skew of nonzero streamflows is needed to guide the choice of a limited but representative set of values for simulating the potential effect of runoff on receiving waters. To this end, the nonparametric rank correlation coefficient (Spearman’s rho) was calculated among these statistics for each of the three streamflow datasets (table 12). Correlations between drainage area and the geometric mean streamflow in cubic feet per second are very strong (greater than or equal to 0.97) because drainage area is the controlling variable. However, the residual correlation between drainage area and the normalized geometric mean streamflow in cubic feet per second per square mile is weak (less than or equal to 0.34). This indicates that the selected values of the normalized geometric mean streamflow may be applied across drainage areas. Any other correlations to drainage area are only moderate to weak indicating that the geometric standard deviation and skew of the logarithms also may be applied across drainage areas. The correlations between the normalized geometric mean streamflow and geometric standard deviation are moderately strong (−0.74 to −0.82), indicating that these statistics covary. Similarly, the correlations between the geometric standard deviation and the skew of the logarithms of streamflow are moderately strong (−0.66 to −0.73), indicating that these statistics also covary. The standard deviation is calculated by using the mean and the skew is calculated by using the mean and the standard deviation, so these correlations can be used to assess the potential strength of a regression relation but should not be used for statistical inference.

Based on these correlations (table 12), regression relations were developed to select representative statistics for the geometric standard deviation and skew of the logarithms of streamflow values by using each of the three datasets (table 13). Because the potential effects of high-leverage outliers in datasets ranging over several orders-of-magnitude on regression relations can be large, the Kendall-Theil robust line method (Granato, 2006) was used to develop these equations. The objective was to produce the best estimate of the geometric standard deviation and skew of the logarithms of streamflow values given a selected value for the geometric mean of streamflows. Because the direct correlation between the geometric mean and the skew was weak, a fourth equation, which was calculated by algebraic combination of the regression equations of the mean to the standard deviation and the standard deviation to the skew, was developed and tested for each dataset. Despite the weak correlation, the direct regression equation between the geometric mean and skew was slightly more predictive than the algebraic combination. The equations in table 13 indicate that if the geometric mean flow was increased from 0.5 to 2 cubic feet per square mile over a series of simulations, then the associated geometric standard deviation would decrease from about 5.5 to about 1.8, and the associated skew would increase from about −0.7 to 0.3. At the median of geometric mean flows of about 1.03, the associated geometric standard deviation and the associated skew would be about 3.1 and −0.2, respectively.

SELDM uses the specified fraction of zero flows to simulate the effect of runoff on ephemeral or intermittent streams by using conditional-probability methods. When the prestorm streamflow value is zero, the runoff or BMP discharge volumes are likely to be a larger fraction of downstream flows than for similar runoff events with nonzero prestorm streamflows. The upstream flow also will depend on the upstream area and lagtime, which are deterministic variables in SELDM, and the upstream runoff coefficients and recession ratios, which are stochastic variables in SELDM. An ephemeral stream has no baseflow; it flows only in response to runoff. Therefore, in theory, the maximum fraction of zero prestorm streamflow for an ephemeral stream is 1. Because the EPA definition of a runoff-generating event has an interevent dry time of 6 hours, multiple runoff events could take place within one day. Given an average number of runoff-generating events per year equal to 52 (table 9), the fraction of days with zero flows is at a minimum 0.8576 for an ephemeral stream. Perennial and intermittent streams both are defined as having flows that are sustained by groundwater discharge between runoff events. Distinctions between these categories, however, are operationally defined and such definitions are not consistent from State to State in southern New England. If the Massachusetts definition of streamflows less than 0.01 cubic feet per second at the 99-percent flow duration (Massachusetts Department of Environmental Protection, 2002) is used, then the risk, expressed as a probability fraction, of zero flows for intermittent streams would range from about 0.8576 to 0.01 (the 99-percent flow duration), and the risk of zero flows for perennial streams would range from about 0.01 to 0. Based on precipitation statistics (table 10), SELDM simulations for southern New England commonly result in about 30 years of runoff-generating events; algebraically, 1 zero-flow day during a 30-year period would have a risk of zero flows equal to about 0.0001. Stochastically, in theory, a risk value less than 0.0001 would not produce any prestorm flow values equal to zero; conversely, because SELDM is a Monte Carlo model, it can produce one or more prestorm flows that are equal to zero if the proportion of prestorm flows is specified as any number greater than zero.

Although the fraction of zero flow is commonly thought to be a function of drainage area, it also depends on physiography, geography, and water use. Streamflow statistics in the datasets selected to represent conditions in southern New England indicate that zero flows occur across a wide range of drainage areas (table 14). Spearman’s rank correlation coefficients between drainage area and the fraction of zero flows are −0.32 for the SELDM dataset, −0.46 for the 1901–2015 dataset, and −0.58 for the index dataset, which indicates that drainage area may be important but it is not the only variable of interest. Comparison of the means, standard deviations, and skews of streamflow for streamgages with zero flows to streamgages without zero flows indicates that these statistics are different between these groups. The geometric means are lower, the variability (standard deviations) is higher, and the skews are more negative for streamgages with zero flows than for streamgages without zero flows (table 11). These difference in streamflow statistics may be attributed to the effects of physiography, geography, and water use.

Regression equations were developed to estimate the fraction of zero flows from the average of the logarithms of flow (table 13; eq. 5) to guide the choice of simulation values at ungaged sites. The average of the logarithms of flow (in cubic feet per second) was used as the predictor variable rather than the normalized average (in cubic feet per second per square mile) to capture the variation in drainage area and normalized mean flow of the sites with one or more zero flows in each dataset. The regression equation developed by using the SELDM dataset (table 13) should not be used for small basins (less than about 10 square miles) because there are only six streamgages in this dataset with one or more zero flows and the smallest drainage area among these streamgages is 10.6 square miles. Given the geometric mean streamflow of 0.87 cubic feet per second square mile for the SELDM dataset of southern New England streamgages (table 11), the fraction of zero flows calculated by using this equation is greater than or equal to 1 (100 percent of flows) if the drainage area is less than or equal to about 1.31 square miles. In comparison, using the geometric means of 0.796 and 0.809 cubic feet per second per square mile for the zero-flow streamgages in the 1901–2015 and Index datasets (table 11) and solving the associated regression equations in table 13 for a ratio of 1 results in drainage-area estimates of about 0.010 and 0.039 square miles for each dataset, respectively. These areas are much less then the drainage areas of 0.35 square miles for the 1901–2015 dataset streamgages and 0.49 square miles for the Index dataset streamgages with one or more zero flows (table 11). Given the geometric mean flows (table 11) and zero-fraction equations (table 13) for streamgages with one or more zero flows, the estimates of the fraction of zero flows for a 1-square-mile basin would be about 0.023 and 0.006 for the 1901–2015 and Index streamgage datasets, respectively. Similarly, the estimates for the fraction of zero flows in a 10-square-mile basin would be about 0.0034 and 0.00014 for the 1901–2015 and Index streamgage datasets, respectively. Given the large standard error of the estimates of the equations for the fraction of zero flows in table 13, these estimates of the fraction of zero flows among the streamgage datasets may be substantially different but are not significantly different at the 95-percent confidence limit.

Although SELDM will simulate nonzero prestorm streamflows below the commonly used USGS minimum reported streamflow measurement of 0.01 cubic foot per second (ft³/s; Rantz, 1982; Granato, 2010, 2013), the USGS streamflow records used to calculate the fraction of zero flows and other statistics are censored at this flow rate. SELDM uses the frequency factor method to generate a population of nonzero streamflows; that method also can be used to estimate the probability of nonzero flows below the reporting limit given the statistics for nonzero flows. If a value of 0.01 ft³/s is substituted into the frequency factor equation and it is rearranged to solve for the lognormal (if the skew is near 0) or log Pearson Type III (if the skew is substantially different from 0) frequency factor, then the resulting equation for the zero-streamflow reporting limit of 0.01 ft³/s (Rantz, 1982) is the following:

The frequency-factor value (K) calculated by using equation 4 can be converted to a probability value by using probability tables or algebraic transfer functions (Haan, 1977; Natural Resources Conservation Service, 1998; Granato, 2010, 2013). Using the average and standard deviation of the logarithms of normalized streamflow among streamgages with one or more zero flows for the 1901–2015 dataset in table 11 results in K values of −1.339, −2.824, and −4.31 for drainage areas equal to 0.1, 1.0, and 10 square miles, respectively. Using probabilities for these K values associated with a skew of −0.775 (table 11) results in zero-flow-fraction estimates of about 0.0997, 0.0111, and 0.0009 for drainage areas of 0.1, 1.0, and 10 square miles, respectively. The equivalent estimates made using the Index dataset statistics in table 11 are 0.0914, 0.0067, and 0.0003 for drainage areas of 0.1, 1.0, and 10 square miles, respectively. Alternatively, the FrequencyFactors2021 program (Granato and others, 2022) can be used to generate a table of exceedance values from the statistics of the logarithms of streamflow; these values can be used to estimate the risks for flows less than 0.01 cubic feet per second. Calculations show that using the medians of the average, standard deviation, and skew of all streamgages in the 1901–2015 dataset and Index dataset results in the fraction of zero flows of less than 0.0001 for basins greater than one square mile. The different dataset estimates are substantially different from each other and are different from the regression-based estimates. Also, substantial percentages of streamgages with drainage areas less than 30 square miles have zero-flow fractions greater than the 0.0001 (30 year) threshold (table 14). Therefore, an average of estimates calculated by using different methods and different datasets may inform professional judgement for estimating the percentage of zero flows at an ungaged site or a site with a very short record.

Runoff Coefficient Statistics

SELDM simulates runoff from precipitation by using stochastic runoff coefficients, which are the ratio of the volume of runoff in watershed inches to the volume of basin-average precipitation (in inches) during each storm event (Granato, 2013). SELDM simulates runoff coefficients from the site of interest and the upstream basin by using the Pearson type III distribution, which is defined by the input average, standard deviation, and skew of the runoff coefficients. The effects of antecedent conditions on upstream runoff coefficients are simulated by using the rank correlation to prestorm streamflow. Wetter antecedent conditions, which tend to increase runoff coefficients, commonly are associated with higher prestorm streamflows (Granato, 2010, 2013). The user-entered rank correlation coefficients between prestorm flows and upstream runoff coefficients commonly are about 0.75 for high quality datasets (Granato, 2010, 2013). The effect of antecedent conditions on runoff coefficients for the site of interest are simulated indirectly by using correlations between runoff coefficients for the upstream basin and the site of interest. These correlations, which are calculated within SELDM, are at a maximum if the imperviousness of the upstream basin and the site of interest are equal and are reduced as the impervious fractions diverge (Granato, 2010, 2013).

In SELDM, runoff coefficient statistics can be calculated as a function of the total impervious area of the site of interest and the upstream basin by using regression equations (table 15) or by entering user-defined values. The regression equations developed to calculate the average, standard deviation, and skew of the runoff coefficients for highway sites were developed with rainfall-runoff data from 58 highway basins across the country, and the regression equations for the upstream basins (or nonhighway sites of interest) were developed with data from 167 basins across the country with various nonhighway land uses (Granato, 2010, 2013). The average, standard deviation, and skew of the runoff coefficients calculated by using regression equations for completely impervious highway areas (TIA equal to 1.0) are 0.785, 0.1917, and −1.19, respectively (table 15). The average, standard deviation, and skew of the runoff coefficients of completely impervious nonhighway areas calculated by using regression equations are 0.769, 0.114, and −0.51, respectively. The average runoff coefficient for the completely impervious highway sites is higher than the average for the completely impervious nonhighway sites; this may be caused by random sampling of different sites or the difference may be caused by highway-engineering design practices to rapidly drain runoff from the roadway and efficiently convey runoff from the road to stormwater discharge locations to maintain safe driving conditions (Brown and others, 2009). Both of these average runoff coefficient values are lower than the commonly used values, which can be as high as 0.96 for completely impervious areas (Granato, 2010, 2013). Field studies show that evaporation and infiltration from paved surfaces commonly reduce the average runoff from such areas by 20 to 30 percent (Ragab and others, 2003; Mansell and Rollet, 2006; Ramier and others, 2006; Wiles and Sharp, 2008; Wanielista and others, 2010; Redfern and others, 2016; Timm and others, 2018; Salt and Kjeldsen, 2019; Rammal and Berthier, 2020). Therefore, the default average values used by SELDM (0.785 or 0.769) are more representative of measured runoff results than higher average runoff coefficient values commonly used in the literature (Granato, 2010, 2013). Although the stochastic runoff-generation algorithms in SELDM produce many events with simulated runoff coefficients at or near a value of 1 for highly impervious sites, a substantial number of events with lower simulated runoff coefficients will show the precipitation losses evident in high-quality runoff monitoring datasets (Granato, 2013).

Hydrograph Recession-Ratio Statistics

The timing of runoff from the upstream basin is defined by using the basin lagtime and the hydrograph recession ratio, which is the ratio of the duration of the falling limb to the rising limb (or time to peak) of the hydrograph (Granato, 2010, 2012, 2013). In SELDM, a triangular hydrograph is used to simulate the timing of runoff stochastically. The duration of the runoff-generating precipitation event is simulated as a random variable. The basin lagtime, which is the time between the centroid of precipitation to the centroid of runoff, is simulated as a constant value based on the main-channel length, slope, and imperviousness of the upstream basin. The hydrograph recession ratio is used to calculate the time to peak, the recession time, and therefore the duration of the highway and upstream hydrographs. The rational-method hydrograph recession ratio for highly impervious basins, which is equal to 1, is used in SELDM to simulate runoff for the highway site. The hydrograph recession ratio for the upstream basin, however, is simulated as a stochastic variable by using a triangular distribution. The triangular distribution of ratios is parametrized by using the minimum, most probable value, and maximum of ratios (Granato, 2012, 2013). The upstream hydrograph is used with the highway-runoff and BMP-discharge durations to calculate the proportion of the total upstream stormflow that is used in the mass balance and dilution calculations (Granato, 2013).

Granato (2012) calculated hydrograph recession-ratio statistics by using least-squares optimization techniques with measured runoff hydrographs from multiple storms from 41 streamgages across the United States. This dataset included 30 stream basins in Massachusetts, 1 stream basin in Connecticut, and 1 in Rhode Island. In the current study, methods developed by Granato (2012) were used to calculate hydrograph recession-ratio statistics from 20 or more runoff events for an additional 13 basins in Connecticut and 6 basins in Rhode Island (Granato and others, 2022). These analyses were done to build a 51-streamgage dataset to simulate runoff events in southern New England. These recession ratios are presented with selected basin properties in table 16. The minimums of recession ratios in the combined dataset ranged from 1.00 to 2.05 with a median of 1.10 and an average of 1.20. The most probable values of recession ratios in the combined dataset ranged from 1.00 to 2.87 with a median of 1.50 and an average of 1.67. The maximums of recession ratios in the combined dataset ranged from 2.49 to 9.67 with a median of 4.42 and an average of 4.80. In comparison, the default values for the minimum, most probable value, and maximum recession ratios within SELDM are 1.00, 1.85, and 4.40, respectively (table 17). In comparison, the median recession ratios for non-New England gages calculated by Granato (2012) were 1.0 for the minimum, 1.16 for the most probable value, and 4.05 for the maximum of the triangular distribution. Similarly, Weaver and others (2019) analyzed hydrographs from multiple storms from 30 sites in North Carolina and calculated median recession ratio values of 1.0 for the minimum, 1.07 for the most probable value, and 4.72 for the maximum of the triangular distribution. Stonewall and others (2019) analyzed hydrographs from multiple storms from 13 sites in Oregon and calculated median recession ratio values of 1.0 for the minimum, 2.22 for the most probable value, and 4.37 for the maximum of the triangular distribution. The similarities in these statistics from hydrographs in the different areas of the country indicate that, like physiography, common hydrologic processes result in similar outcomes.

Spearman’s rank correlation analyses were done in an attempt to provide guidance on the selection of hydrograph recession-ratio statistics by using basin properties. The rank correlations among recession-ratio statistics were about 0.193 between the minimum and the most probable value, about 0.189 between the minimum and maximum, and about 0.228 between the most probable value and the maximum. Correlations between the three triangular-distribution statistics and the drainage area, main-channel length, main-channel slope, basin-lag factor, and imperviousness were of mixed sign and had absolute values ranging from 0.00027 and 0.292. These low correlations indicate that basin properties cannot be used to quantitatively select recession-ratio statistics. These results are similar to the results of correlation analyses done by Granato (2012), which also included 15 other variables that included land cover, wetlands, impoundments, and other hydrologic variables. Therefore, the hydrograph recession-ratio statistics are random variables with respect to each other and to basin properties, and the median values are robust estimates for unmonitored sites.

Highway Runoff

Highway-runoff quality statistics for commonly measured properties and constituents (tables 18, 19) were calculated by using version 1.1.0b of the Highway-Runoff Database (HRDB; Granato and Cazenas, 2009; Granato, 2019a; Granato and Friesz, 2021b). All highway-runoff concentrations were simulated as random variables by using the frequency-factor method with the average, standard deviation, and skew of the transformed (logarithmic) values (Granato, 2013). Dependent relations were not used to simulate highway-runoff quality, but two bacterial constituents (p31625 and p50569) were simulated by using statistics for other equivalent bacterial constituents, and one bacterial constituent (p31673) was simulated by using statistics for urban-runoff data (tables 18, 19).

Rank correlation (rho) analysis using Spearman’s rho was used to evaluate the correlation between the average and standard deviation of the logarithms of concentrations to determine whether the values used for simulation could be selected independently. Rank correlation values for highway-runoff constituents ranged from −0.8 to 0.393 with a median correlation of −0.18 (table 19), but none of these values were statistically different from 0 at the 95-percent confidence limit (Haan, 1977). Therefore, the average and standard deviation of the logarithms of concentrations could be selected independently to simulate highway runoff quality.

Concentrations of highway runoff for all constituents were simulated by using a skew value of 0. This value was selected because the percentage of datasets with skew values that were not significantly different from 0 at the 95-percent confidence limit (Haan, 1977) was large. For highway runoff, the percentage of constituent datasets with calculated skews that were not significantly different from 0 ranged from 75 to 100 percent with an average of 93.8 percent (table 19). Skew values of 0 were used because most constituents could be characterized as lognormal and because use of a skew equal to 0 would reduce the risks for generating improbable extreme outliers in the simulated urban-runoff concentration populations (Risley and Granato, 2014).

Rank correlation analysis using Spearman’s rho also was used to examine relations between the geometric mean concentrations of constituents and the AADT reported for each highway monitoring site. These rank correlation values ranged from −0.8 to 0.93 with a median correlation of 0.445 (table 19). Only 9 of 20 highway-runoff constituents had rank correlation values that were statistically different from 0 at the 95-percent confidence limit (Haan, 1977). Regression equations between the logarithms of AADT and geometric mean concentration were developed for these 9 constituents, which included total nitrogen, total phosphorus (p00665), trace elements (except mercury), and total PAHs (table 20). Although some of the rank correlations are strong, there is considerable uncertainty in geometric-mean predictions for a given AADT value. Based on the root mean square error in logarithmic space (Helsel and Hirsch, 2002), the 95-percent prediction interval ranges by a factor of 1.4 for total nickel (p01067) and as high as a factor of about 6.86 for total phosphorus (p00665).

Although regression equations are provided, the results will not be used for simulations in this study because there are many complications for application of specific AADT estimates. AADT may not be the primary causal variable. Increased AADT is associated with increases in the imperviousness of land covers within a mile radius of highway-runoff monitoring sites (Smith and Granato, 2010; Granato and Friesz, 2021a, b). Smith and Granato (2010) determined that surrounding-area imperviousness may have a greater effect on runoff quality than AADT. Similarly, Wagner and others (2011) found only weak relations between concentrations and AADT but found that bridge deck runoff concentrations were higher in urban areas than in rural areas. Prediction uncertainties in regression results are calculated assuming that the value of the predictor variable (in this case AADT) is known (Helsel and Hirsch, 2002). Uncertainty in point AADT estimates from short-term monitoring stations commonly is on the order of plus or minus 20 percent and as high as plus or minus 50 percent for roads with less than 1,000 vehicles per day (Krile and others, 2015). Uncertainties in AADT may be much greater in basin-wide areas than at single road-stream crossings because traffic counts can change dramatically from route to route and as a road crosses each intersection. Because the application of the AADT-based estimates of the geometric mean concentrations may be highly uncertain, the median of geometric mean concentrations is the robust choice for simulating runoff quality at unmonitored sites (Stonewall and others, 2019; Weaver and others, 2019; Jeznach and Granato, 2020; Granato and Friesz, 2021a). Two other estimates, labeled the low- and high-concentration values in table 19, represent the 15th and 85th percentile of at-site geometric mean concentrations. These low- and high-geometric-mean concentration estimates may be used in rural low-traffic or urban high-traffic areas, respectively.

Urban Runoff

National urban-runoff quality statistics for commonly measured properties and constituents (table 21) were calculated by using version 1.2.0 of the Best Management Practices Statistical Estimator (BMPSE), which includes urban runoff data from the December 2019 version of the International BMP database (Granato, 2021a). Although this report is using the term “urban runoff” to describe the stormwater quality from developed areas and many of the stormwater monitoring sites in the BMPSE are fully impervious, these sites may exist outside U.S. Census Bureau (1994) defined urban areas. These urban-runoff statistics were calculated by using data from available monitoring sites, which included commercial, industrial, mixed use, parking, residential, roadway, and open-space land-cover areas. All urban-runoff concentrations were simulated as random variables by using the frequency-factor method with the average, standard deviation, and skew of the transformed (logarithmic) values (Granato, 2013). Dependent relations were not used to simulate urban-runoff quality, but two bacteria constituents (p31625 and p50569) were simulated by using statistics for other equivalent bacterial constituents (tables 18, 21).

Rank correlation analysis using Spearman’s rho was used to evaluate the cross-correlation between the average and standard deviation of the logarithms of concentrations to determine whether the values used for simulation could be selected independently. Rank correlation (rho) values for urban-runoff constituents ranged from −0.4 to 0.66 with a median correlation of −0.040 (table 21). The only seemingly strong correlation (0.66) was for Mercury (p71900), but because of the small number of sites (8) for this constituent (table 18), this correlation value is not statistically different from 0 at the 95-percent confidence limit (Haan, 1977). Therefore, the medians of the average and standard deviation of the logarithms of concentration were used to simulate urban-runoff concentrations with a skew of 0 (table 21). Because urban-runoff simulations are not the focus of this study, alternate runoff statistics were not used for these simulations.

Concentrations of urban runoff for all constituents were simulated by using a skew value of 0. This value was selected because the percentage of datasets with skew values that were not significantly different from 0 at the 95-percent confidence limit (Haan, 1977) was large. For urban runoff, the percentage of datasets with calculated skews that were not significantly different from 0 ranged from 50 to 100 percent with an average of 88 percent (table 21). Skew values of 0 were used because most constituents could be characterized as lognormal and because use of a skew equal to 0 would reduce the risks for generating improbable extreme outliers in the simulated urban-runoff concentration populations (Risley and Granato, 2014).

Risk-Based Analyses

The risk assessment process is the foundation of the regulatory framework for numeric and narrative water-quality criteria (U.S. Environmental Protection Agency, 1991, 1998, 2002; National Research Council, 2009a; National Academies of Sciences, Engineering, and Medicine, 2013). Numeric criteria are based on a concentration, frequency of occurrence, and exposure duration from which an aquatic ecosystem can recover. Based on ecological research (Niemi and others, 1990; U.S. Environmental Protection Agency, 1991, 1998), the U.S. Environmental Protection Agency (EPA) has specified three years as the acceptable risk-based frequency of occurrence for water-quality excursions. In southern New England, where the simulated long-term average number of runoff-generating events per year can range from about 55 to 59, the risk of one event in three years calculated by the Cunnane plotting position formula is about 0.559 to 0.595 percent. In many of the comparisons, however, an approximate value of 0.5 percent is used to approximate a 3-year risk. Narrative water-quality criteria commonly are statements of an objective for one or more intended uses for the waterbody; although they do not assign a numeric value, these criteria are employed with causal presumptions that trigger targeted load reductions in the watershed. In such cases, risk assessments can be used to examine the potential load reductions from different areas and the necessary margin of safety to meet various objectives.

Available stormwater data can provide information about the distribution of event mean concentrations (EMCs) to estimate the potential for exceeding a specified concentration standard or an assigned load limit, but statistics calculated from available data must be extrapolated to estimate long-term exceedance probabilities. Version 1.1.0 of the HRDB contains 106,441 concentration values with data for 414 different water-quality constituents (Granato, 2019a). However, large datasets for individual water-quality constituents are not common. For example, the HRDB has suspended solids (table 18, p00530) data for only 216 sites across the country and total phosphorus (table 18, p00665) data for only 201 sites across the country, even though these constituents are among the most commonly measured constituents. For suspended solids, the most commonly measured constituent in the HRDB, the number of EMCs per site range from 1 to 127, with a median of 16 EMCs and an average of 18 EMCs. For total phosphorus, also a commonly measured constituent, the number of EMCs per site range from 1 to 53, with a median of 15 EMCs and an average of 16 EMCs. Similarly, among instream water-quality monitoring sites with sufficient total phosphorus data to estimate transport curves (table 22), the number of samples per site ranged from 10 to 712, with a median of 24 and an average of 79. The long-term average number of events per year in southern New England ranges from about 41 to 67 events by site (table 10). Nationwide, only 10 sites in the HRDB have more than 40 suspended-sediment EMCs, and only 5 sites have more than 40 total phosphorus EMCs. Only 30 of the 69 instream sampling sites had more than 40 total-phosphorus concentration values. Using available statistics, SELDM simulates a population of long-term values based on available data.