Model-Based Analysis of Increased Loads on the Performance of Activated Sludge and Waste Stabilization Ponds

: In a way to counter criticism on low cost-effective conventional activated sludge (AS) technology, waste stabilization ponds (WSPs) offer a valid alternative for wastewater treatment due to their simple and inexpensive operation. To evaluate this alternative with respect to its robustness and resilience capacity, we perform in silico experiments of different peak-load scenarios in two mathematical models representing the two systems. A systematic process of quality assurance for these virtual experiments is implemented, including sensitivity and identiﬁability analysis, with non-linear error propagation. Moreover, model calibration of a 210-day real experiment with 31 days of increased load was added to the evaluation. Generally speaking, increased-load scenarios run in silico showed that WSP systems are more resilient towards intermediate disturbances, hence, are suitable to treat not only municipal wastewater, but also industrial wastewater, such as poultry wastewater, and paperboard wastewater. However, when disturbances are extreme (over 7000 mg COD · L − 1 ), the common design of the natural system fails to perform better than AS. Besides, the application of sensitivity analysis reveals the most inﬂuential parameters on the performance of the two systems. In the AS system, parameters related to autotrophic bacteria have the highest inﬂuence on the dynamics of particulate organic matter, while nitrogen removal is largely driven by nitriﬁcation and denitriﬁcation. Conversely, with an insigniﬁcant contribution of heterotrophs, the nutrient removal in the pond system is mostly done by algal assimilation. Furthermore, this systematic model-based analysis proved to be a suitable means for investigating the maximum load of wastewater treatment systems, and from that avoiding environmental problems and high economic costs for cleaning surface waters after severe overload events.


Introduction
Conventional activated sludge (AS) systems, the most common application for sewage treatment, have recently been criticized due to their low cost-effectiveness, with high energy demand and limited recovery potential [1]. While the applicability of advanced technologies, such as membrane bioreactors, sand filtration, and aerobic granulation, is still being questioned in developing countries due to the barrier of affordability, waste stabilization ponds (WSPs) appeared to be an inexpensive, but effective, alternative. In fact, thanks to low cost and minimal operation and maintenance (O&M) the systems, with an average flow of around three L d −1 , for 210 days. The recipe of the artificial wastewater was based on the OECD [17] guideline, resulting in a COD of 275 g O 2 ·m −3 , total nitrogen (TN) of 40 g N·m −3 , and total phosphorus (TP) of 7 g P·m −3 . A specific configuration was selected for each treatment type, i.e., the Wuhrmann process for AS systems, and conventional WSPs, including three compartments in series, an anaerobic (AP), a facultative (FP), and a maturation (MP) pond. These configurations were chosen because of their basic, conventional, and common settings for removing organic matters and nutrients. The overview of the experimental setup of two treatment systems is illustrated in Figure S1 (Supplementary Material A) and its detailed description can be found in Ho, Van Echelpoel, Charalambous, Gordillo, Thas and Goethals [3]. To ensure the stability of the systems, their start-up period was maintained for 179 days, with samples being collected and analyzed two times per week. After the stabilization period, the peak load scenario was implemented in three phases. Standard artificial wastewater was fed to the systems for eight days of the first phase. Subsequently, the influent pollutant concentrations were tripled and then kept for five days of the second phase. The recovery of the systems was followed for the next 18 days with the initial wastewater in the third phase. The data collected from this peak load period were used for model calibration before the two models were applied in the scenario analysis.

Model Description
Two models representing the two systems were developed in the software, AQUASIM 2.1 [18]. Specifically, two compartments of an AS model were simulated within the completely mix module and connected via an advective link. The DO concentration was kept at 4 g O 2 ·m −3 in the first aerated compartment while there was no aeration supply in the second compartment. The AS model was developed, based on ASM2d, with 20 processes, including different microorganisms and microbial transformations. Sludge recirculation in the system was simulated as the biomass in the influent. A WSP model with three different compartments was simulated within the plug flow module, with each compartment containing different microorganisms, processes, and variables, hence increasing significantly the model complexity compared to the AS model. Within the plug-flow module, its assumption of a homogenous depth profile of variables is valid since the photosynthesis activity in facultative (FPs) and maturation ponds (MPs) can extend down to a depth of 20-30 cm, creating a homogenous aerobic condition in these oxidation ponds whose height is 20 cm [19,20]. Regarding the processes in the WSP model, following the Constructed Wetland Model No.1 of Langergraber, et al. [21], the anaerobic processes in the anaerobic compartment included the hydrolysis process of slowly biodegradable COD and the metabolism of anaerobic bacteria. FPs and MPs were modeled following the model of Sah, et al. [22], without pathogen removal and anaerobic bacteria due to relatively no sludge accumulation. The removal of phosphorus was assumed to have occurred mainly by chemical precipitation, hence, biological phosphorous removal processes were not considered. In addition, the gas exchange at the pond surface was included with reaeration of oxygen and ammonia volatilization. Light attenuation was accounted in the model to describe the exponential decrease of light intensity with depth according to Beer's Law. More importantly, the cycle of day and night of light intensity are also included in the model to illustrate the fluctuation of algal photosynthesis and respiration activities. The interactions between all involved microbial groups are summarized in Supplementary Material B. The initial and influent conditions of the model variables are calculated and shown in Supplementary Material C. The details of stoichiometric matrix, kinetic rate expressions, and their values in the two models are shown in Supplementary Material D.

Screening for Important Parameters
Sensitivity analysis (SA) evaluates the degree to which model inputs affect model output, and from that the universality and robustness of these parameters can be further investigated. The following technique of SA was proposed in Brun, et al. [23]. Firstly, the model is defined as described above. Subsequently, the prior uncertainties of model parameters and inputs (θ j ) are estimated, based on the literature [9,24,25]. The scale factors of model outputs, which are used to make the results for the various model endpoints comparable, are calculated, based on their mean concentrations in the pond system [26]. After these steps, Gaussian error propagation is applied to compute the sensitivity function, s ij , for each state variable, y i , against the changes of any parameter, θ j , as shown in Equation (1): where ∆θ j represents uncertainty ranges (∆) of the parameters and inputs were divided into three categories: Precisely known parameters (class 1), ∆θ j = 5%; poorly known parameters (class 3), ∆θ j = 50%; and moderately known parameters (class 2), ∆θ j = 20%; sc i is a scale factor of the state variable, y i . The values of ∆θ j and sc i are listed in Supplementary Material E. The importance ranking was then determined based on the sensitivity measure, δ msqr j , in Equation (2): A high δ msqr j means that the value of the parameter, θ j , has an important influence on the simulation result while the sensitivity measure of zero indicates that the parameter has no effects on the model outputs.

Identifiability Assessment of Parameter Subsets
After the determination of δ msqr j , the identifiability of the parameter subset, K, is estimated to avoid the compensation effect of the changes in the parameter values, θ j , on the model output, y i . This identifiability analysis of Brun, Reichert and Kunsch [23] is based on two different measures, collinearity index, γ, and determinant measure, ρ. The collinearity index represents the compensability of the parameter subset, K, which can be calculated as shown in Equation (3): where S K is an n × K submatrix of the normalized matrix, S = s ij , with a normalized value of the sensitivity function, s iK = s iK / s iK ; β = (β 1 , . . . , β K ) T is a vector of coefficient of the length, k, with the constraint, β = 1; EV is the eigenvalue of [ S T K S K ]. γ K quantifies the degree of approximate linear dependence of the s ij of the parameters. A value of γ over 10 indicates a poor identifiability of the parameter subsets [25]. The second criterion, the determinant measure, ρ, is defined as shown in Equation (4): The determinant measure, ρ, is proposed by Weijers and Vanrolleghem [27], which combines the information provided by δ msqr j and γ. A high value of ρ K indicates a low value of γ and a high value of δ msqr j , hence a good "conditional identifiability" of a parameter subset [26].

Model Calibration
After SA and IA, a subset of identifiable and influential parameters is chosen for calibration with the data collected from the shock-load experiment. Generally, to objectify the calibration process, a function representing the agreement between the model and data is defined to demonstrate the desire to fit the model to the data [28]. In this case, the weighted sum of squares (WSS) of the residuals is minimized, from which selected parameters are calibrated using the simplex algorithm in the AQUASIM 2.1 software [18]. To calculate the WSS, each residual is divided by a scale factor (sc i ) of the corresponding variable, y i , from that the residuals become non-dimensional as sc i and the model output have the same dimension [25].
In this equation, y meas,i is the i-th measurement, y i θ j is the calculated value of the model output corresponding to the i-th measurement and evaluated at the time and location of this measurement, and n is the number of data points.

Scenario Analysis
To evaluate the impact of increased loads on the performance of the two systems, several prospective scenarios were simulated using the UNCSIM package [29]. These theoretical simplified scenarios include quantitative changes in the pollutant concentrations, which are expected to occur in real WWTPs. According to Joseph [30], four baseline scenarios representing two intermediate and two extreme conditions are a good starting point for the development of further policy scenarios. As such, we implemented four scenarios of different wastewater strengths, i.e., 2, 5, 10, and 25 times higher than the standard OECD domestic wastewater. These scenarios allow investigation of the increased load threshold of which the two systems can still be able to recover to the initial conditions, hence, the conclusion of their robustness and recoverability can be withdrawn. More importantly, to ascertain the water management decision process, the Monte Carlo simulation was applied to compute the model output uncertainty as a result of the prior uncertainty of model parameters and inputs (listed in Supplementary Material E). 500 sets of samples were generated using the Latin Hypercube Sampling (LHS) technique, which provides a sufficient coverage of the parameter space with optimal computation time [31]. Scenarios caused by hydraulic overloading were not considered in this study as their practical solution is normally based on flow management with the manipulation of the proportion of by-passed water flow and the design of control structure.

Sensitivity Analysis
Sensitivity analysis (SA) is designed as a tool to identify the most influential model parameters for the variability of the state variables. For the overparameterized models, SA is considered very useful as the model output is often strongly influenced by few key inputs [16]. In this case study, we investigate particularly the degree to which model inputs affected two groups of model output, i.e., organic matter and nutrient removal.

Activated Sludge Model
The Most Influential Parameters for Organic Matter Removal Organic matter (OM) content in wastewater can be measured via COD, which includes slowly biodegradable particulate COD (X S ), fermentable and readily biodegradable soluble COD (S F ), fermentation products as acetate (S A ), inert soluble, and particulate COD (S I and X I ), and COD from bacterial biomass. However, not all of these components are of equal importance [32]. Particularly, the COD fraction from microorganisms and the inert COD are also not in our interest because of the marginal variation of their values. To identify the most influential model parameters for the first three COD fractions, i.e., X S , S F , and S A , we calculate the proportions of sensitivity function (s ij ) for each of the three state variables against the changes of each parameter, θ j , over the total sum of sensitivity measures. As seen in Figure 1, the most significant parameters influencing the concentration of X S is related to autotrophic bacteria, i.e., b A , µ A , K A P , K A O2 , and Y A . These parameters are responsible for more than 90% of the total variance of the concentration of X S . The main degradation process of the particulate COD is hydrolysis, which is sensitive to the availability of dissolved oxygen in the system. On the one hand, the presence of dissolved oxygen is needed for the aerobic hydrolysis, but on the other hand, the prevalence of O 2 inhibits the process rate of hydrolysis under anoxic and anaerobic conditions. Therefore, the parameters related to autotrophs, which are one of the main oxygen consumers, can affect considerably the removal of particulate COD. The influence of oxygen on OM removal in the AS system is also indicated in the case of S F in which the saturation/inhibition coefficient of heterotrophs for oxygen (K H O2 ) explains 20% of its total variance. Also noteworthy is the substantial contribution of saturation coefficient for PHA storage in PAOs (K P PHA ) in the total variances of soluble COD, i.e., around 50% for both S F and S A . The high importance ranking of this kinetic parameter is also recorded in several studies on the parameter identification of biological wastewater model in the literature [26,33,34]. conditions. Therefore, the parameters related to autotrophs, which are one of the main oxygen consumers, can affect considerably the removal of particulate COD. The influence of oxygen on OM removal in the AS system is also indicated in the case of SF in which the saturation/inhibition coefficient of heterotrophs for oxygen (K ) explains 20% of its total variance. Also noteworthy is the substantial contribution of saturation coefficient for PHA storage in PAOs ( K ) in the total variances of soluble COD, i.e., around 50% for both SF and SA. The high importance ranking of this kinetic parameter is also recorded in several studies on the parameter identification of biological wastewater model in the literature [26,33,34].

The Most Influential Parameters for Nutrient Removal
Nutrient content in municipal wastewater mainly includes two fractions of nitrogen and phosphorus. Particularly, the total nitrogen concentration contains particulate and soluble Kjeldahl nitrogen (TKN), and nitrate-and nitrite-nitrogen (SNO) [32]. Particulate Kjeldahl nitrogen, as the sum of nitrogen bound to all organic particulate fractions, is not investigated due to its marginal values compared to the value of soluble Kjeldahl nitrogen, which is dominated by ammonium-nitrogen (SNH), likewise with particulate phosphorus. As such, the most influential model parameters for nutrient removal represented by the variance of SNH, SNO, and SPO4 are shown in Figure 2. Regarding the first state variable, it appears that nitrification is a main process of ammonium removal in the AS Figure 1. Ten most significant parameters influencing the variability of slowly biodegradable particulate COD, fermentable and readily biodegradable soluble COD, and fermentation products as acetate in the activated sludge (AS) systems. The results in the radar graphs are the proportion of sensitivity function (s ij ) for the state variable against the changes of each model parameter over the total sum of sensitivity measures. The description of the influential parameters in the graph can be found in the Supplementary Material F.

The Most Influential Parameters for Nutrient Removal
Nutrient content in municipal wastewater mainly includes two fractions of nitrogen and phosphorus. Particularly, the total nitrogen concentration contains particulate and soluble Kjeldahl nitrogen (TKN), and nitrate-and nitrite-nitrogen (S NO ) [32]. Particulate Kjeldahl nitrogen, as the sum of nitrogen bound to all organic particulate fractions, is not investigated due to its marginal values compared to the value of soluble Kjeldahl nitrogen, which is dominated by ammonium-nitrogen (S NH ), likewise with particulate phosphorus. As such, the most influential model parameters for nutrient removal represented by the variance of S NH , S NO , and S PO4 are shown in Figure 2. Regarding the first state variable, it appears that nitrification is a main process of ammonium removal in the AS system, which is indicated via the growth of ammonium and nitrite oxidizers as Y A , µ A , and K A P and are responsible for up to 70% of the total variance. On the other hand, parameters related to denitrification, i.e., Y H , K A P , and b H , contributes around 50% of the total variance of S NO , indicating the significant role of denitrification in nitrite and nitrate removals. Besides, hydrolysis also appears to be important to the variability of S NO , which can be explained by the fact that, apart from the influence, hydrolysis is the main source of S F, the only electron donor of the denitrification process in this case. As such, the changes in values of K h , K X , and n hy can affect considerably the availability of S NO in the AS systems. Note that the process of bacterial assimilation plays an insignificant role in both nitrogen and phosphorus removal. The latter process is mainly carried out via PAOs, the group of organisms that have the ability to accumulate phosphorus in excess of normal metabolic requirements [35]. As such, from the radar graph, most of the influential parameters for S PO4 are related to the metabolic process of PAO.
Water 2018, 10, x FOR PEER REVIEW 7 of 17 system, which is indicated via the growth of ammonium and nitrite oxidizers as YA, µA, and K and are responsible for up to 70% of the total variance. On the other hand, parameters related to denitrification, i.e., Y , K , and bH, contributes around 50% of the total variance of SNO, indicating the significant role of denitrification in nitrite and nitrate removals. Besides, hydrolysis also appears to be important to the variability of SNO, which can be explained by the fact that, apart from the influence, hydrolysis is the main source of SF, the only electron donor of the denitrification process in this case. As such, the changes in values of Kh, KX, and nhy can affect considerably the availability of SNO in the AS systems. Note that the process of bacterial assimilation plays an insignificant role in both nitrogen and phosphorus removal. The latter process is mainly carried out via PAOs, the group of organisms that have the ability to accumulate phosphorus in excess of normal metabolic requirements [35]. As such, from the radar graph, most of the influential parameters for SPO4 are related to the metabolic process of PAO.

The Most Important Parameters Driving Model Outputs
After calculating the , the further step of SA, following the procedure of Brun, Reichert and Kunsch [23], is the computation of importance rankings of the parameters ( ) to detect the parameters driving the variability most in the model outputs. Having the highest influence on model outputs and significant contribution to model uncertainty, 20 parameters with the highest value of are presented in Table S16 (Supplementary Material G). The ranking shows that the PAO-

The Most Important Parameters Driving Model Outputs
After calculating the s ij , the further step of SA, following the procedure of Brun, Reichert and Kunsch [23], is the computation of importance rankings of the parameters (δ msqr ) to detect the parameters driving the variability most in the model outputs. Having the highest influence on model outputs and significant contribution to model uncertainty, 20 parameters with the highest value of δ msqr are presented in Table S16 (Supplementary Material G). The ranking shows that the PAO-related parameters, accounting for eight out of the 20 most important parameters, drive the variability in most of the model outputs, which is in line with previous studies on parameter identifiability of ASM2d, i.e., Brun, Kuhni, Siegrist, Gujer and Reichert [26], Ferrero, Chai, Diez, Amrani and Lie [33], and Liau, Shoji, Ong, Chua, Yeoh and Ho [34]. However, in contrast to the low importance of autotroph-related parameters in these enhanced biological phosphorus removal (EBPR) systems, high sensitivity of the model output toward nitrification is indicated via the position of µ A and Y A in this post-denitrification system.

Waste Stabilization Pond Model
The Most Influential Parameters for Organic Matter Removal As can be seen from Figure 3, the following parameters were found to have a significant impact on the variability of X S (listed in the decreasing order of importance): K HMB H2 , K AMB O2 , b A , and b ALG . The first two coefficients suggest the important role of methanogenic bacteria in anaerobic digestion, which is the main removal process of X S . While the low concentration of H 2 as the electron donor in pond systems is the limiting factor of the growth of HMB, a very low value of K AMB O2 , 0.0002 g O 2 .m −3 , indicates the high sensitivity of this bacteria toward O 2 . The influence of the oxygen level in the system on OM removal is displayed in the presence of the decay rate of algae, the main oxygen producer, and autotrophs, the main oxygen consumer, whose changes also significantly affect the total variance of S F and S A . In fact, as algal photosynthesis is the only source of O 2 for heterotrophs to degrade S A in the FPs and MPs, the decay rate of algae appears to be a main contributor to its total variability, together with other photosynthesis-related parameters, including the light extinction coefficient and light saturation constant. Also noteworthy is that AMB proves to have a more important role in removing S F compared to HMB in anaerobic digestion as b AMB contributes to 35% of the total variance of S F , which can be explained by the fact that their electron donor, S A , is more available than H 2 .
Water 2018, 10, x FOR PEER REVIEW 8 of 17 related parameters, accounting for eight out of the 20 most important parameters, drive the variability in most of the model outputs, which is in line with previous studies on parameter identifiability of ASM2d, i.e., Brun, Kuhni, Siegrist, Gujer and Reichert [26], Ferrero, Chai, Diez, Amrani and Lie [33], and Liau, Shoji, Ong, Chua, Yeoh and Ho [34]. However, in contrast to the low importance of autotroph-related parameters in these enhanced biological phosphorus removal (EBPR) systems, high sensitivity of the model output toward nitrification is indicated via the position of µA and YA in this post-denitrification system.

Waste Stabilization Pond Model
The Most Influential Parameters for Organic Matter Removal As can be seen from Figure 3, the following parameters were found to have a significant impact on the variability of XS (listed in the decreasing order of importance): K , K , bA, and bALG. The first two coefficients suggest the important role of methanogenic bacteria in anaerobic digestion, which is the main removal process of XS. While the low concentration of H2 as the electron donor in pond systems is the limiting factor of the growth of HMB, a very low value of K , 0.0002 g O2.m −3 , indicates the high sensitivity of this bacteria toward O2. The influence of the oxygen level in the system on OM removal is displayed in the presence of the decay rate of algae, the main oxygen producer, and autotrophs, the main oxygen consumer, whose changes also significantly affect the total variance of SF and SA. In fact, as algal photosynthesis is the only source of O2 for heterotrophs to degrade SA in the FPs and MPs, the decay rate of algae appears to be a main contributor to its total variability, together with other photosynthesis-related parameters, including the light extinction coefficient and light saturation constant. Also noteworthy is that AMB proves to have a more important role in removing SF compared to HMB in anaerobic digestion as bAMB contributes to 35% of the total variance of SF, which can be explained by the fact that their electron donor, SA, is more available than H2.

The Most Influential Parameters for Nutrient Removal
As shown in Figure 4, the variability of nutrient variables was affected significantly by the rate of growth and decay of algae. Particularly, around 40% of the total variance of S NH within the pond systems was caused by the metabolisms of algae, suggesting the important role of both algal assimilation and oxygen supply to nitrification of autotrophs in ammonium removal. While the decay of bacteria is a minor source of S NH in the system, which is reflected by a small contribution percentage of decay rates of other bacteria, these processes contributed around 25% of the variability of S NO and S PO4 . These results suggest that as a result of long hydraulic retention time, bacteria and algae in the pond systems can release nutrients into the water body via their decay process. This reflects on the list of influential parameters for S PO4 containing high contribution percentage of phosphorus fraction in different components, which are related to the decay process of algae and different bacteria. In contrast to this contribution, the disappearance of parameters related to heterotrophic bacteria in nitrogen removal and the chemical precipitation process in phosphorus removal suggest their marginal role in pond systems.

The Most Influential Parameters for Nutrient Removal
As shown in Figure 4, the variability of nutrient variables was affected significantly by the rate of growth and decay of algae. Particularly, around 40% of the total variance of SNH within the pond systems was caused by the metabolisms of algae, suggesting the important role of both algal assimilation and oxygen supply to nitrification of autotrophs in ammonium removal. While the decay of bacteria is a minor source of SNH in the system, which is reflected by a small contribution percentage of decay rates of other bacteria, these processes contributed around 25% of the variability of SNO and SPO4. These results suggest that as a result of long hydraulic retention time, bacteria and algae in the pond systems can release nutrients into the water body via their decay process. This reflects on the list of influential parameters for SPO4 containing high contribution percentage of phosphorus fraction in different components, which are related to the decay process of algae and different bacteria. In contrast to this contribution, the disappearance of parameters related to heterotrophic bacteria in nitrogen removal and the chemical precipitation process in phosphorus removal suggest their marginal role in pond systems.

The Most Important Parameters Driving Model Outputs
Based on the values of δ msqr , the 20 parameters contributing the most on the variability of the outputs and uncertainty of the WSP model are presented in Table S17 (Supplementary material G). Noticeably, the five highest ranking parameters driving the most output variance are the parameters that can affect the process rate of the algal photosynthesis. Indeed, besides µ ALG , light availability appears to be highly influential on the performance of both algal activity and the pond system. High sensitivity of the model output toward other physical parameters representing the temperature influence, i.e., θ Tw and T w , is also observed. Indeed, temperature and other climatic conditions, such as solar radiation, are very important to the pond performance [36]. Also noteworthy is the high importance of the flow rate on the performance of the WSP system. Finally, two stoichiometric parameters (Y FB and Y HMB ) among the top ten parameters suggest the important role of anaerobic bacteria in the first compartment of the WSP system.

Activated Sludge Model
The possibility of compensating a change in model output caused by a change of one parameter by the others is evaluated via two measures, i.e., the collinearity index, γ, and determinant measure, ρ. The top 20 parameters of importance ranked in Table S12 are divided into three functional groups, including hydrolysis, autotroph, and PAO-related process, to investigate their correlations via γ and ρ (see Table 1). Turning to parameters related to PAO, it is noticeable that the top four parameters show a relatively small γ of 4.98 and a high ρ of 24.92, which makes them identifiable from the available data. On the other hand, we observe a strong interdependence between the PHA saturation coefficient (K P PHA ) and these four parameters, which could be a result of the high influence of K P PHA on both the growth and polyphosphate storage processes of PAO [37]. A strong correlation between µ A and b A is recorded, which infers the simultaneous estimation of these two parameters is likely to fail. Considering the determinant measure, a higher value in the PAO group indicates its more influential role on the system performance compared to the other two groups. After inspecting the correlation among parameters of all subsets, γ and ρ are calculated for the possible subsets of size 10 and 11. They can be combined by four hydrolysis parameters, three linked to autotrophs, and five PAO-related parameters. Noticeably, most of the subsets with size 11 have γ higher than the threshold of 10 while their value range of ρ fluctuates between 10 and 13. Hence, after assessing all of the combination group of size 10, the subset in bold is selected based on its low value of γ and high value of ρ compared to other subsets.   Table 2 shows the collinearity index (γ) and determinant measure (ρ) of the top 20 influential parameters, which are divided into five functional groups, i.e., physical, anaerobic, algal activity, autotrophic, and heterotrophic processes, to investigate their correlations. As seen in Table 2, the value range of former criterion stays low, between 1 and 5, which suggests that the proposed threshold of Brun, Kuhni, Siegrist, Gujer and Reichert [26] can be irrelevant for classifying non-identifiable parameters in this case. As this subjective threshold can be varied in different case studies, the combination with another criterion, such as the determinant measure, is necessary [38]. According to the value of ρ, the physical and algal-related parameters appear to have a high influence on the model outputs while parameters linked to autotrophic bacteria can generate less impact on the performance of the pond system. It is noticeable that when the subset size increases in all of the groups, except for the first group, there is a significant drop of the value of ρ. Hence, the chosen combined subsets have the size of 12 and 13, including four physical parameters and two parameters of each of other groups. From the obtained results, it is revealed that there is a small difference in the value of γ and ρ between the combined subsets of size 12 and 13. Noticeably, the higher value of the second criterion belongs to subsets that include more physical parameters. As such, with the highest value of ρ, the parameter subset of size 13 in bold, including all five physical parameters, is selected for calibration. Table 2. Collinearity index (γ) and (ρ) of selected parameters subsets. Based on their function, the parameters are categorized into different process groups, i.e., physical, anaerobic, algal activity, autotrophic, and heterotrophic processes. The subset in bold is selected for parameter estimation based on values of γ and ρ.

Model Calibration
After the establishment of a subset of 10 selected parameters for the AS model, an automatic calibration was performed on the basis of the collected data until the convergence in simplex algorithm is achieved in AQUASIM [18]. The original and calibrated results are shown in Table 3. The first conclusion that can be drawn is that the calibrated model evidences the increase of the presence of heterotrophs in the AS systems, which is reflected by the increase of µ H (2.55%) and the significant decrease of b H (45.27%). A higher heterotrophic biomass means a higher removal efficiency of organic matters and nitrogen, which is demonstrated via the considerable drop of the WSS in Table 5. From the calibrated results, the removal of TN is expedited due to the higher process rate of nitrification in the calibrated model as a result of the increase of µ A and Y A , i.e., 8.83% and 15.19%, respectively. The changes in the value of PAO-related parameters, i.e., q PP , b PAO , Y PAO , and Y PHA , appears to have little impact on the removal of TP as the value of WSS stays similar after the calibration. In general, a better predictive performance is obtained in the calibrated model in estimating the removal efficiencies of the AS systems as the total WSS drops by around 60% of its original value, demonstrating the improvement of the calibrated model.  Table 4. In general, there are few changes in the value of selected parameters after the calibration. The most changes are located in parameters involved in the processes of algal photosynthesis and autotrophic bacteria, which also explains the small difference of the error measure before and after the calibration in Table 5. The calibrated model shows a higher degree of photosynthetic activities via its higher values of µ ALG and pH, with an increase of 15.03% and 6.77%, respectively. It is also observed that a higher amount of autotrophic bacteria is present in the calibrated model as a result of a lower decay rate (b A ) and higher yield coefficient (Y A ), 9.73% and 12.5%, respectively. Greater consumption of ammonium in the growth of algae and the expansion of nitrification by autotrophs can explain the better fit of the model in predicting the nitrogen removal. Conversely, insignificant roles of other bacteria and physical parameters barely improve the model's goodness of fit regarding OM and phosphorus.
Comparing the two models, it can be drawn from Table 5 that by having much lower WSS of model outputs, the AS model displays higher accuracy than the WSP model. This fact can be associated with the higher complexity of the natural treatment system when its performance is highly dependent on external factors, such as climatic conditions, leading to significant spatial and temporal variations of their removal efficiency. Especially noteworthy is that the WSP model shows its difficulty in predicting the nitrogen removal of the system, which, theoretically, can be caused by many processes, including microbial assimilation, ammonia volatilization, adsorption, and nitrification/denitrification [39]. Moreover, after the calibration, a significant improvement in the model accuracy can be found in the AS model, but not in the other. In fact, as a result of very few changes in the parameters, high values of error measurement are still obtained in the calibrated WSP model, which is not the case in the AS model. This contradictory output can be explained by the overparameterization of the large mechanistic model representing the natural system with numerous parameters and inputs, which makes the number of data collected during the peak load experiment insufficient for a proper calibration.

Scenario Analysis
The robustness of the two systems is analyzed via in silico experiments where different shock-load scenarios are simulated. These virtual experiments allow the provision of numerical responses of environmental systems to possible events in cooperation with a statistical technique for nonlinear error propagation. Figure 5 shows the performance and their uncertainty ranges of the AS and WSP systems in four scenarios varying from moderate to extreme influent concentrations. Generally speaking, the WSP can endure the peak load better than the AS system, except for the most extreme case of 25 times higher strength of wastewater. It is interesting that the natural system can produce comparable results of relatively low effluent concentrations when the wastewater strength increases from double to five times higher while the effluent quality of the AS system reduces by two times. This result highlights the robustness of WSPs in treating not only municipal wastewater with high strength, but also industrial wastewater with composition under 1500 mg COD·L −1 . However, when the wastewater strength increases higher, their durability is reduced. From the wastewater containing around 2700 mg COD·L −1 , WSPs show relatively similar removal efficiencies to AS systems. Especially at the last extreme scenario, their organic removal efficiency during the peak load is only half of that of the conventional biological treatment system. This deterioration reveals that with an extreme organic load of more than 7000 mg COD·L −1 , the conventional design of WSPs appears to be ineffective with very long recovery periods of up to 50 days. It is also noteworthy that the higher degree of complexity in the pond model compared to the AS model causes its broader uncertainty range, illustrating the intrinsic property of large model predictions. In fact, this overparameterization is also reflected in the calibrated value of model parameters, with few deviations compared to initial values as a result of the significantly lower number of available data for parameter estimation compared to the model complexity.

•
We performed in silico experiments of four different shock-load scenarios in two sophisticated mechanistic models representing the two systems, i.e., AS and WSP. A systematic procedure of quality assurance for these virtual experiments was implemented to assess their uncertainty outputs, including sensitivity and uncertainty analysis with non-linear error propagation, and, more importantly, model calibration with a 210-day real experiment with 31 days of an increased load scenario. The simulation outputs highlight that the WSP can generally endure the increased load better than AS system, except with extremely high strength wastewater (over 7000 mg COD·L −1 ), where a specific design focusing on the primary anaerobic pond is needed. From this result, the robustness of WSPs is proved suitable in treating not only municipal wastewater with high strength, but also industrial wastewater, such as poultry wastewater and paperboard wastewater. For further research, different characterizations of these types of wastewater could be applied in the two models to simulate their performance, and from that a concrete conclusion of preferential choice can be withdrawn. Besides removal performance, other factors related to plant footprint, operational and maintenance costs, energy efficiency, and greenhouse emissions should also be considered in this pre-selection process.

Conclusions
• We performed in silico experiments of four different shock-load scenarios in two sophisticated mechanistic models representing the two systems, i.e., AS and WSP. A systematic procedure of quality assurance for these virtual experiments was implemented to assess their uncertainty outputs, including sensitivity and uncertainty analysis with non-linear error propagation, and, more importantly, model calibration with a 210-day real experiment with 31 days of an increased load scenario. The simulation outputs highlight that the WSP can generally endure the increased load better than AS system, except with extremely high strength wastewater (over 7000 mg COD·L −1 ), where a specific design focusing on the primary anaerobic pond is needed. From this result, the robustness of WSPs is proved suitable in treating not only municipal wastewater with high strength, but also industrial wastewater, such as poultry wastewater and paperboard wastewater. For further research, different characterizations of these types of wastewater could be applied in the two models to simulate their performance, and from that a concrete conclusion of preferential choice can be withdrawn. Besides removal performance, other factors related to plant footprint, operational and maintenance costs, energy efficiency, and greenhouse emissions should also be considered in this pre-selection process.

•
The practical sensitivity analysis casts light on the most influential parameters on the performance of the conventional AS and pond systems. Particularly, as the AS system's behavior is strongly dependent on the variability of oxygen, parameters related to autotrophic bacteria, the main oxygen consumer, initiate the most variability of particulate organic matter. PAOs emerges to be a main user of phosphorus whereas nitrogen removal is largely driven by nitrification and denitrification in the AS system. In contrast, the nutrient removal in the pond system is mostly done by algal assimilation while the absence of heterotrophs-related parameters indicates the insignificant role of the denitrification process. Also noteworthy is that the five top parameters in the importance-ranking list are all related to photosynthetic activity, which displays its crucial role in the pond performance. • Model calibration displays a significant improvement in the prediction performance of the AS model, but not the WSP model. This contradictory result can be explained by the overparameterization of the large mechanistic model representing the natural system with numerous parameters and inputs, leading to a high requirement of both the quality and quantity of available data for proper calibration. • The systematic model-based analysis proved to be a suitable mean for assessing the maximum load of wastewater treatment systems, thus avoiding environmental problems and high economic costs for cleaning surface waters after severe overload events. Moreover, these virtual experiments can be also a handy tool to find a proper solution for system overload, which is currently one of the main challenges of pond treatment technology.