Comparison of SELDM Simulated Total-Phosphorus Concentrations with Ecological Impervious-Area Criteria

The objective of this study was to simulate TP event mean concentrations (EMCs) from developed impervious areas and undeveloped areas in concurrent downstream stormflows over a range of impervious percentages to examine relations between the percentage of IC and the risk for exceeding commonly used water-quality criteria. These relations may provide planning-level estimates of phosphorus EMCs as an explanatory variable for ecosystem changes commonly observed as IC increases in developing drainage basins. Although many water-quality constituents in runoff from impervious areas have adverse effects on receiving waters (Walsh et al. 2005; Brown et al. 2009; National Research Council 2009), TP was selected because data are available for runoff, receiving waters, and BMP performance.

Stormwater runoff in a hypothetical stream basin that represents hydrologic and physiographic basin properties in southern New England was simulated using the Stochastic Empirical Loading and Dilution Model (SELDM) to conduct a numerical experiment designed to explore relations between IC and receiving-water quality. These simulations included a range of IC from 0.1% to 30%. A numerical experiment was needed because long-term characterization data for runoff quality, BMP treatment, and background stormwater quality in hydrologically similar basins with varying impervious fractions were not available. The potential effectiveness of structural stormwater control BMPs for decreasing instream concentrations by hydrograph extension, flow reduction, and water-quality treatment were simulated with statistics that are representative of actual BMP performance (Granato 2014). Simulations were done using input values that are representative of southeastern New England because this area is highly developed and is within the area where IC-based TMDLs are being used.

SELDM simulates stormflows from runoff-generating events using statistics for prestorm streamflows, precipitation, and runoff coefficients (Granato 2013). Individual prestorm streamflows, precipitation event characteristics, and runoff coefficient values are simulated by using the log-Pearson Type III distribution, the two-parameter exponential distribution, and the Pearson Type III distribution, respectively. Furthermore, the effects of antecedent conditions on runoff coefficients are simulated by incorporating the rank correlation to prestorm streamflow. The SELDM user can select regional statistics, statistics from sites that represent a site of interest or user-defined statistics. Statistics available within SELDM for USEPA Level III Ecoregion 59, the Northeastern Coastal Zone, were selected to develop planning-level estimates of prestorm streamflows and precipitation in this study because values for this ecoregion are representative of hydrologic conditions in southeastern New England. The average of streamflow statistics from 79 stream gauges and 31 precipitation stations in the ecoregion were used in the simulations. The proportion of zero prestorm streamflows was 0.00057. The retransformed average, standard deviation, and skew of the logarithms of normalized streamflow (in liters per second per square kilometer) were 11.36, 3.0004, and $-0.1358$, respectively. The average precipitation volume, duration, and time between storm-event midpoints were 17.8 mm, 9.69 h, and 155 h, respectively. The average, standard deviation, and skew of runoff coefficients for the pervious area were 0.129, 0.099, and 1.08, respectively. The average, standard deviation, and skew of runoff coefficients for the impervious area were 0.769, 0.114, and $-0.51$, respectively.

SELDM uses five basin-property variables to simulate the volume and timing of runoff from a site of interest and the upstream basin (Granato 2013). These variables include the drainage area, main channel length, main channel slope, and imperviousness (Granato 2012, 2013). The drainage area, the percentage IC, and the precipitation volumes determine the runoff volumes from each drainage area. A total drainage area of approximately $26{\mathrm{km}}^2$ ($10{\mathrm{mi}}^2$) was used for all the simulations because this is a moderately sized basin for Ecoregion 59. The impervious area was increased and the pervious area was decreased proportionally to simulate impervious percentages of 0.1% (to simulate essentially natural conditions) 1%, 2.5%, 5%, 10%, 15%, 20%, and 30% to simulate various levels of development. The timing of runoff from the impervious and pervious areas is determined by the basin lag time, which is simulated using the main channel length, main channel slope, and imperviousness (Granato 2013). Increasing main channel length will increase basin lag time; increasing main channel slope or imperviousness will decrease basin lag time (Granato 2012).

Although the pervious and impervious areas were varied for the different simulations, the sum of areas remained at $26{\mathrm{km}}^2$ and the “upstream” length and slope were held constant, as would be the case for a developing basin. The main channel length and slope of the pervious area were estimated using available physiographic information for stream gauges in the northeastern United States (Granato 2012). The main channel length had a strong correlation to drainage area, with a Spearman’s rho value of 0.958; a length value of approximately 10.17 km was calculated as the median channel length for $26\text{-}{\mathrm{km}}^2$ basins in this area. The main channel slope had a moderately strong correlation to drainage area, with a Spearman’s rho value of $-0.75$; a slope value of approximately $6.17\mathrm{m}/\mathrm{km}$ was calculated as the median channel slope for $26\text{-}{\mathrm{km}}^2$ basins in this area. These values resulted in a basin lag time of approximately 4 h for the upstream pervious area. Runoff from the pervious area was simulated using a triangular hydrograph with standard stochastic hydrograph-recession ratio statistics (Granato 2012, 2013).

Although the entire impervious area is simulated as a single discharge at the point of interest, the average main-channel length and slope were simulated to represent the timing of runoff from many small impervious areas distributed along the tributaries of the stream above the site of interest. The length of overland flow for contributing areas commonly has been estimated as half the length between stream channels (Carlston 1963). Using this assumption, the flow length has been estimated as half of the reciprocal of drainage density. Bent et al. (2014) calculated drainage density values for 197 streamflow measurement sites in southern New England. In this data set, the drainage areas ranged from 0.728 to $761{\mathrm{km}}^2$ and the drainage density ranged from 0.347 to $2.35\mathrm{km}/{\mathrm{km}}^2$. Analysis of these data indicate that, with a rank correlation coefficient (Spearman’s rho) of only 0.129, drainage density is not a function of drainage area at these sites in southern New England. The average drainage density for these sites is $1.26\mathrm{km}/{\mathrm{km}}^2$, and the reciprocal is approximately 0.8 km. Therefore, the Horton half-distance would be approximately 402 m (Horton 1945). Assuming that the storm sewers are designed to maintain a flow velocity greater than approximately $0.76\mathrm{m}/\mathrm{s}$ to avoid deposition of sand and sediment, an average drainage-system slope of approximately $1.9\mathrm{m}/\mathrm{km}$ is assumed (Peavy et al. 1985). Using these basin properties for the IC areas results in a basin lag time of approximately 6.5 min (Granato 2012). Although the actual flow distance from different impervious areas may range from several meters to a kilometer, the differences in the estimated basin lag times for these distances vary by orders of a few minutes for completely impervious areas (Granato 2012); these differences are small in comparison to the average storm duration (approximately 10 h) in the Northeastern Coastal Zone ecoregion. Runoff from the impervious area was simulated using a hydrograph with an equal rising and falling limb.

The quality of runoff from the developed areas was simulated using EMC data for unfiltered TP from the September 2016 version of the International Stormwater BMP Database. EMC data collected from BMP inflow-monitoring points represent stormwater from developed areas that is routed to the BMPs for treatment. Queries were written to identify BMP-inflow monitoring sites, obtain stormwater-quality data, and calculate sample statistics for each monitoring site to represent urban-runoff quality for nonhighway land uses. A simple replacement method was used to estimate statistics for sites with one or more values below the reported detection limits (censored data). The International Stormwater BMP Database does not provide the means to calculate censored values using statistical methods but instead has a concentration field that has either an uncensored value or half the reported detection limit. Antweiler and Taylor (2008) indicated that substituting concentrations equal to half the detection limits was sufficient for developing planning-level estimates for water-quality statistics. Croghan and Egeghy (2003) demonstrated that substituting concentrations equal to the detection-limit concentration divided by the square root of 2 produced unbiased estimates up to censoring levels of approximately 50%. Therefore, the Croghan and Egeghy (2003) approach was used and only data sets with censoring levels that were less than or equal to 50% of concentrations were included for estimating sample statistics. Approximately 18% of the TP data sets obtained from the International Stormwater BMP Database have one or more censored values.

SELDM uses the average, standard deviation, and skew of the logarithms of EMCs to simulate a population of runoff concentrations from the developed area. There are wide ranges in the site statistics for TP in the International Stormwater BMP Database (Table 1, Fig 1). For example, the maximum geometric mean is 107 times the minimum geometric mean for TP. EMC statistics were grouped by land-use categories, but there are no clear relations between land-use category and the magnitude of the geometric means (Fig. 1). The statistics listed in Table 1 were ranked and selected independently; the averages, standard deviations, and skew values on each row in the table may be from different study sites. To evaluate the validity of this approach, a rank-correlation analysis was done to assess potential relations among the average, standard deviation, and skew of the logarithms of EMC values. The rank correlation coefficient between the average and standard deviation and the average and skew of the logarithms of TP EMCs were approximately $-0.17$ and $-0.18$, respectively. The rank correlation coefficient between the standard deviation and skew of the logarithms of TP concentrations was approximately 0.1. Analysis of the site statistics indicated that the average, standard deviation, and skew of these EMC data are not correlated. Therefore, the developed-area runoff quality was simulated using the median of the average ($-0.7988$) and standard deviation (0.2953) of the logarithms of concentrations (Table 1).

Only 24% of logarithmic skew values were outside the 95% confidence limits of a zero skew for TP and the statistically significant skew values had both positive and negative signs, so the developed-area runoff quality was simulated as a lognormal variable with a skew of zero. Use of zero skew also may be warranted because the logarithmic skew values were only weakly correlated with the average and standard deviation of the logarithms of EMCs. This approach also is supported by the literature because the lognormal distribution commonly is used to characterize and simulate urban-runoff quality (Di Toro 1984; Novotny 2004; National Research Council 2009).

The upstream water quality was simulated using a water-quality transport curve with stochastic variations above and below the regression line. A water-quality transport curve is a regression relation between streamflow volumes and constituent concentrations (O’Connor 1976; Glysson 1987; Granato 2006, 2013). A two-segment water-quality transport curve (Fig. 2) was developed using the logarithms of paired streamflow and concentration data with the Kendall Theil Robust Line program (Granato 2006). The transport curve was developed with 159 paired streamflow and TP concentration values measured at USGS Station 01170100 Green River near Colrain, Massachusetts, during the period from March 25, 1993, through August 2, 2016. These data are available in the USGS National Water Information System database (USGS 2019). Data from this stream gauge were selected because the basin is rural (approximately 91% forest and 0.28% impervious based on the USGS StreamStats delineation), and paired values of streamflow and concentration were available over a large range of flows. Among these paired measurements, streamflow ranged from 2.657 to $\mathrm{1,431}\mathrm{L}/\mathrm{s}/{\mathrm{km}}^2$ and TP concentrations ranged from 0.002 to $0.125\mathrm{mg}/\mathrm{L}$. The equation for the first segment isand is used when streamflow is less than $73.851\mathrm{L}/\text{s}/{\mathrm{km}}^2$ and the equation for the second segment isand is used otherwise. In these equations MAD is the median absolute deviation of the regression residuals (1.33 for the first segment and 1.61 for the second segment) and ${\mathit{K}}_{\mathit{random}}$ is the random normal variate, which is used to simulate the scatter of concentrations above and below the regression line. MAD is used to simulate variability rather than the standard deviation of residuals because the transport curve is being used to simulate EMCs rather than instantaneous concentrations (Granato 2013). The transport curve is being used as a planning-level estimate of the central tendency of EMCs for a given streamflow value. The zero slope of the first segment indicates that concentration varies randomly with a lognormal distribution with a geometric mean of $0.004\mathrm{mg}/\mathrm{L}$ when streamflow is below $73.851\mathrm{L}/\mathrm{s}/{\mathrm{km}}^2$. Concentration varies stochastically with flow and lognormal random variation above this threshold. Although the maximum flow measured concurrently with TP sampling is approximately $\mathrm{1,431}\mathrm{L}/\mathrm{s}/{\mathrm{km}}^2$, the regression equation produces TP concentrations that are well within observed TP concentrations in minimally developed basins in southern New England when much higher flow rates are simulated.

SELDM can simulate three stormwater-treatment mechanisms, including runoff hydrograph extension, volume reduction, and water-quality treatment, using a trapezoidal distribution with correlation to the inflow values (Granato 2013, 2014; Granato and Jones 2014). Three types of generic BMPs were simulated: one with hydrograph extension only, one with flow reduction only, and one comprehensive BMP with hydrograph extension, flow reduction, and water-quality treatment. BMPs with hydrograph extension lengthen the runoff-discharge hydrograph, increasing the amount of concurrent upstream stormflow by retention, thereby increasing the dilution of runoff into a greater volume of receiving-water stormflow. Hydrograph extension may reduce the instantaneous erosivity of stormflow but may increase the fully mixed mass-balance concentration at the mixing point if upstream concentrations exceed the runoff concentrations. BMPs with volume reduction decrease stormflow by infiltration or evapotranspiration. Comprehensive BMPs with water-quality treatment use biological and physicochemical processes to decrease outflow concentrations. The SELDM BMP module also simulates the irreducible minimum concentration by substituting the specified value when simulated BMP effluent concentrations are below the specified value. In this study, a concentration of $0.007\mathrm{mg}/\mathrm{L}$ TP was used as the irreducible minimum concentration (Granato 2014). In SELDM, the flow and water-quality treatment statistics used can result in a net increase in flow or concentrations for some storms, which is consistent with actual BMP monitoring results. In this study, the median of treatment values for the bioretention, composite, detention basin, biofilter (swale), infiltration basin, manufactured device, media filter, retention pond, wetland basin, and wetland channel BMPs from Granato (2014) were used (Table 2).

Table 2 also shows the theoretical median and average for the trapezoidal distribution of each treatment statistic because the average and median are more familiar than the trapezoidal values (Granato 2014). These measures of central tendency are not applied to each event; the actual simulated treatment statistics for each event depend on the correlation to inflow values and, for concentrations, the effect of the irreducible minimum concentration value.

The potential effects of untreated runoff were examined by comparing the risks of exceeding example water-quality criteria of 0.02375 and $0.1\mathrm{mg}/\mathrm{L}$ TP (Fig. 3). Four slices of the simulated TP populations were used to examine the potential effects of increasing IC on TP EMCs. The first slice (horizontal dashed line in Fig. 3) represents the USEPA reference condition (which reflects pristine or minimally impacted waters) for TP of $0.02375\mathrm{mg}/\mathrm{L}$ based on the 25th percentile of nutrient data for Ecoregion 59 (USEPA 2000). This criterion is designed to represent conditions of surface waters that are minimally impacted by human activities. The second slice, also shown as a horizontal dashed line in Fig. 3, represents a fixed concentration of $0.1\mathrm{mg}/\mathrm{L}$ (0.1 ppm) TP, which is a TP criterion for preventing nuisance plant growth in streams (USEPA 2017). This fixed concentration was used to examine how the risk of exceeding this threshold changed with variations in the IC percentage. The third and fourth slices of the population (vertical dashed lines on Fig. 3) represent fixed exceedance frequency values of 0.58% and 50%. The exceedance frequency of 0.58% is the once-in-3-years exceedance frequency, which represents the USEPA return period selected as a protective measure to provide for ecological recovery from periods of severe stress (USEPA 2017). The 3-year percentage may be slightly different in this ecoregion depending on the number of storms and annual-load accounting years that are simulated using different random seeds (Granato 2013).

In Fig. 3, the percentage of simulated events that exceed a TP EMC of $0.02375\mathrm{mg}/\mathrm{L}$ ranges from approximately 3.33% of events with an IC of 0.1% to approximately 99.12% of events with an IC of 30%. The percentage of events that exceed a TP EMC of $0.1\mathrm{mg}/\mathrm{L}$ ranges from approximately 0.16% of events with an IC of 0.1% to approximately 59.74% of events with an IC of 30%. The TP EMCs associated with the once-in-3-years exceedance risk range from 0.05 to $0.68\mathrm{mg}/\mathrm{L}$ for the range of ICs from 0.1% to 30%, respectively. The TP EMCs associated with the median exceedance risk range from 0.01 to $0.12\mathrm{mg}/\mathrm{L}$ for ICs ranging from 0.1% to 30%, respectively.

The potential effects of treated runoff were examined by comparing the risks of exceeding the 0.02375 and $0.1\mathrm{mg}/\mathrm{L}$ criteria (Fig. 4). In these simulations, all runoff from the IC area was simulated as being treated, which may overrepresent the results that may be achieved using structural BMPs in such basins. In Fig. 4, four different scenarios with and without BMPs are shown to illustrate the effects on downstream EMCs in a drainage area with 10% IC. The BMPs that were simulated included a comprehensive BMP (with flow reduction, hydrograph extension, and water-quality treatment) or BMPs with only flow reduction or hydrograph extension alone. The distribution of “natural” or upstream EMCs shown in the figure illustrates the performance of the BMPs relative to upstream water quality. Based on simulation results, BMPs have the greatest effect on decreasing downstream EMCs (compared with no BMP) when downstream EMCs are greater in higher-probability events. For example, with a comprehensive BMP, simulated TP EMCs below the 50th percentile are on average 6.2 times greater than the simulated natural background EMCs compared to EMCs above the 50th percentile that are on average 4.9 times the natural background concentration. Hydrograph extension and flow reduction alone decreased downstream EMCs compared to the simulations without BMP treatment, but not as effectively as with a comprehensive BMP. The BMPs with only hydrograph extension were slightly more effective than BMPs with only flow reduction, decreasing EMCs to an average of approximately 10.4 times the natural upstream EMCs. Flow reduction decreased EMCs to an average of approximately 11.2 times the upstream EMCs.

The simulations also indicate the potential for load reduction. Hydrograph extension BMPs alone do not effectively reduce downstream loads of TP because hydrograph extension alters the timing of flow but not the total mass load of urban-runoff constituents to the stream. BMPs with flow reduction decrease the total amount of flow to the stream and therefore decrease the long-term average annual mass load of TP to the receiving waters from the impervious area by 31%. With a comprehensive BMP (one that includes hydrograph extension, flow reduction, and water-quality treatment), long-term average annual runoff loads from the impervious area were decreased by 55% compared to no BMP.

Fig. 5 shows the simulated effects of BMPs on TP downstream EMCs at the 0.58% (once-in-3-years event) and 50% exceedance frequencies with increasing IC. The simulated BMP results shown in this figure include the effects of the hydrograph extension, flow reduction, and water-quality treatment statistics described in Table 2. The downstream concentration at the 0.58% exceedance frequency is 46% lower in the BMP simulations than the untreated simulations. BMPs have a similar effect on downstream EMCs at the 50% exceedance frequency, decreasing downstream EMCs on average by 49%. Generally, the greater the IC, the greater the positive effect of the BMP on decreasing downstream EMCs.