Mathematical Modeling of a Domestic Wastewater Treatment System Combining a Septic Tank, an Up Flow Anaerobic Filter, and a Constructed Wetland

Systems combining anaerobic bioreactors with constructed wetlands (CW) have proven to be adequate and efficient for wastewater treatment. Detailed knowledge of removal dynamics of contaminants can ensure positive results for engineering and design. Mathematical modeling is a useful approach to studying the dynamics of contaminant removal in wastewater. In this study, water quality monitoring was performed in a system composed of a septic tank (ST), an up flow anaerobic filter (UAF), and a horizontal flow constructed wetland (HFCW). Biological oxygen demand (BOD5), chemical oxygen demand (COD), total Kjeldahl nitrogen (TKN), NH3, organic nitrogen (ON), total suspended solids (TSS), NO2−, and NO3− were measured biweekly during a 3-month period. First-order kinetics, multiple linear regression, and mass balance models were applied for data adjustment. First-order models were useful to predict the outlet concentration of pollutants (R2 > 0.87). Relevant multiple linear regression models were found, which could be applied to facilitate the system’s monitoring and provide valuable information to control and improve biological and physical processes necessary for wastewater treatment. Finally, the values of important parameters (μmax, Ks, and Yx/s) in mass-balance models were determined with the aid of a differential neural network (DNN) and an optimization algorithm. The estimated parameters indicated the high robustness of the treatment system since performance stability was found despite variations in wastewater composition.


Introduction
The development of adequate technologies for wastewater treatment in developing countries is urgent to minimize environmental degradation and satisfy water supply necessities [1]. The Sustainable Development Goal number six of the United Nations is dedicated to clean water and sanitation. It seeks to substantially improve water quality by increasing the safe reuse and recycling of water systems and lowering the release of pollutants and toxins into wastewater by the year 2030 [2]. Wastewater treatment facilities may be part of centralized or decentralized clustered systems. The management of centralized systems is usually the responsibility of public agencies, as they treat vast amounts of The objective of this work was to characterize the pollutant removal processes occurring in a system consisting of an ST, a UAF, and an HFCW with the aid of mathematical models-first-order kinetics, multiple linear regression, and mass balance models. Furthermore, the use of a DNN and an optimization algorithm were proposed to determine the parameters of the Monod´s kinetic model. These models were developed for a system that has proven its effectiveness treating wastewater with a high organic and nutrient load, requires a smaller surface area per inhabitant, and produces effluents in compliance with environmental regulations [3]. To the best of our knowledge, there are few detailed studies regarding passive treatment systems in tropical or subtropical countries where climatic conditions, especially temperature, enhance the performance and stability of both anaerobic units and wetlands. In addition, several modeling studies have focused on a sole treatment unit. In this study, the complete system combining an ST, a UAF, and an HFCW was assessed with the aid of the mathematical models demonstrating the contribution of each unit to global performance. These models provided a description of the dynamics of pollutants removal within the treatment system.

Site Description
The study site was an on-site pilot wastewater treatment plant consisting of an ST, a UAF, and an HFCW planted with Agapanthus africanus (Figure 1). The system was designed to treat heavily polluted domestic wastewater produced in a research and development (R&D) center, looking to reach zero energy consumption and to eliminate the use of chemical additives. For 5 years, this plant has received the wastewater generated by approximately 150 people with a daily water consumption estimated at 50 L/person/day [16].
Water 2020, 12, x FOR PEER REVIEW 3 of 20 The objective of this work was to characterize the pollutant removal processes occurring in a system consisting of an ST, a UAF, and an HFCW with the aid of mathematical models-first-order kinetics, multiple linear regression, and mass balance models. Furthermore, the use of a DNN and an optimization algorithm were proposed to determine the parameters of the Monod´s kinetic model. These models were developed for a system that has proven its effectiveness treating wastewater with a high organic and nutrient load, requires a smaller surface area per inhabitant, and produces effluents in compliance with environmental regulations [3]. To the best of our knowledge, there are few detailed studies regarding passive treatment systems in tropical or subtropical countries where climatic conditions, especially temperature, enhance the performance and stability of both anaerobic units and wetlands. In addition, several modeling studies have focused on a sole treatment unit. In this study, the complete system combining an ST, a UAF, and an HFCW was assessed with the aid of the mathematical models demonstrating the contribution of each unit to global performance. These models provided a description of the dynamics of pollutants removal within the treatment system.

Site Description
The study site was an on-site pilot wastewater treatment plant consisting of an ST, a UAF, and an HFCW planted with Agapanthus africanus (Figure 1). The system was designed to treat heavily polluted domestic wastewater produced in a research and development (R&D) center, looking to reach zero energy consumption and to eliminate the use of chemical additives. For 5 years, this plant has received the wastewater generated by approximately 150 people with a daily water consumption estimated at 50 L/person/day [16]. Wastewater is received in the pump sump (PS) at a rate of approximately 7.5 m 3 /day and is then pumped to the ST (Figure 2). This stage requires intermittent energy, which is supplied by solar panels, to pump the wastewater to the treatment system. Energy is also required periodically for sludge removal from the anaerobic stages (ST and UAF) every six to twelve months [16,27]. This procedure regularly requires the use of external suction systems to remove the excess sludge from the ST chambers and from the bottom of the UAF. Additional energy may be required if motor devices are used to keep the HFCW clean and to control the growth of plants. Wastewater is received in the pump sump (PS) at a rate of approximately 7.5 m 3 /day and is then pumped to the ST (Figure 2). This stage requires intermittent energy, which is supplied by solar panels, to pump the wastewater to the treatment system. Energy is also required periodically for sludge removal from the anaerobic stages (ST and UAF) every six to twelve months [16,27]. This procedure regularly requires the use of external suction systems to remove the excess sludge from the ST chambers and from the bottom of the UAF. Additional energy may be required if motor devices are used to keep the HFCW clean and to control the growth of plants.
The ST is divided into two communicating chambers with respective volumes of 11.3 m 3 and 7.1 m 3 . Inside the septic tank, the water flow is slow enough to allow the sedimentation of solid particles; in addition, fats are retained in scum formed on the surface of water inside the tank. The outflow of the ST is received in the UAF. The UAF has a volume of 69.5 m 3 , and it is filled with porous volcanic rock (tezontle), which serves as an efficient support for the growth of methanogenic bacteria, creating a fixed-bed biofilm system [28,29]. pumped to the ST (Figure 2). This stage requires intermittent energy, which is supplied by solar panels, to pump the wastewater to the treatment system. Energy is also required periodically for sludge removal from the anaerobic stages (ST and UAF) every six to twelve months [16,27]. This procedure regularly requires the use of external suction systems to remove the excess sludge from the ST chambers and from the bottom of the UAF. Additional energy may be required if motor devices are used to keep the HFCW clean and to control the growth of plants. The UAF is composed of two chambers. In the first one, the water flow is received and directed toward the bottom of the second chamber, in which the water flows upward. Direct contact with bacteria attached to tezontle allows for the degradation of organic matter present in wastewater through anaerobic biochemical reactions. The final stage of an HFCW is a pond with 336 m 2 of surface and a depth of 0.7 m filled with porous rock (tezontle), which serves as a support for the roots of macrophyte plants (Agapanthus africanus) planted with a density of about three plants per square meter. The water flows below the surface at a 0.6 m hydraulic depth; this system promotes the removal of organic matter and nutrients (nitrogen and phosphorus) by the interaction between the growth media tezontle, plant roots, and microorganisms [16].
The mean hydraulic retention times (HRT) for this system were previously reported as 2.45, 6.4, and 11.75 days for the ST, the UAF, and the HFCW, respectively, adding up to 20.6 days for the total mean HRT of a sample [16].

Water Quality Sampling
Water quality monitoring and sampling was carried out fortnightly for 3 months (January 16th and 31st, February 12th and 26th, March 11th and 25th, 2020) at four sampling points within the treatment system: the sump pump (SP1), the UAF inlet (SP2), the UAF outlet (SP3), and the HFCW outlet (SP4).
Grab sampling was used because steady-state conditions were reported previously for the removal of pollutants within the wastewater treatment plant [16]. Most activities at the R&D center happen between 9:00 a.m. to 5:00 p.m. on Monday through Friday; consequently, the highest flow is received during these hours, and it is highly reduced during the weekends and holidays. This is the reason why the samples were collected during normal working hours.
Temperature, pH, conductivity (CT), and dissolved oxygen (DO) were measured in situ using a multi-parameter probe (HANNA, HI 9828). In addition, water samples were analyzed by the Analytical and Metrologic Services Unit (USAM) of the Centro de Investigación y Asistencia Tecnológica del Estado de Jalisco (CIATEJ) according to Federation, W. E., & American Public Health Association [30] to determine biological oxygen demand (BOD 5 ), chemical oxygen demand (COD), total Kjeldahl nitrogen (TKN), ammoniacal nitrogen (NH 3 ), organic nitrogen (ON), and total suspended solids (TSS). Water samples were received by the laboratory within the first four hours after sampling. Nitrates (NO 3 − ) and nitrites (NO 2 − ) were measured within the first 3 h after sampling using multiparametric kits (HACH, TNT 835 and TNT 839) by spectrophotometry (HACH, DR 5000) following the manufacturer's protocol. All concentrations were measured in triplicate and means, and standard deviations are presented. To compare the effect of different operational conditions on the performance of the wetland, the removal efficiency and mass removal rate were calculated. Removal efficiencies were calculated considering the concentration data of each pollutant at the inlet and outlet of each treatment stage as follows: ST (SP1-SP2), UAF (SP2-SP3), HFCW (SP3-SP4), and for the overall system (SP1-SP4).

First-Order Kinetic Model
The first-order kinetic model is a nonlinear relationship that has been used to predict pollutant removal in wastewater treatment systems [31]. The first-order kinetics for the removal of contaminants is described as: The monitoring data were adjusted to generate equations of the following form: where C out is the concentration of the pollutant at the system outlet in mg/L, C in is the concentration of the pollutant at the system inlet in mg/L, k is the removal constant in days −1 , and HRT is the duration of hydraulic retention of the system expressed in days, which were previously reported for this treatment system by [16]. Data adjustment was made for the entire treatment system (black box approach) and for each step (ST, UAF, and HFCW) by nonlinear least squares. The statistical significance of each parameter modeled was measured by p-value calculation, and R 2 was estimated to measure the goodness of fit. First-order kinetic models were used because they provide a parameter (k) related to the rate of removal of contaminants that can be compared with the removal efficiency and the mass reduction.

Multiple Liner Regression Model
Multiple linear regression is a statistical model with two or more independent variables (predictors) [32]. In this model, the different water quality parameters measured were used as "predictors", and a specific water quality parameter outlet concentration was considered as the dependent variable (response). The linear regression equation is as follows: where Y is the response variable, and X is the predictor variable. The parameters considered for model development were DO, CT, pH, temperature, BOD 5 , COD, TSS, TKN, NH 3 , ON, and NO 2 − .
Models were considered to be relevant when R 2 was found to be greater than 0.99 and p-values of predictor variables were found to be lower than 0.1. Additionally, obvious relations such as COD-BOD 5 , and TKN = NH 3 + ON were not reported.

Mass Balance Model
Mass balance equations were used to model the removal of BOD 5 , COD, TKN, NH 3 , and ON in each treatment stage (ST, UAF, and HFCW), because these parameters are degraded primarily by microbial action [26]. Monod's equation Equation (4) was used to express the specific growth rate of microorganisms as a function of substrate concentration S (BOD 5 , COD, TKN, NH 3 , and ON in this study). Similar approaches have been reported in the literature modeling microbial growth as a function of COD and NH 3 (as substrates) in an anaerobic reactor and a CW, respectively [20,21,25].
where S in is the substrate concentration at the inlet (mg/L), µ max is the maximum growth rate (h −1 ), and K s is substrate affinity constant (mg/L). The mass balance for every stage and for every contaminant was written as: Water 2020, 12, x FOR PEER REVIEW 6 of 20 The accumulation and transformation terms are dependent on the volume of the reactor expressed in m 3 , whereas the inlet and outlet depend on the hydraulic load of wastewater (m 3 /day). Transformation processes that occur in biological treatment units involve substrate utilization for the formation of new cells. The ratio of biomass formed to substrate consumed is defined as / ; mg of biomass/mg of substrate. The overall rate of utilization ( ) can be defined as [33]: where ( ) is the substrate concentration at time , and is the biomass concentration in the biological reactor (overall system (OS), ST, UAF, or HFCW). An average value of 800 mg/L was assumed as previously suggested by [34]. The final model for data adjustment was: where and are the inlet concentration of substrate (mg/L). The nonparametric identification of the derivative term ( ) ℝ (each component of this vector corresponds to the substrate concentration ( ) = BOD , ( ) = COD, ( ) = TNK, ( ) = NH − N, ( ) = ON was accomplished using a differential neural network (DNN) as previously described by [35]. In the case of our experimental setup, the DNN identifier had the following structure: where ( ) ℝ is the identified state vector, and the current concentration ( ) at time given was used for the determination of the set of parameters ( , , / ). ℝ is the Hurwitz matrix and corresponds to the linear part of the model. The adaptable weight matrix ( ) ℝ provides the adaptation of the nominal dynamics, and ( ) ℝ is the activation vector function [35,36]. In order to obtain the set of parameters ( , , and / ), the following optimization problem was set: * = arg min , , where * is the vector of estimated parameters. The Nelder-Mead algorithm was applied to solve such optimization problems. This algorithm minimizes a nonlinear function of variables and can be applied to solve parameter estimation problems [37,38]. MATLAB R2020a was used to fit the concentration curve for each water quality parameter (using the Curve Fitting toolbox) as well as to develop the DNN algorithm (using the Simulink toolbox) and to estimate Monod´s model parameters (using the fminsearch function). Statistical analysis was performed using R version 3.6.2 and RStudio version 1.2.5033 (using the RcmdrMisc, dplyr, ggplot2 and tidyverse packages). A personal computer (Intel i7-4500 CPU at 2.4 GHz, RAM 8 GB) was used.
The accumulation and transformation terms are dependent on the volume V of the reactor expressed in m 3 , whereas the inlet and outlet depend on the hydraulic load of wastewater Q (m 3 /day). Transformation processes that occur in biological treatment units involve substrate utilization for the formation of new cells. The ratio of biomass formed to substrate consumed is defined as Y x/s ; mg of biomass/mg of substrate. The overall rate of utilization r(t) can be defined as [33]: where S(t) is the substrate concentration at time t, and X is the biomass concentration in the biological reactor (overall system (OS), ST, UAF, or HFCW). An average value of 800 mg/L was assumed as previously suggested by [34]. The final model for data adjustment was: where V and S 0 are the inlet concentration of substrate (mg/L). The nonparametric identification of the derivative term dS(t) dt R 5 (each component of this vector corresponds to the substrate concentration S 1 (t) = BOD 5 , S 2 (t) = COD, S 3 (t) = TNK, S 4 (t) = NH 3 − N, S 5 (t) = ON was accomplished using a differential neural network (DNN) as previously described by [35]. In the case of our experimental setup, the DNN identifier had the following structure: whereŜ(t) R 5 is the identified state vector, and the current concentrationŜ(t) at time t given was used for the determination of the set of parameters (µ max , K s , Y x/s ). A R 5x5 is the Hurwitz matrix and corresponds to the linear part of the model. The adaptable weight matrix W 1 (t) R 5x5 provides the adaptation of the nominal dynamics, and σ(t) R 5 is the activation vector function [35,36]. In order to obtain the set of parameters (µ max , K s , and Y x/s ), the following optimization problem was set: where k * is the vector of estimated parameters. The Nelder-Mead algorithm was applied to solve such optimization problems. This algorithm minimizes a nonlinear function of n variables and can be applied to solve parameter estimation problems [37,38]. MATLAB R2020a was used to fit the concentration curve for each water quality parameter (using the Curve Fitting toolbox) as well as to develop the DNN algorithm (using the Simulink toolbox) and to estimate Monod´s model parameters (using the fminsearch function). Statistical analysis was performed using R version 3.6.2 and RStudio version 1.2.5033 (using the RcmdrMisc, dplyr, ggplot2 and tidyverse packages). A personal computer (Intel i7-4500 CPU at 2.4 GHz, RAM 8 GB) was used.

System Performance
The sampling procedure yielded a 24 by 11 data matrix (24 observations for each of 11 parameters), resulting in a total of 264 values. Out of the 24 observations gathered for each parameter, six corresponded to each of the sampling points. Table 1 shows the average and standard deviation of the measurements corresponding to each sampling point and provides an overview of the overall performance of each unit of the system. It can be observed that mean concentrations decrease as wastewater advances through the sampling points. Table 1. Average values of water quality parameters measured at each sampling point.  As seen in Table 2, the system efficiently removed organic matter with average removal efficiencies of 90 ± 5% for BOD 5 and 90 ± 5% for COD. Each treatment stage (ST, UAF, and HFCW) contributed to the reduction of organic matter in the wastewater; however, most of the organic load was eliminated in the first two stages (the anaerobic stages, ST and UAF) of the system with an average mass reduction for BOD 5 of 201.6 ± 60.2 mg/L and 243.2 ± 105.1 mg/L, respectively, and of 330.8 ± 92.8 mg/L and 345.2 ± 136.9 mg/L, respectively, for COD. These results are consistent with the optimal levels of DO and temperature [39] in the ST and UAF shown in Table 3 (the mean DO levels and temperature were 0.26 ± 0.13 mg/L and 0.72 ± 0.53 mg/L and 21.1 ± 2.0 • C and 20.9 ± 2.3 • C for the ST and the UAF, respectively). In addition, as shown in Table 2, there was a significant degradation of organic matter in the HFCW, with an average mass removal of 103.3 ± 22.1 mg/L for BOD 5 and of 171.0 ± 37.0 mg/L for COD. This degradation could be attributed to different anaerobic pathways, due to the low oxygen levels, as observed in Table 3 (0.58 ± 0.67 mg/L of DO was measured at the outlet of the system). Additional DO measurements taken in nine secondary points within the HFCW were even lower (0.28 ± 0.15 mg/L) than those reported at the outlet of the system. Anaerobic pathways such as methanogenesis, denitrification, fermentation, or sulfate reduction can be occurring within the HFCW as previously reported for CW where low DO levels occur [40]. The organic load entering the HFCW was considered adequate: approximately 3.64 g of BOD/m 2 day, which does not exceed the recommended value 6 g of BOD/m 2 day for horizontal flow CW [41]. However, this low organic load can limit nitrogen removal, as a ratio of 5 COD/N is recommended for the effective removal of nitrogen in subsurface CW [42]. In the studied system, the average COD/N ratio was around 1. For TSS, the overall removal efficiency was high, reaching an average of 97 ± 0% during the entire sampling period ( Table 2). As shown in Table 2, the highest performance in terms of solids removal was found in the ST, which eliminated 257.7 ± 65.2 mg/L of TSS, while UAF and HFCW eliminated 38.8 ± 15.0 mg/L and 42.5 ± 11.4 mg/L, respectively. These results prove that the ST efficiently fulfills its function of preventing the accumulation of solid particles in the UAF and the HFCW, thus avoiding obstruction (clogging), which is one of the main problems present in this type of system after long periods of operation [43]. After 4 years of operation, this treatment plant presents no evidence of clogging events. Table 2 shows an average removal efficiency of 51% for TKN and 45% for NH 3 for the overall system, which corresponds to an overall mass reduction of 173 ± 20 mg/L of TN and 69 ± 7 mg/L of NH 3 (Table 2). de Anda et al. [16] obtained similar results for the TKN, reporting a removal efficiency of 48% for the same system; however, the mass reduction increased by 80% (from 96.02 mg/L).
In this study, an initial average load of 336.92 ± 12.17 mg/L TKN (Table 1) was found, and this concentration can be compared with industrial loads [33]. The wastewater that is treated in the system comes from the sanitary services, the cafeteria, and the laboratories of the R&D center, where research on industrial and plant biotechnology as well as food technology is performed, producing residual effluents that can be rich in nitrogen. Furthermore, the excess presence of various substrates (such as ammonia/ammonium and nitrite) in wastewater can inhibit the growth of Anammox bacteria when influent concentration is high [44]. Since in this study, a very atypical high load of TNK was observed, and the removal efficiency of TKN could be affected by substrate inhibition. As shown in Table 2, most nitrogen removal occurs in the HFCW, with a 95 ± 26 mg/L mass reduction of TKN compared with mass reductions of 42 ± 38 mg/L and 36 ± 8 mg/L that occurred in the ST and the UAF, respectively. In accordance with the DO levels observed at the HFCW inlet and outlet (Table 3), it can be assumed that the main nitrogen removal mechanism in the wetland is denitrification [45].
As shown in Table 2, the nitrite and nitrate removal rates were found to be 91% and 95%, respectively, corresponding to 2.69 ± 1.68 mg/L and 0.25 ± 0.07 mg/L of global mass reduction. According to Table 2, the greatest mass reduction for both nitrogenous compounds occurred in the ST, with removal rates of 1.90 ± 1.82 mg/L for NO 2 − and 0.12 ± 0.05 mg/L for NO 3 − . The concentrations of nitrates and nitrites entering the treatment system were low, as can be expected for domestic wastewater, for which the presence of these pollutants is commonly close to zero [33]. However, the low levels of NO 2 − and NO 3 − at all sampling points, along with the high mass removal rates of TKN and NH 3 , indicate efficient nitrogen degradation in the treatment process through denitrification, especially since there is no evidence of accumulation of intermediate components of this degradation pathway (nitrates and nitrites) [40].

First-Order Models
First-order kinetic models were used to model the removal kinetics of each measured contaminant for the overall treatment system and for each stage (ST, UAF, and HFCW). The resulting models are summarized in Table 4. These models assume that outlet contaminant values are proportional to inlet values and exponentially related to the hydraulic retention time (HRT). The kinetic parameters are important for evaluating the performance of the modeled process, C 0 is an average value of the inlet load, and k values are interpreted as decay rate constants (day −1 ) and are indicative of the removal efficiency. p-values were used to determine the statistical significance of each adjusted kinetic parameter. Higher values represent better removal, and lower volumes are thus required in treatment units to accomplish sufficient wastewater treatment [46]. The goodness of fit of the models developed for the overall treatment system was found to be good for BOD 5 , COD, and TSS, with R 2 values equal or greater than 0.90. A decent fit was obtained for TKN, NH 3 , ON, and NO 2 − values with an R 2 ranging between 0.83 and 0.89. Finally, a poor adjustment was found for NO 3 − , with an R 2 equal to 0.61.  Furthermore, the p-values for the C 0 and k terms were below 0.05. The models developed for the ST exhibited a decent fit only for the cases of TSS with R 2 of 0.87; however, the p-values for C 0 and k were below 0.05, except for NO 3 − . The better fit obtained for TSS is advantageous, since a reduction in this parameter at this stage is a priority [16]. In the case of the UAF, a good fit was obtained for BOD 5 and COD with R 2 values equal to 0.81 and 0.85, respectively; in addition, the p-values of C 0 and k for BOD 5 , COD, TSS, NO 2 − , and NO 3 − were below 0.05. Finally, in the case of the models developed for HFCW, p-values of C 0 and k were below 0.05 for all pollutants. However, decent R 2 (0.82-0.90) values were only obtained for the TSS, NH 3 , and NO 3 − models.
In Figure 3, it is easy to observe the differences in decay rates between stages. Observing the adjusted first-order kinetic curve (red line), it can be seen that in the first two stages (ST and UAF), the removal rates of organic matter (COD and BOD 5 ) are greater than that of the HFCW, as indicated by the curve inclination. In the case of TSS, it is clear that the highest removal rates occur in the ST, with a lower contribution of the UAF. Almost constant removal rates were observed for TKN, NH 3 , and ON at the UAF and the HFCW; however, ST presented important variations. These variations can be related to the ammonification process, which converts dissolved organic nitrogen into ammonia [47]. As mentioned previously, the R&D produces residual effluents that can be rich in dissolved organic nitrogen compounds as proteins, amino acids, nucleic acids, and urea, among others. The ammonification of these compounds can cause a variable rise on NH 3 concentrations, as there are simultaneous reactions of degradation and generation of NH 3 occurring in the ST [47], thus affecting the adjustment of the first-order kinetic models. A higher removal rate was also observed in the ST and the UAF, in comparison with the HFCW, for NO3 − and NO2 − .
For COD, the models developed for the ST and the UAF yielded values that were at least 10 times higher than those reported previously [25,48] for anaerobic reactors (specifically, up flow anaerobic sludge blanket (UASB) and continuously stirred tank reactors (CSTR)). The reported COD For COD, the models developed for the ST and the UAF yielded k values that were at least 10 times higher than those reported previously [25,48] for anaerobic reactors (specifically, up flow anaerobic sludge blanket (UASB) and continuously stirred tank reactors (CSTR)). The reported COD decay rate constants for the anaerobic treatment of wastewater in fixed film reactors range from 0.01 to 0.1 [49], which are lower than those obtained here for the OS, the ST, and the UAF, but similar in range to those obtained for the HFCW. Furthermore, the k values obtained for the OS, the ST, and the UAF were within the range of those reported as first-order hydrolysis rate constants for COD [50]. This is important, as hydrolysis is the first step of anaerobic organic degradation by methanogenesis [14]. Based on this information, we can suggest a highly efficient organic removal in the system, with the major contributions occurring in ST and UAF. In addition, we can expect a determinant contribution for organic matter hydrolysis occurring in the ST and the UAF. In the case of the HFCW, reported k values for COD, TSS, and TN are generally between 40 and 60 times higher than those obtained in this work [31]. Using the equation reported by [31] that correlates the volume-based decay rate constant (days −1 ) with the area-based decay rate constant (m days −1 ), we can compare our findings with reported area-based decay rate constants [46,51], for which values are between 5 and 10 times greater than those obtained in this study, except for k values for TSS models for ST and OS. In previous studies where CW were the sole system units, this is attributable to the designs of different systems [31,46,51]. Additionally, in those cases, lower contaminant loads were received in the influent to the treatment plants, which were around 5-15 times lower than those received this study, for organic matter, suspended solids, and nitrogen. This information reinforces the statement that wastewater with high pollutant loads requires a multistage system combining different treatment technologies such as anaerobic reactors and wetlands to achieve efficient treatment; however, for waters with lighter pollution loads, the use of a treatment system based only on CW is possible. Finally, the exceptionally high k values obtained for TSS in the ST can be explained by the maximum removal of this contaminant accomplished in this specific stage, reducing the contribution of the next two steps. Thus, the high removal rate of TSS is obtained mainly by the ST, highlighting the important contribution of this stage to overall system performance. Table 5 summarizes the most relevant equations generated in the case of the multiple linear regression models. Models with an R 2 greater than 0.99 and p-values of predictor variables lower than 0.1 were considered relevant. Obvious correlations such as COD-BOD 5 and TKN = NH 3 + ON were disregarded. For the complete treatment system, relevant models relating to nitrogen output (TKN and NH 3 ) with DO and pH measurements were obtained. These models are relevant, since they would allow an estimation of the nitrogen concentration that exits the system from parameters that are easy to measure in situ with sensors. Models relating variables that require analytical processes to be measured with those that can be measured through sensors are useful, as this facilitates the monitoring of the performance of the treatment system. Eight models were found with this characteristic and are marked in Table 5. DO is an important parameter that determines the specific pathway and rate of nitrogen removal [52]. In addition, pH is also relevant, as it affects the performance of the microbial communities involved in wastewater treatment [53].  These types of models were reported previously but with lower goodness of fit (R 2 ranging from 0.67 to 0.74) [21,54]. It is useful to highlight one of the models for the complete system that indicates that there is an important relationship between the suspended solids and ammoniacal nitrogen that are removed in the system. This relationship has been previously reported, with R 2 ranging between 0.64 and 0.67 [55]. For all stages (ST, UAF, and HFCW), a consistent relationship was found between nitrogen forms and organic matter loads. Models relating nitrogen forms with organic matter (BOD 5 and COD) can be explained by removal biological mechanisms, in which multiple groups of bacteria are involved, whose metabolic pathways can be interconnected [56,57]. Carbon availability can become a limiting factor for nitrogen removal processes such as heterotrophic denitrification; furthermore, in the studied wetland, adequate temperature, DO concentrations, and pH conditions exist for this process to occur [58][59][60]. As this relationship was constant for models in every stage of treatment, we can suggest a relevant interaction between different bacterial species responsible for the removal of nitrogen and organic matter contributing to the removal performance of the overall treatment system. Regression models relating nitrogen and organic matter were reported previously for CW, with R 2 values of approximately 0.9 [55,61]. In the case of the ST, two models were obtained associating TSS with conductivity and organic matter. Associations between TSS, conductivity, and organic matter could be explained by sedimentation processes, where a carryover of salts and organic matter occurs [62]. Additionally, methanogenic activity could also occur inside the tank within flocs that contribute to the removal of solids, organic matter, and even nitrogen [63,64]. Relations between solids and organic matter were also reported in previous regression models with R 2 values of 0.6-0.7 [55,61]. For the UAF, the dependency of NH 3 removal on temperature is expected, since temperature can be determinant for removal rates and performance [65,66]. For the UAF, a relevant model showed the dependence of NH 3 on temperature. For HFCW, models were generated relating organic matter with CT and DO. A relationship between the removal of organic matter (BOD 5 and COD) and oxygen levels was expected, since DO determines the biological pathway of organic removal (aerobic or anaerobic) and, to a large extent, the system's performance [67][68][69]. Regression models presenting associations between organic matter and oxygen levels presented R 2 values of 0.78 [54]. Conductivity is related to the content of salts in wastewater, which can interfere with the microbial processes occurring within the wetland [70]. Additionally, the dependence of TSS removal on temperature was found to be relevant. Finally, a model exhibiting the dependency of TSS removal on temperature within the HFCW is expected as a result of the effect of temperature on physical and biological processes involved in the removal of solids, such as adsorption and microbial degradation [71,72].

Mass Balance Models
The results for nonparametric identification by the DNN for COD removal are shown in Figure 4 as an example. Figure 4a shows the comparison between the original data and the identified system. The concentration curves for each treatment stage are shown. It can be noted that the model curve (which was obtained through an interpolation of experimental data) and that obtained through the DNN are consistent. Figure 4b shows the derivative of the substrate dS(t) dt (S= COD in this case) estimated by the neural network. It should be mentioned that it is always negative but decreases progressively in magnitude. It coincides with the interpretation of the concentration graph (Figure 4a). At the beginning, the removal rate is high, but it decreases until it reaches practically zero at the end of the process. Figure 5 shows the values of the Monod parameters (µ max , K s , and Y x/s ) for the mass balance of COD in HFCW. These parameters were obtained using the Nelder-Mead algorithm, and r(t) was obtained as presented by Equation (6). For K s and Y x/s , it must be observed that the variation is very low; the points indicate the values obtained by the algorithm for each t. DNN are consistent. Figure 4b shows the derivative of the substrate ( ) ( = COD in this case) estimated by the neural network. It should be mentioned that it is always negative but decreases progressively in magnitude. It coincides with the interpretation of the concentration graph ( Figure  4a). At the beginning, the removal rate is high, but it decreases until it reaches practically zero at the end of the process.  Figure 5 shows the values of the Monod parameters ( , , and / ) for the mass balance of COD in HFCW. These parameters were obtained using the Nelder-Mead algorithm, and ( ) was obtained as presented by Equation (6). For and / , it must be observed that the variation is very low; the points indicate the values obtained by the algorithm for each . The average values and standard deviation for , , and / obtained from the mass balance model are detailed in Table 6. The mean values obtained for / and were similar for all three stages. / is the rate at which biomass is generated by substrate consumption, and is a measure of microbial affinity to the substrate [33]. In contrast, values of presented important variations between stages. The is related to specific characteristics of the microbial groups involved in the removal processes, as it is their maximum growth rate [20]. Consequently, the combined interpretation of these parameters indicates that the variations between stages are more related to the composition and characteristics of the bacterial groups within each reactor than to the composition of the substrate, in this case, wastewater. This indicates a high robustness of the treatment system, as the variations in wastewater composition do not cause variations in performance.  The average values and standard deviation for µ max , K s , and Y x/s obtained from the mass balance model are detailed in Table 6. The mean values obtained for Y x/s and K s were similar for all three stages. Y x/s is the rate at which biomass is generated by substrate consumption, and K s is a measure of microbial affinity to the substrate [33]. In contrast, values of µ max presented important variations between stages. The µ max is related to specific characteristics of the microbial groups involved in the removal processes, as it is their maximum growth rate [20]. Consequently, the combined interpretation of these parameters indicates that the variations between stages are more related to the composition and characteristics of the bacterial groups within each reactor than to the composition of the substrate, in this case, wastewater. This indicates a high robustness of the treatment system, as the variations in wastewater composition do not cause variations in performance. For BOD 5 , Y x/s, and K s values ranged between 0.65-0.67 mg of biomass/mg of substrate and 57.62-64.47 mg/L for the three stages, respectively, representing the highest values when compared with the other water quality parameters modeled. These parameter's values are within the range of those reported previously for BOD 5 for the biological treatment of domestic wastewater [73]. The estimated parameters for COD are also in accordance with those reported previously for anaerobic reactors [25,74,75]. This is true even for HFCW, probably because the OD level in wetlands are very low, causing an anaerobic degradation of organic matter in the system. These parameters indicate efficient organic matter removal in each stage; however, the µ max values decrease from the initial to the final stage for both cases BOD 5 and COD. This is expected, since the volume of the reactor decreases in that same order. For NH 3 , TKN, and ON, Y x/s values were around 0.32, 0.56, and 0.56 mg/mg, respectively. K s values were around 4.6, 5.75, and 3.45 mg/L for the same parameters, respectively, with standard deviation values below 0.001. For all nitrogenous parameters (TKN, NH 3 , and ON), the results were very similar in range and comparable to those reported by other studies in CW, ANNAMOX processes, and anaerobic reactors [17,26,[62][63][64]. These results prove the efficiency of the system for nitrogen removal not only in the HFCW but also in the anaerobic reactors ST and UAF.

Conclusions
Mathematical modeling provides relevant and valuable information about the performance and specific processes of wastewater treatment. In this study, first-order kinetic models were useful to predict outlet pollutant concentrations based on inlet concentrations with significant precision. These models require the estimation of only two parameters (C in and k). In particular, the parameter k is of direct interpretation for the performance of the system. Multiple linear regression models proved to be efficient to correlate different physicochemical parameters in accordance with the logical relations explained by biological and physical processes of contaminant removal. The applicability of linear regression models was highlighted, as these relate variables that require analytical processes to be measured with those that can be measured through sensors. These models can reduce the cost and time required to monitor the performance of the treatment system. Furthermore, multiple linear regression models provide an interpretation of the extent to which the independent variables (predictors) affect the independent variable (response). In the case of mass balance models, the differential neural network (DNN) allowed for the nonparametric identification of the derivative of the concentration of each water quality parameter measured, while the optimization algorithm allowed for the estimation of the parameters that express the rate of microbiological consumption as proposed by Monod´s model, thus providing a description of the dynamics of pollutants removal within the treatment system. Y x/s and K s values obtained from the mass balance models were similar between stages; in contrast, the µ max values varied significantly from one stage to the next. This fact indicates that performance variations could be related to characteristics of the specific microbial composition in each stage (ST, UAF, and HFCW); however, variations in inlet wastewater composition are not a cause of performance variability. This information also allows an improvement of the design of this wastewater treatment system composed of two anaerobic units (ST and UAF) and an HFCW.