Appraisal of Steady-State Stormwater Control Measure Pollutant Removal Models within a Dynamic Stormwater Routing Framework

Three different SCM pollutant removal models were evaluated in this study. The models were selected primarily due to their prevalent use in other SCM modeling studies and recommendations as SCM pollutant removal algorithms by the Water Environment Research Foundation (Leisenring et al. 2013).

The Fair and Geyer model [Eq. (1), Fair and Geyer 1954] is suggested by USEPA (1986) as an appropriate model for simulating particle removal in stormwater detention basins under flow-through conditionswhere $R$ = the fraction of particles removed; $v$ = particle settling velocity ($\mathrm{m}/\mathrm{s}$); $A$ = basin surface area (${\mathrm{m}}^2$); $Q$ = steady-state flowrate through basin (${\mathrm{m}}^3/\mathrm{s}$); and $n$ represents a hydraulic efficiency factor.

The model was originally developed for application to water/wastewater treatment facilities where steady-state conditions apply, and the existence of hydraulic dead zones (i.e., non-ideal settling) within settling basins reduced particle settling efficiency from that predicted using Camp’s (1946) concept of ideal basins. Although such conditions rarely exist in stormwater basins due to the intermittent and highly variable nature of rainfall/runoff, basins that fill over a short period of time and empty over an extended period of time (i.e., extended detention basins, EDB) can be reasonably assumed to be operating at steady-state over the duration of an entire event.

Here, a modified version of the Fair and Geyer model (MFG) is used that replaces the removal term ($R$) with influent ($C_i$) and effluent ($C_e$) TSS concentrations and simulates removal of particles with different settling velocities (Olson et al. 2020)where $C_i$ = influent TSS concentration ($\mathrm{mg}/\mathrm{L}$); $C_e$ = effluent TSS concentration ($\mathrm{mg}/\mathrm{L}$); and $ps{f}_k$ = the fraction of particles in particle size bin $k$.

Incorporating multiple particle bins into the model allows for a more representative analysis, considering the large distribution of particle sizes that are found in stormwater runoff (Selbig and Bannerman 2011; Roseen et al. 2011; Kim and Sansalone 2008; Greb and Bannerman 1997; USEPA 1986). We follow recent recommendations by Leisenring et al. (2013) to use five particle bins representing particle sizes in the following ranges; $2–10\mathrm{\mu m}$, $11–30\mathrm{\mu m}$, $31–60\mathrm{\mu m}$, $61–100\mathrm{\mu m}$, and $>100\mathrm{\mu m}$.

The $\mathrm{k}\text{-}{\mathrm{C}}^{*}$ model was first proposed by Kadlec and Knight (1996) to model pollutant removal in wastewater treatment wetlands. The model uses a first-order decay coefficient ($k$ in $\mathrm{m}/\mathrm{day}$) to reduce the influent concentration ($C_i$ in $\mathrm{mg}/\mathrm{L}$) towards a background or irreducible concentration ($C^{*}$ in $\mathrm{mg}/\mathrm{L}$) as a slug of pollutants moves through a treatment device. The model assumes steady-state and plug-flow conditions exist with the treatment device. The effluent concentration ($C_e$) is estimated as:where $q^{\prime }=Q/A$ denotes the hydraulic loading rate ($\mathrm{m}/\mathrm{day}$), in which $Q$ is the average inflow rate (and/or discharge rate, assuming steady-state conditions) in ${\mathrm{m}}^3/\mathrm{s}$.

Wong et al. (2006) recommended the $\mathrm{k}\text{-}{\mathrm{C}}^{*}$ model be adopted as a “unified approach” for simulating pollutant removal in SCMs, and demonstrated the model’s ability to be calibrated to measured data from several different types of SCMs for several different pollutants. It should be noted, however, that the calibration datasets in Wong et al. (2006) included spatially-variable measurements of the pollutants within the SCM. Park et al. (2011) calibrated the $\mathrm{k}\text{-}{\mathrm{C}}^{*}$ model to simulate TSS effluent discharged from EDBs using data obtained from the International Stormwater Best Management Practices (BMP) Database (www.bmpdatabase.org). The authors assumed a ${\mathrm{C}}^{*}$ value based on the smallest TSS effluent EMC in the dataset and calibrated $k$ assuming steady-state conditions. Leisenring et al. (2013) included the $\mathrm{k}\text{-}{\mathrm{C}}^{*}$ model as a potential model for simulating pollutant removal in SCMs, but noted that the estimation of $k$ and ${\mathrm{C}}^{*}$ parameters can be difficult using published SCM data (e.g., the BMP Database).

For a large number of SCMs studies, the only information available for assessing SCM performance are measured values of influent and effluent pollutant EMCs (Loáiciga et al. 2015). Recognizing the limitations of evaluating SCM performance using the percent removal metric (Strecker et al. 2001), Barrett (2005) suggested a new methodology for evaluating SCM pollutant removal performance by using linear regression of measured influent and effluent EMCs:where $m$ = the regression slope; $b$ = the regression intercept; and $\epsilon$ denotes residuals between the simulated and measured values of $C_e$. Barrett (2005) discussed several useful attributes of this model for evaluating SCM pollutant removal. One is that the intercept value provides a reasonable approximation of the “irreducible minimum effluent concentration” that is frequently observed in SCM monitoring datasets. Another is that hypothesis testing can indicate whether the slope value is significantly different than zero. If it is not, then one may conclude that effluent concentrations are independent of influent concentrations for a particular BMP. Another benefit of this methodology can be readily applied to data contained in the BMP Database for any type of BMP and pollutant (Barrett 2008).

Table 1 presents the mean parameter values and their uncertainty used in the pollutant removal models and UA. All the parameter values and estimates of their uncertainty were obtained using data from the BMP Database, either as part of this study or others. Parameter values for the linear regression model were determined as part of this study using 137 pairs of influent/effluent TSS EMCs obtained from the BMP Database. Those data and results of the linear regression analysis are presented in Fig. S1 and Table S1. Parameter values for the MFG model are adopted from the results of Olson et al. (2020), which used the same BMP Database dataset as the linear regression analysis described previously. The parameter values for the $\mathrm{k}\text{-}{\mathrm{C}}^{*}$ model are adopted from the results of Park et al. (2011), which used a smaller subset of the BMP Database dataset (approximately 46 pairs of influent/effluent TSS EMCs), but which were also included in the larger dataset used for linear regression and MFG model parameterization.

For dynamic modeling, we assume pollutants are routed through the SCM using a variable-volume, continuously stirred tank reactor (CSTR) model. This reactor model assumes that the pollutant concentration within the reactor is equal throughout (vertically and horizontally) and that the pollutant concentration discharged from the reactor is equal to the pollutant concentration within the reactor. The mass balance equation for this type of reactor can be written as (Chapra 1997):where $V$ = volume of runoff within the SCM (${\mathrm{m}}^3$); $C$ = concentration of pollutant within the SCM ($\mathrm{mg}/\mathrm{L}$) (and also discharging from the SCM due to the assumption of a CSTR); $Q_i$ = rate of runoff entering the SCM (${\mathrm{m}}^3/\mathrm{s}$); $Q_e$ = rate of runoff discharging from the SCM (${\mathrm{m}}^3/\mathrm{s}$); $\partial V/\partial t$ = rate of change of runoff stored within the SCM (${\mathrm{m}}^3/\mathrm{s}$); and $\partial C/\partial t$ = rate of change of pollutant concentration within the SCM in $\mathrm{mg}/\mathrm{L}$ per second. In typical reactor applications, Eq. (5) will also include additional terms describing the time-rate of pollutant removal within the reactor; however, such rates cannot be quantified for SCMs without intra-event pollutant concentration data. Instead, we apply the steady-state pollutant removal equations [Eqs. (2)–(4)] to runoff as it enters the reactor, such that the values of the $C_e$ terms in those equations become the $C_i$ term in Eq. (5). In order to apply Eqs. (2)–(4) to this dynamic model, the values for $A$, $Q$, and $q^{\prime }$ in those equations are computed at the timestep ($t$) that runoff enters the SCM, i.e., $A=A(t)$, $Q=Q(t)$, and $q^{\prime }=q^{\prime }(t)$.

Eq. (5) can be rewritten as a finite-difference approximation over a timestep interval $\mathrm{\Delta }t$ :

Eq. (6) can further be simplified if the dynamic model is operated over relatively small increments of $\mathrm{\Delta }t$ (i.e., on the order of minutes), where inflows and outflows are assumed to not change significantly over the timestep:

This assumption is justified under the conditions of this study because (1) the precipitation data used to generate the inflow hydrographs are aggregated on 1-h increments, so runoff rates generated from the watershed model assume that precipitation is equal within each 1-h increment; and (2) the inflow pollutographs are assumed constant throughout the entire storm event, as is typically performed for planning-level studies. More discussion on these assumptions is provided in the following sections.

Finally, Eq. (7) can be rearranged to solve for $C_{t+1}$, which is the pollutant concentration within the SCM at the end of the simulation timestep and the concentration of pollutant discharging from the SCM during the next timestep:

Dynamic runoff routing through the SCM is based on the continuity equation [Eq. (9)], where the value for $Q_{e,t}$ is estimated using the storage-indication method (Viessman and Lewis 1996) and a priori knowledge of the stage-volume-discharge relationship of the SCM, characterized as:

Inputs for the SCM models were generated using SWMM. SWMM is a widely used dynamic stormwater model capable of simulating runoff and pollutant concentrations from urban areas (Rossman 2010). Runoff quality and quantity were simulated from a hypothetical 12.1 ha (30 acre) watershed in Fort Collins, CO (Table 2) for three different precipitation events. All events had a total precipitation depth of 12.7 mm (0.5 in.), but different durations were used (Table 3). These represent actual events obtained from the hourly National Climatic Data Center (www.ncdc.noaa.gov) record of Gage 053005 recorded in Fort Collins (Table 3). A total precipitation depth of 12.7 mm was selected because the stormwater drainage criteria for Fort Collins requires that SCMs be designed to capture and treat the runoff that occurs from approximately 12.7 mm of precipitation (City of Fort Collins Colorado 2018). The three different storm durations were selected to evaluate if storm duration affects the performance of the dynamic algorithms.

Pollutant concentrations were simulated using the EMC method in SWMM. The SWMM EMC method applies a constant pollutant concentration to runoff. While it is possible to calibrate the SWMM buildup/washoff algorithms to generate intra-event variability of pollutant concentrations in runoff, it is more common to apply the EMC method for planning-level studies, due to lack of available intra-event calibration data. In this study, we applied TSS EMCs of $50\mathrm{mg}/\mathrm{L}$, $125\mathrm{mg}/\mathrm{L}$, and $200\mathrm{mg}/\mathrm{L}$.

Fig. 1 shows the hydrographs and pollutographs generated from SWMM for the three different precipitation events. These hydrographs and pollutographs were used as inputs to the dynamic SCM models.

The objectives of this study were evaluated using an EDB SCM. The simulated EDB was designed according to design criteria published by the Colorado Department of Transportation Urban Drainage and Flood Control District (UDFCD 2015), which describe methods for calculating the water quality capture volume (WQCV) and drawdown time (40 h) of EDBs implemented along the Front Range of Colorado. We assumed the EDB has vertical side slopes and a WQCV depth of 0.91 m. Table 4 presents the EDB design parameters used in the SCM simulations.

SWMM was used to generate a stage-surface area-discharge table for the EDB, assuming a single 102 mm (4 in.) circular orifice is used to control the drawdown of the WQCV. The stage-surface area-discharge curve (Fig. 2) was used as input to the dynamic SCM model.

The dynamic model generates estimates of effluent flows, pollutant concentrations, and uncertainty of the pollutant concentrations at every timestep of the simulation. To compare the dynamic results to steady-state model results, the TSS effluent EMC and variance of the EMC were computed using Eqs. (10) and (11):

Estimates of TSS effluent EMC uncertainty were generated using three different UA methods: MC, prediction interval for linear regression, and FOVE. The parameters for which uncertainty was propagated are listed in Table 1.

Parameter uncertainty for the MFG and $\mathrm{k}\text{-}{\mathrm{C}}^{*}$ models was propagated using 1,000 simulations of the steady-state and dynamic models. (Initial testing of higher numbers of simulations showed that 1,000 simulations of both models were sufficient to generate stable results.) Random parameter combinations were generated using the Latin Hypercube Sampling (LHS) method (McKay et al. 1979). The LHS method is a stratified sampling method that divides the multi-variate space into a specified number of equally probable intervals, from which one sample is generated randomly. The primary benefit of using the LHS method is that it provides similar estimates of uncertainty to traditional random sampling methods, but uses fewer model simulations. Fifth and 95th percentile values were obtained for each scenario by sorting the 1,000 model outputs lowest to highest and retrieving the value of the 50- and 950-ranked values.

The FOVE method was also used to generate uncertainty estimates for all three pollutant removal models applied to the dynamic algorithms. The FOVE method propagates the variance of random parameters using a Taylor-series expansion of the model function to generate an estimate of the variance of the model output. Consider a function $Y=g(X)$, where $X$ is a series of random variables. $X=\{X_1,X_2,\dots ,X_k\}$ and $X_o=\{{\overline{x}}_1,{\overline{x}}_2,\dots ,{\overline{x}}_k\}$ denotes the set of mean parameter values. The value of the function at $X_i=\{x_1,x_2,\dots ,x_k\}$ can be approximated using a first-order Taylor series approximation (Morgan and Henrion 1990; Tung et al. 2006):

According to properties of variance, the variance of function $Y=g(X)$ can be expressed as:

The partial derivative terms for each pollutant removal model response with respect to each model parameter were determined analytically and are provided in Table 5. The probability distribution function, expected value, and variance of each random parameter assumed in this study are provided in Table 1. Additionally, it was assumed that all random parameters were independent and uncorrelated; hence, the second term in the right-hand side of Eq. (13) was dropped from the final analysis.

Recalling that Eq. (8) is the finite-difference approximation of effluent pollutant concentrations for the dynamic algorithm, the values of both $C_t$ and $C_{i,t}$ are considered random parameters that affect the variance of $C_{t+1}$. The FOVE derivation of the dynamic pollutant routing model iswhere $a=(V_t-Q_{e,t}\mathrm{\Delta }t/V_{t+1})$; $b=(Q_{i,t}\mathrm{\Delta }t/V_{t+1})$; $\mathrm{COV}(C_t,C_{i,t})=r{\sigma }_{C_t}{\sigma }_{C_{i,t}}$; and r is the correlation coefficient between $C_t$ and $C_{i,t}$). It is reasonable to assume that $C_t$ and $C_{i,t}$ are perfectly and positively correlated (i.e., $\mathrm{r}=1$) because if $C_{i,t}$ is greater than the expected value (at all timesteps of the simulation) then $C_t$ will also be greater for all timesteps of the simulation. The 5th and 95th percentiles of $C_e$ were estimated assuming that $C_e$ was lognormally distributed (Williams et al. 2014).

The prediction interval method (Walpole et al. 1998) is used to estimate uncertainty intervals for the linear regression model. The prediction interval provides the uncertainty of the dependent variable (i.e., effluent EMC) of a single, unknown event as a function of the independent variable (i.e., influent EMC), and represents how well the best-fit linear regression parameters fit the data. For an individual storm event with influent concentration $C_i$, the uncertainty of the predicted effluent concentration $C_e$ can be estimated using Eq. (15), which is the $(1-\alpha )100\%$ prediction interval (Walpole et al. 1998):where $t_{\alpha /2,n-2}$ denotes the $t$-distribution value at the $\alpha$ significance level, with ($n-2$) degrees of freedom; $n$ = the number of data point pairs; $s$ represents the standard error of the linear regression model; $\overline{C_i}$ denotes the mean of the influent concentration data points; and $C_{i,j}$ denotes the individual influent concentration data point. The 5th and 95th percentiles of $C_e$ were obtained using $\alpha =0.1$. A flowchart depicting the steps involved in the methodology is provided in Fig. S2.

The results of deterministically modeling SCM effluent EMCs using the mean value of pollutant removal model parameters are shown in Table 6. The linear regression model produced approximately the same ($<1\%$ difference) effluent EMCs for both the steady-state and dynamic models. This was expected, because the linear regression equation is only a function of the influent concentration, and the influent concentration was constant throughout the duration of the storm events. These results demonstrate that linear regression models calibrated to EMC data can be applied to SCM models with dynamic routing without introducing additional error to the model results. It is likely that the results may be different if the influent concentration varies throughout the duration of the event; however, calibrating pollutant buildup/washoff models at the planning level is generally not possible due to lack of sufficient data.

The results of applying the MFG model showed that effluent EMCs are higher when applied to a dynamic model compared to a steady-state model. The dynamic model produced effluent EMCs that were approximately 50%, 38%, and 22% higher than the steady-state model for storm durations of 1 h, 6 h, and 24 h, respectively. The effect of influent concentration on the difference between steady-state and dynamic effluent EMCs was considerably smaller (1%–4%) than the effect of storm duration.

The dynamic $\mathrm{k}\text{-}{\mathrm{C}}^{*}$ model also produced much larger effluent EMCs compared to the steady-state model, with differences ranging from 27% to 92% higher. The largest difference (92%) was produced by applying an influent concentration of $200\mathrm{mg}/\mathrm{L}$ during a 1-h storm duration, while the smallest difference (27%) was produced by applying an influent concentration of $50\mathrm{mg}/\mathrm{L}$ during a 24-h storm duration. Overall, the differences between dynamic and steady-state model results increased with larger influent concentrations and decreased with longer storm durations.

Because the EDB was assumed to have constant surface area and influent concentration, the only difference between the steady-state and dynamic models is the value of the discharge rate ($Q$). All other parameters being equal, a higher value of $Q$ will result in a higher effluent concentrations for both the MFG and $\mathrm{k}\text{-}{\mathrm{C}}^{*}$ models. For the steady-state model, the value of $Q$ was assumed to be constant and equal to the WQCV divided by the design drawdown time (40 h), an assumption also used by Park et al. (2011). The dynamic models assumed that the discharge rate in an EDB is not constant, as the discharge rate increases with increasing depth of storage due to hydrostatic head on the outlet orifice. Fig. 3 shows the discharge rate generated from the dynamic model compared to the assumed discharge rate for the steady-state model for all three precipitation events tested. Clearly, the discharge rates simulated using the dynamic model are higher than the average discharge rate during a considerable portion of all storm events. Most importantly, the discharge rates are generally higher than the average discharge rate when the discharge hydrograph is rising, which represents conditions when runoff is entering the EDB. The dynamic models apply the pollutant removal equations instantaneously at the time that runoff enters the EDB, using the discharge rate computed for that timestep. When the instantaneous discharge rate is larger than the average discharge rate when runoff enters the EDB, the MFG and $\mathrm{k}\text{-}{\mathrm{C}}^{*}$ models generate higher effluent concentrations of pollutants. Fig. 3 also shows that the difference between the instantaneous discharge rates and the average discharge rate is much lower for the 24-h storm event than for the 1-h and 6-h storm events. This explains why the difference between the dynamic and steady-state model results are lower for the 24-h storm duration compared to the results of the 1-h and 6-h storm durations.

The 95% prediction interval (PI) TSS effluent concentrations are presented in Fig. 4 for the MFG model, Fig. 5 for the $\mathrm{k}\text{-}{\mathrm{C}}^{*}$ model, and Fig. 6 for the linear regression model. The marker on each figure represents the mean effluent TSS EMC, while lines show the PI. Each panel represents the results for a specific influent TSS concentration (e.g., $50\mathrm{mg}/\mathrm{L}$, $125\mathrm{mg}/\mathrm{L}$, or $200\mathrm{mg}/\mathrm{L}$). For the MFG and $\mathrm{k}\text{-}{\mathrm{C}}^{*}$ models, three different UA results are presented for each storm duration (e.g., 1 h, 6 h, or $24\mathrm{h}$). The first interval for each storm duration category (circle marker) represents the results of 1,000 Monte Carlo simulations of the respective steady-state models using the random parameter distributions shown in Table 1. The second interval (square marker) represents the results of 1,000 Monte Carlo simulations of the respective dynamic model. The third interval (diamond marker) represents the results of applying the FOVE method to the respective dynamic models.

For the Monte Carlo scenarios, the percentiles were obtained directly from the 1,000 simulation outputs while the percentiles for the FOVE scenario were obtained assuming a lognormal distribution of the outputs. For the linear regression model, the left interval (circle marker) was generated using the steady-state model, with uncertainty estimated using the prediction interval method and the right interval (square marker) was the result of applying the FOVE method to the dynamic model.

The widths of the uncertainty intervals generated using Monte Carlo simulations were approximately the same for both the steady-state and dynamic MFG models. The primary difference between the results of those two scenarios is that the dynamic mean EMC values and uncertainty intervals were consistently higher. As discussed in the previous section, this shift in outputs from the dynamic model is likely due to the use of the instantaneous discharge rate to compute pollutant removal in the dynamic model. The difference between the steady-state and dynamic model results was smaller for larger duration storms, but the overall PI remained mostly unchanged.

Application of the FOVE method to the dynamic MFG model produced a slightly narrower PI compared to the MC method under all combinations of influent concentration and storm duration. The 5th percentile values generated using the FOVE were all larger than those generated using the MC simulations, and the 95th percentile values generated using the FOVE were all smaller than those generated using the MC simulations. The maximum absolute difference between FOVE- and MC-generated 95th percentile values was approximately $24\mathrm{mg}/\mathrm{L}$, which occurred with influent concentration at $200\mathrm{mg}/\mathrm{L}$ and storm duration of 24 h. The maximum relative difference between those values was approximately 21%, which occurred with influent concentration at $125\mathrm{mg}/\mathrm{L}$, and storm duration of 24 h. The width of PI generated by the FOVE method tended to decrease with increasing storm duration.

The PI values for the FOVE method were narrower than the steady-state MC results for all storm durations and it appears that the narrower intervals generated using the FOVE method somewhat compensated for the dynamic model’s tendency to predict higher effluent concentrations compared to the steady-state model. However, the difference between those values increased with longer storm durations, with the FOVE method generating upper PI values approximately 20% smaller for a 24-h storm duration.

MC simulations of the dynamic $\mathrm{k}\text{-}{\mathrm{C}}^{*}$ model produced noticeably narrower uncertainty intervals compared to MC simulations of the steady-state $\mathrm{k}\text{-}{\mathrm{C}}^{*}$ model. However, most of the difference in uncertainty results between those applications is in the estimation of the lower percentiles, as the 95th percentile values are similar under all conditions. The maximum absolute difference between 95th percentile values is approximately $13\mathrm{mg}/\mathrm{L}$, which occurs with influent concentration at $200\mathrm{mg}/\mathrm{L}$ and storm duration of 1 h. The maximum relative difference between those values is approximately 7%, which occurs with influent concentration at $125\mathrm{mg}/\mathrm{L}$, and storm duration of 1 h. This is an important finding when considered within the context of UA, as decision makers may be interested in planning for worst-case scenarios to reduce the probability that discharged pollutants exceed a regulatory threshold. The 95th percentile value, for example, has a 5% chance of being exceeded under the simulated conditions. Thus, one may conclude that applying the dynamic $\mathrm{k}\text{-}{\mathrm{C}}^{*}$ model with MC methods produces estimates of higher percentile values that are reasonably close to the steady-state model with MC methods.

The FOVE method generally produces narrower uncertainty intervals compared to the MC method, when applied to the dynamic model. One exception is for the case of influent concentration of $50\mathrm{mg}/\mathrm{L}$ and storm duration of 24 h, where the FOVE uncertainty estimates are slightly wider. Similar to the results generated using the MFG model, the FOVE method produced lower values of 95th percentiles and higher values of 5th percentiles. The maximum absolute difference between 95th percentile values is approximately $33\mathrm{mg}/\mathrm{L}$ and the maximum relative difference between those values is approximately 20%, which both occur with influent concentration of $200\mathrm{mg}/\mathrm{L}$ and storm duration of 24 h. In general, the FOVE uncertainty intervals tend to increase in width with increasing storm duration.

Overall, the FOVE method produces narrower uncertainty intervals than the prediction interval for all combinations of influent concentration and storm durations; however, the primary effect of the smaller uncertainty intervals is the estimate of the 5th percentile value. The 95th percentile estimates are very similar for both methods, with the maximum absolute difference between 95th percentile values being approximately $4\mathrm{mg}/\mathrm{L}$ (for influent concentration of $200\mathrm{mg}/\mathrm{L}$, regardless of storm duration) and the maximum relative difference being approximately 5% (for storm durations of 1 h, regardless of influent concentration). As discussed in the previous section, the estimates of higher percentile values may be more important to decision makers as they consider how to plan and implement pollutant reduction strategies with low probability of failure; thus the application of the FOVE method within a dynamic modeling framework can likely be used without significantly affecting decision making under uncertainty.