A mathematical model for a granular biofilm reactor for leachate treatment was validated by long-term measured data to investigate the mechanisms and drivers influencing biological nitrogen removal and microbial consortia dynamics. The proposed model, based on Activated Sludge Model (ASM1), included anaerobic ammonium oxidation (anammox), nitrifying and heterotrophic denitrifying bacteria which can attach and grow on granular activated carbon (GAC) particles. Two kinetic descriptions for the model were proposed: with and without soluble microbial products (SMP) and extracellular polymeric substance (EPS). The model accuracy was checked using recorded total inorganic nitrogen concentrations in the effluent and estimated relative abundance of active bacteria using quantitative fluorescence in-situ hybridization (qFISH). Results suggested that the model with EPS kinetics fits better for the relative abundance of anammox bacteria and nitrifying bacteria compared to the model without EPS. The model with EPS and SMP confirms that the growth and existence of heterotrophs in anammox biofilm systems slightly increased due to including the kinetics of SMP production in the model. During the one-year simulation period, the fractions of autotrophs and EPS in the biomass were almost stable but the fraction of heterotrophs decreased which is correlated with the reduction in nitrogen surface loading on the biofilm.
Mathematical modelling of microbial community dynamics and interactions in mixed consortia is of great significance not only in ecological simulation but also for activated sludge models in wastewater treatment (Klimenko et al. 2016; Arashiro et al. 2017). In biological nitrogen removal processes, simulation of microbial population can help to develop support tools to investigate the interactions between microorganisms. This step is also beneficial to improve calculation of substrate removal in bioreactors.
The anammox-based technologies are potential energy-positive alternatives to conventional nitrification-denitrification, especially for ammonia-rich wastewater treatment (Cao et al. 2016; Ma et al. 2016). Anammox bacteria belong to the Planctomycetes phylum, with six known genera (Azari et al. 2017a). The bacteria consume ammonium as their energy source (electron donor) and dissolved CO2 and HCO3− is used for cell biosynthesis (Suneethi et al. 2014; Ma et al. 2017). The anammox process has brought advantages for the wastewater treatment industry, since it does not require any extra organic carbon sources, up to 60% energy for aeration is saved and 90% less excess sludge is produced compared to conventional nitrification-denitrification (Azari et al. 2017a). Due to the slow growth of anammox bacteria, improvement of solids retention time (SRT) is required and therefore the application of anammox in biofilm reactors has been suggested as a possible approach (Lackner et al. 2014). Mathematical modelling of this type of anammox biofilm reactor could be an advantageous tool for the simulation, prediction and understanding of microbial community structure and biofilm characteristics. However, this is controversial due to present research gaps. First, current implications from mathematical modelling of biofilm reactors mainly focus on simulation of the substrates removal but little is known about modelling the active and inactive microbial groups regarding to dominance, abundance and biocoenosis. Also, in terms of dynamic modelling of extracellular polymeric substances (EPS) and soluble microbial products (SMP) with a focus on nitrogen-converting bacteria only few studies are available (Ni et al. 2012; Liu et al. 2016) while other studies have concluded that production of EPS components is linked to the formation of biofilms and granular sludge and the stability of bioreactor (Tan et al. 2017). In addition, model verification requires a systematic algorithm which has not been adequately addressed in current studies (Zhu et al. 2016). For example, based on the sensitivity analysis of previous works (Azari et al. 2017b), the model output results of such complex biofilm models are prone to parameter collinearity, numerical instability and minor changes in kinetic and stoichiometric parameters and sometimes initial parameters. Hence these values must be tuned and identified properly. Finally, among studies explaining models of anammox plants (Dapena-Mora et al. 2004; Ni et al. 2009, 2014; Liu et al. 2017), few checked the model against good long-term experimental data, not only for verification of nitrogen components and organic carbon, but also for calibration and validation of microbial community dynamics (Ni et al. 2012). However, low-quality long-term experimental data can hamper the verification of biokinetic models on a statistical basis (Kovárová-Kovar & Egli 1998). To overcome the claimed challenges, this study investigates the modelling of simultaneous partial nitrification, anammox and denitrification (SNAD) in a full-scale biofilm reactor assisted by granular activated carbon (GAC) particles under anoxic conditions.
The main aim of the work is to develop and evaluate a comprehensive and novel mechanistic model including SMP and EPS definitions to explain the relative abundance of independent microbial groups and to further describe the dynamics of microbial and non-microbial groups in the biofilm. Specific objectives include: (1) modelling of nitrogen removal and microbial community structure, (2) application of quantitative fluorescence in-situ hybridization (qFISH) to verify the simulated results of the microbial abundances, (3) comparing the model performance with and without EPS and SMP kinetics.
Operation of biofilm reactor
A full-scale anammox plant for municipal landfill leachate treatment is operated by LAMBDA GmbH (Herten, Germany) consisting of three main stages, including biological treatment with activated sludge, ultrafiltration and activated carbon biofilm process (Azari et al. 2017a) (Figure 1). The entire processes are very stable to sudden variations in flow and nitrogen concentration, and it can convert nearly 97% of the ammonium to nitrogen gas. In this study, the last part of the plant before the outflow containing GAC-assisted biofilm reactors (essential for final polishing of remaining amount of ammonium, nitrite and nitrate after the activated sludge process) was modelled. The anammox organisms grow on the surface of activated carbon with silica gel as the adsorbent medium. The GAC biofilm also comprises nitrifying and heterotrophic denitrifying bacteria (Ke et al. 2015). Daily concentrations of NH4-N, NO2-N and NO3-N were recorded photometrically using German standards (Hach, Germany). Total organic nitrogen (TON) is negligible in the influent and effluent of the biofilm reactor and the pH value is controlled between 7 and 7.4. Total inorganic nitrogen (TIN) was calculated as the sum of NH4-N, NO2-N and NO3-N.
Quantitative fluorescence in-situ hybridization (qFISH) and digital image analysis
A series of experiments was designed to develop a protocol adapted from Azari et al. (2017a, 2017b) for qFISH, in-situ microscopy and image analysis for quantification of relative abundance of dominant microorganisms (Figure A1 of Supplementary Materials, available with the online version of this paper). The biomass was frequently taken from the steady-state operation of the GAC biofilm reactor in five monthly timeframes. Three groups of bacteria were targeted, including anammox (AMX), ammonia-oxidizing bacteria (AOB) and nitrite-oxidizing bacteria (NOB). In total, five specific probes mixed with a general probe were chosen (Table 1). Probe Amx820 and Sca1309 target most anammox bacteria, Nso190 targets most AOB, NIT3 and Ntspa662 target most NOB, and EUB338 is a general probe to target most domain bacteria (more than 90%). 20 qFISH images from randomized field of observations (FOVs) were obtained by an epifluorescence microscope (Axio imager 2, Carl Zeiss AG, Jena, Germany). The relative abundances of AMX, AOB and NOB were estimated by calculating the percentage of their respective areas taken up by fluorescence targeted cells from the specific probes compared to the total area taken up by the EUB 338 general probe using ImageJ and daime software (Daims et al. 2006; Almstrand et al. 2014). The relative abundance using qFISH will be further used for model evaluation.
Model definition and evaluation
Activated Sludge Model No. 1 (ASM1) based on the biofilm modelling platform in AQUASIM 2.1 (EAWAG, Switzerland) was modified for enhanced nitrogen elimination including four main microbial groups i.e. AOB, NOB, AMX and heterotrophic denitrifiers (Table 2) according to the Monod kinetics and the revised stoichiometric matrix of the anammox process (Lotti et al. 2014). As previously discussed in the introduction, it is essential to establish a robust protocol for model evaluation. Therefore, this work formulated a protocol adapted from Zhu et al. (2016) for experiment design, data collection, model development, sensitivity analysis, calibration, validation and uncertainty analysis (Figure 2). The new contributions to this model evaluation algorithm are (1) adding the experimental design procedure for qFISH and (2) defining statistical criteria for model calibration and validation for TIN and the relative abundance of microbial consortia. An additional definition of SMP and EPS was added to the model structure to compare the model performance with and without EPS and SMP kinetics. This structure can model SMP production and uptake, EPS production and EPS hydrolysis rates according to Liu et al. (2016). SMPs were classified into two groups, utilization-associated products (UAP) from the natural product of bacterial growth and biomass-associated products (BAP) from the hydrolysis process of the active biomass (Laspidou & Rittmann 2002). Figure 3(a) and 3(b) illustrate how two model structures with and without EPS kinetics were developed. Since GAC has the potential to adsorb chemical oxygen demand (COD) and NH4-N, a simple adsorption isotherm equation (Kd approach) as a linear form of the Freundlich equation was also added to both models to define adsorption kinetics (Ranjbar & Jalali 2015). Bioprocess kinetic rates equations including those for EPS, UAP and BAP are presented in Table 3. The stoichiometric matrix and definition of biofilm matrix parameters are provided in Tables A1 and A2, respectively, of the Supplementary Materials (available online).
RESULTS AND DISCUSSION
Simulation of nitrogen removal
The general behaviour of the biofilm reactor's influent is characterized by significant fluctuations in NH4-N and COD with a COD:TN ratio of 4 to 20. The one-year performance of the anammox biofilm reactor was calibrated for two models using measured data of NH4-N, NO2-N and NO3-N. The regression analysis between the simulated and measured concentration of TIN and NH4-N implies a higher R2 value for the model with EPS and SMP. For the model with EPS and SMP, R2 values for TIN and NH4-N are equal to 0.65 and 0.66 respectively while R2 values of 0.64 and 0.61 were calculated for the model without EPS kinetics. Calibrated values for model parameters and regression analysis graphs for TIN and NH4-N are provided in Tables A3 and A4 and Figure A4 of the Supplementary Materials (available with the online version of this paper). After the model calibration, validation results of both models for TIN over 120 days of simulation showed that the model with EPS and without EPS predict the fate of remaining nitrogen in the effluent with a similar accuracy. Although the model with EPS undermined some measured data in the middle of validation period, the model without EPS overestimated some data points in the beginning of validation period (Figure A5 of Supplementary Materials, available online). Nevertheless, after regression analysis both models gave the identical R2 values of 0.65 for the validation period.
Model evaluation for relative abundance of microbial consortia
The relative abundances of AMX, AOB and NOB were visualized (see Figure A2 of the Supplementary Materials for anammox bacteria and Figure A3 for AOB and NOB, available online) and estimated by qFISH. As shown in Figure 4, the averaged relative abundance of anammox bacteria as the dominant species to the total domain bacteria over one year was 52.2% ± 5.4 of total microbial groups while the simulated results give the averaged relative abundance of 45.8% for the model without EPS (Figure 5(a)) and 53.7% for the model with EPS (Figure 5(b)). For AOBs, the average estimated relative abundance of bacteria to the total domain bacteria throughout the year was 5.7% of total microbial groups while the modelled results give a relative abundance of 15.8% for the model without EPS (Figure 5(a)) and 8.5% for the model with EPS (Figure 5(b)). The NOB was the least dominant bacterial group, at only 2.6 ± 3.2% from qFISH results (Figure 4). Simulation of NOB relative abundance was in good agreement with observations: the model with and without EPS predicted 2.6 ± 1.8% and 3.3 ± 2.8% respectively (Figure 5(a) and 5(b)). Comparing the simulated relative abundance with the observed values for AMX, AOB and NOB, it could be seen that statistically the model with EPS fits better than model without EPS, in which averaged relative abundance was closer to the measured data and the trends and variations were better predicted (Figure 5(b)). Moreover, throughout the simulation period, the model without EPS overestimated the AOB relative abundance while the relative abundance of heterotrophs decreased continuously to less than 10%. This could have been due to the underestimate of heterotrophic growth rates due to neglecting heterotrophic synthesis from the utilization of UAP and BAP as extra organic carbon sources (Figure 3(b)).
With regard to the kinetics of BAP and UAP formation and uptake, UAP and BAP have distinct kinetics parameters, with UAP more readily biodegradable for heterotrophs. Considering the set parameters in our model for UAP and BAP, the specific growth rates, and the yield coefficient of readily soluble COD (Ss) is higher than UAP and for UAP is higher than BAP. On the other hand, the affinity constant for BAP is higher than UAP and for UAP is higher than readily soluble COD (Ss) (see Table A4, available online). These reasons will lead to the favourable consumption of UAP than BAP by heterotrophs. This finding follows earlier theoretical hypotheses to the most recent model-based findings claiming that UAP is more important than biomass-associated SMP (BAP) (Namkung & Rittmann 1986; Liu et al. 2016).
Correlation analysis of nitrogen removal efficiency and biomass fractionations
The model predicts total biofilm fractionations in terms of percentage of active (microbial) groups e.g. heterotrophic denitrifiers, anammox, AOB, NOB and other microbial groups not involved in nitrogen removal, as well as inactive (non-microbial) groups e.g. EPS, inert, and slowly biodegradable matter (Figure 6). To link the fractions of microbial and non-microbial groups in the biofilm with system performance, the observed and modelled nitrogen removal efficiency (NRE) (average of 10 days' values) are sketched. Inert matter contributed 37.5–44% of biomass and EPS abundance varied between more than 6% to 10%.
The simulation results indicated that a period of 15 days for both models is the time necessary to reach a steady-state. During the steady-state period, the fraction of autotrophs and EPS was stable, but the fraction of heterotrophs decreased and the fraction of inert matter rose while the NRE fluctuated between 37% and 82% (Figure 6(a) and 6(b), black connected line). The highest removal efficiency of the system was found between days 90 and 250, which can be predicted by models as well. Nevertheless, no specific link can be found to correlate the system performance in terms of efficiency and stability with fractions of microbial and non-microbial groups in the biofilm.
Correlation analysis of the nitrogen surface loading and biomass fractionations
The nitrogen surface loading (NSL) was defined as an average of TIN loaded on a square metre of biofilm active surface. In the first 150 days of the simulation, the NSL was higher, estimated to be at 59–105 kg N m−2 d−1 for the model with EPS and 44.1–85.0 kg N m−2 d−1 for the model without EPS (Figure 6, yellow connected line). But after day 150, NSL declines to less than 20 kg N m−2 d−1 due to an increase in biofilm thickness for both models. Consequently, the concentration of heterotrophs, autotrophs (AMX, AOB and NOB) and EPS decreased slightly when the biofilm thickness increased from 0.62 to 1.2 mm for the model without EPS and from 0.57 to 1.01 mm for the model with EPS. Nevertheless, over the one-year simulation period, the fractions of autotrophs and EPS were approximately stable, but the fraction of heterotrophs decreased significantly which is correlated with reduction in NSL of the biofilm. The fraction of inert matter in the biofilm correlates inversely with NSL values.
Our findings are compatible with a study from Ni et al. (2012), which demonstrated how the NSL impacts on the biomass fractions of anammox biofilm and it indicated that during a long-term simulation, the NSL increase will lead to an increase in the fractions of heterotrophs, EPS, and anammox and a reduction in the fraction of inert matter. After a while, these fractions will level off even though the nitrogen loading rate is increasing (Ni et al. 2012).
This kinetic study proposed a novel model structure for biological nitrogen removal by including EPS, UAP and BAP components and new process rates, which made ASM1 more complete. The model is verified for GAC biofilm reactors and results elucidated that the model can simulate inorganic nitrogen components in the bulk liquid and the relative abundance of targeted bacterial groups with sufficient accuracy. The key findings are as follows:
The model was improved after the addition of EPS and SMP kinetics.
The model proved that the microbial and non-microbial composition of a stable biofilm reactor is not correlated with the NRE, but it might be more influenced by any change in operation conditions such as flow rate, hydraulic retention time, dissolved oxygen, temperature, salinity and pH level.
The biomass fractionations of microbial and non-microbial groups predicted by the biofilm model can be correlated with a factor called NSL.
The model without EPS and SMP kinetics overestimated the relative abundance of AOB.
The model without EPS underestimated heterotrophic growth compared with the model with EPS, because the model with EPS and SMP incorporates UAP and BAP for the support of heterotrophic growth.
UAP is more readily biodegradable for heterotrophs than BAP.
The real-time production of EPS and SMP can be predicted as well.
This work suggests that the growth and existence of heterotrophs in anammox biofilm systems increases due to incorporating the production of SMP so that SMP supports the growth of heterotrophic denitrifiers. Hence, application of the model with EPS is an advantage when we want a better insight on microbial consortia. With a reliable prediction of SMP release, we can better calculate and control soluble COD and inorganic nitrogen concentrations in the effluent of wastewater treatment plants.
We acknowledge the support of the German Academic Exchange Service DAAD and remarkable supports from Zentrum für Wasser- und Umweltforschung (ZWU) at the University of Duisburg-Essen. We also appreciate the collaboration with the AGR Group and LAMBDA Gesellschaft für Gastechnik mbH for their technical assistance. We are also grateful to the Biofilm Centre, Faculty of Chemistry, University of Duisburg-Essen and Dr Frösler and Dr. rer. nat. Jost Wingender for providing advanced microscopy facilities.
Contributed equally to this work.