Unlocking the potential of sidestream EBPR: exploring the coexistence of PAO, GAO and DGAO for effective phosphorus and nitrogen removal Open Access

The discharge of phosphorus and nitrogen into the environment is a serious concern due to the impacts of eutrophication (Di Capua et al. 2022). Untreated municipal wastewater typically contains 20–40 mgN L⁻¹ total nitrogen and 3–6 mgP L⁻¹ total phosphorus (Metcalf & Eddy et al. 2013). To reduce these excess nutrient loadings to aquatic environments, wastewater resource recovery facilities (WRRF) use biological methods, namely enhanced biological phosphorus removal (EBPR) and nitrification–denitrification, to sequester phosphorus into biomass and convert bioavailable nitrogen into inert nitrogen. However, EBPR can experience process upsets if adequate influent carbon is not available or if anaerobic zones are undersized or unable to maintain proper anaerobic conditions (Sabba et al. 2023). In these cases, a sidestream EBPR (S2EBPR) process is introduced to improve the reliability of biological phosphorus removal in the presence of weak raw influents. An anaerobic reactor, ‘fermenter’, receives a portion of biomass that will undergo fermentation and produce volatile fatty acids (VFA), essential for biological phosphorus removal. Based on the source of biomass and fermenter operation, S2EBPR exists under four configurations: sidestream return activated sludge (RAS) fermentation (SSR), sidestream RAS fermentation with supplemental carbon addition (SSR + C), sidestream mixed liquor suspended solids (MLSS) fermentation (SSM) and unmixed in-line MLSS fermentation (UMIF) (Izadi & Andalib 2023).

In the SSR or SSR + C configurations, a higher anaerobic mass fraction can be attained in a smaller volume with more stable and lower oxidation-reduction potential (ORP) conditions (Barnard et al. 2017). This provides an ideal environment to supplement carbon for EBPR without competition from other electron acceptors while stored carbon production through fermentation will offset some of the carbon required as a function of the RAS fermenter inventory (WRF 2023). By removing the anaerobic zone from the main flow path through the secondary process and by requiring only a portion of RAS, numerous economical options to upgrade existing facilities to EBPR using existing tankage may be presented by considering S2EBPR concepts (Onnis-Hayden et al. 2018).

EBPR requires the alternation of anaerobic and aerobic or anaerobic and anoxic conditions to enrich phosphorus-accumulating organisms (PAO). Under anaerobic conditions, PAO uptake VFA or other forms of readily biodegradable chemical oxygen demand (rbCOD) and store them internally as polyhydroxyalkanoates (PHA), with the energy for this carbon storage process derived from glycogen and intracellular stored polyphosphate (polyP) (Coats et al. 2017). Under aerobic conditions, PAO uptake orthophosphate (OP) to gain back the intracellular polyP levels via the oxidation of the stored PHA. Some PAO can also perform polyP uptake using nitrate or nitrite as the electron acceptor, referred to as denitrifying PAO (DPAO). PAO and DPAO may compete for organic carbon with ordinary heterotrophic organisms (OHO) and glycogen-accumulating organisms (GAO) (Yang et al. 2010). The GAO cycle internal carbon storage polymers under alternating redox conditions like PAO but do not take up polyP (Roots et al. 2020). While GAO are often perceived as undesirable in EBPR due to their lack of phosphorus removal benefits, GAO may exert minimal impacts on phosphorus removal performance (Stokholm-Bjerregaard et al. 2017), and denitrifying GAO (DGAO) may contribute to total nitrogen removal (Li et al. 2023).

The Calumet Water Reclamation Plant (WRP) is a large facility in the Chicago region (a design average flow of 1,340,000 m³ d⁻¹ and a design maximum flow of 1,650,000 m³ d⁻¹) currently operating as a single-stage aerated process. Given the size of the facility, high-quality conventional anaerobic zones would be expensive to implement. Furthermore, the primary effluent is very carbon-limited (median rbCOD:P ∼ 6:1 during the study period), and internal studies have demonstrated that carbon dosing is necessary for reliable EBPR. Previous work at Calumet has also shown that with aeration turned off at the head of the process, oxygen back-mixing and loading lead to poor anaerobic conditions which made it challenging to sustain EBPR without high carbon supplementation (Andalib et al. 2015). Andalib et al. (2015) showed that by feeding 11,320 L d⁻¹ of MicroC 2000, a glycerin-based carbon source, to the unaerated head of the process, 2–3 mg L⁻¹ of additional phosphorus was removed relative to the control. This translates to ∼30 mg rbCOD mg P⁻¹, deviating from the ideal values (10–15 mg rbCOD mg P⁻¹) due to oxygen and competition for rbCOD. Considering these operational issues, a full-scale S2EBPR demonstration was tested at Calumet WRP and evaluated for efficacy, cost-effectiveness and reliability of phosphorus and nitrogen removal.

DPAO, DGAO and denitrifying OHO all play crucial roles in the combined removal of phosphorus and nitrogen. However, these microorganisms compete for limited amounts of organic carbon, and it is unclear how DPAO and DGAO maintain their desired functions without competing for the same carbon sources (Li et al. 2023). In previous work at this facility, Sabba et al. (2023) used process modeling to show that the presence of GAO did not impact EBPR performance, and GAO accumulation was not predicted by the model. However, Farmer et al. (2023) used 16S rRNA sequencing to show that GAO flourished throughout the entire demonstration, particularly when supplemental carbon dosing was eliminated, and showed that carbon dosing improved phosphorus removal while denitrification was achieved regardless of carbon dosing. Given that GAO bloomed when carbon dosing was removed, and that denitrification was also maintained without carbon dosing, there is a possibility that DGAO were key community members for denitrification. However, no study has attempted to model full-scale data and determine the metabolic potential for DGAO populations in S2EBPR. In this study, wastewater process modeling and molecular data were used to investigate nitrogen removal mechanisms and microbial community trends in the Calumet S2EBPR demonstration. The objectives of this work are as follows:

• – Capture OP and nitrate-N plus nitrite-N (NO_x-N) trends via modeling.

• – Predict trends for a DGAO population in the S2EBPR demonstration.

• – Explore the effects of external carbon addition on the modeled microbial populations.

• – Profile denitrification pathways of PAO and GAO in the S2EBPR process.

The key operational changes for this study are summarized in Table 1. The study period consisted of four phases. Phase I consisted of complete mixing at the beginning and end of fermenters, low mixing/daily bumping for the middle mixers and no carbon in the fermenter, while Phases II–IV consisted of complete mixing at the beginning and end of the fermenter and low mixing/daily bumping for the middle mixers. Phases II–IV were used to test the impact of operational strategies, including carbon dosage and RAS diversion, on P removal performance. Phase II was characterized by the addition of external carbon with a carbon feed of 9,071 kg COD d⁻¹. The feed was then reduced to 6,033 kg COD d⁻¹ during Phase III, and the RAS diversion was reduced from 36,368 to 15,141 m³ d⁻¹ during Phase IV. Additional details on the operation of the fermenter, RAS diversion and external carbon addition can be found in Sabba et al. (2023).

Weekly data were used to create OP and NO_x-N concentration profiles. This information, along with other data, including total suspended solids (TSS), volatile suspended solids (VSS) and influent/effluent parameters, were used to evaluate pilot performance and make operational changes. Wet chemistry data was collected using methods as described in Table 2.

where sCODt_x is the soluble COD concentration (1.2 μm filter) at t = x, mg L⁻¹; OPt_x is the orthophosphate concentration at t = x, mg L⁻¹; NO₃t_x is the nitrate-nitrogen concentration at t = x, mg L⁻¹; sbCOD for P is the sbCOD required per phosphorus and is constant (2.3 gCOD gP⁻¹); COD for denit, is the COD required for denitrification and is constant (4.0 gCOD gN⁻¹) and VSS are the volatile suspended solids, mg L⁻¹.

A parallel configuration was developed in the model for the control basin, and it was built in a similar layout but without a sidestream fermenter and anaerobic zones with the same data sources used for model process control. Modeling of this reactor resulted in close matches for effluent NO_x-N and OP, indicating that the primary effluent loading of COD, nitrogen and phosphorus and the degree of assimilation of nutrients were all accurately captured by the model.

Two models were used in this study. The first model, called the ‘Sumo2 model’, represents the default platform model with all the parameters left at default values. The second model, called ‘modified Sumo2’, included parameter adjustments and simplified model structure to better capture the fermentation rates, effluent OP and NO_x-N trends, in the basin nutrient profiles, and the relative fraction of PAO and GAO (discussed more in the following section). In this model, the reduction factor for anaerobic hydrolysis was reduced substantially from 0.50 to 0.05 to approximately match the apparent fermentation rates observed in bench-scale testing at the start of the pilot. This corresponds to an anaerobic hydrolysis rate of 0.10 d⁻¹ which is at the low end of the range observed for biological nutrient removal (BNR) systems (WRF 2023).

The default Sumo2 model predicted almost complete out-selection of GAO due to most VFA uptake occurring in the RAS fermenter at low ORP conditions. During the periods of excess carbon addition, some GAO activity was predicted when VFA bled through to the higher ORP unaerated portion of the mainstream reactors. To capture GAO accumulation and activity, the ORP/temperature-based GAO uptake activity switch was overridden with a manual input (0–100%). The manual override was tuned to match the average competition apparent from the molecular data (presented below) in which there were slightly more GAO than PAO on average. Additionally, the uptake rate of VFA by GAO (qGAO,GLY) was adjusted to match the uptake rate of VFA by PAO (qPAO,PHA). This was done for simplicity so that the manual activity switch value would represent the actual relative rates. The ORP-based PAO fermentative growth activity switch was overridden with a manual input (0–100%) which was set to 100%. Additionally, the ratio of phosphorus released to VFA stored was decreased from 0.65 to 0.55 g XPP g SVFA⁻¹. This resulted in a better match to OP release in the fermenter and improved the overall efficiency of EBPR by resulting in less OP taking up the same amount of PHA stored. This helped to offset the reduction in carbon stored by PAO when the manual GAO activity switch was introduced to increase GAO abundance to observed levels.

The intent of the modeling component was to calibrate a model of the facility to match the observed performance while predicting the relative PAO and GAO abundance observed in the molecular analysis below. This helps to understand how theoretical stoichiometry from process models used in industry can be paired with molecular data to elucidate how observed performance can be achieved concurrently with the high proliferation of GAO observed. The general approach of using the manual GAO–PAO competition switch as opposed to the ORP-based switch allowed for an empirical tuning that could be helpful in modeling scenarios with uncertainty around PAO–GAO competition or metabolism, or to ease calibration of operating plants to observe performance. However, this should not be taken as a calibrated model, and site-specific PAO–GAO competition should be considered.

16S rRNA gene amplicon sequencing and shotgun metagenomic sequencing were used to profile microbial communities in the S2EBPR demonstration. All DNA were extracted using the FastDNA SPIN Kit for Soil following the manufacturer's protocols. The 16S rRNA gene sequencing approach has been described previously (Farmer et al. 2023). In brief, extracted DNA from S2EBPR-activated sludge underwent polymerase chain reaction (PCR) for amplification using primers 515F and 926R (Parada et al. 2016). The primer sequences are reported in Supplementary Table S2. PCR products were sent to Rush University Genomics and Microbiome Core Facility for sequencing. Raw paired-end reads were imported into QIIME2 (2020.11) for processing into amplicon sequence variants (ASVs) and taxonomic classification against the MiDAS v4.8.1 database (Dueholm et al. 2022).

For metagenomic sequencing, DNA was extracted from biomass samples from the S2EBPR pilot collected on days 33 and 124. Extracted DNA was sent to the NUSeq Core Facility for library prep and sequencing on an Illumina HiSeq 4000 for 2 × 150 bp reads. Raw reads were filtered and trimmed using fastp. Clean reads were used as input to HUMAnN3 for taxonomic and functional profiling (Beghini et al. 2021) and were mapped using DIAMOND (Buchfink et al. 2015) against denitrification genes of key PAO and GAO identified with 16S rRNA sequencing obtained from UniRef100.

Both models performed comparably for NO_x removal prediction, with the modified Sumo2 model showing slight improvement in NO_x removal performance (Figure 3(b)) during the carbon addition starting in Phase II. While the model predictions for NO_x are similar qualitatively, the modified Sumo2 model had lower mean squared error (MSE) values by about 5% across phases (Figure 5(a)). This is also true when comparing the percentage of error of NO_x and OP for the two models (Supplementary Figure S1). Of relevance is Phase II, where the NO_x model prediction for the modified Sumo2 model was closer to the field data, as shown in Figure 3(a). This moderate improvement could be attributed to the change in GAO accumulation and activity. Previous modeling efforts along with the molecular supporting data have highlighted the possibility of GAO bacteria coexisting with and being present at the same or higher concentrations than PAO (Li et al. 2023). Initial model runs in this study indicated that denitrification mediated with storage products may be playing a role in the nitrogen removal observed, and it was theorized that this may represent a GAO ecological niche enabling the high GAO presence. However, refined modeling with more granular operational data showed that this was not necessary to explain the degree of N removal in general and so model changes to specifically have GAO emphasize denitrification-based growth were not necessary.

Furthermore, there were more noticeable differences between the two models comparing the predicted effluent OP. During the no carbon addition phase (Phase I), the Sumo2 model predicted lower EBPR stability which did not consistently match the observed effluent OP. The difference in the models' OP performance became more significant during the carbon addition phase (Phase II), where the Sumo2 model predicted substantial short-term upsets to effluent OP which were not observed, while the modified Sumo2 model generally predicted more accurate OP effluent values. Both models predicted that the transition from no carbon to high carbon (i.e., the Phase I/Phase II transition) happens rapidly with the models showing an almost immediate adjustment, while the full-scale data suggested the potential metabolic acclimation took approximately 1 week. The modified model showed a slight response to an upset occurring in Phase IV due to low RAS diversion to the RAS fermenter while the default Sumo2 model does not show this effect.

The combination of changes improved the model prediction of effluent OP and OP profiles while better matching PAO and GAO molecular data. In the default Sumo2 model, without GAO competition, there was excess PAO uptake capacity much of the time (indicated by the polyphosphate to PAO ratio, PP:PAO, being well below model inhibitory levels). In the modified model, the incorporation of a substantial carbon sink in the form of GAO along with other model adjustments seems to have improved the overall model fit to performance and molecular data. Incorporating the GAO competition helped with model calibration, but the ORP/temperature-based switch was challenging to tune to predict this competition.

To quantify the statistical error in the two models, the MSE was estimated for both models. MSE evaluates the average squared difference between the predicted (model) and observed (field) data. A model with zero error would have an MSE value of zero. Conversely, as the error in the model increases, the MSE value also increases. The MSE was used to determine the quality of fit based on NO_x (Figure 3(a)) and OP (Figure 3(b)) performance for both models for the entire study, as well as for each phase (Figure 5). Figure 5(b) compares the MSE based on OP performance in both models. In general, the modified Sumo2 model predicted OP trends closer to the field data, with all phases showing lower MSE values. Additionally, as observed for the NO_x performance (Figure 5(a)), the modified model had the largest discrepancy from the default Sumo2 model during Phase II. It is noteworthy that Phase II was characterized by the highest carbon addition of the entire study, and the improved performance of both OP and NO_x removal suggested augmented bioP removal and denitrifying activity. To further explore this hypothesis, the modified model was used to explore GAO and PAO population dynamics throughout the four phases.

16S rRNA gene sequencing was used to track PAO and GAO populations as a point of comparison to the process modeling approach (Farmer et al. 2023). PAO Ca. Accumulibacter, Ca. Phosphoribacter, Ca. Dechloromonas phosphorivorans and Tetrasphaera with greater than 0% median relative abundance were detected. GAO Ca. Competibacter and Defluviicoccus were detected with greater than 0% median relative abundance, though Ca. Competibacter had greater abundance (1.14% median relative abundance) than Defluviicoccus (0.06% relative abundance) over the study period. Overall, the relative abundance of the GAO Ca. Competibacter was greater than all detected PAO, which was reflected in the predicted biomass concentrations from the modified model (Figure 6(a)).

An important deviation was observed between the modeling estimates (Figure 6(a)) and 16S results (Figure 6(b)), namely that the modeling estimate did not capture the bloom of Ca. Competibacter observed in Phase I as there was not adequate carbon available during this period to account for this relative abundance of GAO compared with other periods. The contributing factors that led to this bloom are not immediately clear from the reactor operation. One factor could be warmer temperatures during the summer; Qiu et al. (2022) showed that Ca. Competibacter were elevated in abundance over Ca. Accumulibacter and Dechloromonas (also referred to as Azonexus) for most of their study in laboratory-scale EBPR reactors operated at 30 and 35 °C. Peces et al. (2022) also observed that multiple Ca. Competibacter strains are subject to seasonal dynamics in Danish full-scale EBPR systems, with increases in abundance during warmer summer months. However, if temperature was the primary factor influencing Ca. Competibacter enrichment, a longer enrichment period throughout the summer would have been expected. Another factor could be facility-wide selection unrelated to the S2EBPR operation, such as the influent composition. During the study, a completely aerated nitrification basin with separate clarifiers and solids recycle was operated in parallel with the S2EBPR basin. Ca. Competibacter abundances were lower in the nitrification basin than the S2EBPR basin, with 0.11 and 1.14% median relative abundance, respectively, which was expected given the lack of advantageous conditions for Ca. Competibacter growth (i.e., no anaerobic zone). However, the abundance of Ca. Competibacter in the nitrification and S2EBPR basins had a strong positive correlation (Spearman correlation coefficient of 0.71, p = 0.005). Molecular samples were not collected from the nitrification basin during Phase I, but the positive correlation between the dissimilar basins suggests that facility-wide selection factors, such as influent immigration or cross-basin migration from solids dewatering, may have exerted an influence. Finally, it should be noted that the Sumo2 modeling estimates of PAO and GAO populations represent active biomass in each functional group with uniform reaction rates and kinetics rather than an estimate of dynamic bacterial populations.

Overall, nitrate reductases were in greater abundance than the rest of the reductases in the pathway. Differential abundances of denitrification pathway genes have been observed previously in activated sludge systems. For example, in their study of three full-scale activated sludge systems, Xia et al. (2014) found copies of nitrate reductase narG and nitrite reductases nirS and nirK at approximately 10 times higher abundances than nitric oxide reductase norB and approximately three times higher abundances than nitrous oxide reductase nosZ. Additionally, Ca. Competibacter had a greater abundance of N₂O reductase genes compared to all PAO. Partial denitrification by PAO is a well-known phenomenon and has been evaluated extensively to enable combined phosphorus and shortcut nitrogen removal (Roots et al. 2020; Jia et al. 2023). Partial denitrification by GAO has also been documented, albeit to a lesser degree than that of PAO. Denitrification kinetics in DPAO and DGAO are generally driven by available electron acceptors and specific strains of bacteria present in the biomass. For example, Ribera-Guardia et al. (2016) found that using as the primary electron acceptor led to the accumulation of N₂O in DPAO and DGAO enrichment reactors. Considering specific taxa, it was previously demonstrated that Ca. Accumulibacter were primarily reducers, while Ca. Competibacter were primarily NO and N₂O reducers in a reactor fed with as the electron acceptor (Wang et al. 2021).

In this study, the modified model improved predictions of bioP performance with moderate improvements to NO_x removal through changes to PAO and GAO activity and hydrolysis kinetics. These simple model modifications also resulted in higher modeled populations of GAO compared with the default model, which reflects recent microbial ecology findings of higher GAO populations in S2EBPR relative to conventional EBPR configurations (Farmer et al. 2023; Zhang et al. 2023). This study also identified putative DPAO and DGAO in the S2EBPR reactor. In future work, the relative contributions of DPAO and DGAO for NO_x transformation and accumulation of denitrification intermediates in S2EBPR should be investigated through modeling and direct measurement. Sumo and other commercial modeling packages typically do not include partial denitrification kinetics in default models and measurement techniques for denitrification intermediate NO and N₂O are not widely adopted in the field, leading to a gap between known bacterial phenotypes and performance predictions in process models.

The purpose of this study was to compare the performance of two models in forecasting NO_x and OP removal from a full-scale S2EBPR reactor based on the field data for the four operational phases. A default Sumo2 and modified Sumo2 model had similar NOx removal predictions, with the modified model having slightly better performance in predicting NO_x removal, with lower MSE values overall. However, the models showed significant differences in predicting effluent OP, with the modified Sumo2 model accurately predicting values during the carbon addition phase. Incorporating GAO competition helped with model calibration, though generally, the degree of GAO selection is not well explained by conventional metrics for assessing PAO and GAO competition. For example, the pH was near neutral throughout the process and both the temperature and COD:P ratio were not excessively high, yet there was a high GAO fraction through seasonal periods and phases with stable and unstable EBPR. This suggests that incorporating GAO into models when they occur in a substantial fraction is useful, but the prediction of GAO competition may still present a challenge in many instances.

The study also explored the microbial populations and denitrification pathways in a full-scale S2EBPR system using 16S rRNA and shotgun metagenomic sequencing. The results showed the presence of both PAO and GAO populations, with Ca. Competibacter is the dominant GAO. Evidence of both DPAO and DGAO with unevenly distributed denitrification reductases was found. While the modeling approach did not capture the bloom of Ca. Competibacter observed in Phase I, neither the process modeling approach nor marker gene microbiome profiling capture the complete phenotype of the organisms in situ. Furthermore, current process models are limited in their representation of denitrification as a single-step process. Accurate modeling of denitrification intermediates is crucial to quantifying and reducing N₂O emissions, and aligning microbial community function with process prediction. Overall, these findings provide a better understanding of the microbial ecology and process modeling in S2EBPR systems, which can aid in optimizing and improving their performance for efficient nutrient removal.

The authors thank Calumet Water Reclamation Plant operations and the Metropolitan Water Reclamation District of Greater Chicago for on-site testing facilities and assistance. The authors express their gratitude to Black & Veatch's Innovation Platform for their support. These results were presented at the 9th IWA Water Resource Recovery Modelling Seminar 2024, University of Notre Dame, USA (WRRmod2024), and the fruitful discussions are kindly acknowledged.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.