Optimization of the Anaerobic-Anoxic-Oxic Process by Integrating ASM2d with Pareto Analysis of Variance and Response Surface Methodology

Abstract

Wastewater treatment plants (WWTPs) are high-energy-consuming units. Reasonable operation strategies can enable WWTPs to meet discharge standards while reducing the operating cost. In this study, the activated sludge model 2d (ASM2d), Pareto analysis of variance (ANOVA), and response surface methodology (RSM) were jointly used to simulate and optimize the operation of a lab-scale anaerobic-anoxic-oxic (AAO) reactor. The optimization objective was to determine the optimal design and operational parameters (DOPs) that could enhance both pollutant removal and energy saving. The DOPs that had significant influence on the optimization objective, such as sludge retention time (SRT), dissolved oxygen (DO), and the ratio of biodegradable chemical oxygen demand to total nitrogen (BCOD/TN), were identified by Pareto ANOVA. The optimal DOPs with SRT of 15 days, DO concentration of 0.5 mg/L, and BCOD/TN of 5.21 were determined by RSM. Under the optimal conditions, the removal efficiencies of NH₄⁺-N, total nitrogen (TN), and total phosphorus (TP) were 96.2%, 76.8%, and 92.8%, respectively, and the annual operating cost was $26.4. Furthermore, this combination of DOPs was validated using a pilot-scale AAO system. The TN and TP removal efficiencies were improved by 11.0% and 5.0%, respectively, and the annual operating cost could be reduced by 15.0%. Overall, this study confirmed that the method integrating ASM2d with Pareto ANOVA and RSM was effective in optimizing wastewater treatment processes.

1. Introduction

In recent years, water environment protection has received increasing attention in China. Many wastewater treatment plants (WWTPs) need to upgrade their treatment processes or configurations to meet stricter wastewater discharge standards [1]. Wastewater treatment is an energy-intensive industry, and energy-saving technologies can enable WWTPs to meet strict discharge standards while reducing the operating cost. For example, the aeration energy consumption of a biological treatment unit accounts for 30.0–70.0% of the total energy consumption of the WWTPs [2]. Sean et al. [3] found that every 1.0 mg/L decrease of dissolved oxygen (DO) concentration in the aeration tank could reduce the energy consumption by 5.0–7.0% with the effluent still meeting the wastewater discharge standards.

However, the activated sludge system is complex and uncertain. It has the characteristics of nonlinear dynamics, multi-time scale in the internal processes as well as a multivariable structure [4]. Optimizing the design and operational parameters (DOPs) is time-consuming and labor-intensive. In the past few decades, many advanced mathematical models of the activated sludge process have been developed to evaluate and predict the treatment performance of various WWTPs, such as the artificial neural network model [5,6], the Neuro-Fuzzy model based on the statistical method [7], and activated sludge models (ASMs) based on biokinetics [8]. These models are effective tools for simulating the complex biochemical wastewater treatment processes. They can economically and effectively guide the operation of the biological treatment processes to achieve an excellent treatment effect and reduce energy consumption. A series of ASMs (ASM1, ASM2, ASM2d, ASM3) developed by the International Water Association (IWA) has been used to simulate complex activated sludge processes and predict biological treatment efficiency under dynamic conditions. Wu et al. [9] used the ASM3 to simulate and optimize the operation of a coking WWTP, reducing the operating cost from 7.2 $/m³ to 6.4 $/m³. Many factors, including DO concentration [10], the ratio of carbon to nitrogen [11,12], the volume ratio of the anoxic tank to the aerobic tank (V_anx/V_aer ratio) [13], the sludge retention time (SRT), and the hydraulic retention time (HRT) [14], can significantly affect the removal efficiencies of pollutants in biological treatment systems. Therefore, it is a multi-factor and multi-objective optimization task to adjust multiple parameters to achieve effluent quality compliance and energy saving at the same time. Xie et al. [15] combined the ASM2d, support vector regression and accelerated genetic algorithm to simulate and optimize the operation of a Carrousel oxidation ditch. The results showed that better effluent quality was achieved by reducing the total HRT, SRT and internal recirculation ratio (IRR) while extending the aerobic HRT. Cao et al. [16] combined the GPS-X and response surface methodology (RSM) for determining the optimal DOPs to improve the total nitrogen (TN) removal efficiency; by optimizing DO concentration and SRT; the removal efficiencies of NH₄⁺-N and TN were increased to 97.1% and 85.3%, respectively. RSM is an experimental design method that can be used to analyze the relationships between explanatory variables and one or more responses. While reducing the number of experiments, RSM can help researchers build models to evaluate the effects of multiple factors, thus achieving optimal conditions for ideal responses [17]. However, selection of the parameter sets used for optimization is a key step in RSM. Analysis of variance (ANOVA) is a commonly used statistical method, which can identify important sources of uncertainty and quantify the individual and interactive effects of factors [18]. Pareto ANOVA based on the Pareto principle is a simplified ANOVA method [19], which can quickly and easily screen out important DOPs for the subsequent optimization step using RSM. Although RSM has already been used to optimize contaminant removal and resource recovery in wastewater treatment processes [20,21,22], Pareto ANOVA and RSM have not been coupled with ASMs to optimize wastewater treatment processes.

In this study, the combination of Pareto ANOVA and RSM for multi-objective optimization was applied to optimize the lab-scale anaerobic-anoxic-oxic (AAO) system based on the ASM2d. Eleven stoichiometric and kinetic parameters were calibrated and validated at different phases according to the operation data of the system, and then the validated model was used to determine the optimal DOPs for the removal of pollutants as well as for energy saving. Finally, the optimal DOPs were tested and validated in a pilot-scale AAO system. The results indicated that the combination of ASM2d, Pareto ANOVA and RSM could optimize the key DOPs accurately and effectively. This research provides a useful tool for optimizing the operation of activated sludge processes.

2. Materials and Methods

2.1. Operation of the AAO Reactor

The lab-scale AAO system used for model calibration and validation consisted of an AAO reactor and a sedimentation tank, with working volumes of 15.3 L and 4.0 L, respectively. The AAO reactor was divided into three compartments (anaerobic zone, anoxic zone, and aerobic zone) by baffles. The working volume of the anaerobic zone was 2.0 L, and the total volume of the anoxic and aerobic zones was 13.3 L. The V_anx/V_aer ratio could be changed by adjusting the position of the baffle. Agitators were set up in the anaerobic and anoxic zones, and the aerobic zone was aerated by microporous aerators. Real wastewater and inoculated sludge were collected from a municipal WWTP in Shanghai, China. The flowrate of the influent wastewater was 20 L/day. The external recirculation ratio (ERR) and the IRR were 300% and 100%, respectively. The SRT was controlled at 15 days.

Three different operation phases were applied to the lab-scale AAO system, namely the calibration phase, validation phase I, and validation phase II. The DOPs of each phase are shown in

Table 1. DOPs of the lab-scale AAO system.

	Calibration Phase	Validation Phase I	Validation Phase II
Vanx/Vaer ratio	1:3	1:3	3:1
Addition of sodium acetate (COD mg/L)	0	120	150
DO concentration (mg/L)	2.0	2.0	2.0
Average temperature (°C)	22.0	20.0	18.5
Influent flowrate (L/day)	20	20	20
Duration (day)	28	24	35

. To examine whether the ASM2d could accurately predict the system performance under different operational conditions, the V_anx/V_aer ratio and the amount of external carbon source affecting the ratio of biodegradable chemical oxygen demand to total nitrogen (BCOD/TN) were set differently in the two validation phases. In validation phase I, the V_anx/V_aer ratio was 1:3, and about 150 mg/L sodium acetate (providing approximately 120 mg/L of chemical oxygen demand (COD)) was added to the anoxic zone. In validation phase II, the V_anx/V_aer ratio was 3:1, and about 190 mg/L sodium acetate (providing approximately 150 mg/L of COD) was added to the anoxic zone. The average concentration of COD in the influent of the lab-scale AAO system was 209 ± 20 mg/L. The physical-chemical method was used to fractionate the influent COD and the result is shown in Figure 1 [23]. The specific influent quality of the model input is shown in Table S1.

The pilot-scale AAO system was established in a municipal WWTP in Anhui Province, China. It also consisted of an AAO reactor and a sedimentation tank made of steel (Figure S1). The working volumes of the anaerobic zone, the anoxic zone, and the aerobic zone were 4.3 m³, 9.4 m³, and 18.7 m³, respectively. The working volume of the sedimentation tank was 8.4 m³. Agitators were set up in the anaerobic and anoxic zones, and the aerobic zone was aerated by microporous aerators. The wastewater and inoculated sludge were collected from the grit chamber effluent and the excess sludge, respectively, in the WWTP. The influent flowrate of the pilot-scale AAO system was 50 m³/day. The ERR and the IRR were 300% and 100%, respectively. The SRT was controlled at 15 days, which was the same with the SRT of the AAO system of the WWTP. The influent COD fractionation of the pilot-scale AAO system is shown in Figure S2, and the DOPs and specific influent quality of the model input are shown in Tables S2 and S3, respectively.

The removal efficiency of pollutants was calculated using Equation (1):

$$\mathrm{Removal}\mathrm{efficiency}=\frac{{\mathrm{M}}_{\inf }-{\mathrm{M}}_{\mathrm{eff}}}{{\mathrm{M}}_{\inf }}\times 100\%$$

(1)

where ${\mathrm{M}}_{\inf }$ is the influent concentration of pollutants; ${\mathrm{M}}_{\mathrm{eff}}$ is the effluent concentration of pollutants.

2.2. Sensitivity Analysis and Simulation Evaluation

Thirty-six kinetic and 6 stoichiometric parameters related to autotrophs, heterotrophs, and phosphorus accumulating organisms (PAOs) were selected for local sensitivity analysis (LSA). The central relative sensitivity (CRS) was used to express the influence of each parameter on the simulation of the lab-scale and pilot-scale AAO systems. The larger the CRS, the greater was the influence of the parameter on the simulation. CRS was calculated as follows [24]:

$${\mathrm{FAS}}_{\mathrm{i},\mathrm{j}}(\mathrm{t})=\frac{\partial {\mathrm{y}}_{\mathrm{i}}}{\partial {\mathsf{\theta }}_{\mathrm{j}+}}=\frac{{\mathrm{y}}_{\mathrm{i}}(\mathrm{t},{\mathsf{\theta }}_{\mathrm{j}}+\Delta {\mathsf{\theta }}_{\mathrm{j}})-{\mathrm{y}}_{\mathrm{i}}(\mathrm{t},{\mathsf{\theta }}_{\mathrm{j}})}{\Delta {\mathsf{\theta }}_{\mathrm{j}}}$$

(2)

$${\mathrm{BAS}}_{\mathrm{i},\mathrm{j}}(\mathrm{t})=\frac{\partial {\mathrm{y}}_{\mathrm{i}}}{\partial {\mathsf{\theta }}_{\mathrm{j}-}}=\frac{{\mathrm{y}}_{\mathrm{i}}(\mathrm{t},{\mathsf{\theta }}_{\mathrm{j}})-{\mathrm{y}}_{\mathrm{i}}(\mathrm{t},{\mathsf{\theta }}_{\mathrm{j}}-\Delta {\mathsf{\theta }}_{\mathrm{j}})}{\Delta {\mathsf{\theta }}_{\mathrm{j}}}$$

(3)

$${\mathrm{CAS}}_{\mathrm{i},\mathrm{j}}(\mathrm{t})=\frac{{\mathrm{FAS}}_{\mathrm{i},\mathrm{j}}(\mathrm{t})+{\mathrm{BAS}}_{\mathrm{i},\mathrm{j}}(\mathrm{t})}{2}$$

(4)

$${\mathrm{CRS}}_{\mathrm{i},\mathrm{j}}(\mathrm{t})=\frac{{\mathrm{CAS}}_{\mathrm{i},\mathrm{j}}(\mathrm{t})\times {\mathsf{\theta }}_{\mathrm{j}}}{{\mathrm{y}}_{\mathrm{i}}(\mathrm{t},{\mathsf{\theta }}_{\mathrm{j}})}$$

(5)

where FAS is the forward absolute sensitivity; BAS is the backward absolute sensitivity; CAS is the central absolute sensitivity; ${\mathrm{y}}_{\mathrm{i}}$ is the model output function; ${\mathsf{\theta }}_{\mathrm{j}}$ is the jth parameter; $\Delta {\mathsf{\theta }}_{\mathrm{j}}$ is the change of the jth parameter; $\frac{\partial {\mathrm{y}}_{\mathrm{i}}}{\partial {\mathsf{\theta }}_{\mathrm{j}}}$ is the absolute sensitivity of the ith model output function to the jth parameter.

When one parameter was changed, the remaining parameters kept the default values, and the absolute values of CRS of all parameters were compared. The steady-state LSA has a good correlation with the dynamic LSA, so the result of the steady-state LSA can be used for dynamic calibration [25].

Sensitivity analysis was used to obtain the parameters which had significant influence on the simulation, and the simulated values which fitted well with the measured values were obtained by adjusting the values of these parameters. The fitting degree between the measured values and the simulated values was characterized by the Theil inequality coefficient (TIC). When TIC is less than 0.3, it can be considered that there is a good fitting degree between the measured values and the simulated values [26]. The TIC calculation is shown in Equation (6):

$$\mathrm{TIC}=\frac{\sqrt{{\displaystyle \sum {({\mathrm{Y}}_{\mathrm{sim}}-{\mathrm{Y}}_{\mathrm{real}})}^2}}}{\sqrt{{\displaystyle \sum {\mathrm{Y}}_{\mathrm{sim}}^2}}+\sqrt{{\displaystyle \sum {\mathrm{Y}}_{\mathrm{real}}^2}}}$$

(6)

where ${\mathrm{Y}}_{\mathrm{sim}}$ is the simulated value, ${\mathrm{Y}}_{\mathrm{real}}$ is the measured value.

The simulation and the sensitivity analysis of the model in this study were carried out on the WEST^® platform. The default values of the stoichiometric and kinetic parameters were used in the first run of the simulation, and some of the parameters were adjusted based on the results of the sensitivity analysis to reduce the differences between the simulated and measured values.

2.3. Screening for Important DOPs

To screen the DOPs that had significant influence on the dependent variables (removal efficiencies of NH₄⁺-N, TN, total phosphorus (TP), and the operating cost), ANOVA was performed with RStudio (version 0.94.102) to obtain the interaction between the DOPs (BCOD/TN, DO concentration, SRT, V_anx/V_aer ratio), and the dependent variables (the principle and operation files of ANOVA can be found in Supplementary Materials). Furthermore, the Pareto analysis was applied based on the results of ANOVA to screen the important DOPs. The contribution rate was used in the Pareto analysis to numerically evaluate the degree of the influence of DOPs on the dependent variables [27]. The F-value conforming to the F-distribution was used to calculate the contribution rate of each DOP. The contribution rates of DOPs were the proportion of the F-value of each DOP to the sum of the F-values of all DOPs [28]. The calculation process of the contribution rate is summarized in Supplementary Material S1. The contribution rate of each DOP was sorted, and the cumulative contribution rate was calculated. A combination of DOPs whose cumulative contribution rate exceeded 80% was considered to have a major impact on the dependent variables. In this study, Pareto ANOVA was performed based on the results of an orthogonal experiment. The orthogonal experiment was simulated on the ASM2d in the lab-scale system and the experimental design table is shown in Table S4. In the orthogonal experiments and optimization experiments of DOPs, influent quality different from that in the model calibration and validation stages was selected (Table S5).

2.4. DOPs Optimization Based on ASM2d

The results of Pareto ANOVA were used to determine the DOPs to be optimized in the RSM. The RSM software used in this research was Design Expert^® 8.0. Box–Behnken design, which is a common second-order experimental design method in RSM [29], was adopted for the experimental design. The optimization experiment was carried out based on the ASM2d in a lab-scale system. RSM was performed on the simulated results, and the influence of each DOP on the second-order linear model was evaluated by ANOVA (details are shown in Supplementary Material S2). DOPs that did not contribute to the final second-order linear model were eliminated and a three-dimensional response surface map was drawn according to the establishment results of the second-order linear model. Multiple responses were optimized simultaneously using Design Expert^® 8.0 to determine the optimal DOPs for both pollutant removal and energy saving.

2.5. Global Sensitivity Analysis

The optimal DOPs obtained from the lab-scale experiment were tested in the pilot-scale AAO system aformentioned. To optimize the operation of the pilot-scale AAO system, global sensitivity analysis (GSA) was used to investigate how the model outputs were influenced by the parameters that changed simultaneously. GSA can compare the influence of different parameters based on comprehensive modeling using all the parameters [30]. The principle of GSA is as previously described by Sin et al. [31]. The parameters in the ASM2d can be divided into two categories: (i) DOPs, such as different biochemical zone volumes, ERR and IRR, etc. (ii) Model parameters (including stoichiometric and kinetic parameters) related to sludge properties in the ASM2d. The probability density function (PDF) and the value range of the parameters involved in GSA are shown in Table S6.

2.6. Analytical Methods

All the samples were filtered through a 0.45 µm membrane before analysis. COD, mixed liquor suspended solids (MLSS), mixed liquor volatile suspended solids (MLVSS), NH₄⁺-N, NO₂⁻-N, NO₃⁻-N, TN, PO₄³⁻-P, TP, and the biochemical oxygen demand for 5 days (BOD₅) were analyzed according to the standard methods [32]. The concentrations of NH₄⁺-N, NO₂⁻-N, NO₃⁻-N, TN, PO₄³⁻-P, and TP were analyzed using an ultraviolet-visible spectrophotometer (L6S, Shanghai INESA, Shanghai, China). The pH was measured using a pH meter (MODEL6010, JENCO, Shanghai, China). The DO concentration was measured using a DO meter (HQ30d, Hach, Loveland, CO, USA). A gas chromatography (6890N, Agilent, Santa Clara, CA, USA) with flame ionization detector and automatic sampler was used to quantitatively determine volatile fatty acids. Molecular nitrogen gas was used as the carrier gas, and DBWA × 125–7332 (30 m × 30 mm × 0.25 mm) was used as the capillary column. The temperature of the sampler and the detector was set at 220 °C and the furnace temperature was run at 55 °C for 1 min, and then raised to 110 °C at 30 °C/min and kept for 1 min. After that, the temperature was raised to 220 °C at 30 °C/min and kept for 1 min. All chemicals with analytical grade in this study were bought from Sinopharm Chemical Reagent Co., Ltd. of China (Shanghai, China).

3. Results and Discussion

3.1. Calibration and Validation of ASM2d in the Lab-Scale AAO System

By ranking the absolute values of CRS of the stoichiometric and kinetic parameters in the ASM2d in terms of effluent quality indexes (NH₄⁺-N, NO₃⁻-N, TN, PO₄³⁻-P, TP, and COD), 11 parameters that had great impacts on the effluent quality were obtained (Table S7).

The hydrolysis rate constant (k_h) had a great influence on the effluent quality. That is because the biodegradation of the slowly biodegradable organic matter (X_S) in wastewater starts from hydrolysis, which is usually slower than the growth of heterotrophic bacteria. Therefore, the hydrolysis process is considered to be the rate-limiting step for the removal of organic matter and nutrients [33]. The yield for heterotrophic biomass (Y_H), the yield for autotrophic biomass (Y_AUT) and the decay coefficient for autotrophic biomass (b_AUT) also had great influences on all effluent quality indexes, while the half saturation coefficient for oxygen (K_O) and the half saturation coefficient for oxygen of autotrophic biomass (K_O,AUT) had great influences on the processes related to nitrogen removal. Except for poly-phosphate requirement per polyhydroxyalkanoate stored (Y_PO), other stoichiometric and kinetic parameters closely related to PAOs in ASM2d had little influence on the phosphorus removal process. This might have been caused by the low content of BCOD in the influent of the reactor, which could not meet the growth of PAOs [34]. Therefore, compared with the heterotrophic bacteria and autotrophic nitrifying bacteria, the contribution of PAOs to the pollutant removal was negligible.

Y_H and the maximum specific growth rate for autotrophic biomass (μ_AUT) were measured [35], which were 0.67 ± 0.01 gCOD/gCOD and 0.71 ± 0.14 day⁻¹, respectively. The measured values of Y_H and μ_AUT were within the reported range of literature and slightly altered in the calibration [36,37]. Other parameters were adjusted according to the result of the sensitivity analysis by comparing the simulated results and the experimentally measured results of the calibration phase. The values of the adjusted parameters are shown in

Table 2. Default and adjusted values of 11 key parameters in the ASM2d based on the lab-scale AAO system.

Parameters	Characterization	Unit	Default Value [8]	Values Used in This Study	General Range
YH	Yield for heterotrophic biomass	gCOD/gCOD	0.625	0.7	0.5–0.67 [36]
bAUT	Decay coefficient for autotrophic biomass	day−1	0.15	0.14	0.05–0.2 [38]
µH	Maximum specific growth rate for heterotrophic biomass	day−1	6	5	2–6.25 [36]
KNH,AUT	Half saturation coefficient for ammonium (substrate) of autotrophic biomass	gSNH/m3	1	0.68	0.1–4 [39]
KNO	Nitrate half saturation coefficient	gSNO/m3	0.5	0.55	0.1–1 [39]
KO,AUT	Half saturation coefficient for oxygen of autotrophic biomass	gO2/m3	0.5	0.4	0.1–4 [39]
ηNO,het	Anoxic growth reduction coefficient of heterotrophic biomass	-	0.8	0.55	0.55–0.85 [8]
QPP	Rate constant for storage of polyphosphate	day−1	1.5	1.275	1–1.5 [40]
ηNO,hyd	Anoxic hydrolysis reduction factor	-	0.6	0.55	0.55~0.9 [8]
bH	Rate constant for lysis and decay of heterotrophic biomass	day−1	0.4	0.45	0.02–1.44 [38]
µAUT	Maximum specific growth rate for autotrophic biomass	day−1	1	0.79	0.2–1.2 [37]

and the values of the remaining parameters in ASM2d were set according to those reported by IWA [8]. Using the adjusted parameters, the simulated values of NH₄⁺-N, TN, PO₄³⁻-P, TP, COD in the effluent and MLSS of the aerobic zone were well fitted to the measured ones (details of simulated and experimentally measured results are shown in Figure S3). The TIC values of each index in the calibration phase were all less than 0.3 (

Table 3. The TIC values concerning effluent quality indexes and MLSS of the aerobic zone of ASM2d in the lab-scale AAO system.

Phase	NH4+-N	TN	PO43−-P	TP	COD	MLSS
Calibration	0.08	0.12	0.06	0.05	0.17	0.10
Validation I	0.07	0.26	0.19	0.15	0.08	0.10
Validation II	0.17	0.14	0.44	0.24	0.20	0.07

), indicating that the adjustment of these 11 parameters was appropriate and effective in reducing the difference between the simulated and measured values.

The calibrated model was further validated using the operation data of the two validation phases (Figure S3). Except that the TIC of the PO₄³⁻-P exceeded 0.3 in the validation phase II, the TIC values of the other effluent quality indexes in each phase were all less than 0.3 (Table 3). This indicated that the model was successfully validated and could be used for further optimization.

3.2. Effects of DOPs on Pollutant Removal and the Operating Cost

To select suitable ranges for DOPs, it was necessary to study how DOPs affected pollutant removal and the operating cost (calculation formulas of the operating cost are shown in Table S8). Previous studies have shown that DO concentration [10], BCOD/TN [11,41], V_anx/V_aer ratio [13], and SRT [14] are the main DOPs affecting the pollutant removal in biological processes. Therefore, the variations of effluent NH₄⁺-N, TN, TP, COD, and the operating cost at different DO concentrations, BCOD/TN, V_anx/V_aer ratio, and SRT were investigated in the calibrated ASM2d.

The simulated results of pollutant removal and the operating cost changing with different DO concentration are shown in Figure 2a. When the DO concentration was in the range of 0.1–0.3 mg/L, the removal efficiency of NH₄⁺-N was only about 3.0% and the removal efficiency of TN was about 20.0%. This is because the nitrifying bacteria in the system were significantly inhibited. When the DO concentration was greater than 0.5 mg/L, the nitrification capacity of the system was restored, and the removal efficiency of NH₄⁺-N gradually increased with the increase of DO concentration. However, when the DO concentration rose to 1.0 mg/L, further increase of DO concentration had little effect on the improvement of NH₄⁺-N removal efficiency, which was consistent with the phenomenon observed by Liu et al. [42]. Excessive DO concentration would also adversely affect the denitrification process and cause a decrease in TN removal efficiency. The TP removal efficiency decreased significantly with the increase of DO concentration. The main reason was that the nitrate concentration in the aerobic zone increased along with the DO concentration. This may have resulted in competition for the carbon source between PAOs and denitrifying bacteria, which led to the poor phosphorus removal [43]. Since most of the COD had been utilized in the anaerobic and anoxic zones, less COD reached the aerobic zone. Consequently, the change of DO concentration had little effect on the COD removal efficiency. The operating cost of WWTPs mainly includes the aeration cost and the sludge disposal cost, and the amount of aeration is closely related to the DO concentration. Thus, the operating cost increased with the increase of the DO concentration.

The simulated results of pollutant removal and the operating cost with different SRT are shown in Figure 2b. When the DO concentration in the aerobic zone was kept at around 1.0 mg/L, the removal efficiency of NH₄⁺-N was almost not affected by SRT (longer than 15 days) and was always around 98.0%. However, the TN removal efficiency was only 65.0%. This may be because the mixed liquor recirculating from the aerobic zone brought oxygen into the anoxic zone and affected the denitrification process. When the SRT was longer than 25 days, the longer the SRT, the lower was the TP removal efficiency. This may be because the decay ratio of activated sludge in the system became higher with the extension of SRT. Thus, the observed biomass yield and the abundance of PAOs decreased [44,45], which led to the poor phosphorus removal performance of the system. In general, microorganisms can make good use of the biodegradable carbon source in municipal wastewater. Therefore, the COD removal efficiency had little change with the increase of SRT, and the effluent COD was always less than 50.0 mg/L, meeting the WWTP effluent discharge standard of China (EDSC) [46]. The operating cost decreased slightly with the increase of SRT. This is because the increase of SRT reduced the daily excess sludge production and the corresponding sludge disposal cost.

Although the removal efficiency of NH₄⁺-N was not affected by the change of BCOD/TN, the increase of BCOD/TN could significantly improve the removal efficiency of TN. When the BCOD/TN increased from 4.07 to 6.16, the TN removal efficiency increased by 30.3%. When BCOD/TN was further increased from 6.16 to 8.25, the TN removal efficiency only increased by 8.6% (Figure 2c). In addition, when the BCOD/TN was 6.16, the effluent quality met with the requirements of the EDSC, and the TN removal efficiency was 72.0%. It can be seen from Figure 2c that, for a system unable to effectively remove phosphorus due to insufficient carbon source, the addition of a carbon source could significantly improve the TP removal efficiency [47]. Meanwhile, when BCOD/TN was further increased from 6.16 to 8.25, little increase in the phosphorus removal efficiency was observed. Moreover, the addition of a carbon source would increase the operating cost of the wastewater treatment unit. When the BCOD/TN increased from 4.07 to 6.16, the operating cost increased by about 33.0%.

According to Figure 2d, the larger the V_anx/V_aer ratio, the higher were the removal efficiencies of TN and TP. The increase in the volume of the anoxic zone improved the anoxic HRT, which was conducive to the enrichment of denitrifying organisms and denitrifying phosphorous accumulating organisms, resulting in higher TN and TP removal efficiencies. The increase in the volume of the aerobic zone would have increased the retention time of nitrifiers in the aerobic zone, which was conducive to the removal of NH₄⁺-N. However, an increase in the volume of the aerobic zone would also increase the required airflow, leading to an increase in the operating cost. The removal efficiency of COD did not change with the V_anx/V_aer ratio and was stable at 89.0%.

3.3. Pareto ANOVA Analysis

Orthogonal experiments were designed and simulated on the ASM2d in the lab-scale system, and Pareto ANOVA was performed based on the simulated results. The orthogonal experimental design table and simulated results are shown in Table S4.

The result of ANOVA analysis for the removal efficiencies of NH₄⁺-N, TN, TP, COD, and the operating cost is shown in Tables S9–S13. Figure 3 shows the Pareto diagram of each objective. DO concentration, SRT, and BCOD/TN had significant impacts on pollutant removal and the operating cost, while the V_anx/V_aer ratio only had a greater effect on phosphorus removal. As most WWTPs use post-chemical methods to remove phosphorus, DO concentration, SRT, and BCOD/TN were selected for subsequent optimization by RSM.

3.4. DOP Optimization by RSM Based on the ASM2d

The value ranges of DO concentration, BCOD/TN, and SRT in the optimization experiments were determined according to single factor experiments (

Table 4. Value range of DOPs in optimization experiments.

DOPs	DO Concentration (mg/L)	BCOD/TN	SRT (Day)
Value range	0.5–2.0	4.07–6.16	10–25

). It can be seen from Figure 2 that when DOPs were within the value ranges in Table 4, the variation of COD removal efficiency was within 3.0%. Therefore, only the removal efficiencies of NH₄⁺-N, TN, TP, and the operating cost were included as objectives (responses) in the subsequent optimization experments.

To avoid the tedious and time-consuming adjustment of the selected DOPs, RSM was used to effectively optimize the multi-variable DOPs at different levels. Three levels of three factors were selected for response surface analysis, and 20 different scenarios were designed based on RSM software (i.e., 20 sets of DOPs were designed for optimization). First, each scenario was simulated based on the ASM2d in the lab-scale system to obtain the simulated value for the optimization objective, and the detailed information of the optimization experiment is shown in Table S14. ANOVA analysis with a confidence level of 0.05 was then used to determine the significance of all variables. The factors that had no significant influence on the response were removed to obtain the final quadratic model for each response, as shown in Table S15. Table S16 summarizes the ANOVA analysis of each regression quadratic equation. The graph of the regression equation was presented as the response surface graph, and the corresponding contour map was plotted using Design-Expert^® 8.0. The 3D response surfaces and contour maps of removal efficiencies of NH₄⁺-N, TN, TP, and the operating cost are shown in Figure 4.

The response surface graphs obtained from the regression equations clearly showed the relationship between factors and responses. For example, the maximum removal efficiency of NH₄⁺-N occurred at the intersection of the maximum values of SRT and DO concentration. Both long SRT and high DO concentration were beneficial to the nitrification process. However, the operating cost rose rapidly with the increase of DO concentration. The decrease of BCOD/TN would lead to a decrease of the removal efficiencies of TN and TP, while the increase of the carbon source would lead to an increase of the operating cost.

Based on the regression quadratic equations, multiple responses were optimized simultaneously by using Design-Expert^® 8.0 to determine the optimal DOPs that could simultaneously achieve high NH₄⁺-N, TN, and TP removal efficiencies as well as low operating cost. The function used for optimization was the expected function. The optimal DOPs obtained by RSM were as follows: SRT = 25 days, DO concentration = 0.5 mg/L, BCOD/TN = 5.22. The results obtained by running the optimal DOPs in the model of the lab-scale system are summarized in Plan A in

Table 5. Comparison between the predicted results by RSM and ASM2d on the WEST^® platform (Plan A and Plan B) as well as the simulated results of ASM2d in different scales under the optimal conditions.

DOPs				Operation Results
Plan A
	SRT (day)	DO concentration (mg/L)	BCOD/TN	RENH4	RETN	RETP	Operating cost ($/m3/day)
RSM	25	0.5	5.22	97.7%	78.2%	93.6%	3.54
WEST	25	0.5	5.22	97.7%	78.6%	89.8%	3.56
RD	-	-	-	0.05%	0.52%	4.19%	0.56%
Plan B
RSM	15	0.5	5.21	95.2%	76.1%	93.4%	3.60
WEST	15	0.5	5.21	96.2%	76.8%	92.8%	3.61
RD	-	-	-	1.03%	0.90%	0.68%	0.28%
WEST (in the lab-scale system)	15	0.5	5.21	96.2%	76.8%	92.8%	3.61
WEST (in the pilot-scale system)	15	0.5	5.21	94.6%	66.1%	53.4%	3.34
RD	-	-	-	1.66%	16.23%	73.92%	8.08%

. The accuracy of RSM predictions was measured by relative deviation (RD). RD was calculated as follow:

$$\mathrm{RD}=\frac{\left|\mathrm{m}-\mathrm{p}\right|}{\mathrm{m}}\times 100\%$$

(7)

where, m is the simulated result of the ASM2d in the lab-scale system, and p is the predicted result of RSM.

It can be seen from Plan A in Table 5 that the prediction of TN removal efficiency and the operating cost under the optimal conditions obtained by RSM was very close to the results obtained by the ASM2d calibrated in the lab-scale system, but the predictions of TP removal efficiency were slightly different. This is because the quadratic term of SRT was considered to have less influence on TP removal efficiency due to its large p-value in ANOVA analysis when establishing the quadratic term equation related to TP removal efficiency, and therefore it was removed from the equation. Although the extension of SRT can affect the TP removal efficiency [45], it is less significant than the influence of DO concentration and BCOD/TN. The SRT of the pilot-scale AAO system was 15 days. Accordingly, the second expectation function was established and the SRT was set as 15 days. The obtained results are shown in Plan B in Table 5. In Plan B, the RD between the predicted values of RSM and the simulated values of the ASM2d in the lab-scale system was very small, indicating that the optimal DOPs obtained from the quadratic polynomial equation were reliable. Compared with plan A, the removal efficiencies of NH₄⁺-N and TN decreased slightly. However, when the NH₄⁺-N removal efficiency was greater than 84.0% and the TN removal efficiency was greater than 72.0%, the NH₄⁺-N concentration (1.1 mg/L) and the TN concentration (7.0 mg/L) in the effluent could meet the EDSC. Therefore, Plan B could meet the discharge requirements and improve the TP removal efficiency at the same time. In terms of the operating cost, Plan B increased by about 1.4% compared with Plan A, and the increase was not significant. Therefore, Plan B can be selected as the optimal condition that can simultaneously meet the requirements of pollutant removal and energy saving.

3.5. Validation of the Optimal DOPs Based on the Pilot-Scale AAO System

To validate the optimal DOPs obtained by the ASM2d based on the lab-scale system and RSM, Plan B was tested using the ASM2d implemented on the pilot-scale reactor. Table S17 shows the default and calibrated values of parameters in the ASM2d in the pilot-scale system and Figure S4 shows the calibration and validation results of the pilot-scale system. The results of the validation were compared with the simulated results of the ASM2d in the lab-scale system (Table 5). In the pilot-scale modeling, when the DO concentration was reduced to 0.5 mg/L, although the NH₄⁺-N removal efficiency was slightly decreased, the effluent NH₄⁺-N concentration ranged from 0.4 to 1.7 mg/L, and it still met the EDSC. In addition, after increasing BCOD/TN, the effluent TN concentration was less than 10.0 mg/L. The operating cost of the WWTP in Anhui province was 3.93 $/m³/day, and it was reduced significantly in the pilot-scale modeling. Changes in aeration cost, sludge disposal cost, mixed cost, and chemicals cost before and after the process improvement of the pilot-scale system are shown in Figure S5. The reduction of the operating cost was mainly due to the reduction of the airflow rate. In this case, the operating cost can be reduced by 15.0% in a one-year perspective. Moreover, it can be seen from Table 5 that, in addition to ammonia removal efficiency and the operating cost, the simulated results of the TN and TP removal efficiencies in the pilot-scale system with optimal DOPs were significantly different from those of the lab-scale system. This may be due to the different influent quality. Sodium acetate as the external carbon source was added to the lab-scale system, resulting in more easily degradable carbon source in the influent, which was conducive to the removal of TN and TP.

GSA analysis was also used to investigate the effect of DOPs and model parameters related to sludge properties on model prediction. The parameters examined included DOPs and model parameters related to sludge properties. In the t-value analysis of standard regression coefficients (SRC), it was found that the absolute values of the t-value of all parameters for the three outputs SRC were greater than 1.96, indicating that all parameters had significant effects on the three outputs. This study only focused on the absolute values of the t-value of the top 10 parameters that had a significant influence on each output variable (Table S18). Overall, the model parameters related to sludge properties had a great influence on the three outputs. Among the top ten parameters influencing the three outputs, only the volume of the anoxic zone (Anoxic.Vol) was a DOP, and it had a negative influence on all three outputs. In addition, the difference in treatment scale can affect the dispersion of sludge flocs in the system. Therefore, the influence of Anoxic.Vol could be regarded as the difference in the degree of sludge dispersion. Among the model parameters related to sludge properties adjusted in the lab-scale and pilot-scale modeling, only k_h was in the top 10 influential parameters (Table 2 and Table S17). In the model in the lab-scale system, the default value was used for k_h. However, in the model in the pilot-scale system, the value of k_h was slightly lower than the default value to reduce the hydrolysis rate of X_S. The hydrolysis of X_S is considered to be the rate-limiting step for the removal of organic matter as well as nitrogen and phosphorus [33]. Therefore, different k_h may cause different predictions of the TN removal efficiency.

3.6. Application Limitation and Perspective

In this study, the method integrating ASM2d with Pareto ANOVA and RSM was used to successfully obtain the optimal DOPs of the lab-scale AAO system. Pareto ANOVA can screen out DOPs to be investigated for optimization by RSM, and this method has been rarely reported. Moreover, the combination of ASMs, Pareto ANOVA, and RSM has not been used to optimize wastewater treatment processes. This optimization method is simple, effective, and easy-to-use, which is suitable for practical applications.

However, it should be noted that the simulation experiments in this study were carried out under the condition of a constant flow. Although the steady-state simulation is not as practical as the dynamic simulation for the actual operation of WWTPs, the results of the steady-state simulation can also be used as a reference. Hvala et al. [48] used the steady-state simulation to choose the upgrade process for an existing WWTP. Koch et al. [49] used the steady-state model to compare the configurations of different biological phosphorus removal processes and preliminarily estimate the nutrient removal efficiency. The steady-state simulation allows the plant performance and key parameters to be determined very quickly. In addition, the buffer tank in some WWTPs can control the excessive inflow of wastewater in the event of peak flow and reduce the diurnal variability of wastewater strength [50,51]. Meanwhile, some WWTPs are designed with long HRT, and the HRT required by the biological treatment unit can be met, despite fluctuations in the influent flowrate. Due to the use of a buffer tank and long HRT, the operation of WWTPs may be relatively stable. In this case, the steady-state simulation may be significant and can provide a theoretical basis for the subsequent dynamic simulation and the operation strategy optimization. However, the practicability of this novel method in the dynamic simulation and whether the optimal DOPs obtained in the dynamic simulation are the same as those obtained in the steady-state simulation, also need to be studied in the future.

The optimization process of DOPs in this study was based on the ASM2d model. However, some restrictive conditions of the ASM2d model make it somewhat inconsistent with the real situation. For example, some complex biological processes are simplified into one step (e.g., nitrification and denitrification), and sludge settleability cannot be accurately predicted by the settling models in ASM2d [52]. However, with appropriate calibration, the difference between the simulation results and the real situation can be reduced to an acceptable range (e.g., the TIC of the simulated value is less than 0.3). In addition, it can be seen from Table S18 that the composition ratios of nitrogen and total suspended solids in the sludge had a significant impact on nitrogen removal and the operating cost. In future studies, the transformation of nitrogen and total suspended solid components during the AAO process should be further studied to improve the accuracy of the model.

4. Conclusions

Based on the ASM2d, the effluent quality and operating cost of the lab-scale AAO reactor were simulated and optimized using Pareto ANOVA and RSM. SRT, DO concentration and BCOD/TN were the DOPs that had great influences on the optimization objectives. Under the conditions of DO concentration = 0.5 mg/L, BCOD/TN = 5.21, and SRT = 15 days, satisfactory effluent quality could be obtained (concentrations of NH₄⁺-N, TN, and TP in effluent were 1.1 mg/L, 7.0 mg/L, and 0.2 mg/L, respectively), and the operating cost could be reduced by 5.0%. A pilot-scale experiment was used to validate the optimal DOPs. The TN and TP removal efficiencies were improved by 11.0% and 5.0%, respectively, and the annual operating cost could be reduced by 15.0%. In addition, GSA analysis revealed that the changes in the degree of sludge dispersion and sludge hydrolytic parameters would lead to uncertainty in the simulation work when expanding the experimental scale. This study provides a useful method combing ASM2d and statistical analysis for the optimization of the design and operation of WWTPs. Further evaluation of this method under dynamic simulation of full-scale WWTPs is still needed.