Numerical Models of Subsurface Flow Constructed Wetlands: Review and Future Development

: Numerical model as a simulation tool was used to describe the pollutants transformation and degradation process in constructed wetlands (CWs). It can help provide insight into the “black box” and increase the understanding of the complex processes in CWs. In the last few decades, several process-based numerical models were developed to depict the pollutants removal processes in CWs, which include biochemical model, hydraulic model, reactive-transport model, plants model, clogging model, and coupling model combining two or more sub-models. However, there was a long way to go before fully understanding the decontamination mechanisms of CWs. On the one hand, single or a composite model coupling a small number of sub-models cannot fully reveal the decontamination processes. On the other hand, a comprehensive model including all sub-models of current cognition involves numerous parameters, most of which are interaction and cannot quantitatively determined, thus making the model complex and leading to di ﬀ use interaction. Therefore, in order to describe the reaction processes in CWs more accurately, it is expected that all parameters should be quantiﬁed as far as possible in the future model. This study aims to provide a review of the numerical models of CWs and to reveal mechanism of decontamination. Based on the advantages and disadvantages of existing models, the study presented the improvement method and future research direction: (1) new detection / monitoring technique or computing method to quantitatively assess the parameters in CWs models, (2) correcting the simulation errors caused by the assumption of Activated Sludge Models (ASMs) and developing a complete bioﬁlm reaction sub-model, (3) simpliﬁcation of the comprehensive model, and (4) need of emerging pollutants modeling.


Introduction
Subsurface flow constructed wetlands (SSCWs) are artificial systems planted with aquatic macrophytes, which are currently recognized as an environmentally friendly and sustainable environmental biotechnology for various wastewater treatments due to their cost-energy efficiency, operation convenience, and good performance [1,2]. Over the last decades, constructed wetlands (CWs) have been wildly used and implemented for removing pollutants from a wild range of wastewaters [3]. The types of wastewater treated include domestic wastewater [4,5], rural wastewater [6], industrial wastewater [7], swine wastewater [8], agricultural wastewater [9], winery wastewater [10], and even landfill leachate [11]. Flores et al. [10] compared six scenarios and three common methods to purify winery wastewater and found that the CW scenarios were the most environmentally friendly alternatives. Zhao et al. [12] used dewatered alum sludge cakes as the main wetland substrate to treat In recent years, various sub-models were developed to describe the processes corresponding to different reaction mechanisms. The main equations in process-based models used for calibration and validation of process-based models are derived from processes occurring in wetlands rather than experimental data.

Hydraulic Sub-Model
Hydraulic sub-model characterizes the transportation and transformation paths of contaminants in CWs. Therefore, It is crucial to understand the hydraulic behavior in CWs [24]. Different hydraulic sub-models were used for describing the hydraulic properties in different type of CWs. Neither horizontal subsurface flow (HF) nor vertical subsurface flow (VF) CWs has free water level. HF CWs are generally simulated under saturated water flow conditions with Darcy's equation. As to VF CWs, which was highly dynamic for intermittent load, Richard's equation was used to simulate the water flow under transient variably saturated flow [19].
Samsó et al. [18] described the hydraulic properties of HF CWs by using Darcy's law. The flow through a multidimensional heterogeneous and anisotropic porous media was considered as saturated flow. The simulation results indicated that the coupling degree of simulation results and tracer curve can be greatly improved by dividing the main layer into different regions according to the hydraulic characteristics.
Richard's equation was usually used to simulate two-dimensional variably saturated water flow in SSF-CWs with intermittent loading [15,25]. By solving Richard's equation for the gas and water phases simultaneously in MATLAB, Forquet et al. [26] improved the hydrodynamic description in CWs and provided an accurate estimation of the gaseous fluxes in porous media. A large number of experiments have proved that the sub-model with Richard's equation can well simulate the hydrodynamic characteristics of VF CWs [19,27,28].

Reactive-Transport Sub-Model
The physical and chemical mechanisms of pollutants removal depend on the interaction of contaminants and substrate, including filtration, interception, absorption-desorption, convective-diffusion, complexation, redox reaction, and ion exchange, etc. All those processes did not occur one by one in CWs but in simultaneous manner and having interactive effects. Therefore, it is almost impossible to simulate the contribution of a single process to pollutant removal.
Reactive-transport sub-model generally includes adsorption and desorption processes between the solid and liquid phases, convective-dispersive transport in the liquid phase, and diffusion in the gaseous phase [19]. Researchers usually describe convective-dispersive transport processes by hydraulic equations and multicomponent chemical equilibrium equations, which were based on the convection-dispersion theory of solute transport in porous media [18,20,21]. The adsorption and desorption processes between the solid and liquid phases were often described as a source sink term in a complete CWs model [29][30][31].
Llorens et al. [21] described the dynamics of the contaminants reactive transport by Fick's law and calculated the dispersive and diffusive fluxes. However, the adsorption and desorption processes between mobile zones and immobile zones were not considered in this sub-model. However, Samsó et al. [18] presented a new reactive-transport sub-model via adding reaction, attachment, and detachment rates to Fick's law to successfully depict the reactive and transported pollutants in porous media.

Biochemical Sub-Model
It has been proved that microbes play a vital role in degradation of various contaminants in CWs [24,29,32]. The main mechanisms of removal organic compounds in CWs were degradation and assimilation by special microorganisms under the conditions of aerobic, anoxic, and anaerobic. Over the past decade, a large number of numerical models have been studied to simulate biochemical processes and pollutants degradation in CWs. Most of those biokinetic models are based on Activated Sludge Model 1 (ASM1) [1,17]. Among these models, CW2D [15] and CWM1 [17] are the most advanced and accepted biokinetic models that reflect the processes of microbial degradation of pollutants in CWs. A brief recall of CW2D and CWM1 is as follows.
CW2D was developed to simulate the biochemical transformation and degradation processes for organics, nitrogen, and phosphorus under aerobic or anoxic condition in CWs. There are 3 functional bacteria groups, 9 processes, and 12 components being contained in CW2D [15], which has been implemented in various software platform by lots of researchers [27,28].
CWM1 contains most of the biochemical processes that occurred both in VF and HF CWs, such as biochemical transformation and degradation processes for organics, nitrogen, and sulphur under aerobic, anaerobic, and anoxic conditions. CWM1 also inherits the formulation of the ASM1 series and is composed of 6 functional bacteria groups, 17 processes, and 16 components [17]. Comparing CW2D, the prediction of sulphur and the biochemical reaction process under anaerobic condition was added in CWM1 [33]. The biochemical reaction and transformation process rates in the model were modeled by Monod-type rate expressions. In addition, two kinds of switching function were added to process rates to describe catalysis and inhibition factors, respectively. Therefore, the biochemical reaction process rates of the biochemical sub-model were generally composed of the switch function and Monod equation.

Plants Sub-Model
Plant sub-model describes uptake of substances and roots oxygen release. Root growth influences the water flow pattern in the subsurface and the water balance through evapotranspiration, which should also be considered [19]. Researchers believe that plants in CWs played an important role in reducing flow velocity and in changing the flow pattern [34], which is helpful for pollutants sediment and nutrient uptake by plants.
Jin and Ji [35] used the Lake Okeechobee Environment Model (LOEM) to simulate the hydraulic and transport-reaction behavior in CWs. The model offered detailed information regarding flow velocity, circulation patterns, sediment resuspension and deposition, and plants distribution by adding the vegetation stress terms in Navier-Stokes equations [36]. After calibration and validation with 6 years of long-term measured data, the results indicated that the modeled flow velocity in CWs without plants was 10 cm/s larger than that with plants and that the CWs without plants showed strong circulation patterns.
Plants provide oxygen to the rhizosphere through passive or active oxygen transport, forming a microcosmic aerobic zone, which was advanced for degradation of pollutants [37]. However, many researchers considered that the oxygen releasing from plants root was too little and could be ignored in a lot of completed CWs models [29]. Therefore, there is little literature on roots oxygen release [18,38]. Dong et al. [39] calculated the amount and diurnal fluctuations of oxygen released by roots based on mass balance method. The results showed the oxygen release rate from root being ranged from 20.3 to 58.3 g O 2 /m 2 ·d and was greatly influenced by light intensity. However, the respiratory oxygen consumption of the roots almost equaled the oxygen release, leaving little for organics degradation and nitrogen-nitrification.
The ability of plants to uptake contaminants is related to the plants characteristics, the properties of the pollutants, as well as evapotranspiration [30,40]. Langergraber [41] explored the plant uptake models from HYDRUS-2D. In the uptake models, nutrient (nitrogen and phosphorus) uptake was coupled with water uptake and indicated that it was feasible to simulate the nutrient uptake by plant in high load. However, the simulation result overestimated the situation for low-strength wastewater. Tang et al. [25] reported that the ability of plants for taking up pesticides is relatively weak. It is well-known that evapotranspiration has a positive effect on the migration of pollutants in plants [40,42]. However, the detailed dynamics of evapotranspiration and pollutant uptake are still in their infancy.
Overall, relatively, plants have little influence on the simulation results of the comprehensive model. Due to the complexity of pollutant removal processes and kinetic mechanisms by plants, the plant sub-model was generally ignored [16,29] or described as a source and sink term in published comprehensive models [18,27].

Clogging Sub-Model
The substrate clogging in subsurface flow constructed wetlands (SSF CWs) is a most serious problem in CW practice and is caused by several physical, chemical, and biological factors [43]. Clogging has varied effects on the hydrodynamic characteristics of CWs by cutting down the porosity and permeability of porous media. Those influences brought some serious operation problems associated with CWs, such as short-circuiting and reduction of hydraulic conductivity and hydraulic retention time, and finally affected the pollutants removal function of SSF CWs [44]. Clogging model was constructed to evaluate the life and long-term performance of SSF CWs by describing the processes underlying the clogging. It should be noted that the clogging model is not strictly a model of pollutant removal mechanisms but that it influences the predictions of other sub-models. On one hand, the clogging model is a response to the other sub-models; on the other hand, the clogging process affected the parameters of the other four sub-models. It is generally believed that accumulation of inert particles is the main reason for clogging, although the growth of microorganism and plant root also influenced the clogging process [45].
Rajabzadeh et al. [43] used the Carman Kozenys (C-K) equation to simulated the spatially varying hydraulic conductivity during clogging process. Since biological and chemical effects are not considered, the simulate value of hydraulic conductivity was less compared with the measured data. Giraldi et al. [46] simulated the clogging process by computing the cumulative of biofilm and six particulate components. The various porosity and permeability coefficients of substrate were also described by C-K equation. However, the biochemical reaction of particulate components was not considered in the mass conservation equation, and only hydraulic module was calibrated and validated while the clogging sub-model was not considered.
Zhao et al. [47] established a HSSF CW clogging model by coupling one-dimensional mass transfer with multicomponent biochemical reactions. The new model used total particulate concentration within media space as an indicator of clogging to dodge the complex relations among porosity, hydraulic conductivity, and retention capability of particles. The degree of clogging was estimated by comparing total particulate concentration with the greatest destiny of solid accumulation. By assuming the particulate release ratio (R j ) was constant, the clogging module simulated the particulate migrated to the next compute unit with fluids through overflow while the clogging occurred in the former compute unit that reflected the real performance of the HSSF CWs. However, the sub-model did not consider the variation of fluid flow and the influence on parameters related to biochemical degradation caused by overflow.

Comprehensive Model
A number of sub-models were required for a complete process-based model of CWs. In order to truly reflect the degradation process of pollutants in wetlands and to obtain accurate simulation results, the biokinetic models have been coupled to hydrodynamic models or other sub-models and carried out in software to solve the differential equations for dynamic simulations [22]. With the help of the development of computer technology, most of the comprehensive models are implemented in variety software platforms by coupling two or more sub-models, such as RTD/GPS-X [48], Diph_M [49], FITOVERT [46], HYDRUS-CW2D [27], HYDRUS-CWM1 [38], CWM1-RETRASO [21], CFD Model [43], and BIO_PORE [18]. Several representative coupling models published recently and their implementation are summarized in Table 1.  (9) aerobic, anaerobic, and anoxic processes (17) aerobic, anaerobic, and anoxic processes (19) aerobic, anaerobic, and anoxic processes (19) aerobic, anaerobic, and anoxic processes (8) aerobic, anaerobic, and anoxic processes (8)

HYDRUS Wetlands Module
HYDRUS is a software package for simulating water, heat, and solute movement in two-and three-dimensional variably saturated media, which simulates variably-saturated water flow, transport of contaminants, and influence of plants by solving systems of partial differential equations in three dimensions. Both CW2D and CWM1 are implemented in the HYDRUS software and obtained better simulation results [27]. Due to the fact that the microbial anaerobic degradation process was included in HYDRUS-CWM1, the simulation results were better than HYDRUS-CW2D. Except for the biokinetic one, the HYDRUS wetlands module included reactive-transport sub-model and plants sub-model. However, the clogging model was not involved in either model. The flow in the HYDRUS wetlands module was described by Richard's equation, so it can be used to simulate both VF and HF CWs. The advantage of HYDRUS is that the model can simulate the fixed biofilm, which is very important for the actual prediction of the treatment efficiency of CWs.

BIO_PROE
BIO_PORE was also a comprehensive model which was implemented in COMAOL Multiphysics, a finite elements simulation platform. The objective of BIO_PORE model is to simulate the processes of pollutants removal in CWs through detailed description of the hydraulics and hydrodynamics, reactive-transport of solute and biochemical reaction kinetic, as well as the influence of plants. Different from original CWM1, the biokinetic sub-model in BIO_PORE model was modified [18]. Two growth-limiting functions were added to the original formulation of process rates in CWM1 to prevent the unrealistic growth of bacteria within high substrate concentrations, which was one of the advantages of the model.

CWM1-RETRASO Model
The RetrasoCodeBright (RCB), a potent modelling tool, has been used to study hydrogeologic and simulate reactive transport of dissolved inorganic and gaseous species in non-isothermal saturated or unsaturated problems by finite elements [50,51]. Ojeda et al. [16] modified the original RCB code by coupling four microbial reactions: aerobic respiration, denitrification, sulphate reduction, and methanogenesis. The new model simulated the contribution of different microbial reactions to organic matter removal. The simulation results indicated that anaerobic reactions contributed to the highest COD removal rate (60%-70%). The study also indicated that acetotrophic methanogenic bacteria and acetotrophic sulphate reducing bacteria were more widespread in HSSF CWs, which was consistent with previous studies [22,52,53].
The CWM1 was firstly implemented in RCB by adding a reaction term to RCB code, called CWM1-RETRASO, to simulate hydraulics; reactive-transport; and the main biochemical reaction for organics, nitrogen, and sulfur biodegradation; and transformation in HF CWs [21]. The new mass balance equation of this model consisted of three parts: the advective flux, diffusive-dispersive fluxes, and biochemical reaction term. The first two items were calculated using Darcy's and Fick's laws, meanwhile the reaction rate was computed by the kinetic rate law formulation describing in CWM1.

FITOVERT
Different numerical models have been developed to describe VF CWs for its high dynamic [15]. However, even highly complex dynamic models are difficult to accurately describe the behavior of wetland systems, especially to provide reliable parameter estimates. Giraldi et al. [46] developed a dynamic multicomponent reactive transport model for variably saturated conditions with volume filtration of particulates on MATLAB, referred to as FITOVERT. The model can maintain the complexity of the model to an acceptable level so as to guide design and operation optimization. One advantage of FITOVERT over other published models is that it can simulate the porosity reduction and the variation of hydraulic conditions caused by microbial growth and particle composition accumulation. Since FITOVERT can accurately describe the saturated and unsaturated flows in the system, it can simulate the high dynamic VF CWS well. The mixed form of the Richard's equation and the transport equation for the dissolved components and the mass conservation equations for the particulate components were discretized in space using the Galerkin finite element (FE) method and the Petrov-Galerkin finite element scheme in FITOVERT.

CFD Model
Different from the activated sludge method, microorganisms in the CW system attach to and grow on porous media to form the biofilm, while the loss of microorganisms is corresponding to the loss of biofilms instead of discharging residual sludge. The growth and loss of biofilms affect the hydraulic properties of porous media, which are often ignored in previous models, causing overestimation of clogging time. Rajabzadeh et al. [43] used COMSOL MultiphysicsTM to couple flow dynamics, biokinetics, and biofilm development/detachment into computational fluid dynamics (CFD) model, which is known as the most comprehensive and encompassing model developed in the field.
In the CFD model [43], instead of Fick's law, Brinkman's equation with viscous transport and inertia terms was used to describe the solutes transport and flow behavior in porous media. The model simulates the temporal and spatial distribution of biofilm development in a vertical flow wetland system. Simulation results showed that, after 365 days of operation, the average porosity of wetlands was reduced by 6.6% due to the biofilm development. A dead zone was formed near the area opposite the outflow of the wetland. Due to the small fluid shear stress, the thickness of the biofilm in dead zone was 330 times that in other zones in the CW. The simulation results revealed the realistic bio-clogging processes in CWs system.
For simplicity, Rajabzadeh et al. [43] did not consider the biochemical transformation and removal of dissolved N and P as well as the plant sub-model, which should be improved in future studies. Fortunately, due to the strong inclusiveness of CFD model, these missing sub-models and dynamic models can be easily added to the comprehensive model.

Balance of Prediction Results and Model Structure
Hydrodynamics is the basis of studying solute migration and degradation in porous media. The earliest numerical models of CWs were implemented in hydraulic simulation software, for example, HYDRUS [15,54]. By summarizing the previous papers on process model of CWs, it is obvious to realize that the experimental data would match the simulation results when the hydraulic behavior of CWs was described well [19,32]. This has also been notified by Giraldi et al. [46] that proper calibration of hydraulic module was an essential condition for the stable operation of the whole model.
It is well known that Richard's equation and Darcy's laws are generally used to describe the hydraulic behavior of wetlands [16,27,54,55]. However, Richard's equation is the basic theoretical equation for vertical unsaturated flow, while the equation is not applicable to macropore flows. As for the substrate with high porosity, the dual-porosity model would get better prediction results [56]. Advective flux and dispersive and diffusive fluxes in those models were usually computed in a quasi-steady-state, which made most process-based models only work with constant influent. Llorens et al. [21] considered the matrix as a uniform medium; their simulation results were somewhat different from the measured data, which could be seen from the simulation results of different mesh numbers. Therefore, increasing the number of computational grids can reduce the impact of matrix heterogeneity on model output, but it will also greatly increase the computational load of the model.
As the same as the method of activated sludge wastewater treatment, microorganism played an important role in biochemical transformation and degradation processes for organic matter, nitrogen, and sulphur in subsurface CWs [57]. Lagngergraber et al. [54] calculated the biomass COD of dry weight (DW) soil of each layer of CWs by measuring the content of C and N at different depths of a CW system. The result showed that most of the microbial would be found in the top 10 cm of the main layer while, from the first centimeter to 10 cm, the mean values of biomass ranged from 5100 to 640 mg COD/g DW. The study on the abundance and spatial distribution of bacteria indicated that, under the stable condition of CW operation, the anaerobic bacteria dominated the CWs while the sulphate reducing bacteria, the dominant bacteria, accounted for 46%. The strong relationship between the temporal and spatial distribution of bacteria and the removal rates of pollutants has been observed [58]. According to the previous study [29], anaerobic processes contributed to almost 72%-79% removal rate of COD in HF CWs, which indicated that anaerobic bacteria, such as methanogens and sulfate reducers, play a main role in pollutants removal. All of those studies proved that biochemical sub-model offers the highest contribution to pollutants removal.
Compared with biodegradation, the absorption of pollutants by plant roots was limited. Actually, under high load conditions, the absorption of nutrients by plants can be negligible compared with the microbial degradation [59]. This was consistent with the research results of Samsó et al. [18] and Giraldi et al. [46]. In addition, the release of oxygen from plant roots is considered negligible by many researchers [60].
Due to the complex clogging mechanism of CWs, the clogging model development was lagging behind. Moreover, the clogging models developed in recent years were independent models [31,48,61,62]. Only a few published comprehensive wetland models contain clogging models, which were only coupled in the form of modules but not closely associated with other sub-models [46]. Therefore, the addition of these clogging sub-models does not significantly improve the prediction of the comprehensive model.
With the deeper understanding of the mechanism of pollutant removal, it has been well accepted that a single sub-model cannot meet the accuracy requirement of model prediction. Thus, more single sub-models were embedded to form comprehensive models [18,32,35,63]. However, with the increase of the number of sub-models, the determination of parameters and the solution of governing equations become more complicated. Therefore, the balance of model structure and prediction results is a crucial issue.
It is fair to say that balancing the number of sub-models with the accuracy of simulation is a problem of contradiction and shield. On one hand, the model should contain as many reaction processes as possible; on the other hand, the model parameters should be measurable and the governing equation should be solvable. However, the computation burden of the model increases with the number of reaction processes, while the closer the process rate expression to the actual, the more calculation complex it is. Such models are usually difficult to solve or time-consuming to compute. Therefore, high-quality models which show a reasonable balance should only include the most essential process of the system, should choose the appropriate rate expression, should simplify the calculation as much as possible, but do not cause the prediction result deviation because of simplification.
According to the purpose of modeling, the model structure can be appropriately simplified and optimized. The model of BIO_PORE [18] was developed for scientific purposes, so the more processes considered the better. Different governing equations can be introduced for different reaction processes. For example, the process of P removal can be improved by introducing the matrix adsorption analytical equation and the phosphate chemical reaction equation.
For the purpose of guiding the design, models like RSF_Sim and RTD/GPSX were developed for easy use tools with good balances between simplification and simulation result quality. All kinds of those models were based on user-friendly software and usually easy to use once it has been accurately calibrated.

The Sensitive Factors
Llorens et al. [29] confirmed that the longitudinal dispersion coefficient was the sensitive parameter in hydraulic model through repeated experimental calibration. The study by Langergraber et al. [19] also indicated that the longitudinal dispersion coefficient was the most sensitive parameter and has the highest influence on the simulation results. However, both the longitudinal diffusion coefficient and diffusion coefficient are essentially determined by wetland type and substrate. At the same time, physical and chemical properties of the substrate (particle size, specific surface area, porosity, permeability coefficient, pH, etc.) are the basic parameters of hydraulic and biochemical models. Adsorption of contaminants is another function of the substrate. Previous studies [64,65] have shown that matrix adsorption and precipitation were the main mechanisms of P removal in CWs. However, the actual P removal processes of CWs are very complicated as well.
The complexity of different models was determined by the number of component assumptions [17,18,63,66]. According to the different substrates consumed by microorganisms, the growth and hydrolysis of heterotrophic bacteria and aerobic bacteria in the original CWM1 model were divided into two parts [21]. However, the model did not achieve the expected effect, so the refinement of components could not significantly improve the simulation results. The reason may lie in that these modifications were other ways of representing those processes without changing their meaning. On the contrary, the modification of stoichiometric and kinetic parameters significantly affected the prediction results [18,21]. It is recognized that the growth of microorganisms and metabolic rates are closely related to temperature [2]. Not only biological processes but also temperature sensitive parameters are also included in several other sub-models. According the study of Samsó et al. [18], water temperatures had a great impact on the model output. Obviously, temperature is one of the important parameters in the wetland comprehensive model.
Sensitive parameters or factors usually have a greater influence on the output of the numerical model. Therefore, in future research, the sensitivity of parameters in each sub-model or comprehensive model is an important direction, which is of great significance to the optimizing of numerical models.

Species and Distribution of Bacteria in CWs
In the past decade, many process-based models have been developed. However, it is still a long way from fully revealing the mechanism of pollutants degradation. The development of the CW model was mostly based on ASMs [1,17]. However, many developers ignored the premise of ASMs based on completely mixed. In contrast with ASM, the microorganisms in CW are fixed and the species and quantity of microorganisms in the reactor were not evenly distributed [29,57]. Esther Ojeda et al. [16] used a two-dimensional mechanism model to study the contribution of each microbial population to pollutant removal in CW systems; the results showed that anaerobic bacteria, especially methanogens and sulfate-reducing bacteria, contributed the most to the degradation of the pollutants, which was also confirmed by Roger Samso et al. [57]. However, due to the lack of measured data of microorganism, it is difficult to calibrate the predicted results. Therefore, it is urgent to develop a new measurement technology to quantify the microorganism in CWs to calibrate this part of data.

The Modification of Kinetic Mechanism of Biofilm
The description of substrate consumption rate in the biochemical sub-model is also worth reconsidering. Different from the fully mixed activated sludge process, microorganisms in CWs are attached to the substrate surface in the form of biofilm, which makes microbes at different locations capture substrates and oxygen with different resistance. The microorganisms at the inner layer of the biofilm cannot directly contact nutrients and oxygen in sewage as the microorganisms at the surface, so that a concentration gradient is formed from the inside out. It is obvious that a single constant mass transfer coefficient is insufficient to properly describe the process. Therefore, in future CW model studies, the influence of biofilm diffusion resistance on biochemical reaction rate should be appropriately modified, which was also mentioned by Galvao et al. [67,68].
In addition, compared the development course of ASMs model, the theory of "Death Regeneration" and "Endogenous Respiration" were used in ASM1 and ASM3 to describe the processes of microbial decay, respectively. Most of the current published CWs models still adopted the hypothesis of the theory of death regeneration, which may have little impact on the simulation results. However, compared with the theory of "Death Regeneration", perhaps the "Endogenous Respiration" theory can better explain the fact that the effluent from wetlands still contains nitrate during a long period of time (48 h+) after the aeration system stops, which was observed in the study of Murphy et al. [68]. A study on the storage mechanism of biofilm in CWs indicated that the growth and metabolic characteristics of heterotrophic bacteria were not adequately described in the current biodynamic sub-model [67]. Therefore, it is the time to modify the biochemical sub-model, while the kinetic model of biofilm with different mass transfer coefficients based on "Endogenous Respiration" theory may be able to more realistically reflect the biochemical process in CWs.

Quantitative Parameters and Computer Technology Development
Quantitative techniques of pollutant composition and kinetic parameters are also important factors affecting the prediction results. The accuracy and practicability of the model will be greatly improved by quantifying model empirical parameters [69]. As the platform support of numerical model development and operation, the development of high-performance computers and cloud computing as well as related programming software will promote the development of numerical modeling of CWs and save the computing time of the models.

Simplification
The purpose of simplifying the model is to guide the design. Some parameters that are difficult to be accurately measured can be defined as constants or many related parameters can be modularized within the range of prediction accuracy. For example, in a CWs clogging model [47], by introducing a new index of total particle concentration, the complex relations among porosity, hydraulic conductivity, and retention capability of particles were dodged dexterously. Simultaneously, the porosity was set as a constant in the simulation process, which simplifies the simulation process and saves calculation time.

Emerging Pollutants Modeling in CWs
One issue worth considering in the next generation of water treatment models is emerging organic contaminants (EOCs), such as antibiotics, antibiotics resistance genes, pesticide, and perand polyfluoroalkyl substances (PFASs), etc. [70]. Due to their difficult elimination, EOCs have been highlighted as a worldwide environmental issue, and some of EOCs have been listed in EU legislation [71]. In the last few decades, CWs was also used to treat wastewater with EOCs. According to the study of Christofilopoulos et al. [72], the bisphenol A and ciprofloxacin were significantly removed by the CW, with an overall removal of 76.2% and 93.9%, respectively. However, due to the different properties of EOCs, the removal mechanisms are different, so the study on the removal mechanism of EOCs has become a hot topic in this field [6]. There is no doubt that numerical simulation will provide a reliable technical means for relevant study. Most recently, Ji et al. [70] have reviewed the PFAS removal and have proposed the strategy that CWs may have the potential to curb these kinds of newly emerging pollutants. Accordingly, the development of CW modeling should consider these kinds of compounds in the future development. Nevertheless, traditional numerical models of wastewater treatment do not include EOCs components, so the development of numerical models that include EOCs is of great significance to eliminate EOCs by CWs.
On one hand, the development of quantitative techniques for matrix properties and parameters are the basis for improving model prediction accuracy and for understanding of pollutant removal mechanism. On the other hand, computer technology, software development, and model simplification will guarantee the practical application of the model. Furthermore, the modification of biofilm dynamics model by the theory of "Endogenous Respiration" and the new components introduction will contribute to the innovation of water processing model. Overall, the challenges and directions on future CW model development are shown in Figure 1.

Conclusions
As an effective research tool, numerical simulation can greatly reduce the workload and time of experimental research. In the last few decades, several process-based numerical models were developed to depict the pollutants removal processes in CWs. Both the single process sub-model and the integrated model have achieved good simulation and prediction results. However, there is a long way to go before the models can guide design in CW practice.
As the basis of wetland reaction system, a good match to measured data can be attained when the hydraulic performance of the system is properly depicted. Microorganisms are major contributors to the removal of pollutants in CW systems, so that accurate description of biochemical reaction process is the key point to develop CW models. Different parameters have different sensitivity to the models based on different processes. However, wetland type (hydraulics properties), matrix property, and temperature are the basic factors that affect all sub-models.
From the current review, in future research, the CW model development should include 1) new detection/monitoring technique or computing method to quantitatively determine the parameters in CWs models, 2) correcting the simulation errors caused by the assumption of ASMs and developing a complete biofilm reaction sub-model, 3) simplification of the comprehensive model, and 4) newly emerging pollutants modeling in CWs.

Conclusions
As an effective research tool, numerical simulation can greatly reduce the workload and time of experimental research. In the last few decades, several process-based numerical models were developed to depict the pollutants removal processes in CWs. Both the single process sub-model and the integrated model have achieved good simulation and prediction results. However, there is a long way to go before the models can guide design in CW practice.
As the basis of wetland reaction system, a good match to measured data can be attained when the hydraulic performance of the system is properly depicted. Microorganisms are major contributors to the removal of pollutants in CW systems, so that accurate description of biochemical reaction process is the key point to develop CW models. Different parameters have different sensitivity to the models based on different processes. However, wetland type (hydraulics properties), matrix property, and temperature are the basic factors that affect all sub-models.
From the current review, in future research, the CW model development should include (1) new detection/monitoring technique or computing method to quantitatively determine the parameters in CWs models, (2) correcting the simulation errors caused by the assumption of ASMs and developing a complete biofilm reaction sub-model, (3) simplification of the comprehensive model, and (4) newly emerging pollutants modeling in CWs.