Adapting the Relaxed Tanks-in-Series Model for Stormwater Wetland Water Quality Performance

: Across the globe, water quality standards have been implemented to protect receiving waters from stormwater pollution, motivating regulators (and consequently designers) to develop tools to predict the performance of stormwater control measures such as constructed stormwater wetlands (CSWs). The goal of this study was to determine how well the relaxed tanks-in-series ( P-k-C* ) model described the performance of CSWs in North Carolina. Storm events monitored at 10 CSWs in North Carolina were used for calibrating the model, and statistical evaluations concluded the model could adequately predict the performance for all pollutants except organic nitrogen. Nash–Sutcliff calibration/validation values were determined to be 0.72/0.78, 0.78/0.74, 0.91/0.87, 0.72/0.62, 0.88/0.73, and 0.91/0.63 for total nitrogen, total ammoniacal nitrogen, oxidized nitrogen, organic nitrogen, total phosphorus, and total suspended solids, respectively. Sensitivity analysis revealed only one calibration parameter with strong sensitivity, the Arrhenius coefﬁcient (temperature dependent model coefﬁcient). With this model, CSWs can be optimized to treat watershed-speciﬁc inﬂuent concentrations to meet efﬂuent targets. In general, the current design technique used in North Carolina and many other locations (a ﬁrst ﬂush volume detention method) oversizes CSWs for water quality vis- à -vis the method herein, suggesting improved designs for water quality may be possible through scientiﬁcally-informed methods.


Introduction
Constructed stormwater wetlands (CSWs) have become one of the more commonly utilized stormwater control measures (SCMs) [1][2][3][4]. A well-functioning CSW has an ecosystem with diverse plant and animal species that also removes water-borne pollutants (source). Most regulatory entities set design standards to capture and detain the water quality volume of runoff, which is released over two to five days and subjected to treatment during the inter-event period [5,6]. If these design requirements are met, pollutant removal credits for total nitrogen (TN), total phosphorus (TP), and total suspended solids (TSS) of 40%, 40%, and 85%, respectively, are granted in North Carolina (NC), with similar removal credit systems being in place in other municipalities in the United States [7]. This sizing methodology implies a degree of treatment that may not be accurate, necessitating more mechanistic approaches to accurately size and credit CSWs when water quality reduction is the primary goal.
Research performed on CSWs suggests that substantial water quality improvement may be possible in smaller systems. Hathaway and Hunt (2010) [2] studied the water quality treatment of three CSW cells in series near Charlotte, NC, USA. The authors observed at least 80% of the total concentration reduction occurred within the first wetland cell for all pollutants examined, 1.08 0. 12 25 Notes: 1 Urban Runoff = U, Agricultural Runoff = A, Golf Course = GC; 2 Wetland 1 Only per nomenclature defined by authors; 3 [7]; 4 [6].
A strong body of knowledge is available for analogous installations of free surface water constructed wetlands for wastewater treatment, hereafter referred to as "treatment wetlands". Such installations have proven to be successful since their introduction as a water treatment technology in the 1970s [3,13]. For 40 years, water quality data have been collected from treatment wetlands, with the finding that the rate of concentration reduction for pollutants is well represented by first order (the k-C* model) decay kinetics [13][14][15]. Subsequently, these first order decay functions have been adapted to event-driven, CSW, systems [14,16,17]. Although Kadlec and Knight (1996) [16] reasoned the plug flow assumption associated with the k-C* model might not be representative of the actual hydrodynamics in CSWs, using the model did provide a conservative estimate of treatment performance provided background concentrations were considered. The k-C* model is often used as an interpolator on data sets, but difficulties arise when the model extrapolates outside of the calibrated concentration ranges, or for comparing design configurations [18]. Thus, although the k-C* model is simple and widely used, it fails to adequately characterize the complex processes that occur in treatment wetlands or CSWs [19].
More recently, the P-k-C* reaction rate model (relaxed tanks-in-series, PTIS, reaction rate model) was adapted to stormwater applications by Wong et al. (2006) [20]. This model has been endorsed because it accounts for velocity (i.e., short-circuiting and stagnant zones) and reaction rate heterogeneity (i.e., concentration weathering) in wetlands, while using a minimal number of parameters [3,21]. However, there has not been an effort to investigate the applicability of this model on CSWs utilizing field-based data, a critical need before this more scientifically informed method of CSW design for water quality can be widely utilized. Thus, the goal of this study was to determine the ability of the P-k-C* model to describe the treatment response of CSWs in North Carolina and demonstrate its application as a design optimization tool.

Model Description
During unsteady flows (descriptive of stormwater systems), water quality varies as "parcels" of variably treated water passes through a CSW [3]. To more thoroughly describe CSW treatment, it is necessary to simultaneously consider both pollutant depletion and flow hydrodynamics [3]. The Tanks in Series (TIS) model (Equation (1)) applies a combination of the k-C* model and the continuously stirred tank reactor (CSTR) model to adequately describe both pollutant depletion and flow hydrodynamics [3,20,22]. This model assumes steady state conditions: no infiltration, no evapotranspiration, and constant flow while this parcel of water is in the wetland. Despite this assumption not truly reflecting field conditions, the TIS and k-C* models have been widely and successfully used to predict discharge concentrations in treatment wetlands [21,23,24].
where C out = effluent concentration (mg/L), C in = influent concentration (mg/L), C* = background concentration in the wetland (mg/L), k T = reaction rate constant (m/year) at ambient Temperature, T ( • C), τ = detention time (d), N = tracer test determined number of tanks in series (TIS), and h = wetland free water depth (m). Per previous studies [1,3,19,21,22,[25][26][27], N, k T , and C* model inputs can be relaxed to model fitting parameters: P, k T,app (apparent reaction rate constant), and C*. This forms the relaxed TIS model (PTIS or P-k-C* model). Because the P-k-C* method encompasses the final concentration profiles from a variety of causative distributions, it better represents multiple situations [3,21]. Nitrogen and phosphorus biogeochemical removal mechanisms and their rate constants are temperature sensitive [3,21,28,29]. The modified Arrhenius equation (Equation (2)) describes temperature dependence and is often combined with reaction models for model calibrations [3,21,28]. A θ value = 1.000 indicates temperature does not affect treatment, while values below or above 1.000 signify a negative or positive effect on treatment, respectively, at higher temperatures [3].
where k T = removal rate constant (m/year) at water temperature, T ( • C), k 20 = removal rate constant at 20 • C (m/year), and θ = dimensionless temperature coefficient. The result of combining Equations (1) and (2) utilizing model filling parameters is presented in Equation (3). This model combines the first order decay function with a reaction model and coefficient adjustment to account for varying hydrodynamics (decay function) and temperature (coefficient adjustment) (Equation (3)).

Site Description
For model calibration and validation, influent and effluent concentrations, design depths, and detention times were compiled from ten CSWs studied in North Carolina ( Figure 1). Physical  Table 2.

Monitored Influent (Cin) and Effluent (Cout) Concentrations
The methods used to collect and analyze water samples were similar among studies: flow-paced water quality samples were collected from the inlets and outlets of all CSWs using ISCO automatic samplers such that nutrient and sediment event mean concentrations (EMC) could be determined (Table 3). Analytes from each site included total ammoniacal nitrogen (TAN), nitrite-nitrate nitrogen (NO2,3-N), TP, and TSS. Organic nitrogen (ON) concentrations were calculated by subtracting TAN from total Kjeldahl nitrogen (TKN) measurements for individual storm events. Total nitrogen concentrations were also calculated by adding TKN and NO2,3-N values for each storm event.

Monitored Influent (C in ) and Effluent (C out ) Concentrations
The methods used to collect and analyze water samples were similar among studies: flow-paced water quality samples were collected from the inlets and outlets of all CSWs using ISCO automatic samplers such that nutrient and sediment event mean concentrations (EMC) could be determined (Table 3). Analytes from each site included total ammoniacal nitrogen (TAN), nitrite-nitrate nitrogen (NO 2,3 -N), TP, and TSS. Organic nitrogen (ON) concentrations were calculated by subtracting TAN from total Kjeldahl nitrogen (TKN) measurements for individual storm events. Total nitrogen concentrations were also calculated by adding TKN and NO 2,3 -N values for each storm event.

Background Concentrations (C*)
Background concentrations (C*) were developed by combining data from CSW literature and previous studies that applied first order decay models on treatment wetlands. Treatment wetland modeling studies typically assumed that the background concentrations for TP, NO 2,3 -N, and TAN were zero [3,23]. Conversely, using effluents of CSWs and a reference natural wetland in North Carolina, Moore et al. (2011) [33] determined the background concentration for ON to range from 0.7 to 0.8 mg/L. The minimum observed effluent TSS concentration (2 mg/L) was used for the TSS background concentration. Seasonal variations (temperature) were not considered for the background concentration estimates in this analysis.

CSW Design Characteristics
The detention time, τ, and free water depth, h, are assumed to remain relatively constant for treatment wetlands, but for CSWs, these parameters depend on antecedent conditions and vary from event to event [20,27]. The average free water depth and detention time for the water quality storm event (25 or 38 mm) reported for each CSW (Table 4) were used for this analysis. Since water temperature was not measured for any of the CSWs, data were compiled from U.S. Geological Survey (USGS) stream stations located closest to each CSW [34]. Although stream water temperatures are most likely to be less than those observed in CSWs, this was the best available estimate of this parameter. The historical mean monthly water temperatures ( • C) of each station were calculated; CSW storm events were assigned a water temperature, on a monthly basis, from the nearest USGS station. Temperatures ranged from 3.3 • C in January at the Asheville, NC (Mountains) station to 31.3 • C in July at the Falls, NC (Piedmont) station.

Model Calibration, Validation, and Sensitivity Analysis
Data from the CSWs were parsed into individual storm events and divided into calibration and validation datasets by assigning odd numbered events (per chronological order) to the calibration dataset and even numbered events to the validation dataset. This method provided a seasonal representation for the calibration and validation datasets for each CSW. For situations when there were not an adequate number of events for both model calibration and validation (i.e., n < 8), all events were used for calibration.
Calibration parameters; k 20 , P, and θ, were found by fitting the modeled data to the observed and minimizing the root mean squared error (RMSE in mg/L) between observed effluent concentrations (O i ) and modeled effluent concentrations (P i ) by simultaneous adjustment of k 20 , P, and θ for all events, n [19,28]. In addition to the RMSE, the coefficient of determination, R 2 , and Nash-Sutcliffe efficiency, NSE [35] were used as statistical performance measures for the model. The NSE was selected as a dimensionless goodness-of-fit indicator.
Sensitivity of fitting parameters (P, k 20 , and θ) was explored using a Monte-Carlo approach, with parameter ranges being developed from literature and the calibration outcomes in this study. For each pollutant, 250,000 simulations were conducted using the Nash-Sutcliffe coefficient as the objective function. Parameter distribution plots were generated for each CSW individually and all CSW data collectively based on the acceptable parameter sets identified in the analysis, allowing an understanding of the sensitivity of the model to each parameter.

Individual CSW Calibration and Validation: Resultant Reaction Rate Parameters
Modeling fitting (reaction rate) parameters and statistical evaluations were generated for both calibration and validation datasets for each pollutant and CSW ( Figure 2). The reaction rates for each evaluation are summarized in Table 5.
The apparent number of TIS (P) ranged from 1.4 to 7.6 for each pollutant and CSW (Table 6). This range fell well within the scope of tanks in series values (1-10) derived from conservative tracer experiments compiled by Kadlec and Wallace (2008) [3] and simulations conducted by Persson et al. (1999) [36] for treatment wetlands. The mean and median P values ranged from 3.0 to 4.0 (Table 5), similar to the value (3.0) recommended for initial modeling estimates and for use in treatment wetland design [3,36,37]. Overall, the CSW sites with lower P values (EB, DB, and JSC1), indicative of well-mixed hydrodynamics, had longer detention times relative to the other sites (Table 4). Higher P values suggested more short-circuiting and, typically, corresponded with deeper storage depths (CMS and Bass, Table 4), where interactions with soil and vegetation are minimized and limit flow retardation and treatment potential [4,27,36,38].  The temperature coefficients (θ) for pollutants derived from the field investigations varied for each CSW. The mean values for TN, TAN, TP, and TSS demonstrated temperature effects to be negligible with values of 1.007, 1.002, 1.002, and 0.993, respectively (Table 5). In relatively warm climates, such as North Carolina, no temperature effect on TP performance has typically been observed [3,28,39], and variability in TSS treatment is not usually attributed to temperature [3]. The lack of effect of temperature on TN and TAN treatment was surprising since it has been welldocumented that nitrogen processes are influenced seasonally [21,28,29,39]. However, when the TN and TAN loading to a wetland is substantially less than the growth requirements of the plants and algae (<120 and 108 g/m 2 /year, [3]), the removal of TN and TAN is mediated by the growth and decay of biomass [29,40]. Plant uptake rates peak in the spring and do not correspond to the annual cycle of water temperatures [3,29]. Therefore, in lightly loaded CSWs, water temperature and modified Arrhenius θ values have no effect on TN and TAN removal. The mean θ for NO2,3-N and ON (1.024 and 1.034, respectively) reflected greater treatment at higher temperatures (Table 6). Kadlec (2010) [37] observed θ values ranging from 1.035 to 1.113 for NO2,3-N treatment in Ohio similar to that (1.043) observed in New Hanover County, NC [3]. Stanford et al. (1973) [41] found θ values for ON ranged from 1.07 to 1.08, while Marion and Black (1987) [42] observed θ = 1.08-1.16. The NO2,3-N and ON temperature effects herein were not as prominent as these studies. This was attributed to the smaller range of annual water temperatures observed at the USGS stream station.
The rate constants (k20) derived from the calibration exercises resulted in a large range for each pollutant ( Table 5). The upper range primarily consisted of rate constants found for the NCSU, JSCA, UNCA sites. These sites had higher influent concentrations compared to the sites with lower k20 values (RB and BES). Because these rate constants describe the rate a wetland can treat influent concentrations to selected background concentrations, higher values are expected for higher influent concentration. This trend does not hold when comparing rate constants for CSWs to those of treatment wetlands. The mean rate constants found in this study were higher than those found for treatment wetlands (Table 6). For example, nitrogen values equated to the 91st (phosphorus, the 84th) percentile in the distribution of treatment wetland rate constants compiled by Kadlec and Wallace (2008) [3]. Although treatment wetlands have higher influent concentrations compared to CSWs, treatment wetlands tend to have longer detention times and steady flow conditions compared to the flashy, episodic nature of stormwater systems [27,36]. Constructed stormwater wetlands yield a higher "rate of treatment" given the flashiness of their hydrology, demonstrating that rate constants are also a function of detention time and hydrodynamics (rather than solely on influent concentrations). This supports Kadlec (2010) [37] where NO2,3-N dynamics in pumped riverine wetlands for pulsed vs. steady flow conditions were evaluated. Rate constants observed were much higher for pulsed flow than for steady state conditions (k20 = 107 vs. 37 m/year). The same was observed for TP in Kadlec (2001) [43] with k20 = 212 vs. 88 m/year for pulsed flow and steady state conditions, respectively. Few studies have been conducted to find rate constants for CSWs. Carleton  (Table 5). In relatively warm climates, such as North Carolina, no temperature effect on TP performance has typically been observed [3,28,39], and variability in TSS treatment is not usually attributed to temperature [3]. The lack of effect of temperature on TN and TAN treatment was surprising since it has been well-documented that nitrogen processes are influenced seasonally [21,28,29,39]. However, when the TN and TAN loading to a wetland is substantially less than the growth requirements of the plants and algae (<120 and 108 g/m 2 /year, [3]), the removal of TN and TAN is mediated by the growth and decay of biomass [29,40]. Plant uptake rates peak in the spring and do not correspond to the annual cycle of water temperatures [3,29]. Therefore, in lightly loaded CSWs, water temperature and modified Arrhenius θ values have no effect on TN and TAN removal. The mean θ for NO 2,3 -N and ON (1.024 and 1.034, respectively) reflected greater treatment at higher temperatures (Table 6). Kadlec (2010) [37] observed θ values ranging from 1.035 to 1.113 for NO 2,3 -N treatment in Ohio similar to that (1.043) observed in New Hanover County, NC [3]. Stanford et al. (1973) [41] found θ values for ON ranged from 1.07 to 1.08, while Marion and Black (1987) [42] observed θ = 1.08-1.16. The NO 2,3 -N and ON temperature effects herein were not as prominent as these studies. This was attributed to the smaller range of annual water temperatures observed at the USGS stream station.
The rate constants (k 20 ) derived from the calibration exercises resulted in a large range for each pollutant ( Table 5). The upper range primarily consisted of rate constants found for the NCSU, JSCA, UNCA sites. These sites had higher influent concentrations compared to the sites with lower k 20 values (RB and BES). Because these rate constants describe the rate a wetland can treat influent concentrations to selected background concentrations, higher values are expected for higher influent concentration. This trend does not hold when comparing rate constants for CSWs to those of treatment wetlands. The mean rate constants found in this study were higher than those found for treatment wetlands (Table 6). For example, nitrogen values equated to the 91st (phosphorus, the 84th) percentile in the distribution of treatment wetland rate constants compiled by Kadlec and Wallace (2008) [3]. Although treatment wetlands have higher influent concentrations compared to CSWs, treatment wetlands tend to have longer detention times and steady flow conditions compared to the flashy, episodic nature of stormwater systems [27,36]. Constructed stormwater wetlands yield a higher "rate of treatment" given the flashiness of their hydrology, demonstrating that rate constants are also a function of detention time and hydrodynamics (rather than solely on influent concentrations). This supports Kadlec (2010) [37] where NO 2,3 -N dynamics in pumped riverine wetlands for pulsed vs. steady flow conditions were evaluated. Rate constants observed were much higher for pulsed flow than for steady state conditions (k 20 = 107 vs. 37 m/year). The same was observed for TP in Kadlec (2001) [43] with k 20 = 212 vs. 88 m/year for pulsed flow and steady state conditions, respectively. Few studies have been conducted to find rate constants for CSWs. Carleton et al. (2001) [17] summarized rate constants for TP, TAN, and NO 2,3 -N in 18 studied constructed wetlands receiving urban runoff. Rate constants for TP ranged from 4.9 to 46.7 m/year aligning with that observed herein (38.0 m/year), but the ranges for NO 2,3 -N and TAN (3.6-57.1 and 1.0-24.6 m/year, respectively) were less than that found here (65.5 and 37.5, respectively: Table 6). For wetlands with varying ponded surface areas and volumes, such as CSWs, Carleton et al. (2001) [17] based their calculations on maximum values of these parameters to generate maximally conservative (that is, minimal) estimates of rate constants. This methodology could explain why their reported k 20 values were, overall, less than those derived in this study.

Calibration and Validation Results
Statistical evaluations were compiled for calibration and validation datasets for each pollutant and CSW (Tables 7 and 8, respectively). Wong et al. (2006) [20] reported the ability of a unified model (P-k-C*) to predict stormwater treatment of various SCMs, and although no rate constants were provided, TN, TP, and TSS model fits (RMSE = 0.08, 0.03, and 6.9 mg/L, respectively) supported this conclusion. Similar findings were observed in this study with RMSE averages of 0.21, 0.02, and 6.85 mg/L for TN, TP, and TSS, respectively ( Table 7). The RMSE for TN was higher than that found by Wong et al. (2006) [20], but given the range of influent concentrations and CSW designs, the RMSE was still relatively low. Total nitrogen and TP validation RMSE values were similarly low: 0.20 and 0.04 mg/L, respectively; however, the RMSE for TSS increased to 21.19 mg/L. This increase was attributed to the poor fit of one CSW, NCSU (RMSE = 124.4 mg/L). The median was similar to results found during calibration (6.59 mg/L).
The average R 2 statistics indicated that the model explained 65-77% of the total variance in observed effluent concentrations for all pollutants (Table 8). This range widened for the validation datasets where 24-77% of observed concentration variability was explained by the model (Table 8). Nash-Sutcliffe model efficiency values across all pollutants ranged from 0.46 to 0.73 and −4.32 to 0.31 for calibration (Table 7) and validation (Table 8), respectively, where a negative NSE signified the observed mean effluent concentration was a better predictor of effluent concentrations than the model. Both of these statistical parameters are sensitive to extreme events [44], commonly found in the relatively small sample size of events (i.e., n < 8) for each individual CSW's calibration and validation datasets.
Plots demonstrating observed vs. predicted relationships for pollutant effluent concentrations from model calibration and validation events across all CSWs are compiled in Figures 3a-f and 4a-f, respectively. Essentially, the CSWs were modeled, calibrated, and validated individually before the observed/predicted values were recompiled for an overall calculation of RMSE and NSE. The statistical metrics of all compiled CSW results for each pollutant are summarized in Table 9.    Plots demonstrating observed vs. predicted relationships for pollutant effluent concentrations from model calibration and validation events across all CSWs are compiled in Figure 3a-f and Figure  4a-f, respectively. Essentially, the CSWs were modeled, calibrated, and validated individually before the observed/predicted values were recompiled for an overall calculation of RMSE and NSE. The statistical metrics of all compiled CSW results for each pollutant are summarized in Table 9.   Model prediction strength varied between calibration and validation. While events were randomly assigned to the calibration and validation data sets, those data sets were in some cases surprisingly dissimilar. For example, the TSS calibration ( Figure 3f) and validation (Figure 4f) dataset distributions were different, even though these datasets were chosen at random. Regardless, the P-k-C* model was typically able to describe the observed behavior well. These results: (1) suggest there is an underlying unity of a number of complex processes; and (2) support the model's intended use for conceptual analysis and as a design tool. As an early attempt to calibrate the P-k-C* model for use in many varying CSWs, these initial results are promising.
Model prediction strength varied between calibration and validation. While events were randomly assigned to the calibration and validation data sets, those data sets were in some cases surprisingly dissimilar. For example, the TSS calibration ( Figure 3f) and validation (Figure 4f) dataset distributions were different, even though these datasets were chosen at random. Regardless, the P-k-C* model was typically able to describe the observed behavior well. These results: (1) suggest there is an underlying unity of a number of complex processes; and (2) support the model's intended use for conceptual analysis and as a design tool. As an early attempt to calibrate the P-k-C* model for use in many varying CSWs, these initial results are promising.

Sensitivity Analysis Results
The results of the sensitivity analysis are presented in Figure 5. The distributions found for the Arrhenius coefficient, θ, indicate that the model outputs were sensitive to the change in this parameter, where a θ value of approximately 1.0 was the optimized value, although some wetlands' optimized values were closer to approximately 1.5. The mean calibrated θ values determined for CSWs in NC (Table 5) fall within the optimized range determined by the sensitivity analysis ( Figure  5) and values determined by previous studies. Kadlec and Wallace (2008) [3] cataloged θ values of

Sensitivity Analysis Results
The results of the sensitivity analysis are presented in Figure 5. The distributions found for the Arrhenius coefficient, θ, indicate that the model outputs were sensitive to the change in this parameter, where a θ value of approximately 1.0 was the optimized value, although some wetlands' optimized values were closer to approximately 1.5. The mean calibrated θ values determined for CSWs in NC (Table 5) fall within the optimized range determined by the sensitivity analysis ( Figure 5) and values determined by previous studies. Kadlec and Wallace (2008) [3] cataloged θ values of 1.002, 1.078, 1.082, and 1.043 for TP, ON, TN, and NO 2,3 -N, respectively, for a wetland system studied in New Hanover County, NC. For TAN, a coefficient of 1.040 was used by Stone et al. (2002) [28] for a livestock treatment wetland located in Duplin County, NC. Finally, total suspended solids concentration time series often display some degree of sinusoidal behavior, similar to nitrogen and phosphorus, through the course of a calendar year; however, variability in performance has not been attributed to temperature [3]. Therefore, a temperature coefficient for TSS of 1.000 is appropriate.  The lack of a defined peak in the distribution for both P and k 20 suggests that model outputs were not sensitive to these parameters. Kadlec and Wallace (2008) [3] compiled tracer response curves for several surface flow treatment wetlands and found them all to be reasonably represented by P = 3, in spite of their different geometries. Likewise, a wetland with a meandering low flow channel with full vegetation was best described as P = 3.1 TIS in Melbourne, Australia [36]; as was NO 2,3 -N treatment in event-driven wetlands [37]. Wong and Geiger (1997) [14] suggested rate constants for TN and TSS of 22 and 1000 m/year, respectively, based on 59 treatment wetland systems. Carleton et al. (2001) [17] narrowed this range in a study of removal rate constants (k 20 ) for 19 constructed wetland systems receiving urban runoff. Lack of sensitivity in k 20 was a surprising result. The consistent concentrations and detention times, which were relatively less than typical values for treatment wetlands, could explain why k 20 was not found to be sensitive in the model, which was a surprising result. The mean, calibrated values for P and k 20 determined in this study (Table 5) were similar those observed in these previous studies, supporting the wide use of these mean values due to their lack of sensitivity on model outputs. That is, these parameters can possibly be set by literature values and no longer require calibration.

Discussion
The practice of stormwater wetland design for water quality improvement is in need of more scientifically-informed methods. Fortunately, there is a strong base of literature to rely on from analogous treatment wetland studies which have primarily targeted wastewater. However, most literature and design methods continue to use the plug flow based k-C* model and univariate analysis, where a single number is used to describe a parameter such as a detention time or rate constant. However, treatment and stormwater wetland variables and parameters do not possess single unique values, but rather a distribution of these values with respect to a wetland attribute [1,3,21]. Wetland detention time distributions alone are often responsible for departures from plug flow behavior [18,19,21], thus affecting the transferability of these models outside their treatment wetland calibration envelops or even to the stochastic nature of CSWs. Kadlec (2003) [21] demonstrated the combination (P-k-C* model) of first-order decay kinetics with a CSTR batch model could adequately embody wetland parameter distributions. This was imperative for the transferability of the model among CSWs with variable design characteristics, detention times, and influent concentrations.
For most pollutants investigated, the P-k-C* kinetic model was able to describe the observed behavior quite well, with the exception of ON, which was not necessarily unexpected given internal residuals and wetland return fluxes [29,33]. These results suggested there is an underlying unity of a number of complex processes and supported the model's intended use for conceptual analysis and as a design tool. Wong et al. (2006) [14] drew the same conclusions when utilizing the model across various SCMs with differing pollutant removal mechanisms.
After model calibration to achieve CSW-centric reaction rate parameters (namely k 20 and P), rearrangement of the P-k-C* model lends designers the ability to solve for the detention time or area required to achieve a target effluent concentration. Designer inputs include watershed-specific influent concentrations for pollutants of concern, effluent targets, design storm event size, and the CSW's design free water depth.
Considering this is one of the first attempts to calibrate the P-k-C* model to numerous individual field-scale runoff events across many CSWs, the results were deemed successful, but also provide a focus for future research to enhance the model. The first priority for such research would be to incorporate storm events from CSWs monitored long term. The age range (1-5 years) for calibrated CSWs in this study was relatively small and primarily consisted of "young" systems. It is imperative to test the applicability of the model for not only different catchments and design characteristics but also different ages. Stochastic effects on treatment, especially in CSWs, could mask the seasonal behavior of a wetland system [39]; therefore, analyzing several years of data to develop seasonal trends in performance is important in accurately establishing temperature coefficients. Further research into the hydrodynamics of CSWs by conducting conservative tracer tests would provide insights into detention time distributions. These tests could also establish field-based short-circuiting/mixing parameters (N) instead of finding these values empirically (P).
The complexity of modeling wetland processes has created two distinctly different reactions among practitioners [21]. Wetland scientists argue for greater quantification and utilization of data and details [45,46], but offer no mechanisms for incorporating that detail or data into practical wetland design [21]. Practitioners and regulators condemn the "black box" approach, but reject detailed modeling methods in favor of practicality [21,47]. The tradeoff between model accuracy and calibration data requirements is difficult. The most reliable prediction of pollutant removal for given event in a specific CSW will be provided by a more detailed model, but often at the expense of extensive calibration requirements. Use of the P-k-C* model minimizes the parameters to be calibrated (k 20 , θ, and P), as demonstrated in this study for CSWs in North Carolina, while also offering the ability to embody the weathering behavior of concentrations and deviation from plug flow observed in real wetland systems. Further, the analysis herein suggests that two of these calibration parameters may be reduced to literature values based on their relative lack of sensitivity. This is important if designers are to use this as a design tool across a range of catchments. Greenway (2004) [4] argued the area of a CSW should be a function of catchment size, the volume of runoff, pollutant characteristics of the runoff, the desired level of water treatment, and the extent to which the CSW was also expected to function as a retention basin. Currently, the design method in North Carolina only considers catchment size, volume of runoff, and peak flow mitigation, but the application of the P-k-C* model can assist practitioners in including runoff pollutant characteristics and the desired level of water treatment. Many studies recommend performance metrics other than percent removal should be employed for SCMs [48][49][50]. Strecker et al. (2001) [48] stated that biological and downstream habitat assessment should be explored as an assessment technique, and McNett et al. (2010) [51] reasoned that evaluating SCMs in context with the environment in which they are located may be as valid a way of measuring how well they are functioning as the currently used concentration and load based metrics. Design of CSWs could also follow this same reasoning by optimizing designs to treat watershed-specific influent concentrations to target effluent thresholds based on receiving waterbodies.
Steady state conditions, no gains or losses, is one main underlying assumption of the P-k-C* model. Volume reductions in CSWs are typically negligible [8,12,52]; however, some systems have experienced significant reductions (24-55%) due to seepage and evapotranspiration [8,50,53,54]. Consequently, CSW designs based on this model would provide conservative area recommendations, and further refinements to the model could include water balance calculations to evaluate CSW performance and sizing on a mass loading basis.

Conclusions
The goal of this study was to determine the ability of the P-k-C* model to describe the treatment response of CSWs in North Carolina, and demonstrate its application as an optimization tool for designing CSWs. Such methods will aid in advancing the science of stormwater management beyond a "black box" understanding of internal SCM processes, allowing more scientifically-informed design. Storm events monitored at 10 CSWs in North Carolina were used to determine reaction rate parameters through model calibration exercises. Statistical evaluations concluded that the P-k-C* model could accurately describe the water quality treatment response for all pollutants except organic nitrogen, due to complex internal wetland fluxes of this parameter, with NSE values ranging from 0.63 for TSS to 0.87 for NO 2,3 -N under the validation exercise. Sensitivity analysis showed that the only sensitive calibration parameter was θ, while k 20 and P were mostly insensitive and could potentially be set to literature values in subsequent studies.
Practitioners have the ability to optimize designs with this model to treat watershed-specific influent concentrations and target effluent thresholds based on receiving waterbodies. Overall, the current design technique (CN Method) is oversizing CSWs with respect to meeting target water quality goals, and reduced CSW areas derived utilizing this new methodology could lead to lower cost of property acquisition and construction. Incorporation of long-term event monitoring datasets are needed to expand the applicability of the model, and further refinements could include water balance calculations to evaluate CSW performance and sizing on a mass loading basis. Further, studies in other geographic locations would corroborate the results of the sensitivity analysis. Overall, the model does provide an efficient approach to describing water quality response, and allows a better understanding of what future research is needed to refine and quantify factors influencing CSW performance.