A hybrid constructed wetland for organic-material and nutrient removal from sewage: Process performance and multi-kinetic models

A pilot-scale hybrid constructed wetland with vertical flow and horizontal flow in series was constructed and used to investigate organic material and nutrient removal rate constants for wastewater treatment and establish a practical predictive model for use. For this purpose, the performance of multiple parameters was statistically evaluated during the process and predictive models were suggested. The measurement of the kinetic rate constant was based on the use of the first-order derivation and Monod kinetic derivation (Monod) paired with a plug flow reactor (PFR) and a continuously stirred tank reactor (CSTR). Both the Lindeman, Merenda, and Gold (LMG) analysis and Bayesian model averaging (BMA) method were employed for identifying the relative importance of variables and their optimal multiple regression (MR). The results showed that the first-order–PFR (M2) model did not fit the data (P > 0.05, and R2 < 0.5), whereas the first-order–CSTR (M1) model for the chemical oxygen demand (CODCr) and Monod–CSTR (M3) model for the CODCr and ammonium nitrogen (NH4−N) showed a high correlation with the experimental data (R2 > 0.5). The pollutant removal rates in the case of M1 were 0.19 m/d (CODCr) and those for M3 were 25.2 g/m2∙d for CODCr and 2.63 g/m2∙d for NH4-N. By applying a multi-variable linear regression method, the optimal empirical models were established for predicting the final effluent concentration of five days' biochemical oxygen demand (BOD5) and NH4-N. In general, the hydraulic loading rate was considered an important variable having a high value of relative importance, which appeared in all the optimal predictive models.


Introduction
The processes of eliminating organic matter and nutrients from polluted wastewaters are quite complicated. They include chemical, physical, and biological factors that demand triply synergistic processes for maximum efficiency. These factors of the constructed wetland (CW) affect the removal efficiency and are affected by shifting local and regional conditions (Gholizadeh et al., 2015;Vo et al., 2018;Vo et al., 2017;Wu et al., 2018a;. As a result, an appropriate CW design that compares favorably to those of other studies is more difficult to achieve. The misunderstanding of the contaminant dynamics of this system can lead to design failures that usually cause a degraded efficiency of the targeted pollutants' removal, clogging, and short circulation in the design's fluid dynamics (Davoodi et al., 2016;Samsó et al., 2016;Wu et al., 2018a;. Therefore, the application of the knowledge of pol-lutant removal kinetics to this data could be expected to ensure positive results if and when these conclusions are effectively applied to the engineering challenges of CW design. There are several models and methods for predicting pollution removal, and the first-order model (Chan et al., 2008;Trang et al., 2010) and linear regression (Babatunde et al., 2011b;Gholizadeh et al., 2015;Reed and D., 1995;Sheridan et al., 2013) are popular. Some approaches 2 3 2 3 2 2 (2) 2 3CO 2 H 2 O 2MnO 2 2Mn recommended recently are the use of artificial neural networks, multicomponent reactive transport module, principal component analysis, clustering tree diagrams, classification and regression trees, and redundancy analysis (Akratos et al., 2008;Babatunde et al., 2011b;Dong et al., 2012;Hijosa-Valsero et al., 2011;Huang et al., 2014). However, these may prove to be overly complicated owing to their numerous empirical parameters. The multiple regression (MR) method has been considered useful for simplified description, analysis of CW performance, and developing predictive models (Murray-Gulde et al., 2008;Tomenko et al., 2007). The majority of the previous studies that were focused on the evalua-tion of the efficacy of MR (Babatunde et al., 2011b;Chan et al., 2008;Hijosa-Valsero et al., 2011) used the response variables of the effluent concentrations while avoiding the application of the removal rate.
To accurately compute MR, the relative importance and optimal models of the given predictors are required, in which relative importance is the proportionate contribution each predictor makes to R 2 (Johnson, 2000). For selecting the optimal models, the previous research was based primarily on R 2 (Babatunde et al., 2011b;Hijosa-Valsero et al., 2011). However, their methods did not involve the evaluation of the interaction analysis and multi-collinearity, which can influence the accuracy of these results. When a model has a large number of variables, it results in a higher R 2 value, which might cause over-fitting. Bayesian model averaging (BMA) is considered to be a good tool for overcoming these limits. By performing averaging over several different competing models, BMA integrates model uncertainty layers of various particle sizes. The filter layers comprised various gravel grades with the total height of the VF being 0.7 m and that of the HF being 0.5 m. This corresponded to 2.65, 1.32, and 1.76 days of hydraulic retention times for HLR 1 , HLR 2 , and HLR 3 , respectively. More details of the characteristics of the sewage, HCW, and operating procedures are described in a previous report by Nguyen et al. (2017).
Wastewater samples were collected from the influent of the VF (SP1) and effluent of the VF (SP2) and HF (SP3) as shown in Fig. 1  . The sampling rate was once a week for more than 6 months of operation. Twenty-three water samples of HCW were ob-tained and tested in this experiment. The pH was determined using a multiparameter water quality meter (HQ40d, Hach, USA). The COD Cr , BOD 5 , total suspended solids (TSS), NH 4 -N, nitrate nitrogen (NO 3 -N), phosphate (PO 4 -P), TN, and total coliforms (TCol) were analyzed according to the Standard Methods 5220D, 5210B, 2540D, 4500-NH 3 F, 4500-NO 3 B, 4500-P D, 4500-P J, and 9221 B, respectively (APHA/ WEF/AWWA, 2005). The spectrophotometer (Cary 60, Agilent Technologies Inc., USA) was used to measure the COD Cr , NH 4 -N, NO 3 -N, and PO 4 -P.
Owing to the intermittent flow, the aerobic and anoxic conditions were dominant in the HCW. Oxygen is a terminal electron acceptor that is reduced while electron donors (mainly organic matter and ammonia) are oxidized, and CO 2 and H 2 O are formed as end products (1). In the bottom layers of the HCW, other reactions (2, 3, and 4) might occur in response to the decrease in oxygen and the redox potential (Kadlec and Wallace, 2009).
into the prediction and estimation of the parameters (Fang et al., 2016;Hoeting et al., 1999).
collinearity and obtain the precise R 2 value. The obtained results clarstructed wetland (HCW), and the optimal models that have not been adequately addressed in previous studies were also selected. Furthermore, by applying more response variables such as hydraulic loading rates (HLRs), we developed more effective and accurate predictive models that contributed significantly in the design, management, and maintenance of the HCW system. The main purpose of this study is to clarify the performance of the HCW and explore the predictive models using kinetics and MR. This study is mainly focused on the comparison of the kinetics, weighting correlation, and relative importance and the development of optimal predictive models for HCW performance using the BMA method.

Description of hybrid constructed wetland and data collection
To clarify the adaptation and removal performance of HCW, a pilot system was installed at the sewage treatment plant of Dong Ha city,

Kinetic models
The first-order and Monod kinetic models were used to describe and evaluate the pollution degradation in this study. The hydrodynamic pattern in the HCW was considered to be the same as that of the plug flow reactor (PFR) and continuously stirred tank reactor (CSTR). Therefore, the new models have been developed by integrating the firstorder and Monod models with PFR and CSTR.

First-order k-C* model with PFR (M 1 )
By assuming an exponential removal rate to reflect a non-zero background wetland concentration (C*) (Ali et al., 2018;Kadlec and Knight, 1996), the removal model based on the first-order and C* models is expressed as follows (Eq. (5): Vietnam, and was then tested in over an operating period of 190 days.

1
The sewer at this location collects the municipal wastewater and runoff water from Dong Ha city. The water level of this sewer fluctuates from

First-order model with CSTR (M 2 )
To model the correlation between the influents and effluents of CW for nutrients and organic matter, the integrated model that combined first-order kinetics with CSTR (Saeed and Sun, 2011) was established. The combined model is expressed as shown in Eq. (6): In this work, we used the new technique of BMA to overcome multiified the relative importance of the influent factors of a hybrid con-the VF and denitrification in the HF in order to reduce the nutrient content in the effluent. Wastewater was pumped from the sewage system (twice a day for an interval of 30 min), stored in the storage tank, and was allowed to flow intermittently to the VF at various HLRs of 0.44 m 3 /d (first stage, HLR 1 ), 0.88 m 3 /d (second stage, HLR 2 ), and 0.66 m 3 /d (final stage, HLR 3 ). The effluent of the VF was drained continuously by gravity into the HF tank. The characteristics of the VF and HF are presented in Table 1. Each tank was built with three filter

Monod kinetics with CSTR (M 3 )
This model combined the Monod kinetics, which comprises half the saturation constant of the limiting substrate and effluent concentration, with the CSTR flow pattern (Saeed and Sun, 2011). It is expressed as shown in Eq. (7) Table 1 Characteristics of hybrid constructed wetland.
applied to determine relative importance or ranking, and BMA was applied for the optimal model selection process. Consequently, using BMA, the optimal MR models that were based on the Bayesian information criterion (BIC) and R 2 were selected. In general, the best models have the lowest BIC value with highest adjusted R 2 (Hoeting et al., 1999).

Statistical analyses
All the statistical analyses of this study were performed using the computing environment R (Version 3.4.0). The statistical differences of the experimental results were evaluated through an analysis of varwhere k 1 (m/d), k 2 (m/d), and k 3 (g/m 2 •d) are the first-order removal rate constants, C i is the inlet concentration (mg/L), C o is the outlet concentration (mg/L), C* is the background concentration (mg/L), HLR is in m/d, and C h is half the saturation constant of limiting substrate (mg/L). Based on the inlet concentrations, as suggested by Kadlec and Wallace (2009), C* for BOD 5 , COD, NH 4 -N, and TN were of 6, 10, 0, and 1.5 mg/L, respectively. In addition, the recommended values of C h for NH 4 -N, COD, and BOD 5 are 0.05, 60, and 20 mg/L, respectively (Tchobanoglous et al., 2004;Saeed and Sun, 2011;Vaccari et al., 2006).

Multiple regression
Multiple regression is a statistical model that has two or more independent (predictor) variables (Darajeh et al., 2016;Katz, 2013). This model was used because the performance of the HCW is influenced by several complicated factors such as climate conditions, hydraulic conditions, plantings, dissolved oxygen levels, microbiota, periodicity, and influent concentrations (Chan et al., 2008). In MR, these factors are expressed as "predictor variables," and the effluent concentrations or removed load is considered to be the "dependent (response) variables".
The purpose of MR for simulating the operation process of the HCW is to find the correlation between the removal of pollutants and the sys-tem's inherent controlling factors. Therefore, the MR equation is as follows (Eq. (8): iance, and a post-hoc test, Tukey's honest significance difference (Tukey HSD) test, was used to compare the means at a 95% confidence level. The R-package (relaimpo package function) was used to determine the predictor's relative importance using the LMG measure and boot-strapped confidence intervals. This method of BMA, supplemented with the use of the "BMA package", proved useful for selecting the optimal models.

System performance
The results of the VF-HF system after more than 6 months of operation showed that the removal efficiency and the effluents differed as the HLR changed (Figs. 2 and 3, and Table 2). The effluent strengths fit well with the discharge limits in the first and third stages of the mon-itoring period. The means of the effluent concentrations were 87 ± 10.6 mg/L of TSS; 31 ± 13.1 mg/L of BOD 5 ; 59 ± 17.2 mg/L of COD Cr ; 5.3 ± 3.0 mg/L of NH 4 -N; 8.4 ± 4.3 mg/L of NO 3 -N; 7.1 ± 2.7 mg/L of TN; 0.9 ± 0.3 mg/L of PO 4 -P; and 1485 ± 1184 of MPN (most probable number)/100 mL of TCol. The average treatment efficiencies for TSS, BOD 5 , TN, NH 4 -N, PO 4 -P, and TCol were 28.3 ± 12.2, 74.9 ± 11.5, 79 ± 7, 76.2 ± 12.9, 3.6 ± 43.7, and  82 ± 11.3%, respectively. The reduction in BOD 5 and COD Cr can be attributed to the biodegradable process that converted the organic matter to water and carbon dioxin in the HCW. In comparison with Vietnam's discharge limits for sewage wastewater (50 mg-BOD 5 /L), the average effluent of BOD 5 in this experiment was quite low and suitable for disposal into water bodies that were used for agricultural purposes or the equivalent. It is also clear that TN and NH 4 -N removal were fairly high owing to both the nitrification (the conversion of nitrogen to NO x in high-oxygen conditions) and denitrification (conversion of NO x into N 2 in low-oxygen conditions) that occur in the HCW (Vymazal, 2005). Furthermore, other mechanisms of nitrogen removal can include assimilation by microbial or plant biomass and adsorption onto the filter as explained by Keffala and Ghrabi (2005). The removal rate of BOD 5 was similar in both tanks (VF and HF), which was approximately 50.9%. The HF unit was considered to have less dissolved oxygen, and yet, it achieved 45 ± 13.4% reduction of NH 4 -N and 44 ± 12% re-duction of TN. There were statistically significant differences in the effluent concentrations among the three HLRs (P < 0.05) except for the case of PO 4 -P. Table 2 presents the inputs and removal efficiency of the system in terms of removal loading rate (L r ) and percentage (%). L r reached a value of 14.5 ± 2.7 g/m 2 •d of BOD 5 (30 g/m 2 •d for the VF unit and 6.9 g/m 2 •d for the HF unit). The overall TN loading rate reached 4.0 ± 1.3 g/m 2 •d, with 9.7 g/m 2 •d for the VF unit, and 1.2 g/m 2 •d for the HF unit. The removal rates observed in this experiment were rela-tively high as compared with those of the previous studies on HCW. Vymazal and Kröpfelová (2015)  Although the PO 4 -P removal rate in this study was quite low, the results show that the HCW is a promising technology solution for treating the sewage in order to meet the current Vietnam's standard discharge limits. Based on these influent loading rates and HLR of lower 0.15 m/d which tested, this HCW system was found to be suitable for treating sewage in Dong Ha city or similar types of wastewater. In addition, the results obtained for mass-balance analysis during operation of this HCW system at a HLR 0.15 m/d indicate that the average pollutant load treatments achieved in VF and HF were 14 g/d and 3.58 g/d for TN, 9.11 g/d and 2.58 g/d for NH 4 -N, 0.0 g/d and 0.21 g/d for PO 4 -P, and 39.95 g/d and 23.58 g/d for BOD 5 , respectively. However, while TSS removal was not very high, with TSS load removal in VF and HF of 3.63 g/d and 20.58 g/d, respectively. The nutrients, organic matter, and TSS reduced by HCW were attributable to vegetative nutrient uptake, bacterial metabolism, adsorption, and/or accumulation by various forms of mechanical separation, such as layers of filter media. Further information about the simple mass balance estimation in each working unit of this HCW system can be found in supplementary data. Average ± SD: Average ± standard deviation.  Fig. 4 shows the correlation coefficients (r) and plots of the useful variables. The value of r ranges from −1-1, with 0 presenting a zero correlation and 1 (or −1) indicating perfectly correlated variables. The  CSTR (M ) intersections between the rows and columns are r values (number in the square) and scatter plots of the variable pairs and histograms of the variables appear along the matrix diagonal. Hence, on investigating the relationship between the influents, operation parameters, and effluents, we can observe that the HLR had a relatively high correlation with C on (r = 0.64), C ob (r = 0.7), L rn (r = 0.46), and L rb (r = 0.56). The role of HLR in effluent concentration and removal efficiency was also con- firmed by Kadlec and Knight (1996). A low negative correlation is obtained between C ob and C ib (r = −0.37), which contrasts with the results of the previous studies. Kadlec and Knight (1996) and Babatunde et al. (2011a) stated that C ib had the highest correlation with C ob . Some variables such as C ib and C its had negative correlations with C ob and C on (r value from −0.31 to −0.46). In addition, the pH correlated weakly with all the effluents (L rb , C ob , L rn , and C on ) with a lower r value at 0.1.

Results of correlation analyses
NH 4 -N, and TN in three models (M 1 , M 2 , and M 3 ), and presented in Table 3. The average influent (C i ) and effluent (C o ) for BOD 5 , COD Cr , NH 4 -N, and TN were 127.7 mg/L and 31 mg/L, 186.8 mg/L and 59 mg/ L, 23.1 m/L and 5.3 mg/L, and 33.6 mg/L and 7.1 mg/L, respectively. Table 3 shows that M 2 does not completely fit the data of all the pollutants (P > 0.05). R 2 of M 3 for COD Cr and NH 4 -N were higher than that of the other pollutants and reached values of 0.60 and 0.68, respectively. This means that the correlation was not coincidental. In general, the models M 1 for COD Cr and M 3 for COD Cr and NH 4 -N (R 2 > 0.5), expressed the correlation required for matching the pollution concentrations and operating parameters well.
The pollutant removal rates of model M 1 for COD Cr was 0.19 m/d, which is similar to the range mentioned in several studies such as Kadlec and Wallace (2009) and higher than the result obtained by Trang et al. (2010). For an HLR of 0.086 m/d, Trang et al. (2010) presented k 1 of COD Cr as 0.07 m/d (R 2 = 0.7). For NH 4 -N, the value of k 3 obtained (2.63 g/m 2 •d) is much lower than that (29.8 g/m 2 •d) re-

Results of multiple regression for BOD 5
The output obtained from the LMG analysis (Table 5) indicates that the 64.9% and 71.7% of the total variation in C ob and L rb , respectively, can be attributed to the seven selected predictor variables. It shows that the HLR can be considered as the most important predictor that ex-plained 39.1% and 35.6% of the total variation of C ob and L rb , respec-tively. Thus, C its and C ib played a smaller role while the pH, C in , C ip , and C it had a weak impact on the variation in response variables. These results demonstrate that there are several factors that influenced the effluent BOD 5 of the HCW among which the HLR is an important parameter. Thus, it can be said that the HLR is required to be carefully selected for the design and operation of the HCW system.
Using the seven predictor and response variables of C ob and L rb from Table 5, the HCW data was calculated using the BMA method and summarized in Table 6. Based on the BMA criteria, four optimal models were selected (B -B ), which have low BIC and high R 2 results. In reported by Saeed and Sun (2011

Correlation analyses and predictive models
The values of the independent and response variables are presented in Table 4. COD Cr changed uniformly with BOD 5 , PO 4 -P was primarily influenced by adsorption and desorption, and TN was treated as a lumped category pollutant; therefore, in this study, MR was used for BOD 5 and NH 4 -N. Table 4 Values of independent and response variables.
Hijosa-Valsero value greater than 0.50 can be considered as a validation for linear models.
B 4 obtained the highest value of R 2 (0.7), whereas B 3 had the lowest BIC (−16.8) which due to fewer variables or more simply. Considering the R 2 value, these results are lower than those obtained in the study of Babatunde et al. (2011a) (R 2 = 0.86). It is likely that the HLR was the most important factor that influenced the removal efficiency of the HCW in the case of all the optimal MR models. This result is not in agreement with the study reported by Babatunde et al. (2011a), which concluded that the optimal model for C ob consisted of C ib , C ic , C in , and C ip (not including HLR). Although the R 2 values reached 0.6-0.7, which means that the predictive models explained 60-70% of the variance of L rb , these results prove that B 1 -B 4 are useful models for predicting the removal loading rate of BOD 5 of the HCW. Table 7 shows that HLR and C its had high relative importance to C on , which showed a 37% and 20% contribution to the total R 2 , respectively. In addition, using the regression model, C in and HLR explained 60. The results of the BMA method for NH 4 -N are presented in Table 8. Models B 8 -B 9 had the highest value of R (0.9) and lowest value of BIC concentration) (−40.5 to −43.3). This means that the MR models could explain 90%  Previous instances of the use of HCW for sewage wastewater Table 6 Optimal models obtained by BMA method for BOD5 effluent and removal rate. treatment in Vietnam is still limited, and therefore, there is some doubt regarding the potential of this technology. For the actual conditions of the sewer and local materials, this work demonstrated that the HCW Table 7 Relative importance for NH4-N effluent and removal rate.  tralized activated sludge wastewater treatment plants that are popular in the cities of Vietnam. Moreover, the use of the statistical method to simulate the HCW performance not only facilitates an investigation of the processes in the HCW but also clarifies the activity trends and the influent parameters. The parameters such as HLR, and TSS had a significant influence on the effluents and will thus be carefully controlled in the real operation of HCWs. Engineers and designers may be able to apply these predictive models for retrofitting or creating more effective HCWs with VF and HF in series.

Results of multiple regression for NH 4 -N
The new local plant (Colocasia esculenta) might be the first plant used in such an experiment with HCW for wastewater treatment. Although this study might contain insufficient control tests and longterm observations, the results obtained for the HF and Colocasia esculenta growth demonstrated the great potential of this plant for application to CWs for treating different types of wastewater. Therefore, testing Colocasia esculenta in various types of CW flows and regimes could create more opportunities for enhancing and promoting the B9 Lrn = 8.6 HLR + 0.14 Cin+ 0.01 Cits -0.05 Cib -2.8 0.9 -41.0 performance of HCWs in the real world. By applying the familiar first-order model for predicting the pollu-(R 2 = 0.69). All the models reveal that the variances in the effluent NH 4 -N and the removal rate were best explained by the predictor variables of HLR, C in , C its , and C ib . These predictors are quite different from those of the optimal model of Babatunde et al. (2011a) which included C its , C ic , and temperature.
The results of the above analysis indicated that the MR equations of B 5 -B 9 can be useful for simulating the variances in the NH 4 -N effluents tion removal in the HCW, this work demonstrated that the MR method is robust in the modeling of HCW performance as well as for selecting the predictive equations. This research also suggests that studies on the use of MR with more inputs and out parameters of the HCW may im-prove the understanding of HCWs. Furthermore, it would be interesting to obtain test results for this method over a longer period of operation time in order to develop a superior HCW system.

Conclusions
The results obtained for a high removal rate of organic matter and nutrients indicate that the use of HCW can serve as an efficient method of sewage treatment. With high R 2 values (0.6-0.9) derived from the MR analysis, it is important to note that random noise was not taken into consideration in the investigated models. In all the predictive models, it was observed that the HLR had a dominant impact on the response variables. In general, the MR linear analysis revealed a su-perior ability for prediction of the HCW performance than those of the kinetic rate constant. Models B 1 -B 9 can be used to evaluate the HCW system performance and predict its output.