Simulation of the Hydraulics and Treatment Performance of Horizontal Subsurface Flow Constructed Wetland Treating Greywater

Constructed wetlands (CWs) have evolved as some of reliable wastewater treatment technologies. Various types of CWs differ in their main design characteristics and in processes responsible for pollutant removal. Classification of CWS is based on the type of vegetation used and hydrological parameters involved and can thus be classified as free water surface or subsurface flow systems. Further, subsurface flow systems can be classified according to flow direction as vertical or horizontal. This study considers horizontal subsurface flow constructed wetlands (HSFCWs) which introduces the mechanistic, dynamic compartmental model-Constructed Wetlands 2D (CW2D). The model has successfully been utilized to evaluate the performance of vertical flow constructed wetlands and is being tested on HFCWs. An outdoor pilot scale HSFCW system was established in Nakuru, Kenya. CW2D was calibrated, validated and used to simulate hydraulic performance of HSFCW system. The model was used in predicting effluent concentrations of the main greywater pollutants. In general, the results obtained showed a good match with the measured data. CW2D is an effective tool for evaluating the performance of CWs and can provide insights in treatment problems at an existing CW. The same methodology can be used to optimize existing systems.


Introduction
Constructed wetlands (CWs) are engineered systems in a more controlled environment acting as filters for wastewater and have evolved as some of reliable wastewater treatment technologies. Various types of CWs differ in their main design characteristics as well as in the processes responsible for pollutant removal. The classification of constructed wetlands is based on the type of vegetation used and hydrological parameters involved. Under hydrological variables, CWs can be classified as free water surface or subsurface flow systems. Further, subsurface flow systems can be classified according to flow direction as either vertical or horizontal.
During the treatment process, interaction between water, the granular medium, macrophytes, litter, and microorganisms in a horizontal subsurface flow constructed wetland system (HSFCWs) presents a complex process. Also plants in constructed wetlands are known to be good water filters since they can tolerate high concentrations of trace elements [2]. In general, the performance of constructed wetlands is good for the removal of suspended solids and organic matter (in terms of COD and BOD 5 ) but inconsistent and often poor for nutrient (N&P) reduction [35].
As such, to study different simultaneous microbial reactions involved in the removal of various pollutant is a challenge. Modeling has thus emerged as a powerful tool for understanding the performance of wastewater treatment systems [20]. Mechanistic mathematical models can be used to determine the relationship between the different processes and weigh their relative contributions. The chosen mathematical model is a tool that helps in predicting the quality of effluent based on other process variables. Many designs models are based on rule of the thumb and regression equations. Some are based on first-order reaction rate models and their extensions and the monod-type equations [34].
However, scientists have developed several mechanistic models in order to better understand the hydraulic and reactive processes of HSFCWs [25]. One example of a dynamic mechanistic, compartment simulation model is the Constructed Wetland 2-Dimensional (CW2D) that was developed by Langergraber [18,20]. The CW2D model has successfully been utilized to evaluate the performance of vertical flow constructed wetlands and is being tested on HFCWs. In this study, the dynamic mechanistic mathematical model (CW2D) chosen was evaluated for effectiveness in predicting various pollutant concentrations of actual greywater treatment system being studied. The focus of pollutant concentrations was on some of the standard water quality variables namely: Biological Oxygen Demand (BOD 5 ), nitrogen and phosphorous with the aim of using the model for prediction of effluent quality.

Constructed Wetland Technology
Constructed wetlands are becoming increasingly vital in the treatment of domestic wastewater. They are widely used and recognized as effective and suitable technique for wastewater treatment systems [6,15,29]. However, the interaction between substrate, water, microorganisms and vegetation is treated as a 'black box' and pollutant removal process is not clearly understood. Another limitation is limited information available in the wetland's actual functioning in terms of hydraulic behavior and removal processes.
Several studies have investigated hydraulics [5,24,31,35] and kinetics [19,22] of pollutant removal processes separately in constructed wetland. However, the processes can be improved through numerical methods. As a result, an explosion of increasingly complicated numerical models for simulating water flow and contaminant transport in a subsurface constructed wetland has been seen in the last decade.
Numerical methods are a means to increase the understanding of the processes occurring in the 'black box' constructed wetland. Once reliable models for constructed wetland are available, they can be used for evaluating and improving existing design criteria [30,18]. Knowledge of hydraulic performance and flow dynamics in the removal process is critical to optimizing the design tools. The authors proposed a kinematics model (CW2D) implemented in the variably-saturated water flow and solute transport program HYDRUS-2D [17,27], to model the biochemical transformation and degradation processes taking place at the same time [20].
The HYDRUS-2D/CW2D model describes the flow in a variably-saturated porous media using the Richard's equation. This model is totally mechanistic and does not consider the wetland as a 'black box' like other models [20,22]. The model is based on HYDRUS 2D model, and incorporates a "wetland" extension term that considers biochemical transformation and degradation process. The CW2D implementation in the HYDRUS software reaction model has a matrix notation based on Activated Sludge Models (ASMs) describing carbon, nitrogen and phosphorous transformation processes [17,21]. It was mainly developed for modeling vertical flow systems. As a result, the CW2D includes only anaerobic and anoxic transformations and degradation processes [20]. These processes are described for the main constituents of wastewater i.e. organic matter, nitrogen and phosphorous [2,10,19]. However, horizontal flow systems can be simulated when only saturated water flow conditions are considered in the models.
Many models oversimplify the constructed wetland treatment systems which have extremely complicated physical, chemical and biological processes [33]. Apart from considering the influent-effluent concentrations in evaluating the treatment performance, hydrodynamic conditions such as wetland dimensions, porosity and conductivity of the medium, need to be considered. In these treatment systems, some have recorded efficiency removal of pollutants of 15-99%, with plant presence significantly enhancing the performance [35]. On the other hand, hydrodynamic conditions and hydraulic residence time influence removal efficiency of pollutants of concern and play a critical role in performance evaluation. In the modeling of treatment systems understanding of hydraulic residence time is a fundamental requirement in performance evaluation.
Several studies have attempted to explain the performance of CWs treatment systems. For example, removal efficiencies for various pollutants in HSFCW and free water surface (FWS) constructed wetland systems was studied by [29]. Two types of models namely first-order plug flow and multiple regression were used to evaluate the system performance. The authors found that the first-order plug flow model estimated slightly higher or lower values than observed when compared with other models. However, a major limitation of the first-order kinetics is that CW system is required to keep the same flow rate, pollutant concentration and ideal plug flow. As such, models that describe the distributed flow hydraulics and reactions need to be explored since the wastewaters are characterized by large variations in influent flow rate, composition and concentration. In an attempt to describe the distributed flow hydraulics, Continuous Stirred Tank Reactors (CSTRs) model were used by [12 and 33].
A series of network of CSTRs is most frequently used to describe the hydraulics of these systems and reactions are modeled with various complexities [20]. Though CSTRs are often employed to overcome the limitations imposed by the assumption of plug flow, a short coming observed in CSTRs models often fails to replicate the sharp rise and long flat tails of residence time distributions observed from tracer tests in wetlands [32]. The authors describe the use of Zone of Diminishing Mixing (ZDM) model that aims to overcome the shortcomings observed in the use of the CSTR approach. However, the approach is still closely linked to plug flow though with a dispersion component since the ZDM modeling approach assumes that water flows through the wetland along a main flow path from the inlet to the outlet. Describing CWs pollutant removal process under plug flow conditions is over simplification of the treatment process since it means that the hydraulic residence time of water flowing through the system is described by a single value. However, flow characteristics through the treatment system exhibit spatial variability and can best be described by a distribution rather than a single value.
Though several studies have investigated the kinematics of pollutant removal, most have concentrated on one line studies concerning flow hydraulics of constructed wetlands and applying different hydraulic models [15,32]. However, the kinetic model has to be assisted by a hydraulic model to achieve an optimum model-based design [26]. Even with abundance of well documented models currently available, one major problem often preventing their optimal use is the extensive work required for data preparation, numeric grid design and graphical presentation of the output results [9,30,34]. However, HYDRUS/CW2D is designed to create, manipulate and display large data files. The model also facilitates interactive data management. Interestingly, more recently, first stage anaerobic up-flow and the remaining stages tidal flow with effluent recirculation operation strategy has proved to be an effective approach to promote the capacity of the CW for high-strength wastewater treatment [35].

Experimental Set-up
The HSFCW treatment system serving a small school population with approximately 240 students was established at Crater View Secondary School, Nakuru-Kenya. The system only treats greywater stream of the school's domestic wastewater. The system includes a trough, pre-treatment chamber partitioned into two, a sense DNN light duty water meter (with protected magnetic transmission) and flexible outlet pipe. Pretreated greywater flowing horizontally through artificial bed consisting of graded sand and planted with vetiver grass was studied by laboratory analysis of samples collected from the treatment system. Steady-state flow through HSFCWs was measured at the entrance (influent) and at the discharge of the wetland (effluent). The HSFCWs had a mean depth of 0.86m, length of 2m, width of 1m, a total surface area of 2.0 m 2 and a nominal hydraulic residence time of 48 hours, with all organic load allowed to settle in the pre-treatment chamber. The width was centrally partitioned into two equal parts to give a length to width (L/W) ratio of 1 to 4. Selecting the correct L/W ratio was critical because it has a significant influence on the flow characteristics through the system [12,21,24]. This ratio was used for the rectangular wetland system at the study site to describe the hydraulic efficiency ( λ λ λ λ ), given by the relationship presented in Equation 1.

Data Collection and Pre-Processing
Data was collected over a period of 16 weeks on a biweekly basis from the outdoor pilot scale HSFCW treatment system. In situ parameters (temperature, EC, pH and DO) were also measured. Influent and effluent samples collected from the treatment system were analyzed in the laboratory within two hours of sample collection for various physicochemical parameters that included BOD 5 , NH 4 -N, NO 2 -N, NO 3-N and Phosphorous. The laboratory method of analysis was by using the procedure described in Standard Methods for the Examination of Water and Wastewater [1].
Water flow at the inlet and outlet from DNN light duty water meter (with protected magnetic transmission) were recorded every 15 minutes for 48 hours. Additionally, tracer studies were carried out during the same period. The pulse tracer, sodium chloride (NaCl) solution, was instantaneously injected at the inlet and exit concentration analyzed using electrical conductivity (EC) meter. The effluent Clconcentration was calculated from EC readings of those samples using a conductivity meter. Results from the tracer studies were used to analyze hydraulic residence time distribution of HSFCW which produces a tracer concentration versus time distribution also called the answers curve. This procedure was used to describe the hydraulic behaviors since it involved evaluating the system's response at the exit through sample collection. To further describe the hydraulic behavior, the following parameters were additionally analyzed in the laboratory in the laboratory using standard methods: porosity, bulk density and saturated hydraulic conductivity of the substrate medium.

Model Calibration
Part of the input data and procedures adopted during the pre-processing stage was based on default parameters obtained from literature described by [17,24]. The required model inputs are nominal hydraulic residence time, evaporation losses, treatment system dimensions, flow rate and concentrations of TP, BOD 5 , NH 4 -N, NO 2 -N and NO 3 -N. The model output consists of flow rate and same five concentrations parameters. Calibration process of CW2D model involved estimation of the model parameters to fit a set of experimental effluent data obtained from HSFCW system under study. The data sets were checked for any inconsistencies, split into two by random selection and from this data, half (data set one) were used for training (examining the model parameters) and the rest (data set two) for testing/ verifying the model.
The model hydraulic soil parameters were adjusted one by one until the model output corresponded with the available HFCWs effluent using data set one. This gave a parameter set that result in acceptable predictions of the effluent concentrations using data set two. To study model hydraulic performance, model output of outflow BOD 5 , nitrogen and phosphorous concentrations were compared to actual effluent values (data set two).

Sensitivity Analysis
Sensitivity analysis studies the "sensitivity" of the system outputs to changes in the parameters, inputs or initial conditions. A sensitivity analysis was conducted on the model parameter and operational conditions at the beginning of simulation to determine which sand hydraulic parameters would require most scrutiny in the future simulations. Sensitivity was calculated according to a method presented by [13,23] as the relative change in state variables (x) divided by the relative change in parameter (p) as presented in Equation 2.

/ /
Each parameter was adjusted a total of two times, parameter values were increased by + 10% and decreased by -10% to determine the relative sensitivity (S r ) of each parameter. Sensitivity analysis can determine whether there is direct (+) or inverse (-) relationship between parameters and output [4,13]. It is a main tool for identifiability analysis which attempts to provide insights on the adequate values.
The variably-saturated water flow and solute transport program HYDRUS-2D [28] and its extension, the multicomponent reactive transport module CW2D [17,20,28] were used in a pilot scale HSFCW to simulate the hydraulic flow and the removal efficiency of pollutants at the study site.

Results and Discussion
The treatment of greywater as it passes through substrate medium relies on the physical, chemical and biological processes taking place within the wetland. Inbuilt van Genuchten-Mualem soil model and other transport values were adopted from the HYDRUS-2D/CW2D catalogue for sandy soil that uses pedotransfer functions (PTF) published by [3]. These soil hydraulic parameters affect the shape and bathymetry of the wetland associated with the spatial variability of flow characteristics within a wetland system. The hydraulic residence time of water depends on the path taken by the water as it flows through the system. It further depends on the evapotranspiration of the system. However, for subsurface flow constructed wetlands, evaporation is minimal. In effect, this flow facilitated the physical, chemical and biological treatment processes that occurred within the wetland system.

Calibration Results
Water flow and soil hydraulic parameters were selected to develop HSFCW system dynamics. Measured system data was divided into two; data set one for model calibration and data set two for validation. The calibration started from adjusting various soil hydraulic parameters presented in Table 1 namely: saturated water content (Q s ), saturated hydraulic conductivity (K s ), pore connectivity parameter (α), and shape parameters (n and l). The measured van Genuchten soil hydraulic parameter values are compared to predicted values and results presented in Table 1. These results were further subjected to a statistical analysis by using the root mean square error (RMSE) method and correlation coefficient. The RMSE between measured and predicted soil hydraulic parameters (R 2 ) was 0.997 as shown graphically in Figure 1, while the correlation (R) was 0.999.

Figure 1. Correlation between measured and predicted values of soil hydraulic parameters in model calibration.
This is suggests that in the final agreement there was minimum discrepancy between predicted and measured values. Since the predicted parameters were close to measured values, then there was good agreement between measured and predicted values for inverse solution of the soil hydraulic parameters of the transport domain studied. Hydraulic efficiency of the wetland was calculated from Figure 2, based on the procedure adapted from [14,25], as the ratio of the time taken for a conservative tracer to reach a peak at the outlet to the nominal retention time given as 35/48.

Figure 2. Graph of Predicted and measured tracer studies.
The two graphs have a coinciding peak after 35 hours of the total simulation time of 48 hours and the corresponding coinciding hydraulic efficiency is 73% (35/48). Hydraulic efficiency provides a good measure of the effective volume as well as pollutant residence time distribution within the treatment system, a view supported by [16,25]. Therefore, to improve the hydraulic, volumetric efficiency, short circuiting and improving mixing conditions there is need to include plants in the system [5,6,8,11]. In this study, the hydraulic characteristics with a significant influence on the efficiency of the wetland are similar for measured and predicted conditions of this treatment system. Meaning these results fit well the simulated to measured data presented in this study. This could suggest that, the predicted tracer concentration curve adequately represents the measured values of the hydraulic efficiency of the system. These results further provide additional specific knowledge of the intrinsic processes of the CW2D system under this study.

Simulation
A simulation was run with part of influent greywater concentration from the sampling (data set two) as inputs to CW2D model until it reached a steady-state. The steady state simulation results were compared to values analyzed from effluent concentrations of selected pollutants. Results of steady-state calibration of CW2D model are given in Table 2.  The model predicted values were compared with laboratory measured and analyzed values using the root mean square error (RMSE) method. The correlation (R) and RMSE indicate how close measured values are to the predicted values. In this study, results presented in Table 3, have a close correlation value that range between 0.825 and 0.956 which means the model could be used to make accurate predictions. Details of each parameter are presented in Table 3. Simulated results from CW2D were in agreement with system data between the actual (measured) and predicted (simulated) values for effluent concentrations. The results showed an overall good match between the measured and predicted results. This could partly be attributed to the good calibration of HYDRUS 2D/CW2D model with site specific soil hydraulic properties that gave representative flow dynamics of the system. The treatment performance is also associated with good hydraulic characteristics within the system and this has a significant influence on the efficiency of the wetland treatment system. However, a very accurate calibration and validation of the model in horizontal flow constructed wetlands can only be achieved with long records of data collection from outdoor treatment systems.

Sensitivity Analysis
Five soil hydraulic parameters namely saturated water content (Q s ), saturated hydraulic conductivity (K s ), shape parameters (n and l) and pore connectivity parameter (α) considered most influential were changed to improve the agreement between predicted and measures values of pollutant concentrations. On adjusting Q s , K s and α using data set one, this procedure yielded BOD 5 , NH 4 -N, NO 2 -N and NO 3-N concentrations values of one or more orders of magnitudes lower than those measured from the system. It was observed during the sensitivity analysis that, the shape parameters "n" and "l" were insensitive. These results are in general agreement with previous findings from similar studies reported by [27,31,32]. Results for sensitivity analysis indicated that changes in "n" and "l" values approached zero, suggesting the effects of individual parameters on model prediction to be minimum.

Conclusions
For HSFCW systems, the Dynamic mechanistic compartment model, Constructed Wetland 2-Dimension (CW2D) showed a good agreement between measured and predicted results though the model was mainly developed for modeling vertical subsurface flow constructed wetland systems. Subsequently, the model and its parameters are reasonably representative of the treatment process. This is an indication that model assumptions are compatible with the system behavior and therefore reveals that the hydraulic performance of HSFCWs in treating greywater can be predicted using the CW2D model.