Evaluation of Three Dynamic Models for Aerated Facultative Lagoons

In this research, three existing dynamic mathematical models for aerated lagoons were applied to an aerated facultative lagoon (AFL) plant for municipal wastewater treatment. The models’ ability to describe the behavior of the primary lagoon was evaluated, and the advantages and limitations of the three models were compared. For this purpose, 7-year operational data from the municipal WWTP were collected and processed; other necessary data—like dynamic temperature values— were estimated. A 2-year period was used for model calibration, while the remaining 5 years served as validation period. One of the models showed poor calibration fit values in the effluent concentration description (R2 of 0.242 and RMSE of 16.8 mg/L); however, with some modifications the adjust was enhanced (R2 of 0.409 and RMSE of 14.0 mg/L); a second model displayed a poor to moderate adjust (0.489 and 13.0 mg/L, respectively), and the third model achieved a moderate fit (0.528 and 11.9 mg/L), though it provided an overestimation of effluent concentration, especially in periods of heavy and frequent rain; this model with some adaptations the adjust was enhanced (R2 of 0.575 and RMSE of 11.4 mg/L). The validation fits are even lower, illustrating the inability of these models to properly describe the AFL behavior. The possible causes of the models’ inadequacy are discussed. Finally, it is established that in terms of AFL behavior description more research and model development are needed. Corresponding authors. C. J. Ábrego-Góngora et al.


Introduction
Many different systems and methods of treating municipal wastewater exist; the choice mainly depends on the initial investment, operating and maintenance cost, and the water quality required at the effluent. Within these systems, lagoon processes are widely used in developing countries, and can be classified into stabilization ponds and aerated lagoons; among the latter, aerated facultative lagoons (AFL) are the most widely used. These are partially aerated, partially mixed lagoons combining aerobic and anaerobic degradation of organic matter. AFL treatment trains can have high BOD removal efficiencies (70% -90%) [1]. Therefore, it is interesting to advance the knowledge of the behavior of these treatment systems in order to make more efficient and/or less expensive the municipal wastewater treatment.
Aerated facultative lagoons (AFL) are a good economic option for wastewater treatment, especially for cities/ towns with intermediate land availability, because AFL are more efficient than facultative ponds and, therefore, require significantly lower surfaces; furthermore, their operation and maintenance are simpler than activated sludge systems [2].
A tool increasingly used to design, analyze, optimize and control systems of wastewater treatment is mathematical modeling [3]. Which in the AFL scope, allows investigating the behavior of such systems under different scenarios without compromising the effluent quality at full-scale.
Mathematical modeling strategies for AFL depend on the mathematical complexity the researcher is willing to adopt, on the quantity and quality of meteorological data, on the available characterization of influent and effluent, and on the available plant operating data. This article discusses the full-scale implementation and evaluation of three AFL dynamic models existing in the literature and involving a limited number of variables (aggregated organic matter/oxygen demand or volatile solids) and processes.

Wastewater Treatment Plant Description
The municipal wastewater treatment plant (WWTP) North of San Luis Potosí, SLP, Mexico, (Planta Norte) is located at the northeast of the city of the same name, at 22˚12'24" of north latitude and 100˚56'49" of west longitude, and 1848 m of elevation. It is a train of four aerated lagoons in series, with an average operation flowrate of 23,126 m 3 /d for the studied period (2004)(2005)(2006)(2007)(2008)(2009)(2010). Figure 1 (elaborated in this research) shows a schematic diagram of the Planta Norte with a brief mention of the main equipment. Details of the plant and, specifically, the primary lagoon are given elsewhere [4].

Sampling and Concentration Analysis
Regular sampling and data recording was performed by Proagua by obtaining grab samples of the WWTP influent and effluent every 4 h, starting at 00:00, recording the date, time, point of monitoring and the instant flowrate. Daily composite samples were prepared from the grab samples, in accordance with the NMX-AA-003-1980 standard [5]. Composite samples were analyzed for total suspended solids (TSS), volatile suspended solids (VSS); chemical oxygen demand (COD) and biological oxygen demand (BOD 5 ) in accordance with the NMX-AA-034-SCFI-2001, NMX-AA-030/1-SCFI-2012, NMX-AA-028-SCFI-2001 standards [6]- [8], respectively. In addition, on Monday and Thursday the same parameters were analyzed at the effluent of each lagoon.

Aerated Facultative Lagoon Models
In the scope of AFL dynamic modeling, there are very few studies published on the scientific literature, most focused primarily on hydrodynamic aspects [9]- [17], other models include the wastewater composition and their transformation processes, but in most of these models the complexity is such that they require advanced waste- water characterization as well as specialized software [18]- [27]. Despite their complexity, these models have not been able to describe the complex behavior of organic matter/oxygen demand in the AFL in a precise or proper manner. This problem motivated the evaluation of existing, simpler aerated lagoon models, based on aggregated organic matter or oxygen demand components. The goal is to determine if they are able to provide a reasonable description of the lagoon behavior at a lower modeling cost.
In the next subsections, three AFL dynamic models (in non-steady state) from the literature are presented [25] [26] [28]. They have been chosen for evaluation in this study.

Bemister Model [28]
The model consists of two state variables representing the mass balance of ultimate BOD (BOD U ) and the total sludge mass in an AFL. His model assumes first-order reactions for substrate utilization rates (aerobic and anaerobic). In this model, it is considered that a fraction of the influent substrate instantly disperses in the lagoon, while the rest settles on the bottom layer, where the accumulation and fermentation of sludge occurs; a fraction of the fermented sludge is released as substrate towards the liquid phase of the lagoon. The model considers one CSTR with continuous flow and was applied to the "Portage la Prairie" AFL-based WWTP sited in Manitoba, Canada. The model equations are: = + = S is the effluent BOD U concentration (mg/L); i ρ is the fraction of influent BOD U dispersed within the body of the lagoon (dimensionless); S i is the influent BOD U concentration (mg/L); R is the lagoon retention time (d); K s is the aerobic degradation rate (d −1 ); s ρ is the fraction of BOD U entering lagoon liquid phase after sludge fermentation (dimensionless); k L is the sludge fermentation reaction rate (d −1 ); L t is the total sludge mass (lbs BOD U ); V is the effective lagoon volume (ft 3 ); V 0 is the initial lagoon volume (ft 3 ); G is the mass density of water (lbs/ft 3 ) and Q is the wastewater influent flowrate (ft 3 /d).
It can be observed that the model has dimensional inconsistencies. Some refer to the International System of Units and others to the English System. A detailed analysis was performed and there was the need to make some modifications to the model for the mass balance presented by Bemister (1978) to have dimensional consistency and to present the expressions in BOD terms; such modifications were: where S i and S are the influent and effluent BOD concentration (mg/L); − L t is the total sludge mass (kg BOD U ); f BODu is the BOD U to BOD ratio; V is the effective lagoon volume (m 3 ); V 0 is the initial lagoon volume (m 3 ); Q is the wastewater influent flowrate (m 3 /d). f l is the sludge accumulation index (dimensionless) and γ l is the sludge specific weight (considered as 1050 kg/m 3 ).

Chagnon Model [29]
This model corresponds to a modified version of the stabilization pond model of [30]; and was applied to the AFL system "As-Samra" in Bal'ama, Jordania. The model consists of three differential governing equations that determine the organic carbon, inorganic carbon and fecal coliform concentrations in the lagoon: According to model author 20 OC is the effluent organic carbon concentration (mg/L); IC is the effluent inorganic carbon concentration (mg/L); CO 2S is the carbon dioxide saturation concentration in the lagoon; CO 2 is the actual carbon dioxide concentration in the lagoon (mg/L); FC is the fecal coliforms concentrations per volume unit (not defined); Q is the flowrate (m 3 /d); V is the lagoon volume (m 3 ); R 12 is the transformation rate from organic carbon to inorganic carbon (d −1 ); R 21 is the transformation rate from inorganic carbon to organic carbon (d −1 ); R 20 is the atmospheric CO 2 re-aeration rate (d −1 ); R 1S is the organic carbon net loss rate (d −1 ); K SC is the half-saturation constant for carbon; R 8S is the overall fecal coliform decay rate (d −1 ). Since there were no data available for CO 2 concentration in the lagoon, (CO 2S -CO 2 ) was considered constant and calculated in parameter estimation. Similarly, K SC was managed as an estimated parameter, since its value was not specified in [29]. In this study, the Equation (5) is not considered since is uncoupled to carbon balance Equation (3) and (4). Since no data were available for measured influentorganic and inorganic carbon concentration, OC was considered proportional to BOD, and IC i as a constant fraction of influent BOD, f l , calculated in parameter estimation; for this, Equation (3) and (4) were modified into Equation (3a) and (4a) expressed in BOD terms. Since the primary AFL of the "Planta Norte" is 1.7 times deeper than the "As-Samra" lagoon, then the R 20 value was changed from 23.86 to 41.02 d −1 , according to the [29] estimation procedure.
where S i and S are the influent and effluent BOD concentration (mg/L); f i is the IC to BOD ratio.

Montalvo et al. Model [31]
A fed-batch model of an aerated facultative lagoon at pilot scale is proposed, consisting of three state variables: substrate (expressed as soluble COD), biomass (expressed as VSS) and the lagoon volume; the pilot study was applied to treat wastewater samples from a winery industry of Molina city, Chile. The model equations are: where, according to the model authors: V is the lagoon volume (m 3 ); F is the volumetric flowrate (m 3 /d); μ m is the maximum specific microbial growth rate (d −1 ); S 0 and S are the influent and effluent substrate concentrations (mg sCOD L −1 ); S nb is the non-biodegradable substrate concentration (mg sCOD L −1 ); X is the cellular or biomass concentration (mg VSS L −1 ); Y is the cell yield coefficient (g VSS g −1 COD); K S is the saturation constant (mg sCOD L −1 ). The model was adapted to a completely-mixed lagoon with equal influent and effluent volumetric flowrate. In this case, dV/dt = 0, so Equation (6) can be discarded. The substrate and biomass balances resulted in the same differential equations, i.e., Equation (7) and Equation (8).
Since there were no data available for sCOD, it was considered proportional to the measured BOD, S nb concentration was assumed as a constant fraction of influent BOD, f nb , and calculated in parameter estimation. Based on these considerations S is expressed in BOD terms, of course, expressing the kinetic parameters in terms of BOD also.

Process of Application, Calibration and Validation of Models
All the mentioned models were built and implemented in MATLAB ® /Simulink ® , their calibrations were performed with the Parameter Estimation tool, by the conjugated gradient optimization method, and using the sum of squared errors of the effluent oxygen demand (measured vs. estimated) as cost function. For this purpose, a period of two consecutive data years (2007-2008) was taken as calibration period and the rest of the time series was used as validation period, this calibration period was chosen due to during this, abrupt oscillations in the influent flowrate occurs with less frequency, in comparing to the rest of the time series.

Inputs of the AFL Models
The daily input data for the proposed models, corresponding to Q, and S, were obtained from the Proagua records of the Planta Norte WWTP corresponding to seven years of measurements, covering the period: January 2004-December 2010. The lagoon temperature input data (T w ) were estimated with the models described below.

Lagoon Temperature Estimation
Influent temperature was estimated using the time-series model developed in a former study [4]. Mean daily ambient temperature and rainfall data (2004-2010) from a nearby weather station [32] were used to dynamically estimate influent temperature. This model was recalibrated for these particular station data, using influent, ambient and lagoon temperatures as well as precipitation data from the period 2006-2007.
The lagoon temperature evolution through the 2004-2010 period was estimated using a dynamic version of the steady state thermal model of [33] as described in [34]. The dynamic temperature model was calibrated with data from [26] for the period 2006-2007, ambient temperature data from [32] and operational data (flowrate, lagoon volume, area) from GrupoProagua SA de CV.

Results and Discussion
All models were applied to the operational data of the Planta Norte primary lagoon, starting with the parameter values proposed by the author of each model. In order for the results to be comparable, the period 2007-2008 was established as calibration period for all models. The model fit evaluations for calibrations and validations were performed by obtaining the coefficient of determination (proportion of variation explained by the model, R 2 ) and the root mean squared error (RMSE). In the figures showing the variations of the measured and estimated concentrations at the lagoon effluent, the dots correspond to the measured data while the solid line show the model estimates, the thick frame within the graph corresponds to the calibration period (2007-2008), while the model validation period comprehends the two sub-periods outside the frame.

Bemister Model [28]
Once resolved the dimensional inconsistencies, the model was executed and its parameters were calibrated. The measured and estimated BOD concentrations at the effluent of the lagoon are presented in Figure 2. The model obtains 0.528 of R 2 with RMSE of 11.9 mg/L for calibration, and 0.345 with 17.6 mg/L for validation, respectively. Table 1 shows the model parameter values (initial and adjusted).
In Figure 2 it can be noticed that although the model describes well the average behavior of BOD in the first   Table 1 shows that the adjusted values are very different to the typical ones; this may be dueamong other causes-to the difference in environmental and operation conditions between Manitoba (Canada) and San Luis Potosí (Mexico), or to the fact that the model does not consider the kinetics of biomass responsible for degradation of the substrate, nor the responsible for sludge fermentation. No reference values or interval was provided for the model parameters by the author in the original study. While the Bemister model [28] visibly follows the general effluent BOD trend-especially during the calibration period-the degree of fit can be considered insufficient in terms of R 2 and RMSE. This is especially true for the validation period and-visibly-around the sharp perturbations appearing in several years during summer/ fall, the latter being associated to intense rain episodes or periods. Another fact to consider is that the Bemister model [28] does not include an oxygen balance, neither the effect of oxygen concentration in the substrate kinetics. While complete-mix aerated lagoons have been traditionally designed using substrate and biomass equations assuming no oxygen limitation-e.g. [34] [35]-dissolved oxygen in aerated facultative lagoons can be low or very low. At the studied plant, effluent dissolved oxygen (DO) ranged 0.0 -2.9 mg/L (average 0.49 mg/L) [26]. Since a typical Monod half-saturation constant value for oxygen in ASM1 is K O,H = 0.2 mg/L [36], the average DO at the lagoon fell in the region where the biodegradation kinetics is still affected by the DO level. So, considering the effect of oxygen limitation could improve the model adequacy. Other aspects that could affect the model performance and are not included in the model are: the actual non-completemix flow in each section of the lagoon, the nitrogen dynamics, and unaccounted effects of intense rains and their associated flow peaks on the primary lagoon performance.
The model considers k L,20 = 0.002 and s ρ = 0.4, according to a study by [37]. In this study, the Bemister model [28] was tested considering that those parameters could also be identified according to the differences in climatic and operation conditions of the stabilization pond studied by [37], with respect to the "Planta Norte" primary AFL. The results were: 0.575 of R 2 with RMSE of 11.4 mg/L for calibration, and 0.362 with 17.2 mg/L for validation, respectively; Table 1 displays the parameter values obtained for this modification to the original model (see "adapted, recalibrated" in Table 1). The estimated BOD effluent is showed in the Figure 1 as BODest, 2; in this figure a significant difference between the original model and adapted is not appreciated. However, the degree of fit in the second case (adapted version) itself increased, this verifies that the assumption made (k L,20 and s ρ with different values than originally proposed) is partially correct.

Chagnon Model [29]
The measured and estimated effluent BOD concentrations are presented in Figure 3. This model achieved a determination coefficient of 0.489 and RMSE of 13.0 mg/L in calibration, and 0.290 with 18.9 mg/L for validation period, respectively. Table 2 shows the parameter values (typical and adjusted). The estimated parameters are remarkably similar to those of the original model.
It can be seen, in Figure 3, a similar behavior to that obtained by the Bemister model [28] for the description of the AFL effluent BOD; however its fit index is lower, and the overestimation at the final period of the time series is greater than that obtained with the Bemister model [28]. For both models, the degree of fit can be considered insufficient, especially for the validation period, and the same considerations on lagoon processes not included in the model can be applied to Chagnon's model [29] too.

Montalvo et al. Model [31]
The results of measured and estimated BOD and VSS concentrations for the lagoon effluent are presented in Figure 4. For BOD, this model achieved R 2 of 0.242 with RMSE of 16.8 mg/L for the calibration period, and 0.230 with 19.6 mg/L for the validation, respectively; while for VSS, it obtains R 2 of 0.078 with 35.8 mg/L of RMSE for the calibration period, and 0.236 with 29.7 mg/L for validation, respectively. The model parameters   values (initial and adjusted) are shown in Table 3. In this case, the model consistently underestimates BOD (112.9 mg/L for estimates vs. 119.9 mg/L for measured data, on average) and, especially, VSS (79.8 mg/L for estimates vs. 104.4 mg/L for measured data, on average). This model presents the poorest fit value of the study. In Figure 4, it can be observed that this model is not able to describe the wide, months long fluctuations in the BOD concentration at the AFL effluent, much less is suitable to represent the sharp pulses observed in the measured lagoon behavior, associated to rain episodes. The model VSS performance is even poorer, with the model estimates clearly departing below the measured data. The considerations done above for the other two models are also applicable here, but the assumption of complete mix is critical when modeling a partially settleable component like VSS. While BOD-soluble or particulate-will largely be biodegraded in the lagoon through a combination of hydrolysis/fermentation and soluble substrate uptake, assuming that VSS will not settle will tend to yield too high modeled VSS, resulting in higher modeled biodegradation rates (since VSS also represents biomass) and lower modeled BOD effluent. When minimizing errors in parameter estimation, the procedure will try to reduce the VSS modeled concentration, by increasing K d (0.60 d −1 , very high) and μ m (1.17 d −1 ). This can give some relief to the BOD lack of fit, but results in a very poor representation of VSS.
Another cause in the lack of fit of this model may be the fact that it does not consider the variation of the parameters with respect to the lagoon temperature, being that according to [39], the kinetic parameters of a model based on biological processes, "per se" are temperature dependent. In Table 3, it can be seen that the adjusted parameter values differ from the identified by the model author; this may be due to the type of substrate utilized by the original model (wastewater from a winery industry) and the substrate and operating conditions of the primary AFL of the "Planta Norte", are very different.
The temperature effect on the model parameters was proved using the Arrhenius correction as follows, in accordance with [39]: 20 20 20 ,20 ,20 ,20 1.07 ; 1.053 ; 1.04 The modified model achieved R 2 of 0.409 with RMSE of 14.0 mg/L for the calibration period, and 0.281 with 19.1 mg/L for the validation, respectively; while for VSS, it obtained an R 2 of 0.062 with 21.7 mg/L of RMSE for the calibration period, and 0.273 with 24.9 mg/L for validation, respectively. Its calibrated parameters values are shown in Table  3 (see "modified, recalibrated"). The proposed modifications enhanced the model degree of fit (from 0.242 to 0.409 of R 2 ), mitigating the underestimation of the BOD (120.1 mg/L for estimates and 119.9 mg/L for measured data, on average) and, particularly, for VSS (107.8 mg/L for estimates and 104.4 mg/L for measured data, on average). However, the fit still remains from poor to moderate. The estimated concentrations at the lagoon effluent are also shown in Figure 4 as BODest, 2 and VSSest, 2, where it can be observed that the modified model-though insufficiently-describes the variation of the measured data in a better way than the original model, specially, for the behavior of the effluent VSS. Finally, the three simple AFL models founded in the literature showed a poor to moderate fit when applied to the operational data of the "Planta Norte" (WWTP) in San Luis Potosí. This can be observed in Table 4, which summarizes the degree of fit (expressed as R 2 ) obtained with each model as well as its associated error (as RMSE) for BOD determination and VSS for the Montalvo et al. model [31], and its modified version, for the selected calibration and validation periods. This table shows that the model with the best fit of all is to the adapted version of the Bemister model [28]. However, it presents an overestimation of effluent concentration specifically at periods when the rains are heavy and persistent. In summary none of the three models provided a good representation of the lagoon effluent BOD.
It is worth mentioning that two additional calibration periods (2004-2005 and 2005-2006) for the three models were tested, finding that it is in the selected period (2007)(2008) in which a best fit was obtained in the description of the BOD behavior, this was verified at all three models.
The evaluation and comparisons made in this research allow confirming that oxygen demand behavior description in AFL requires further investigation for the development of better simple models that can be applied at full-scale plants for different objectives such as, effluent quality control or aeration control.

Conclusions
• Most of the existing AFL models focus on hydrodynamic aspects and to a lesser extent to the composition of the wastewater and its transformation processes, and, due to their complexity, advanced wastewater characterization and specialized software are required for their application. • The Bemister model [28] allows a moderate description of the variation in the BOD measured concentration at the lagoon effluent, but the fluctuations in the measured data are not very well represented by the model, especially in the cases of concentration pulses (associated with flowrate peaks), and overestimates the effluent concentrations in the latter part of the time series. • The Chagnon model [29] provides a poor to moderate fit, with a similar behavior to that obtained with the Bemister model [28], but the overestimation at the final of the time series is greater. The model has the need to generate data for some non-readily measured input variables. • The Montalvo et al. model [31] presents a poor fit value, and is not able to describe the fluctuations measured at the AFL effluent concentrations; neither is capable to describe the pulses observed in the measured behavior. The model does not consider the temperature dependency of the kinetic parameters. This was solved with the modification proposed, and the model degree of fit it increased for BOD and VSS description, especially for the behavior of the solids in the lagoon, its description significantly improved with respect to the original model; however, for both concentrations, the degree of fit still remains from poor to moderate. • In all three models, several not included aspects or processes can be responsible for their lack of fit. Some of these processes are: influence of dissolved oxygen and oxygen supply in the process, non-completely-mixed flow, nitrogen dynamics, and the effect of heavy rains in the lagoon dynamics. Integrating them without adding excessive complexity is a challenge in aerated facultative modeling. • The application, evaluation and comparisons made in this research confirm that a proper description of the oxygen demand behavior in aerated facultative lagoons requires further research for the development of better and simple models that can be applied at full-scale for multiple objectives.