Comparison of Four Nitrate Removal Kinetic Models in Two Distinct Wetland Restoration Mesocosm Systems

The objective of the study was to determine the kinetic model that best fit observed nitrate removal rates at the mesocosm scale in order to determine ideal loading rates for two future wetland restorations slated to receive pulse flow agricultural drainage water. Four nitrate removal models were investigated: zero order, first order decay, efficiency loss, and Monod. Wetland mesocosms were constructed using the primary soil type (in triplicate) at each of the future wetland restoration sites. Eighteen mesocosm experiments were conducted over two years across seasons. Simulated drainage water was loaded into wetlands as batches, with target nitrate-N levels typically observed in agricultural drainage water (between 2.5 and 10 mg L−1). Nitrate-N removal observed during the experiments provided the basis for calibration and validation of the models. When the predictive strength of each of the four models was assessed, results indicated that the efficiency loss and first order decay models provided the strongest agreement between predicted and measured NO3-N removal rates, and the fit between the two models were comparable. Since the predictive power of these two models were similar, the less complicated first order decay model appeared to be the best choice in predicting appropriate loading rates for the future full-scale wetland restorations.


Introduction
Over the next decade the United States is expected to continue and intensify its focus on ecological systems and their capacity to store water and improve water quality [1,2]. The impacts of increased pollutant loads from non-point sources continue to threaten important economic and recreational aquatic ecosystems. These growing pollutant loads have the potential to further impair drinking water, as well as continue to threaten important aquatic dependent industries and recreational water-based environments [3].
Many of the environmental pressures in coastal regions have been linked directly to the drainage systems that were implemented to alter areas with wetland characteristics (e.g., high water tables, hydric soils, and hydrophytic vegetation) into highly productive farmland [4][5][6]. Each year an estimated 13 million tons of nitrogen (N) are purchased to incorporate into cropland soils in the United States enabling crop production to sustain the continually growing world population [7]. Fertilizer applications have been reported to be the dominant N input in most sensitive watersheds found in the Southeastern and Mississippi River regions [8]. Furthermore, reactive nitrogen in drainage water entering sensitive freshwater and marine ecosystems has been found to harmfully alter these systems (e.g., acidification, toxic algae blooms, hypoxia, contamination of drinking water aquifers, loss of biodiversity) [3,[9][10][11][12].
Negative ecological impacts from draining wetland environments could be reversed in part by restoring wetlands. Mitigation or conservation programs that invest resources into strategically positioned wetland restoration projects in sensitive watersheds could help lead to a reversal in both declining wetland ecosystem structure and water quality trends. Recent studies have been conducted to improve restored wetland designs to offer additional water quality ecosystem services that may have not been originally a part of the ecosystem. For example, strategically placed wetlands in coastal landscapes receiving pumped agricultural drainage water have been reported to not only enhance wetland hydrology, but provide reduction of nutrient, sediment, and fecal bacteria up to 97% from incoming drainage water [13,14]. However, N removal rates, particularly nitrate-N (NO 3 -N) which is significant in agricultural drainage water, can be variable in wetland systems depending on such factors as wetland to watershed ratio, soil type, wetland ecosystem type, residence time, and loading rates [15][16][17]. To promote the increased application of wetlands to improve watershed health, the development of improved wetland removal process design tools for sensitive watersheds are needed.
A significant body of work has been completed to determine the N removal capacity of agricultural, municipal wastewater, and stormwater treatment wetlands [18][19][20][21]. However, N removal in constructed and restored surface pulse flow wetlands continues to be based primarily on simple zero and first-order decay models with few studies completed to investigate alternative N removal kinetic models [22][23][24]. Some researchers have identified several potential drawbacks to using zero and first order decay approaches, particularly regarding utilizing models originally developed for constant flow municipal applications in intermittent flow agricultural settings [25,26].
Restoration success evaluations are rarely completed because project funds often are predominately earmarked for construction of on-the-ground projects [27]. Therefore, the utilization of wetland mesocosms may be a more cost effective alternative to enhance our understanding of N processes prior to full-scale construction and provide useful predictions of the NO 3 -N assimilation and removal capacity of these systems [25,[28][29][30].
In this study, the potential to reduce NO 3 -N concentrations at two future wetland restoration sites in eastern North Carolina was evaluated. Resultant freshwater pulses from developed and drained agricultural land that contain excess nutrients, sediment, and fecal bacteria, have likely contributed to the lower quality of shellfishing waters in these areas [31]. Restored forested pulse flow wetlands, which are planned for these sites, are one type of wetland that has been found to provide water treatment, flood abatement, biodiversity, and carbon storage in coastal regions [32]. A series of laboratory mesocosm experiments were completed to determine the best kinetic model to predict NO 3 -N removal rates in the two distinct wetland systems. This was vitally important to maximize the hydrologic and nutrient assimilation and removal potential for these future restored wetlands because the amount of drainage water pumped to the wetland area will assist in restoring a more natural wetland hydroperiod and equal the reduction of drainage water currently pumped directly to a nearby sensitive estuary. The study was unique because eighteen NO 3 -N removal studies were completed to capture differences of each wetland's NO 3 -N removal capacity with respect to N loading and season. Findings provided from this study were intended to guide wetland designers as to the best method to predict NO 3 -N removal capacity in wetland systems that will receive an intermittent load of water with relatively low concentrations of NO 3 -N.
The primary objectives of this study were to: (1) Create a dataset of NO 3 -N removal observations over varying seasons, N loadings, and antecedent conditions for two pulse flow wetland environments at the mesocosm scale; (2) Evaluate the fit of four NO 3 -N removal models (zero order, first order decay, efficiency loss, and Monod) with observed daily NO 3 -N removal observations; and (3) Provide predictions for the maximum drainage water volumes that can be pumped into the future restored wetlands based on the kinetic model that best fit and was most practical for observed NO 3 -N removal datasets.

Experimental Setup
Six large wetland mesocosms (3.5 m long × 0.9 m wide × 0.75 m deep) were constructed in a greenhouse near North Carolina State University 16 months prior to the initial experiment ( Figure 1). Details of the mesocosm construction and experiments can be found elsewhere [33]. The soil/water interface is where critical biogeochemical transformations occur in wetlands; therefore, using restoration site soils was a critical component for the success of estimating site specific NO 3 -N removal rates for these wetlands. In summary, soils were excavated directly from the future restoration sites to be used in the study. Three randomized mesocosm replicates were loaded with Scuppernong soil (a poorly drained, organic soil typically associated with Pocosins), while three more were loaded with Deloss soil (a poorly drained, mineral soil typically associated with marine terraces). The Scuppernong and Deloss wetland mesocosms will be referred to as WET-Org and WET-Min, respectively, for the remainder of this article. Soft-stemmed bulrush (Schoenoplectus tabernaemontani) was chosen because other researchers have successfully utilized this plant to study wetland nutrient dynamics at the mesocosm scale [24,34]. Additionally, the above and belowground biomass of S. tabernaemontani can establish quickly with careful attention to the plant requirements. Three smaller mesocosms served as controls for the experiment and had only the simulated drainage water. Each wetland mesocosm was outfitted with a recirculation system to replicate the conditions of drainage water moving slowly through a wetland. Steady state flow was maintained through each wetland mesocosm by collecting water along one end of the mesocosm and pumping the water through a linear PVC apparatus to the other end of the mesocosm. Pumping rates were adjusted such that average turnover time of water in each mesocosm was approximately once a day. Complete details of the mesocosm setup can be found elsewhere [34].

Experimental Setup
Six large wetland mesocosms (3.5 m long × 0.9 m wide × 0.75 m deep) were constructed in a greenhouse near North Carolina State University 16 months prior to the initial experiment ( Figure 1). Details of the mesocosm construction and experiments can be found elsewhere [33]. The soil/water interface is where critical biogeochemical transformations occur in wetlands; therefore, using restoration site soils was a critical component for the success of estimating site specific NO3-N removal rates for these wetlands. In summary, soils were excavated directly from the future restoration sites to be used in the study. Three randomized mesocosm replicates were loaded with Scuppernong soil (a poorly drained, organic soil typically associated with Pocosins), while three more were loaded with Deloss soil (a poorly drained, mineral soil typically associated with marine terraces). The Scuppernong and Deloss wetland mesocosms will be referred to as WET-Org and WET-Min, respectively, for the remainder of this article. Soft-stemmed bulrush (Schoenoplectus tabernaemontani) was chosen because other researchers have successfully utilized this plant to study wetland nutrient dynamics at the mesocosm scale [24,34]. Additionally, the above and belowground biomass of S. tabernaemontani can establish quickly with careful attention to the plant requirements. Three smaller mesocosms served as controls for the experiment and had only the simulated drainage water. Each wetland mesocosm was outfitted with a recirculation system to replicate the conditions of drainage water moving slowly through a wetland. Steady state flow was maintained through each wetland mesocosm by collecting water along one end of the mesocosm and pumping the water through a linear PVC apparatus to the other end of the mesocosm. Pumping rates were adjusted such that average turnover time of water in each mesocosm was approximately once a day. Complete details of the mesocosm setup can be found elsewhere [34].

Wetland Mesocosm Experiments
Eighteen batch experiments were conducted between September 2012 and September 2014 (16 months after mesocosm construction and planting) to capture differences with respect to N loading and season of each wetland's NO3-N removal capacity. The experiments were conducted as batch

Wetland Mesocosm Experiments
Eighteen batch experiments were conducted between September 2012 and September 2014 (16 months after mesocosm construction and planting) to capture differences with respect to N loading and season of each wetland's NO 3 -N removal capacity. The experiments were conducted as batch studies, with each mesocosm loaded with identical hydraulic and nutrient loading rates during each experiment (Table 1). Initial NO 3 -N concentrations were varied during the study at ranges commonly found in agricultural drainage water at each restoration site (2.5-10 mg L −1 ). The base target load was 0.6 g NO 3 -N m −2 . The load represented an 18 cm water depth and 2.5 mg L −1 NO 3 -N concentration, which is the most frequent loading expected at the future wetland restoration sites. Higher target loads that varied between 0.9 g NO 3 -N m −2 and 3.6 g NO 3 -N m −2 and water depths of 18 to 30 cm were also evaluated during the batch studies. Stage gages were permanently placed in each mesocosm to observe average water depth within the mesocosm at the beginning of each experiment and used for calculations discussed below. Each batch experiment lasted 7 to 10 days. More detailed methods for the initiation of each batch run can be found elsewhere [34].

Sampling Plan and Analysis
Grab water quality samples were collected from the recirculation system on days 0, 1, 2, 3, 5, and 7 during all batch runs and in 14 of 18 batch runs on day 10. Samples were then analyzed for NO 3 -N after being filtered through 0.45 µm filters. All water quality analyses of the grab samples were conducted in the Soil Science Environmental and Agricultural Testing Service Lab (SSC-EATS) with a Quikchem 8000 (Lachat, Milwaukee, WI, USA) using the cadmium reduction method (EPA Method 353.2). Evapotranspiration rates were calculated for NO 3 -N concentration adjustments based on the changes in the water depth measured on a stage gage in each mesocosm over the course of each experiment. Additionally, water temperature in each wetland mesocosm was monitored hourly using 8k HOBO pendant temperature sensors (Onset Computer Corporation, Bourne, MA, USA).

Apparent NO 3 -N Removal Calculations
Area-based NO 3 -N removal rates (J NN , mg m −2 d −1 ) were determined on a daily basis during each experiment.
where, C i was the initial NO 3 -N load applied to each mesocosm (mg) at the beginning of each timestep, C t was the NO 3 -N remaining in the mesocosm water at each timestep (mg), A was the surface area of the wetland mesocosm (m 2 ), and t was the timestep or mean residence time of water held within the wetland mesocosm system. Determination of the removal rates for a particular experiment were terminated when C t < 0.05 mg L −1 [35].

NO 3 -N Removal Kinetic Models
Four models were considered to predict NO 3 -N removal: zero order (ZO), first order decay (FO), efficiency loss (EL), and Monod (M). The zero order, first order, and Monod models were chosen because they have been used to model NO 3 -N removal in wetlands in past experiments [25,36]. The efficiency loss model was chosen because to our knowledge it had not been explored prior to this study in wetlands, but has been observed to model NO 3 -N removal well in streams [37][38][39]. Observations from the first nine of the eighteen completed batch run experiments were used to determine removal rate coefficients, which varied in NO 3 -N load and season. Coefficients were calculated using NO 3 -N concentrations 24 h following NO 3 -N loading to allow initial concentrations to stabilize in each mesocosm and final NO 3 -N concentrations for each model. The kinetic model fits were then evaluated using the remaining nine experimental datasets from the batch run experiments [40].

Zero Order (ZO) Model
The most simplistic quantitative model, the zero order model (ZO), assumes contaminant reduction is independent of NO 3 -N concentration. Zero order NO 3 -N models have been used to model NO 3 -N removal in wetlands assuming a constant consumption rate of NO 3 -N [41,42]. Therefore, the J ZO , the area based NO 3 -N removal rate (g m −2 d −1 ), would be constant. These systems can be characterized by a linear NO 3 -N concentration profile, which is often unrealistic in systems with longer residence times that allow NO 3 -N to become limited. The model assumes the system is closed, anoxic, completely or partially mixed, independent of hydraulic loading rates, and insignificantly influenced by other kinetic reactions occurring within the system [43]. The model can be expressed as [44]: where, C 1 was the NO 3 -N concentration in the wetland mesocosm 24 h after simulated drainage was added to allow the system to be well mixed (mg L −1 ), t F−1 was the mean residence time or period the water was held within the wetland mesocosm system after the water was well mixed (24 h) until C t < 0.05 mg L −1 (d) [35], C F was the NO 3 -N concentration at end of the analysis (mg L −1 ), and D was initial average water depth (m). Once the J ZO was estimated using observed NO 3 -N removal datasets from the first 9 of the 18 mesocosm experiments, C t was predicted using the following equation for the remaining 9 batch experiments: where, C Applied is the NO 3 -N concentration in simulated drainage water applied to the wetland mesocosm on Day 0 of the experiments.

First Order Decay (FO) Model
The second model, the first order decay model (FO), assumes NO 3 -N reduction rates are directly proportional to NO 3 -N concentration [45]. Therefore, removal rates (J FO ) increase linearly with NO 3 -N concentration assuming removal efficiency does to not change in relation to NO 3 -N load [46]. The model also assumes that the substrate concentration is significantly smaller than the half-saturation constant (K s , the concentration supporting an uptake rate of 50% of the maximum rate), the system is well mixed, has no significant influences from water losses or gains, and is dependent on only one reactant [14,43,47,48]. The ρ FO , the mass transfer coefficient (cm d −1 ), accounts for the intrinsic ability for the soil to retain NO 3 -N by including depth, which is often ignored by using other removal rate constants (i.e., k (d −1 )) in first order decay model evaluations. The model can be expressed mathematically as [37]: After ρ FO was determined using the first 9 of the 18 mesocosm experiments, daily C t values were predicted using the following first order decay equation for the remaining 9 batch experiments:

Efficiency Loss (EL) Model
The third model, the efficiency loss model (EL), is similar to the first order decay model, but it accounts for the efficiency of the process rate relative to NO 3 -N concentration decline over time [45]. The surficial removal rates are proportional to the NO 3 -N concentration 'to the α', in which α, the rate order, is less than 1. Indeed, it has been theoretically shown in systems that have diffusion only processes and perfectly flat sediment, water column NO 3 -N removal rates are 1/2 order rates [49] The efficiency loss model thus essentially implies the apparent removal rate above an actual uneven sediment, which results in the exchange surface area that is higher than the projected area and must be adjusted to a rate between 0.5 and 1. Similar to the first order decay rate model, the model assumes that the substrate concentration is significantly smaller than K s , the system is well mixed, and has no significant influences from water losses or gains. However, the model assumes a power relationship represented by the coefficient α, in which the order is less than 1. The model can be expressed as: where ρ EL was the mass transfer coefficient (cm d −1 ), α was a unitless constant that varies between 0 and 1. The model was evaluated empirically in R Studio (2015) [50] to solve for α and ρ EL for individual mesocosms using 9 of the 18 batch run datasets. To test the model, daily C t values were then predicted for the remaining 9 batch experiments, using the ρ EL and averaged α values not used to determine removal coefficients. Observed NO 3 -N concentrations were compared to predicted daily NO 3 -N concentrations (C t ) from the mesocosm studies that were calculated using the following equation: The fourth model evaluated was the Monod model (M), also referred to as the Michaelis-Menten model for theoretical considerations, which often represents biologically mediated reactions that display first order decay kinetics at low concentrations and zero order kinetics at higher concentrations, resulting in a hyperbolic relationship between J M and NO 3 -N concentrations. Zero order kinetics has been observed when NO 3 -N concentrations exceeded biological demand of the system, while first order decay kinetics have been observed when NO 3 -N concentrations are below the biological demand of the system [51]. Therefore, the model ultimately interpolates between a zero order and first order decay models. The model assumes there are no intermediate or product inhibitions and the system is at steady state [52]. The model can be expressed as: where J M was the area based NO 3 -N loss rate (mg m −2 d −1 ), J max was the maximum removal rate achieved by the system (mg m −2 d −1 ), and K s (mg L −1 ) was the half-saturation constant. K s and J max were determined graphically using a Lineweaver-Burke plot with observed values from batch experiments [53]. Daily J M values were then calculated using the remaining 9 of 18 batch experiments with Equation 5, the empirically determined constants, and daily NO 3 -N concentrations. More specifically, K s was the concentration at which point the NO 3 -N removal rate was at half of the maximum NO 3 -N removal rate (J max ) for the WET-Min and WET-Org systems. The relationship between J M and C f was evaluated using a Lineweaver-Burke plot to empirically determine the maximum areal removal rate (J max ) and the half-saturation constants (K s ) for the Monod model for each treatment (WET-Min and WET-Org) and replicate.

Temperature Adjustment of Removal Rate Coefficients
The effect of temperature on the estimated removal coefficients of each batch run (J or ρ) was determined using a modified Arrhenius relationship [54]: where X was either J, the area based NO 3 -N loss (mg m −2 d −1 ) at temperature T ( • C), or ρ, the mass transfer coefficient (cm d −1 ) at temperature T ( • C), X 20 was either J 20 , the area based NO 3 -N loss (mg m −2 d −1 ) at 20 • C, or ρ 20 , the mass transfer coefficient (cm d −1 ) at 20 • C, and θ was an empirical temperature coefficient. A linear form was used to estimate J 20 and ρ 20 were determined for each treatment using the first 9 of the 18 batch run datasets.

Statistical Evaluation
The accuracy and the reliability of the four models for predicting NO 3 -N removal rates were compared utilizing three statistical parameters that are described in Table 2 [55,56]: coefficient of determination (R 2 ), relative root mean square error (RRMSE), and model efficiency (MEF). The statistical parameters were used to evaluate the correlation, differences, and variations between the predicted and measured NO 3 -N removal rates only in the 9 batch experiments that were used to validate the removal rate models. Measures the extent of linear correlation between two datasets.
Relative root mean square error: Measures the differences between the predicted and the measured values.
Measures the variations accounted for by the model.

NO 3 -N Removal and Areal Removal Rates
Hydraulic loading and water quality observations during the batch studies are presented in Table 3. Prior to each batch experiment, average daily recirculation rates were adjusted to 0.5 or 1 m 3 d −1 depending on depth to allow drainage water to turnover within each mesocosm once a day. Applied drainage water had NO 3 -N concentrations within 10-20% of the initial target concentrations. Significant NO 3 -N reduction was observed in both the WET-Min and WET-Org wetland systems during all batch runs, with NO 3 -N reductions as high as 99%. The highest removal rates were observed during the growing season, as expected. NO 3 -N removal in both wetland systems were significantly affected during fall and winter months when plants began to die back, average water temperatures were colder, and dissolved oxygen concentrations were higher, all of which likely limited plant and microbial uptake and inhibited denitrification [58,59]. These findings support the importance of developing removal coefficients at a base temperature (J 20 or ρ 20 ) that can be adjusted dependent on average water temperatures to more accurately predict fluctuations in wetland NO 3 -N removal rates [60].
Average and daily areal NO 3 -N removal rates were higher in the WET-Min mesocosms when compared to the WET-Org wetland mesocosms. Average removal rates in the WET-Min system ranged from 62 to 603 mg m −2 d −1 , while the WET-Org wetland systems had removal rates from 59 to 344 mg m −2 d −1 during the growing season when water temperatures were above 12 • C, with the exception of one fall experiment. Average removal rates dramatically decreased once water temperatures fell below 12 • C. Observed NO 3 -N removal rates from this study were similar to removal rates found in other microcosm and mesocosm wetland NO 3 -N removal studies during the growing season. Gebremariam and Beutel reported comparable NO 3 -N removal rates (175-500 mg m −2 d −1 ) in an evaluation of surface flow constructed treatment wetland mesocosms spiked with NO 3 -N at 15 mg L −1 during the summer [61]. Stringfellow et al. investigated three wetland restoration substrates receiving agricultural drainage water during the growing season utilizing microcosms and reported mean areal NO 3 -N removal rates of 142-380 g m −2 d −1 , also consistent with the estimated wetland mesocosm NO 3 -N removal rates in this study [52]. At a larger scale, a constructed surface flow wetland macrocosm study supplied with NO 3 -N concentrations of 8-10 mg L −1 observed NO 3 -N removal rates to be as high as 834 mg N m −2 d −1 [62]. Table 3. Summary of selected water quality parameters during batch experiments and NO 3 -N removal rate means ± standard deviations (SD).

Model Calibrations
Observed daily NO3-N concentrations from the first 9 of 18 batch experiments were used to calibrate models and determine kinetic constants for each model. Potential impacts of water depth and temperature were also investigated during the calibration stage of the models.
The zero and first order decay kinetic coefficients were calculated using data observed from each batch run experiment for individual mesocosms (Equations (2) and (4)). The zero order kinetic model had R 2 values that ranged between 0.78 and 0.99 with an average of 0.97 in the WET-Min system, while the WET-Org systems had R 2 values that ranged from 0.81 to 0.99 with an average of 0.97. In contrast, WET-Min treatments had R 2 values that ranged between 0.70 and 0.99 with an average of 0.95 in the first order decay models, while the WET-Org systems had R 2 values that ranged from 0.73 to 0.99 with an average of 0.94 for individual batch runs. The lower R 2 values were observed during the winter months when average water temperature was below 12 • C. All evaluations during the growing season had R 2 values above 0.90 for both the zero and first order decay fits to observed datasets. The efficiency loss model was evaluated in R Studio to empirically solve for ρ EL and α in each mesocosm during the first 9 batch experiments using daily NO 3 -N concentration observations at each sampling point. The two unknown values (α and ρ EL ) were determined for seven of the nine batch run datasets and resulted in good fits based on outputs from the R Studio model (p > 0.01 for ρ EL and p > 0.2 for α). The WET-Min treatments had an R 2 values that ranged between 0.93 and 0.99 with an average of 0.98, while the WET-Org systems had R 2 values that ranged from 0.95 and 0.99 with an average of 0.98.
R studio was unable to solve for α values in the efficiency loss model in a two batch runs when NO 3 -N removal periods were short or when the NO 3 -N change over time resembled a zero or first order decay relationship, which resulted in an α equal to 0 or 1, respectively. Therefore, α values that represented the individual wetland systems were difficult to determine since empirically determined α values often varied between 0 and 1 within the same treatment and batch run. Correlations between the α values and season, NO 3 -N load, temperature and DO were not observed. However, differences between the α values and water depth were observed particularly in the WET-Min wetland systems. Consequently, α values empirically determined in R Studio were separated into treatment and initial water depths (18 and 30 cm). Average α values were then calculated based on initial water depth for each treatment and used as the initial input of α for the efficiency loss model to fit observed NO 3 -N removal rates in the 9 of 18 batch run experiments used for calibration. The α values were then adjusted until the best fit was found between observed and predicted values.
The Monod model constants (K s and J max ) were estimated graphically using the first 9 of the 18 batches. J max and K s were then averaged across replicates of each treatment. The Monod model had R 2 values that ranged from 0.19 to 0.95 with an average R 2 of 0.65 in the WET-Min systems, while R 2 values ranged from 0.13 to 0.99 with an average R 2 of 0.55 in the WET-Org systems.

Temperature Adjustment for Removal Rate Coefficients
Removal rate coefficients for 20 • C were calculated for each mesocosm and model utilizing the first 9 of 18 batch run datasets and the Arrhenius relationship. Individual mesocosm values were then averaged for each treatment (i.e., each batch θ obtained for WET-Min and WET-Org was an average of the three replicates). Computed values can be found along with published values from past mesocosm and full-scale wetland studies for comparison in Table 4. To our knowledge, no studies have reported using the efficiency loss and zero order models at the mesocosm scale, therefore, comparisons were made with full-scale systems. In general, our results were consistent with observed values from past published wetland studies.

Validation of Kinetic Model Fits
Observed NO 3 -N removal rates from the remaining 9 batch run datasets (not used for calculating removal rate coefficients and other kinetic coefficients) were used to compare the daily predicted values produced by each of the developed kinetic models and evaluate the goodness of fit for each model. Kinetic coefficients were adjusted for average water temperature in each batch run. A regression of the predicted vs. observed removal rates for each of the four models are displayed in Figure 2.

Validation of Kinetic Model Fits
Observed NO3-N removal rates from the remaining 9 batch run datasets (not used for calculating removal rate coefficients and other kinetic coefficients) were used to compare the daily predicted values produced by each of the developed kinetic models and evaluate the goodness of fit for each model. Kinetic coefficients were adjusted for average water temperature in each batch run. A regression of the predicted vs. observed removal rates for each of the four models are displayed in Figure 2.

Prediction Power
The fit of the zero order and Monod models were weak between the predicted and observed values for all statistical parameters (Figure 2). Both models consistently underpredicted the removal rates for all NO3-N loads. The weak fits between the predicted and observed removal rates were not surprising. Based on the assumption for the zero order model, removal rates are not affected by changes in initial NO3-N concentrations, which tend to be observed at high NO3-N concentrations (>200 mg L −1 ) [36,64]. However, these models were evaluated due to little, if any literature, has evaluated these models for pulse flow systems. During every season, as NO3-N loads increased the observed removal rates also increased. Therefore, the assumptions required to use the zero order model were not met for these systems. Additionally, Ks values, empirically calculated as 5.96 mg L −1 for the WET-Min system and 5.63 mg L −1 for the WET-Org system, ultimately resulted in Ks and Co (initial concentrations) being similar.
Zero and Monod kinetic models are often used in wastewater treatment wetland systems where high initial NO3-N concentrations allow for the systems to become saturated resulting in a maximum removal rate (JZO), which in this case would be equivalent to the Jmax of the Monod kinetic model. However, based on observed removal rates, the wetland mesocosm systems never reached a saturation point, which resulted in poor estimates for the removal coefficients for both models. Therefore, based on observations in this study, these predictive models are not recommended for wetlands loaded with low NO3-N concentrations.
Stronger statistical results were observed for the first order decay and efficiency loss models with R 2 ≥ 0.84 (Figure 2). Each of these models are dependent on initial NO3-N loads and do not require the determination of saturation coefficients such as Ks. Additionally, the first order model and the efficiency models had similar fits. The stronger predictive power of the first order decay and efficiency loss models compared to the zero and Monod models was due to chemical and physical based assumptions of the mathematical theories in terms of NO3-N reduction. Chemically the zero order model assumes substrate concentrations at all points within the water column is notably greater than half-saturation concentrations, while the Monod model required NO3-N saturated environments at least one point during the experiments. NO3-N half-saturation constants can only be determined

Prediction Power
The fit of the zero order and Monod models were weak between the predicted and observed values for all statistical parameters (Figure 2). Both models consistently under predicted the removal rates for all NO 3 -N loads. The weak fits between the predicted and observed removal rates were not surprising. Based on the assumption for the zero order model, removal rates are not affected by changes in initial NO 3 -N concentrations, which tend to be observed only when NO 3 -N concentrations are high (>200 mg L −1 ) [36,64]. However, these models were evaluated due to little, if any literature, has evaluated these models for pulse flow systems. During every season, as NO 3 -N loads increased the observed removal rates also increased. Therefore, the assumptions required to use the zero order model were not met for these systems. Additionally, K s values, empirically calculated as 5.96 mg L −1 for the WET-Min system and 5.63 mg L −1 for the WET-Org system, were similar to C o (initial concentrations) and thus did not meet the definition of the half saturation constant required to make the continued use of the Monod model justifiable.
Zero and Monod kinetic models are often used in wastewater treatment wetland systems where high initial NO 3 -N concentrations allow for the systems to become saturated resulting in a maximum removal rate (J ZO ), which in this case would be equivalent to the J max of the Monod kinetic model. However, based on observed removal rates, the wetland mesocosm systems never reached a saturation point, which resulted in poor estimates for the removal coefficients for both models. Therefore, based on observations in this study, these predictive models are not recommended for wetlands loaded with low NO 3 -N concentrations.
Stronger statistical results were observed for the first order decay and efficiency loss models with R 2 ≥ 0.84 (Figure 2). Each of these models are dependent on initial NO 3 -N loads and do not require the determination of saturation coefficients such as K s . Additionally, the first order model and the efficiency models had similar fits. The stronger predictive power of the first order decay and efficiency loss models compared to the zero and Monod models was due to chemical and physical based assumptions of the mathematical theories in terms of NO 3 -N reduction. Chemically the zero order model assumes substrate concentrations at all points within the water column is notably greater than half-saturation concentrations, while the Monod model required NO 3 -N saturated conditions at least one point during the experiments. NO 3 -N half-saturation constants can only be determined once the systems are saturated (meaning the microbial and biomass pools are removing the maximum amount of NO 3 -N) [25], which was never observed during these experiments. In contrast, the first order and efficiency loss models do not require a saturation point. Physically, the water was recirculated slowly within the mesocosms, as expected in the future restored wetlands, which likely enhanced the ability of NO 3 -N to reach denitrifying hot-spots through advection and diffusion in the system [65], which further limited these systems from becoming saturated.
Both the first order decay and efficiency loss models were further assessed to determine if defined parameters, particularly retention time, impacted the prediction power of the models. Residence time of surface water in the full-scale wetland environments must stay within the range of 1 to 10 days to preserve established and recently planted trees [66,67]. Therefore, both predicted and observed values were compared (Figure 3). Figure 3 shows the first order decay model has a tendency to over predict NO 3 -N concentrations at values < 3 to 4 mg L −1 , while the model typically under predicts NO 3 -N concentrations at values > 3 to 4 mg L −1 for both the WET-Min and WET-Org wetland mesocosms systems (Figure 3a). Similarly, the efficiency loss model typically under predicted at higher NO 3 -N concentrations (Figure 3b).
Water 2017, 9,517 14 of 20 once the systems are saturated (meaning the microbial and biomass pools are removing the maximum amount of NO3-N) [25], which was never observed during these experiments. In contrast the first order and efficiency loss models do not require a saturation point. Physically, the water was recirculated slowly within the mesocosms, as expected in the future restored wetlands, and likely enhanced NO3-N to reaching denitrifying hot-spots through advection and diffusion in the system [65], further limiting these systems from becoming saturated. Both the first order decay and efficiency loss models were further assessed to determine if defined parameters, particularly retention time, impacted the prediction power of the models. Residence time of surface water in the full-scale wetland environments must stay within the range of 1 to 10 days to preserve established and recently planted trees [66,67]. Therefore, both predicted and observed values were compared (Figure 3). Figure 3 shows the first order decay model has a tendency to underpredict NO3-N concentrations at values < 3 to 4 mg L −1 , while the model typically underpredicts NO3-N concentrations at values > 3 to 4 mg L −1 for both the WET-Min and WET-Org wetland mesocosms systems (Figure 3a). Similarly, the efficiency loss model typically underpredicted NO3-N concentrations at lower NO3-N concentrations and underpredicted at NO3-N higher concentrations (Figure 3b).   Additionally, to further examine the prediction power of the first order and efficiency loss models, the predicted C 7 values (NO3-N concentrations after 7 days within the wetland mesocosm) for the first order decay and efficiency loss models were examined and compared to observed NO3-N concentrations for the 9 batch run experiments used for validation ( Figure 4).
Additionally, to further examine the prediction power of the first order and efficiency loss models, the predicted C7 values (NO3-N concentrations after 7 days within the wetland mesocosm) for the first order decay and efficiency loss models were examined and compared to observed NO3-N concentrations for the 9 batch run experiments used for validation ( Figure 4). Therefore, based on the assessment, both the first order decay and efficiency loss models had a tendency to predict lower final NO3-N concentrations during the end of a batch run than was observed, which is consistent with model results in Figure 3 for predicted NO3-N.
Model fits between predicted and observed NO3-N removal rates in both the first order decay and efficiency loss models were similar ( Figure 2). However, in the efficiency loss model the α constant must be empirically calculated for several batch runs to be accurately determined for alternative systems than the ones evaluated in this study, which is often not practical for wetland Therefore, based on the assessment, the first order decay model had a tendency to over predict NO 3 -N concentrations at lower NO 3 -N concentrations (< 3 to 4 mg L −1 ), while both models typically under predicted NO 3 -N concentrations at higher NO 3 -N concentrations values (> 3 to 4 mg L −1 ) which was consistent with model results in Figure 3.
Model fits between predicted and observed NO 3 -N removal rates in both the first order decay and efficiency loss models were similar ( Figure 2). However, in the efficiency loss model the α constant must be empirically calculated for several batch runs to be accurately determined for alternative systems than the ones evaluated in this study, which is often not practical for wetland designers. Additionally, the first order decay model is expected to provide slightly more conservative predictions than the efficiency loss model particularly during the growing season when water retention times are most critical for wetland tree survival. Therefore, the first order decay model was determined as the most suitable prediction model for NO 3 -N removal rates to be used in the design of these future wetland restoration systems.

Theoretical Wetland Predictions
To demonstrate the use of the first order decay NO 3 -N removal model determined at the mesocosm scale in this study, an initial estimate of the maximum hydrologic loading capacity for two future 210 ha wetland restorations (with the modeled soil types) were conducted. The sites were expected to receive agricultural drainage with NO 3 -N concentrations on average of 2.5 mg L −1 . Possible target NO 3 -N effluent concentrations considered were 0.1 mg L −1 (96% reduction), 0.5 mg L −1 (80% reduction), 1.0 mg L −1 (60% reduction), and 1.75 mg L −1 (30% reduction). The maximum acceptable hydrologic loading rates for the wetlands to achieve these target effluent concentrations were determined across a range of average water temperatures.
Based on the estimated removal rates, the maximum loading rates of agricultural drainage water that could be applied were determined to meet the range of target outflow concentrations at the outlet of wetland areas (C eff ) utilizing the following equation [24,68]: where L is the maximum loading rate into the wetland (m d −1 ), C eff is the effluent NO 3 -N concentration (mg L −1 ), C in is the influent NO 3 -N concentration (mg L −1 ), T is the temperature ( • C), ρ 20 is the  (Table 4). Maximum hydrologic loading capacities for the two wetland systems for the various NO 3 -N removal goals and average water temperatures can be found in Table 5. Based on an aggressive 96% removal goal, the wetland restorations similar to the WET-Min and WET-Org mesocosm systems were predicted to have the capability to receive 1.5 and 1.2 cm day −1 of drainage water, respectively, assuming the drainage water has NO 3 -N concentrations of 2.5 mg L −1 and an average water temperature of 20 • C. Current local rules developed for watersheds such as the Neuse and Tar-Pamlico in North Carolina require non-point source N pollution to be reduced by only 30%. If this were the treatment goal, and under the same assumptions described earlier, wetland restorations similar to the WET-Min and WET-Org mesocosm systems would have the capacity to receive up to 13.1 and 10.8 cm day −1 of drainage water, respectively, to meet the target NO 3 -N removal threshold. Table 5. Predicted maximum hydrologic loading capacities for the future restored wetlands assuming n of 0.95, and C in of 2.5 mg L −1 using the first order decay kinetic model.

Conclusions
In this study, the potential removal rates of NO 3 -N for two distinct pulse flow wetland environments slated to receive agricultural drainage water in future wetland restoration projects were evaluated at the mesocosm scale. Three kinetic models that are often utilized to predict NO 3 -N removal rates within wetland systems and one that has never been evaluated (efficiency loss) were assessed. Based on the results of this study, the following conclusions were drawn: 1.
The first order decay and efficiency loss kinetic models provided stronger statistical agreement between predicted and measured NO 3 -N removal rates compared to other models. 2.
The first order decay model was determined the most practical model due to its conservative predictions, simplicity, and reasonable fit compared to the efficiency loss model, NO 3 -N removal rates developed at the mesocosm scale have also been reported to be conservative estimates for full-scale wetlands. Therefore, the first order decay model should provide reasonable but conservative predictions of NO 3 -N removal rates, particularly if maximum residence times are limited, such as to ensure tree survivability during the growing season in restored forested wetlands.

3.
ρ 20 values determined for the first order decay model for NO 3 -N removal in the WET-Min and WET-Org wetland mesocosms were 4.9 cm d −1 and 4.1 cm d −1 . θ values were estimated to be 1.15 and 1.09 for the WET-Min and WET-Org systems, respectively. These values can be used to develop informed water management plans for the restoration sites. This is important to ensure the wetlands are not loaded with NO 3 -N that exceeds their assimilation capacity; a situation that would lead to higher NO 3 -N export and possible unintended eutrophication of downstream water bodies.