Unsteady state series CSTR modeling of removal of ammonia nitrogen from domestic wastewater treated in a vertical flow constructed wetland

Received Apr 22, 2019 Revised Oct 4, 2019 Accepted Oct 12, 2019 This work shows simulation results for subsurface vertical flow constructed wetland (VFCW) using a series CSTR model. The VFCW considered received the outflow from a domestic wastewater treatment plant. In addition, it was planted with Cyperus sp. and filter media was unsaturated. The model was based on an unsteady state mass balance for ammonia, nitrites, and nitrates, using one to three series CSTRs. Nitrogen transformation mechanisms considered were ammonification, nitrification, plant uptake and denitrification. The following effects were evaluated: the number of reacting CSTRs from one to three; the occurrence of the reaction in second and third CSTRs for the case that three CSTRs hold; the use of either equal or different values of reaction rate parameters between CSTRs; and the discretization of the reaction rate parameters. The inflow and outflow measurements of ammonium, nitrites, and nitrates were used for model calibration. The estimated parameters included the reaction rate coefficients and reactor water volume. The coefficient of determination (R) evidenced a satisfactory capability of simulating outlet pollutant concentrations. Two and three reacting CSTRs achieved similar R value (0.54-0.55), whereas one reacting CSTR achieved an R of 0.39, and three CSTRs with reaction only in the first tank achieved an R of 0.42. Discretization of the nitrification rate for the case of two reacting CSTRs led to an R of 0.94. The parameter sensitivity analysis revealed a significant effect of model parameters on the R value.

the effect of the following model features: the number of reacting CSTRs from one to three; the occurrence of reaction in second and third CSTRs for the case that three CSTRs hold; the use of either equal or different values of reaction rate parameters between CSTRs; and finally the discretization of the reaction rate parameters.

A brief review of pollutant removal in SSFCWs
Nitrogen removal in CWs encompasses a number of complex processes, including ammonification, nitrification-denitrification, adsorption, uptake by plants and uptake by microorganisms [2]. Nutrient removal in SSFCWs comprises a limited effect of plant absorption in comparison to biological degradation. Indeed, in the case of high nutrient loadings, plant uptake is negligible in comparison to microbially-mediated removal [19]. The amount of nitrogen plant uptake depends on the plant's growth requirements rather than nitrogen loading [23,24].
The nitrogen removal in SSFCWs depends on hydraulic retention time (HRT), temperature, vegetation type and the filter media [25]. Biological nitrification and denitrification are regarded as the main processes for the removal of ammonia in subsurface flow constructed wetlands (SSFCWs) [2]. It has been noted that nitrification is the limiting step for nitrogen removal. Nitrification is strictly aerobic, and it only occurs if DO is available [26]. It is dependent on dissolved oxygen, pH and temperature. The nitrification rate is significantly affected by DO and pH, whereas temperature has a lower effect [16]. Aeration significantly enhances nitrification, whereas low oxygen concentration and temperature have a decreasing effect [27,28]. The optimum pH range is 7.0 -9.0; whereas pH below 7.0 has an adverse effect [16]. A high content of organics leads to nitrification inhibition because heterotrophic organic degradation depletes the dissolved oxygen [29,30].
In [25] nitrogen and organics removal efficiencies were reported for VFCWs in tropical and subtropical regions. Therein, NH4 + -N removal percentages range from 45% to 99% (average of 66%), and outflow NH4 + -N is below 52.0 mg/L with average of 18.0 mg/L, whereas BOD removal efficiencies average 88%. In [31] nitrification rates in VFCWs were reported, which ranged from 0.01 to 2.15 g N/(m 2 day), with an average of 0.048 g N/(m 2 day). It has been observed that background concentration for ammonia removal is nearly zero [24].

System description
In this study, the VFCW experiment reported by [32] is considered. The system is located at Pereira (Colombia) and planted with Cyperus sp. It receives the outflow from a domestic wastewater treatment plant, which includes a septic tank and an upflow anaerobic filter. The length, width, and depth of the VFCW are 8.65 m, 5.0 m, and 0.80 m, respectively, resulting in a surface area of 43.25m 2 . The porous media consists of medium and thick size gravel. The data were taken on a weekly basis over ten weeks, starting in March 2012. The inflow was intermittent downflow, with three pulses per hour. The average inflow discharge was 314 L/h, and the hydraulic loading rate was 17.4 cm/d. The inflow pH ranged from 7.4 to 7.6 (average of 7.5), whereas the outflow pH ranged from 7.46 to 7.61 (average of 7.5). The inflow water temperature ranged from 24.0 to 26.0°C (average of 25°C), and the outflow temperature ranged from 25.0 to 27.0°C (average 26°C). The inflow COD ranged from 208.7 to 318.8 mg/L (average of 249.1 mg/L) and the inflow NO3ranged from 3.9 to 7.9 mg/L (average 5.0 mg/L). The outflow NO3ranged from 3.9 to 48.8 mg/L (average 27.8 mg/L).
The high outflow concentrations of nitrite and nitrates are due to the low denitrification capability, what is typical of properly aerated VFCWs. The time courses of outflow ammonium and (NO2 -+NO3) involved unexpected effects, mainly abrupt changes. Outflow ammonia exhibited an initial increase and an abrupt decrease after ca. 35 days, whereas outflow (NO2 -+NO3 -) exhibited a significant increase from 5 mg/L at the initial time, to ca. 60 mg/L in the last ten days. These abrupt changes are attributable to the adaptation of the wetland to the incoming wastewater during the start-up process. This is concluded because wetland adaptation to wastewater lasted only one month [32], and similar effects are exhibited by the HFCW start-up described by [24]. CWs start-up involves the establishment of sorption process and adaptation of the microbial community to the loading and aeration. These effects can occur for many months, thus affecting the removal percentage [23,24]. An example of nitrogen performance in HFCW start-up is shown by [24]. The time course of ammonia outflow concentration involves an increase, a smooth decrease, and an abrupt decrease, such that the corresponding removal percentage comprises a decrease, an increase, and an abrupt increase. In contrast, the outflow nitrate concentration departs from zero, follows a smooth bell-shaped curve along the timeline and vanishes at the end.

Mass balance model
Consider ammonification, nitrification, plant uptake and denitrification as being the primary nitrogen removal and formation pathways, and assume: i) ammonium is converted to nitrite and nitrate in one step; ii) plant growth rate is constant. Thus, the mass balance for nitrogen concentrations across a single CSTR gives [16,33].
Where ON is the organic nitrogen concentration, DO is the dissolved oxygen concentration, Q is the flow rate, V is the effective volume, which is the product of volume and porosity of each CSTR; NH4,in is the inflow NH4 + -N concentration; ra is the ammonification rate, rn is the nitrification rate, rp is the plant uptake rate, rd is the denitrification rate. The model parameters are ka, kn, KDO, kAN, kd, kp. CpH is a linear function of pH and CT is an exponential function of temperature. Both Monod and first-order expressions can be used to describe ammonia removal in wetlands [24]. When Monod kinetics exhibit values of the halfsaturation constant significantly higher than pollutant concentration (C<<K), a first order expression is more suitable [34]. Thus, in this study both first order and Monod kinetics in terms of NH4 + -N concentration are evaluated for the nitrification rate. It is worth mentioning that NO2and NO3are produced by nitrification and are consumed by denitrification [16]. In this study, experimental data evidenced a significant formation of (NO2 -+NO3 -), meaning that the denitrification rate (rd) is small with respect to the nitrification rate (rn). Therefore, the (NO2 -+NO3 -) mass balance is considered in this study, in order to simulate its formation.
The following features of the VFCW are considered: i) available VFCW measurements included neither concentrations of organic nitrogen nor dissolved oxygen; ii) the water temperature remained approximately constant, wherein the outflow values ranged from 25.0°C to 27.0°C and averaged 26.0°C; iii) pH variations were negligible, as outflow pH ranged from 7.46 to 7.61 and averaged 7.5. Therefore, the following simplifications are made: i) the temperature and concentrations of organic nitrogen and DO are assumed constant; ii) the term (ra-rp) is assumed to be constant in the NH4 + -N mass balance (ra-rp=kap); iii) the term rn is simplified as Incorporating these simplifications and arranging the notation, gives Where NHin and NH are the inflow and outflow NH4 + -N concentrations, respectively, whereas NOin and NO are the inflow and outflow (NO2 -+NO3 -) concentrations, respectively. Now, we consider three series CSTRs, to account for non-ideal flow mixing through the VFCW. Parameter values are considered to be different between CSTRs, in order to account for the fact that removal efficiencies and reactions rates are different between filter media layers. Also, different CSTR volumes are used in order to account for different media porosities. This is in accordance with [35], who found that NH4 + -N and NO3 --N removal in VFCWs is different along the filter bed, and with [16], who used different parameter values between CSTRs, where each CSTR corresponded to one filter media layer. The subscript I={1,2,3} indicates the i-th CSTR. The mass balance for the concentrations of NH4 + and (NO2 -+NO3 -) across the three reacting CSTRs gives Where NH1, NH2, NH3 are the NH4 + -N outflow concentrations in first, second and third CSTRs, respectively, whereas NO1, NO2, and NO3 are the outflow concentration of (NO2 -+NO3 -) in first, second and third CSTRs, respectively. The general modeling structure is presented in Figure 1.

Model calibration and simulation
Model calibration used the dataset corresponding to the VFCW experiments reported by [32] and described above. The overall flowrate and the concentrations of NH4 + , (NO2 -+NO3 -) at VFCW outflow and inflow are available, but concentrations at intermediate bed locations were not. Therefore, the model output of last CSTR is compared to the VFCW outflow, but comparison of first and intermediate CSTRs to intermediate bed measurements is not possible. In order to simulate the VFCW performance, the proposed model is calibrated by minimization of the sum of the squared errors between the model output and the experimental data, using the least squares criterion (see Appendix A). Matlab 2014 software (the MathWorks Inc., Natick, Mass.) with the 'fmincon' subroutine was used to perform minimization, with the ode113 command for integration of differential equations. The aforementioned staged behavior of ammonium nitrogen removal hampers accurate simulation; hence discrete parameter values are used. Reaction rates can be modeled by using coefficients that take on discrete values over time, in order to account for time-dependent changes [36]. Examples of this type of parameters for batch biological reaction are presented by [37,38]. Thus, in this study an additional model calibration case is performed, which involves discrete values for the ammonium nitrogen model parameters.
The prediction capability of the model was assessed by the coefficient of determination (R 2 ) (see Appendix B). The coefficient of determination was evaluated for the effect of the following model features: number of reacting CSTRs from one to three; the occurrence of reaction in second and third CSTRs for the case that three CSTRs hold; the use of either equal or different values of reaction rate parameters between CSTRs; and the discretization of the reaction rate parameters as shown in Table 1. In summary, the main features of simulation and parameter estimation are:  Model calibration compared pollutant concentrations at the last CSTR with concentrations at VFCW outflow.  All model coefficients were estimated by sum-of-squares minimization, and none of them was borrowed from other studies; the volumes of the three CSTR were different, and they were estimated by sum-of-squares minimization.  The coefficient of determination was evaluated for the effect of several model features, including number of reacting CSTRs, reaction occurrence, equal/different model parameters between CSTRs, as shown in Table 1.  An additional model calibration was performed which involved discretized values of the ammonium nitrogen model parameters. Also, the capability of the used models for simulating VFCW outflow concentrations of NH4 + and (NO2 -+NO3 -) was assessed by drawing the predicted concentrations as a function of observed VFCW outflow concentrations. The precision of the estimated parameters was assessed by the sensitivity analysis used by [39]. It comprised the effect of each estimate on the coefficient of determination, with only one estimate varied at a time, and a -50 to 50% variation range.

Model fitting
The cases of two and three reacting CSTRs (models A1 and A2) achieved similar fitting quality, being higher than in the following cases: one reacting CSTR followed by a non-reacting CSTR (model B1); a single reacting CSTR (model A3); two reacting CSTRs, with model parameters being the same between CSTRs (model C1). These results are in agreement with the VFCW study of [16], where three reacting series CSTRs were used with model parameters being different between CSTRs.
The coefficient of determination for NH4 + simulations is shown in Table 1. Taking the case of two reacting CSTRs (model A2) as the reference model (R 2 of 0.54), the following is concluded: i) two series CSTRs, with reaction only in the first one (model B1) yields a slight decrease of R 2 , from 0.54 to 0.51; ii) two series reacting CSTRs, with the same model parameters excepting volume (model C1) yielded a significant decrease of R 2 , from 0.54 to 0.33. Furthermore, the discrete estimation of nitrification coefficients (model A2b) resulted in an increase of the coefficient of determination from 0.54 to 0.94. Hence, the model comprising two series reacting CSTRs with different model parameters and discrete parameters is suitable. The observed VFCW outflow concentrations of NH4 + and the A2 and A2b model simulations are shown in Figure 2 to Figure 5, confirming its correlation.   Estimated parameters of models A2, A2b, and B1 for NH4 + simulation are shown in Table 2. For model A2b, only parameter kap exhibited time-varying behavior. The VFCW outflow concentrations of (NO2 -+NO3 -) and the A2 and A2b model simulations are shown in Figure 6 to Figure 8, confirming its correlation.   Table 3. Using time-varying coefficients of rn and rap yielded a slight decrease in the coefficient of determination for (NO2 -+NO3 -) from 0.79 to 0.70.

Sensitivity of estimated parameters
Sensitivity analysis for model A2 parameters is shown in Figure 9 (refer Appendix C), consisting of the effect of parameters on the coefficient of determination related to NH4 + simulation. The order of parameters, from higher to lower effect on the coefficient of determination is kns1, knmx1, V2, kap2, V1, knmx2, kap1. Hence, the estimated kap1 has the lowest accuracy, whereas estimated kns1 has the highest accuracy.Sensitivity analysis for model A2b parameters and BOD simulation is shown in Figure 10 (refer Appendix C). The order of parameters, from higher to lower effect on the coefficient of determination is kns1, knmx1, kap1, kap2, V2, V1, kn2. Hence, the estimated kn2 has the lowest accuracy, whereas estimated kns1 has the highest accuracy.

CONCLUSION
The coefficient of determination between VFCW outflow concentrations of NH4 + and (NO2 -+NO3 -) and simulations confirms the CSTR-like behavior of VFCWs stated by [8,18,13]. Also, it enabled the identification of the model features that lead to improved simulations. Indeed, it was found that the model comprising two reacting series CSTRs with parameter values that are different among CSTRs and take on discrete values over time is suitable. The cases of two and three reacting CSTRs have similar fitting quality and are more suitable than a single reacting CSTR. The higher fitting quality obtained for reaction parameters with different values between CSTRs in comparison to the case with the same values is in agreement with the modeling study performed by [16], where different parameter values between CSTRs were used because of the different features of filter media layers. The lack of pollutant measurements at intermediate bed locations did not avoid using series CSTRs modeling, but posed uncertainty on the setting of the initial time value of pollutant concentrations. Nevertheless, setting this value as the VFCW outflow pollutant concentration at the initial time resulted in proper fitting quality. The estimated coefficients of the nitrification rate indicate that a first-order expression is more suitable than Monod expression. Indeed, the values of Monod half-saturation constant are significantly higher than NH4 + concentration. The discretization of the coefficients of the reaction rates of the ammonium nitrogen model yielded a significant improvement of fitting quality. This strategy enabled addressing the staged increase of ammonium nitrogen removal.
The parameter sensitivity analysis confirms that the estimates correspond to an optimum R 2 value, which is attained by using the sum of squares minimization, and indicates that most model parameters have a significant effect on the R 2 value. To use the modeling results for sizing of VFCws, further experiments with different hydraulic loadings must be conducted. Furthermore, the large scale VFCW must have the same features of the VFCW experiment (for instance bed porosity, vegetation, filter media depth, frequency of intermittent loading), and input data must include inflow rate, inflow concentrations of BOD and NH4 + , and the corresponding target limit concentrations. Where θ is the vector of model parameters; N is the number of data points; NH4|exp is the observed concentration of NH4 + at VFCW outflow, NH4|calc is the NH4 + model simulation; NO|exp is the (NO2 -+NO3 -) observed concentration at VFCW outflow, NO|c is the (NO2 -+NO3 -) model simulation.

Appendix B: Coefficient of determination
The definition of the coefficient of determination (R 2 ) provided by [38] is considered in this study. The expression of R 2 for the simulation of NH4 + and (NO2 -+NO3 -) at VFCW outflow is Where N is the number of data points, is the mean of NH4|exp, and is the mean of NO|exp. Further discussion of the coefficient of determination is given by [38].
Appendix C: Figure 9. Sensitivity analysis for model A2 parameters and NH4 + simulation (Model A2 involves reaction occurring in the first and second CSTRs; Monod and first-order kinetics used for nitrification in first and second CSTRs, respectively) Figure 10. Sensitivity analysis for model A2b parameters and NH4 + simulation (Model A2b involves reaction occurring in first and second CSTRs; nitrification parameters are time-varying; Monod and first-order kinetics are used for nitrification in first and second CSTR, respectively)