Modeling nitrogen removal in membrane aerated biofilm reactors: the role of nitritation, denitritation, and anammox nitrogen removal

1.   Introduction

Nitrogen removal from wastewater is a critical process in environmental engineering, particularly due to the adverse effects of nitrogenous compounds on aquatic ecosystems and human health. Traditional methods, such as nitrification and denitrification, have been widely employed to mitigate nitrogen pollution. However, these processes often require significant energy inputs and can lead to the production of greenhouse gases, such as nitrous oxide ($N_2O$). In recent years, the anammox (anaerobic ammonium oxidation) process has gained attention as a more sustainable alternative, allowing for the direct conversion of ammonium and nitrite into nitrogen gas under anaerobic conditions, thereby reducing the overall energy demand and greenhouse gas emissions associated with nitrogen removal ^[1,2,3,4,5]. Although this innovative biological strategy holds significant technological and economic promise, its application is challenging primarily due to the slow growth rate of anammox bacteria and their sensitivity to various external factors. Combining nitritation and anammox in biofilm or granular reactors has been suggested as a strategy to overcome these challenges ^[6,7].

Membrane aerated biofilm reactors (MABRs) have emerged as a new type of biological wastewater treatment process for enhancing nitrogen removal efficiency. They utilize a gas-permeable membrane to supply oxygen directly to the biofilm from the membrane substratum against the concentration gradients of carbon and nitrogen substrate within the biofilm. The counter-diffusion mass transfer characteristics allow MABRs to establish an aerobic zone adjacent to the membrane substratum and an anoxic/aerobic region close to the bulk liquid/biofilm interface to achieve simultaneous nitrification and nitrogen removal within the biofilm. The typical nitrogen removal pathways include conventional nitrification and denitrification (CND) ($NH_4^+-N\rightarrow NO_3^–N \longrightarrow N_2$), the shortcut nitrogen removal or more concisely the nitritation-denitritation (ND) ($NH_4^+-N\rightarrow NO_2^–N\rightarrow N_2$), and partial nitritation-anammox (PN/A) ($NH_4^+-N +NO_2^–N\rightarrow N_2$). Compared to CND, ND can save 25% of oxygen for nitrification and 40% of organic carbon substrate for denitrification, while PN/A can reduce oxygen consumption by 62% and organic carbon substrate consumtion by 100% ^[8,9]. Studies have shown that MABRs can achieve 30% to 97% nitrogen removal for influent total nitrogen concentrations ranging from $15 \; gm^{-3}$ to $2500 \; gm^{-3}$, Chemical Oxygen Demand (COD) (the amount of oxygen required to chemically oxidize organic and inorganic substances in a water sample) concentrations from $7 \; gm^{-3}$ to $5620 \; gm^{-3}$, MABR aeration pressure from $0.3$ to $152 \; kPa$ via pathways including CND, ND, and PN/A ^[8]. Furthermore, a variety of studies have demonstrated that with proper control of the oxygen supply and carbon to nitrogen ratio, MABRs can achieve effective nitrogen removal via ND and anammox ^{[10,11,12,13]}. While the counter-diffusion characteristics of MABRs have been shown to effectively support the growth of bacteria involved in nitrogen removal, controlling the distribution of oxygen and substrates within the biofilm in practical applications remains challenging. This control is crucial for stimulating ecological niches that promote ND and PN/A nitrogen removal functions within the biofilm. For instance, surpassing the growth of nitrite-oxidizing bacteria (NOB) is necessary to foster efficient ND and PN/A nitrogen removal. On the other hand, the enrichment of NOB, influenced by substrate and oxygen supply as well as biofilm thickness, can shift nitrogen removal pathways from ND and PN/A to CND, reducing overall nitrogen removal efficiency. Therefore, to optimize energy-efficient ND and PN/A nitrogen removal with MABRs, a comprehensive understanding of the interactions between oxygen and substrate supply, the distribution of oxygen and nitrogenous intermediates within the biofilm, and the underlying nitrogen removal mechanisms is essential.

Since it is challenging to experimentally characterize the distributions of oxygen, substrates, microbial populations, and their physiological interactions within the biofilm, along with biofilm thickness, mathematical modeling is pivotal for studying the potential of the anammox process and understanding the mechanisms responsible for maintaining high anammox activity. Existing theoretical studies frequently utilize one-dimensional biofilm models within the framework of the 1D Wanner–Gujer model ^[14], which operates under the assumption that biomass volume fractions sum to unity (a reasonable assumption for well-developed, homogeneously stratified biofilms with co-diffusional substrate supply), or they employ kinetic-based models ^{[11,13,15,16,17,18,19]}. Among these works, Li et al. ^[11] developed a one-dimensional multispecies model that included anammox bacteria, AOB, NOB, and aerobic heterotrophic bacteria (HB) to study the effects of temperature, aeration pressure, and COD/N ratio on nitrogen removal mechanisms in MABRs. Their results demonstrated that increasing the C/N ratio and controlling aeration pressure could balance the competition between AOB, NOB, and aerobic heterotrophs for oxygen, ultimately enhancing ND nitrogen removal. Hao et al. ^[16] developed a model to evaluate the influence of variable temperature, ammonium loading rates, and dissolved oxygen concentrations on a full autotrophic N-removal process in an aerated biofilm reactor. Their results indicated that temperature and inflow variations negatively influence the N-removal efficiency due to limited activity of anammox bacteria and an unbalance between applied ammonium surface load and required oxygen concentration, respectively. Moreover, they showed that at a constant temperature and a defined ammonium surface load, there is an optimal biofilm thickness to achieve the maximal N-removal efficiency. Recently, Liu et al. ^[18] proposed a one-dimensional model that considered the impact of anammox and nitrite/nitrate-dependent anaerobic methane oxidation ($n$-DAMO) on nitrogen removal in MABRs. Their findings revealed that the distribution of anammox and $n$-DAMO within the biofilm, along with controlled methane supply under regulated partial pressure, could enhance nitrogen removal efficiency.

The assumption in the 1D Wanner-Gujer model is often valid, particularly for well-developed biofilms in nutrient-rich environments. Studies have shown that under these conditions, biofilm models without this assumption still predict biomass volume fractions close to unity ^[20]. Nevertheless, it is well-known that the biomass in biofilms does not uniformly reach the same volume fraction throughout, particularly in cases of cell lysis or nutrient limitation, which lead to the development of internal regions with significantly reduced cell density. Models that assume the total biomass volume fraction is constant and sums to unity at all locations are inherently unable to capture this phenomenon and will instead predict biofilm shrinkage or contraction ^[21,22]. Moreover, studies have revealed that extracellular polymeric substances (EPS) which play a crucial role in determining biofilm density are not uniformly distributed across the biofilm thickness and, additionally, the chemical composition of EPS can be influenced by the oxygen loading rates ^[23]. Furthermore, imposing the constraint that the sum of all biomass volume fractions must equal unity significantly restricts the model's ability to accurately predict the spatial distribution and interactions of different microbial species within the biofilm ^[24]. The limitation of kinetic-based models lies in their insufficient spatial resolution, which hinders their ability to accurately describe the heterogeneity of biofilms. Therefore, implementing a high-dimensional, comprehensive mathematical model, without the limitations of the Wanner-Gujer model, is essential to capture biofilm heterogeneity driven by counter-diffusion mechanisms. This approach will provide insights into the fundamental processes governing the behavior of MABRs and facilitate the investigation of anammox's role in nitrogen removal under diverse environmental conditions. Considering these factors and addressing the existing gap in the literature, where, to the best of our knowledge, no comprehensive study has been conducted on nitrogen removal and the potential of anammox, this study employs a two-dimensional partial differential equation (PDE) model to investigate the influence of environmental conditions on nitrogen removal in MABRs, with particular emphasis on anammox processes.

2.   Mathematical model

For a comprehensive analysis of nitrogen removal, particularly incorporating the anammox process in MABRs, we extend the model proposed in ^[25] to account for the roles of anammox and NOB in the nitrogen removal process. The model is governed by a system of degenerate reaction-diffusion equations over the computational domain $\Omega\subset \mathbb{R}^2$, which describes biofilm growth and substrate dynamics. To formulate the model, the multi-component biofilm is assumed to consist of inert bacteria (In), ammonium oxidizing bacteria (AOB), aerobic heterotrophic bacteria (HB), anaerobic heterotrophic bacteria (AHB), NOB, and anammox. We also assume that oxygen is supplied through the boundary of the domain on which the biofilm forms, mimicking a membrane there, while the other substrates are introduced from the opposing boundary of the domain indicating that a counter-diffusion mechanism is involved. Since flow field calculations are often significantly more computationally expensive than simulating the actual biofilm growth process, we neglect the hydrodynamic effects and the convective contribution to substrate transport in the aqueous phase. Instead, we consider a hydrostatic environment where nutrients are supplied or removed through the membrane and the concentration boundary layer (introduced in the boundary conditions, Eq (2.2)). Hence, we assume that nutrients are transported to the biofilm from the aqueous phase by a diffusion gradient, more specifically by Fickian diffusion. The dependent variables include the volume fractions of bacterial species: $I$ (In), $B_a$ (AOB), $B_{h1}$ (HB), $B_{h2}$ (AHB utilizing nitrite), $B_{h3}$ (AHB utilizing nitrate), $B_n$ (NOB), and $B_{an}$ (anammox). Additionally, the model tracks the concentrations of dissolved substrates, including ammonium nitrogen $N_1$, nitrite nitrogen $N_2$, nitrate nitrogen $N_3$, oxygen $O$, and organic substrate acetate $A$. The independent variables are time and spatial location which are denoted by $t \geq 0$ and $(x, y) \in \Omega$ respectively. A summary of biofilm components, dissolved substrates, and the corresponding variables is given in Table 1.

Using the above assumptions and notations, the differential mass balances for the dependent variables are obtained as:

where $M$ is the total biomass fraction, i.e.,

The two-dimensional computational domain, denoted as $\Omega$, is composed of two regions:

which characterizes the biofilm with positive biomass density, and

representing the aqueous phase (bulk liquid, channels, and biofilm pores without biomass), c.f. Figure 1. These regions are separated by the biofilm-liquid interface $\Gamma(t) = \partial \bar{\Omega}_1(t) \cap \partial \bar{\Omega}_2(t)$. $\Omega_1$ may consist of several disconnected sub-regions, each representing a distinct biofilm colony, which can merge as they grow. Hence, the regions $\Omega_1$ and $\Omega_2$ are not necessarily connected and $\Gamma$ can consist of several disjoint segments, c.f. Figure 1b. The time-dependence of these regions arises from biofilm structural changes due to growth. Following the previous studies ^{[20,26,27,28,29]}, etc. The diffusion coefficient for all biomass fractions is defined as $D(M) = \delta \frac{M^{4}}{(1 - M)^{4}}\; m^2d^{-1}$, where $\delta\; m^2d^{-1}$ is the biomass motility coefficient, which is positive but significantly smaller than the diffusion coefficients of dissolved substrates in the liquid phase. The function $D(M)$ exhibits two non-linear effects: (ⅰ) a power-law degeneracy similar to the porous medium equation, such that $D(M)$ approaches zero as $M \approx 0$, and (ⅱ) a super-diffusion singularity as $M$ approaches unity. The porous medium degeneracy term, $M^{4}$, ensures that biofilm spreading is minimal when biomass density is low, i.e., $0 < M \ll 1$, and is responsible for the formation of a sharp interface between the biofilm and the surrounding liquid. This guarantees that initial conditions with compact support lead to solutions with compact support. The second effect, (ⅱ), enforces an upper bound on the solution, ensuring that $M$ remains below unity, as demonstrated in ^[30]. This is counteracted by the degeneracy as $M = 0$ at the interface. Consequently, $M$ becomes concentrated within the biofilm region and approaches its maximum value 1. Hence, the interaction between these non-linear diffusion effects and the growth term is essential for accurately describing the spatial spread of biomass ^[31]. The diffusion coefficients for substrates are also dependent on biomass density, however in a less critical manner. These are expressed as

where $D_{(n_1, n_2, n_3, o, a)}^{0} > 0$ represents the diffusion coefficient in water, and $\rho_{(n_1, n_2, n_3, o, a)}$, with $0 < \rho_{(n_1, n_2, n_3, o, a)} \leq 1$, denotes the biofilm-to-water diffusivity ratio. As a result, these diffusion coefficients are bounded from below and above by finite values within one order of magnitude.

The reactions included in our model are:

● $R1$, $R2$: These terms represent biomass loss occurring through two primary mechanisms: lysis and endogenous respiration. Lysis is characterized as a first-order process with a constant loss rate denoted by $d_{(a, h1, h2, h3, n, an)}$, while endogenous respiration follows Monod kinetics with maximum oxidation rates $r_{(a, h, n)}$ and half-saturation constants $\kappa_{a_{o}}$, $\kappa_{h_{o}}$, and $\kappa_{n_{o}}$ for autotrophic, heterotrophic, and nitrite oxidizing bacteria, respectively. Note that shear-induced deformation and detachment have been excluded in this model. Although the concentration boundary layer, which will be incorporated into the boundary conditions for substrates, is qualitatively related to the bulk flow velocity, it will be set to a value that does not simulate rapid flow in subsequent simulations. Therefore, together with the size of the computational domain used in this study, it is reasonable to neglect shear induced deformation and detachment and account for biomass loss only through lysis and endogenous respiration ^[32].

● $R3$: Growth of ammonium oxidizing bacteria ($B_a$) which is regulated by ammonium and oxygen and is characterized by dual Monod kinetics with the maximum growth rate $\mu_a$ and half-saturation constants $\kappa_{a_{o}}$ for oxygen and $\kappa_{n1}$ for ammonium.

● $R4$: Growth of aerobic heterotrophic bacteria ($B_{h1}$) controlled by availability of oxygen and acetate. This is described by dual Monod kinetics with the maximum growth rate $\mu_h$ and half-saturation constants $\kappa_{h_{o}}$ and $\kappa_{a}$ for oxygen and acetate respectively.

● $R5$, $R6$: Growth of anaerobic heterotrophic bacteria ($B_{h2}$, $B_{h3}$) utilizing nitrite and nitrate in the presence of acetate under anoxic condition. These terms are described by oxygen inhibition and acetate, nitrite, and nitrate dependent Monod kinetics with $\eta_d\mu_h$ as the maximum growth rate and $\kappa_{h_{o}}$, $\kappa_{a}$, $\kappa_{h_{n2}}$, and $\kappa_{h_{n3}}$ as oxygen, acetate, nitrite, and nitrate half saturation constants respectively.

● $R7$: Growth of nitrite oxidizing bacteria ($B_n$) regulated by oxygen and nitrite and characterized by dual Monod kinetics with the maximum growth rate $\mu_n$ and half saturation constants $\kappa_{n_{o}}$ for oxygen and $\kappa_{n_{n2}}$ for nitrite.

● $R8$: Growth of anammox that occurs under anoxic condition and in the presence of ammonium and nitrite and is described by oxygen inhibition and ammonium and nitrite dependent Monod kinetics with $\mu_{an}$ as the maximum growth rate and $\kappa_{an_{o}}$, $\kappa_{an_{n1}}$ and $\kappa_{an_{n2}}$ as half saturation constants for oxygen, ammonium and nitrite respectively.

● $R9-1$: Ammonium consumption by AOB and conversion into nitrite with the yield coefficient $Y_a$.

● $R9-2$: Degradation of ammonium via anammox with the yield coefficient $Y_{an}$.

● $R10-1$, $R10-2$: Degradation of nitrite and its conversion into nitrate through NOB and anammox with the yield coefficients $Y_n$ and $Y_{an}$ respectively.

● $R10-3$: Shortcut nitrogen removal (ND) with the yield coefficient $Y_{h2}$.

● $R11$: CND with the yield coefficient $Y_{h3}$.

● $R12$, $R13$: Oxygen and acetate consumption.

To study the mathematical model (2.1) we restrict ourselves to the rectangular domain $\Omega = [0, L]\times[0, H]$. The substratum, i.e., the membrane surface on which biofilm colonies form and through which oxygen is supplied is the bottom boundary, $y = 0$. We assume that the substratum is impermeable to biomass and dissolved substrates ($N_1$, $N_2$, $N_3$, $A$) cannot diffuse out of the domain through this segment. For the remaining three sides of the domain, $y = H$, $x = 0, \; L$, we use the open boundaries conditions discussed in ^[33]. At the lateral boundaries, $x = 0$ and $x = L$, we assume a symmetry boundary condition for biomass species and substrates, which allows us to view the domain as a part of a continuously repeating larger domain. The condition posed at these boundaries is a homogeneous Neumann condition. To mimic the counter-diffusion substrate delivery, non-homogeneous Robin conditions are specified for ammonium and acetate at the top boundary while the homogeneous Neumann condition is imposed at the bottom boundary to reflect a hydrophobic membrane which does not permit these substances to diffuse through the membrane. The Robin conditions at the top boundary quantify mass transfer into the system, c.f. ^[33]. For oxygen, the boundary conditions at the top and bottom are non-homogeneous Robin conditions. Oxygen is added through the bottom and removed through the top boundary. However, we assume that it is not completely washed out at the top segment of the domain. Additionally, it is assumed that nitrite and nitrate are washed out through the top boundary, hence a homogeneous Neumann boundary condition is imposed. Biomass is assumed to be unable to exit the domain through the top and bottom boundaries, thus homogeneous Neumann conditions are defined, as in ^[33]. These conditions introduce the counter-diffusion effect on which membrane aerated biofilm reactors are based. The boundary conditions on domain $\Omega = [0, L]\times[0, H]$ are defined as:

Here, $\lambda_{(\cdot)}\; mm$ represents the concentration boundary layer thickness, which accounts for the influence of bulk hydrodynamics on mass transfer into and out of the simulation domain. Smaller values of $\lambda$ indicate rapid bulk flow, whereas low bulk flow velocity leads to a thicker concentration boundary layer ^[33]. Typically, the concentration boundary layer thickness varies between $y = 0$ and $y = H$. However, for the sake of simplicity, we assume that all substrates share the same concentration boundary layer thickness at both the top and bottom boundaries, denoted as $\lambda_{o} = \lambda_{n1} = \lambda_{n2} = \lambda_{n3} = \lambda_{a} = \lambda$.

3.   Numerical treatment

For the numerical simulation, we first non-dimensionalize the partial differential equation (PDE) system described in Eq (2.1). We use the height of the computational domain and the maximum growth rate of AOB as scaling factors for space and time, respectively, defined as $\tilde{x} = x/H$ and $\tilde{t} = t\mu_{a}$. The substrate concentrations are non-dimensionalized using the following scaling factors:

Note that the biomass volume fractions $I$, $B_a$, $B_{h1}$, $B_{h2}$, $B_{h3}$, $B_{n}$, and $B_{an}$ have already been defined as dimensionless variables relative to the maximum cell density $M_{\infty}$. Subsequently, we introduce a uniform grid for the rectangular domain $[0, 1]\times [0, H/L]$ to discretize the dimensionless PDE model by cell centered finite difference-based finite volume scheme. For time integration of the resulting semi-discrete model, we employ the time-adaptive, error-controlled embedded Rosenbrock-Wanner method ROS3PRL, as introduced by Rang ^[34]. The method requires in each time-step the solution of several sparse, large, nonsymmetirc linear systems, for which we use the stabilized bi-conjugate gradient method. The numerical method has been implemented in Fortran 95. OpenMP was used for the parallelization of selected computationally expensive tasks, such as the linear solver and the formation of the Jacobian matrix within the Rosenbrock-Wanner method, and for the evaluation of non-linear reaction terms. The method's ability to proceed with efficient time-steps for interface moving problems with similar nonlinearities as those in model (2.1) has been extensively studied in ^[22,29,35] and we refer the reader to these references for detailed description of numerical challenges and the applied numerical algorithm.

Remark 3.1. Note that the model represented in Eq (2.1) is a two-dimensional PDE model capable of capturing biofilm heterogeneity arising from an initially sparse inoculation (Figure 1b). Previous studies ^{[25,26,36,37]} have demonstrated the advantages of using high-dimensional models to accurately depict biofilm development from a realistic initial inoculum and to characterize internal biofilm heterogeneity in MABRs. However, as the overall performance of MABRs is generally not highly sensitive to initial inoculation or internal biofilm heterogeneity ^[26,36], we simplified our approach in this study by assuming that biofilm layers homogeneously cover the substratum. This scenario can be effectively captured by a one-dimensional reduced model derived from the proposed PDE system while avoiding the limitations of the Wanner-Gujer model. Although this study primarily focuses on the impact of operational and environmental conditions on nitrogen removal pathways, it lays a strong foundation for future model development that accounts for internal heterogeneity and spatial variations within the biofilm. This will be essential for addressing questions where spatial heterogeneity plays a critical role.

The parameters and their values used in computer simulations are summarized in Table 2.

3.1.   Simulation set-up

It has been shown in ^[26,36] that while variations in the initial biofilm inoculation and internal heterogeneity may influence its internal structure, these factors generally do not have a substantial impact on the overall performance of MABRs. Therefore, in this study, we employ a flat layer coverage of the substratum to investigate nitrogen removal through different pathways as it is a more relevant case from a biofilm reactor performance perspective. Although biofilm structure is highly sensitive to environmental conditions, making the formation of a specific structure challenging, the counter-diffusion mechanism inherent in MABRs allows for layered configurations to be utilized in experimental studies ^[41]. Consequently, we assume that the substratum is initially covered by a flat layer of $B_a$, $B_{h1}$, $B_{h2}$, $B_{h3}$, $B_{n}$, and $B_{an}$, arranged from the membrane surface (bottom segment of the domain) to the top, respectively. The initial thickness of each bacterial layer is set to $10\; \mu m$. Furthermore, since varying the sequence of biofilm layers has minimal impact on nitrogen removal ^[26], we adopted the AOB-HB-AHB-NOB-anammox layer configuration to enhance the clarity of AHB distribution within the biofilm and to better observe the development of biofilm layers over time. For the numerical simulations, the bulk concentration of oxygen, $O_{\infty}$, is varied from $0.25\; gm^{-3}$ to $5\; gm^{-3}$, and the bulk concentrations of ammonium nitrogen and acetate are selected within the ranges $N_{1\infty} = 5-400\; gm^{-3}$ and $A_{\infty} = 6-100\; gm^{-3}$, respectively. The bulk concentration of oxygen is the dissolved oxygen concentration at the interface of the biofilm and membrane substratum. Since a low oxygen supply is essential to stimulate ND and PN/A nitrogen removal, this study focuses on the relatively low $O_{\infty}$ conditions ranging from $0.25\; gm^{-3}$ to $5\; gm^{-3}$, which corresponds to the oxygen partial pressures of $0.57 \; kPa$ and $11.3 \; kPa$ for the Henry's law constant of $0.001378 \; mol/L/atm$. The $N_{1\infty}$ and $A_{\infty}$ are the concentrations of ammonium nitrogen and acetate in the bulk liquid. For a complete mix MABR, these concentrations represent the ones in the reactor effluent. The combinations of $N_{1\infty}$ and $A_{\infty}$ used in the simulations include $A_{\infty} = 6\; gm^{-3}$ for $N_{1\infty} = 5, \; 10, \; 15, \; 100, \; 200, \; 400\; gm^{-3}$; and $N_{1\infty} = 15, \; 100\; gm^{-3}$ for $A_{\infty} = 6, \; 30, \; 100\; gm^{-3}$. These combinations are selected to represent the secondary effluent from biological wastewater treatment and the return liquid from the anaerobic digestate and the municipal wastewater. The secondary effluent and the return liquid normally have low readily biodegradable COD and moderate to high ammonium concentrations, while the municipal wastewater has a moderate readily biodegradable COD and ammonium content. The ammonium nitrogen concentration in typical municipal wastewater and in the effluent from secondary treatment processes that do not include nitrification usually ranges from $10\; gm^{-3}$ to $50 \; gm^{-3}$. In contrast, the return liquid from anaerobic digestate treating wastewater sludge can contain ammonium nitrogen concentrations between $500\; gm^{-3}$ and $1500 \; gm^{-3}$ ^[42]. The concentrations of ammonium and acetate in the effluent discharged from the MABR reactor will depend on the substrate loading rates and the removal efficiency of the reactor. The thickness of the biofilm is a critical factor determining the microbial composition, oxygen transfer efficiency and biological functions of a biofilm in MBfR ^[43]. An overgrown biofilm will increase the resistance of substrate and oxygen diffusion, which could be deleterious to the micro-environment needed for the growth of bacteria within the biofilm. On the other hand, studies have shown that achieving a biofilm thickness between $500\; \mu m$ to $1000\; \mu m$ can facilitate simultaneous nitrification and denitrification ^[12,44,45]. With the biofilm thickness in such a range, complete oxygen penetration through the biofilm can be avoided so that an anoxic zone and an aerobic zone can co-exist within the biofilm to promote the growth of heterotrophic denitrifiers and autotrophic nitrifiers. Thus, in the subsequent numerical studies (except in part of Section 4.5), simulations are terminated when the biofilm occupies 30% of the domain, i.e., when the biofilm thickness reaches $600\; \mu m$.

3.2.   Postprocessing and quantities of interest

For better interpretation of simulation results, we will quantify the influence of environmental conditions on biofilm development and the contributions of various nitrogen removal pathways through the definition of the following lumped output variables:

● Total nitrogen removal ($NR_T$):

where

● Nitrogen removal via PN/A ($NR_1$):

● Nitrogen removal via ND ($NR_2$):

● Nitrogen removal via CND ($NR_3$):

● Nitrification by AOB:

● Conversion of nitrite nitrogen to nitrate nitrogen via NOB

● Total acetate consumption by $B_{h1}$, $B_{h2}$, and $B_{h3}$

● Total oxygen consumption by $B_a$, $B_{h1}$, $B_n$, and endogenous respiration

● Total amount of biofilm species

4.   Results

4.1.   The effect of ammonium and oxygen concentrations on biomass distribution within the biofilm

To examine the influence of environmental factors, particularly the bulk concentrations of oxygen and ammonium, on the growth dynamics of different biomass species, we implement a systematic variation in the values of $N_{1\infty}$ and $O_{\infty}$. The ammonium bulk concentration is set to $N_{1\infty} = 5, \; 10, \; 15, \; 100, \; 200, \; 400\; gm^{-3}$ and oxygen bulk concentration is adjusted to $O_{\infty} = 0.25, \; 0.5, \; 1, \; 2.5, \; 5\; gm^{-3}$ for each ammonium concentration. The acetate bulk concentration is maintained at a constant $A_{\infty} = 6\; gm^{-3}$ and the system is analyzed once the biofilm reaches a thickness of $600\; \mu m$. The results, depicted in Figure 2, indicate that at the lowest oxygen concentration ($O_{\infty} = 0.25\; gm^{-3}$), biofilm growth is limited and it predominantly consists of heterotrophic bacteria ($B_{h1}$), which have a competitive advantage under oxygen-limited conditions over AOB ($B_{a}$). Although increasing ammonium concentration leads to an increase in AOB, their contribution remains minimal relative to other species, with $B_{h1}$ continuing to be the dominant component. Increasing the oxygen concentration from $0.25\; gm^{-3}$ to $0.5\; gm^{-3}$ enhances biofilm growth and AOB production, however the rate of AOB production is still insignificant across all ammonium concentrations, and the biofilm remains largely occupied by other species mainly $B_{h1}$. Further elevation of oxygen to $1\; gm^{-3}$ facilitates formation of more AOB, which, through oxygen depletion and nitrite production, stimulates the growth of $B_{h2}$ and $B_{h3}$. Increasing oxygen to $2.5\; gm^{-3}$ promotes AOB growth, leading to greater nitrite and nitrate production and subsequently supporting the growth of $B_{h2}$ and $B_{h3}$, which utilize acetate and thereby reduce the growth rate of $B_{h1}$. Under ammonium-limited condition ($N_{1\infty} = 5\; gm^{-3}$), further increase in oxygen concentration from $2.5\; gm^{-3}$ results in an overall decrease in total biomass. At higher oxygen concentration ($O_{\infty} = 5\; gm^{-3}$), the activity of $B_a$ increases, resulting in greater $N_2$ production. However, since the growth of $B_a$ is limited due to ammonium shortage, high oxygen inhibits the formation of an anoxic zone and reduces the activity of $B_{h2}$, which competes with $B_{n}$ for $N_2$. This leads to increased formation of $B_{n}$ and, consequently, more $N_3$ production. Nevertheless, under high oxygen concentration, the growth rate of $B_{h3}$, which occurs in the anoxic zone, decreases. This results in the formation of less $B_{h3}$ compared to the lower oxygen condition ($O_{\infty} = 2.5\; gm^{-3}$), and a reduction in the total biomass. Oxygen availability for $B_{h1}$ increases at $O_{\infty} = 5\; gm^{-3}$ which competes with $B_{h2}$ and $B_{h3}$ for acetate. This leads to $B_{h1}$ and $B_{h3}$ becoming the dominant heterotrophic species within the biofilm at $O_{\infty} = 5\; gm^{-3}$. Results in Figure 2 demonstrate that higher ammonium concentrations at $O_{\infty} = 5\; gm^{-3}$ correlate with a monotonic increase in total biomass volume fraction, albeit with varying contributions from each species. Increased oxygen concentrations enhance the formation of AOB and $B_{h3}$, while causing slight modification in the production of $B_{n}$. In summary, at a concentration of $A_{\infty} = 6\; gm^{-3}$, AOB is enriched at $O_{\infty}\geq 1\; gm^{-3}$, shortcut denitrifier ($B_{h2}$) at $O_{\infty}$ between $1 \; gm^{-3}$ and $2.5\; gm^{-3}$, and denitrifier ($B_{h3}$) at $O_{\infty}\geq 2.5\; gm^{-3}$. The fractions of anammox and NOB ($B_n$) remain consistently low across all tested bulk concentrations of ammonium nitrogen and oxygen.

4.2.   The influence of oxygen and ammonium nitrogen bulk concentrations on nitrogen removal

In MABRs, nitrogen is removed through various pathways. The biofilm's sensitivity to environmental conditions, coupled with the counter-diffusion mechanism in MABRs, makes it challenging to predict which pathway will dominate and how the others can be optimized to enhance the reactor's performance. To explore how environmental factors influence total nitrogen removal and the activity of each bacterial species, we vary the bulk concentrations of oxygen and ammonium and compute the total nitrogen degradation/conversion through the pathways incorporated in model (2.1) when the biofilm occupies 30% of the computational domain (biofilm thickness = $600\; \mu m$). Ammonium nitrogen and oxygen concentrations are set to $N_{1\infty} = 5, \; 10, \; 15, \; 100, \; 200, \; 400\; gm^{-3}$ and $O_{\infty} = 0.25, \; 0.5, \; 1, \; 2.5, \; 5\; gm^{-3}$, respectively, while the acetate bulk concentration is held constant at $A_{\infty} = 6\; gm^{-3}$. As illustrated in Figures 3 and 4, the minimum nitrogen removal occurs at $O_{\infty} = 0.25\; gm^{-3}$ across all ammonium concentrations, primarily via the anammox pathway, underscoring the importance of sufficient oxygen availability for effective nitrogen removal. Increasing the oxygen concentration to $O_{\infty} = 0.5\; gm^{-3}$ and ammonium nitrogen concentration enhance the activity of autotrophic bacteria, slightly improving nitrification and making PN/A and ND (via $B_{h2}$) the primary contributors to nitrogen removal. Further increase in oxygen concentration promotes the growth rate of AOB, which enhances nitrification and the metabolism of $B_{h2}$ and anammox (due to increased nitrite nitrogen formation), thereby improving total nitrogen removal mainly through ND and PN/A pathways. As oxygen concentration reaches $O_{\infty} = 2.5\; gm^{-3}$, the metabolic activity of AOB further increases, resulting in higher nitrification rates and more nitrite nitrogen production. This leads to increased nitrogen removal via ND, followed by CND at $N_{1\infty} = 5\; gm^{-3}$ (due to promotion of conversion of $NO_2$ to $NO_3$ via NOB) and PN/A at higher ammonium concentrations. Note that NOB activity decreases as $N_{1\infty}$ increases due to the outgrowth of $B_{h2}$ in competition with other species for nitrite nitrogen and acetate. When the oxygen concentration is further increased to $O_{\infty} = 5\; gm^{-3}$, the formation of an anoxic zone is delayed, thereby reducing the activity of $B_{h2}$ and $B_{h3}$, which lead to a decrease in total nitrogen removal despite increased nitrification. Under ammonium-limited ($N_{1\infty} = 5\; gm^{-3}$) and oxygen-rich conditions, NOB outgrow species competing for nitrite, promoting the growth of $B_{h3}$ and making CND the major pathway for nitrogen removal. As $N_{1\infty}$ increases, NOB activity declines, reducing nitrogen removal via CND due to the decreased growth of $B_{h3}$, and shifting the main nitrogen removal pathways to PN/A and ND. Overall, the simulation results indicate that effective nitrogen removal is only achievable at $O_{\infty} > 1\; gm^{-3}$, with an optimal range for ND and PN/A nitrogen removal observed between $O_{\infty} = 1 \; gm^{-3}$ and $2.5\; gm^{-3}$ under the conditions of $N_{1\infty} = 10- 400\; gm^{-3}$ and $A_{\infty} = 6\; gm^{-3}$.

To facilitate comparison of the results shown in Figures 3 and 4, the numerical data are summarized in Table 3, where $\frac{NR_1}{NR_T}$, $\frac{NR_2}{NR_T}$, and $\frac{NR_3}{NR_T}$ represent the contributions of PN/A, ND, and CND nitrogen removal, respectively, to total nitrogen removal and $\frac{NR_T}{\text{Nitrification}}$ is the ratio of total nitrogen removal to nitrification. Key trends observed from the simulations under varying concentrations of $N_{1\infty}$ include (1) approximately 95% nitrogen removal occurs via the PN/A pathway at $O_{\infty} = 0.25\; gm^{-3}$; (2) PN/A and ND contribute equally at $O_{\infty} = 0.5\; gm^{-3}$; (3) at $O_{\infty} = 1\; gm^{-3}$, 80%–90% of nitrogen removal is achieved via ND, although its contribution decreases as $N_{1\infty}$ increases; (4) at $O_{\infty} = 2.5\; gm^{-3}$, the PN/A pathway contributes 4.1%–22.5%, while ND accounts for 50.8%–74.4% of total nitrogen removal. The contribution of PN/A increases as $N_{1\infty}$ rises, while the effect of ND reduces for $N_{1\infty} > 5\; gm^{-3}$, and (5) at $O_{\infty} = 5\; gm^{-3}$, 99% of nitrogen removal is achieved via CND at $N_{1\infty} = 5\; gm^{-3}$, but its contribution decreases to 82.9% at $N_{1\infty} = 10\; gm^{-3}$, and further to 52.4%, 25.2%, 24.2%, and 23.9% at $N_{1\infty} = 15, \; 100, \; 200, \; 400\; gm^{-3}$, respectively. When nitrogen removal is dominated by the PN/A and ND pathways, the value of $\frac{NR_T}{\text{Nitrification}}$ can be equal to or exceed 1 due to the contribution of the PN/A pathway. Furthermore, it is important to note that the impact of $N_{1\infty}$ on nitrogen removal pathways and $\frac{NR_T}{\text{Nitrification}}$ becomes negligible for $N_{1\infty} > 100\; gm^{-3}$.

4.3.   Oxygen consumption under varying ammonium concentrations

Aeration is the most energy-intensive process in biological nutrient removal wastewater treatment plants, accounting for approximately 50%–90% of the total energy consumption ^[46]. Consequently, oxygen demand emerges as a critical factor influencing the overall energy efficiency of MABRs. Since providing higher concentrations of oxygen results in higher energy usage which is not cost-effective, it is essential to elucidate the impact of oxygen levels on nitrogen removal mechanisms and to develop strategies for enhancing MABR performance without increasing oxygen input. In this section, we investigate oxygen consumption through different oxygen dependent mechanisms within the proposed model to identify the primary contributors under different environmental conditions. Specifically, we vary the bulk concentrations of ammonium and oxygen as $N_{1\infty} = 5, \; 10, \; 15, \; 100, \; 200, \; 400\; gm^{-3}$ and $O_{\infty} = 0.25, \; 0.5, \; 1, \; 2.5, \; 5\; gm^{-3}$, while maintaining the acetate bulk concentration at $A_{\infty} = 6\; gm^{-3}$.

Figure 5 illustrates the total oxygen consumption and the contributions of each oxygen-dependent mechanism at a biofilm thickness of $600\; \mu m$. At $O_{\infty} = 0.25\; gm^{-3}$, oxygen consumption is minimal, with heterotrophic bacteria ($B_{h1}$) being the primary oxygen consumers due to their competitive growth advantage over autotrophic bacteria under oxygen-limited conditions. As the oxygen concentration increases to $O_{\infty} = 0.5\; gm^{-3}$, AOB growth is stimulated, particularly at higher ammonium nitrogen concentrations, leading to increased total oxygen consumption. Further increasing the oxygen concentration to $O_{\infty} = 1\; gm^{-3}$ enhances AOB activity, promoting the growth of acetate-utilizing anaerobic heterotrophic bacteria (AHB), which in turn reduces the abundance of HB and establishes AOB as the dominant oxygen consumers. A subsequent increase in oxygen concentration marginally boosts the growth of nitrifying organisms (NOB), especially at higher ammonium nitrogen levels, resulting in increased total oxygen consumption, nevertheless AOB production continues to be the major pathway for oxygen uptake.

In conclusion, the simulation results indicate that oxygen consumption is insensitive to the increase in $N_{1\infty}$ from $10 \; gm^{-3}$ to $400 \; gm^{-3}$. At $O_{\infty} = 0.25\; gm^{-3}$, where anammox nitrogen removal is the dominant pathway, the oxygen consumption is approximately $0.1\; gO_2/m^2/d$. When nitrogen removal occurs primarily through PN/A and ND pathways at $O_{\infty} = 0.5\; gm^{-3}$, oxygen utilization increases to about $\approx 0.2\; gO_2/m^2/d$. At $O_{\infty} = 1\; gm^{-3}$ and $O_{\infty} = 2.5\; gm^{-3}$, where ND nitrogen removal is dominant, oxygen uptake reaches approximately $0.4$ and $1.1\; gO_2/m^2/d$, respectively. The highest oxygen consumption, approximately $\approx 2.2\; gO_2/m^2/d$, occurs at $O_{\infty} = 5\; gm^{-3}$ when the conventional nitrogen removal is dominant. Noteworthy here is that the contribution of endogenous respiration remains negligible across all scenarios.

4.4.   Impact of organic substrate on nitrogen removal under varying ammonium and oxygen levels

To thoroughly assess the impact of environmental conditions on nitrogen removal, this section explores how variations in acetate bulk concentration influence the efficiency of each nitrogen degradation pathway. Specifically, we calculate nitrogen removal via different mechanisms at $N_{1\infty} = 15, \; 100\; gm^{-3}$, $O_{\infty} = 0.25, \; 0.5, \; 1, \; 2.5, \; 5\; gm^{-3}$, and $A_{\infty} = 6, \; 30, \; 100\; gm^{-3}$. Results obtained at the biofilm thickness $600\; \mu m$ are illustrated in Figures 6 and 7. Moreover, a summary of total nitrogen removal ($NR_T$) for different scenarios, is provided in Table 4. At $A_{\infty} = 6\; gm^{-3}$, the highest nitrogen removal rates are $0.16\; gN/m^2/d$ for $N_{1\infty} = 15 \; gm^{-3}$ and $0.17\; gN/m^2/d$ for $N_{1\infty} = 100 \; gm^{-3}$, both occurring at $O_{\infty} = 2.5 \; gm^{-3}$. When $A_{\infty}$ is increased to $30\; gm^{-3}$, the maximum nitrogen removal rates rise to $0.47\; gN/m^2/d$ and $0.54\; gN/m^2/d$ for $N_{1\infty} = 15 \; gm^{-3}$ and $N_{1\infty} = 100 \; gm^{-3}$, respectively, at $O_{\infty} = 5 \; gm^{-3}$. However, further increases in $A_{\infty}$ to $100 \; gm^{-3}$, lead to significant reductions in maximum nitrogen removal rates, falling to $0.076 gN/m^2/d$ for $N_{1\infty} = 15 \; gm^{-3}$ and $0.089 gN/m^2/d$ for $N_{1\infty} = 100 \; gm^{-3}$ at $O_{\infty} = 5 \; gm^{-3}$. At $A_{\infty} = 6\; gm^{-3}$, nitrogen removal increases by increasing $O_{\infty}$ from $0.25 \; gm^{-3}$ to $2.5 \; gm^{-3}$, driven primarily by ND nitrogen removal and enhanced nitrification. However, at $O_{\infty} = 5\; gm^{-3}$, CND becomes dominant due to increased $NO_2$ to $NO_3$ conversion by NOB ($B_n$) at $N_{1\infty} = 15 \; gm^{-3}$. Conversely, at $N_{1\infty} = 100\; gm^{-3}$, oxygen concentration is reduced by the growth of AOB, which limits NOB activity, making ND the dominant pathway and decreasing total nitrogen removal. The results suggest that at $A_{\infty} = 6\; gm^{-3}$, total nitrogen removal is sensitive to oxygen availability across $N_{1\infty}$ values from $15 \; gm^{-3}$ to $100 \; gm^{-3}$. At $A_{\infty} = 30 \; gm^{-3}$, nitrogen removal rates remain low until $O_{\infty} > 2.5 \; gm^{-3}$, limited by AOB growth and nitrification. At $O_{\infty} = 0.25$ to $1\; gm^{-3}$, increasing $A_{\infty}$ from $6 \; gm^{-3}$ to $30 \; gm^{-3}$ leads to increased growth of aerobic heterotrophic bacteria ($B_{h1}$), which outcompete AOB for oxygen, making nitrification the limiting step in nitrogen removal. Increasing $O_{\infty}$ to $2.5 gm^{-3}$, enhances the nitrification rate due to the proliferation of AOB. However, raising $A_{\infty}$ from $6 \; gm^{-3}$ to $30 \; gm^{-3}$ promotes the growth of aerobic heterotrophic bacteria ($B_{h1}$), which surpass AOB and NOB in oxygen competition. This competition reduces the contribution of CND, thereby decreasing nitrogen removal. Nonetheless, nitrogen removal improves at $O_{\infty} = 5\; gm^{-3}$. At $A_{\infty} = 100\; gm^{-3}$, although nitrification increases with increase in $O_{\infty}$, nitrogen removal rates across $O_{\infty} = 0.25\; gm^{-3}$ to $5\; gm^{-3}$ are significantly lower than those achieved at $A_{\infty} = 6$ and $30\; gm^{-3}$. This reduction is due to low nitrification at low $O_{\infty}$ ($0.25$ to $1\; gm^{-3}$) and inefficient denitrification at $O_{\infty} = 2.5$ and $5\; gm^{-3}$. The high growth rate of aerobic heterotrophic bacteria ($B_{h1}$) surpasses that of anaerobic bacteria ($B_{h2}$ and $B_{h3}$) under conditions where both oxygen and acetate are non-limiting factors for $B_{h1}$. Thus, nitrogen removal by MABR biofilms is highly sensitive to both oxygen ($O_{\infty}$) and acetate ($A_{\infty}$) concentrations. Figure 6 and Table 4 highlight that optimal nitrogen removal occurs at $O_{\infty} = 5\; gm^{-3}$ and $A_{\infty} = 30\; gm^{-3}$ for $N_{1\infty} = 15$ to $100\; gm^{-3}$. For $A_{\infty} = 100\; gm^{-3}$ and higher, the elevated carbon substrate concentration stimulates aerobic heterotrophic growth and reduces the oxygen available for nitrification. This leads to a significant reduction in nitrogen removal across all pathways.

4.5.   Effect of temperature on biofilm composition and nitrogen removal pathways

The reaction rates of biological processes occurring within the membrane significantly influence the performance of MABRs. Therefore, accurately estimating bacterial growth rates is essential for evaluating MABR efficiency. Among the microbial species involved, anammox bacteria are particularly noted for their slow growth rates. To improve the precision of anammox growth rate estimation, Zhang et al. reevaluated the maximum specific growth rate by directly measuring the exponential increase in 16S rRNA gene copy numbers of two distinct anammox species, Ca. Brocadia sinica and Ca. Jettenia caeni ^[47]. Their study, conducted at 37℃, reported the maximum growth rate of $0.33 \pm 0.02 \; d^{-1}$ for Ca. B. sinica and $0.18 \; d^{-1}$ for Ca. J. caeni, both higher than previously recorded values. To assess how these growth rates affect substrate removal in MABRs, we use the average of the reported anammox growth rates and update the growth rates of heterotrophic and autotrophic bacteria using the Hoff-Arrhenius equation as $\mu_{(B_a, B_{h1}, B_{h2}, B_{h3}, B_n)}(T) = \mu_{{(B_a, B_{h1}, B_{h2}, B_{h3}, B_n)}_{20}}\theta^{T-20}$, where $\mu_{{(B_a, B_{h1}, B_{h2}, B_{h3}, B_n)}_{20}}$ represents the rates at 20℃, T[$^\circ$ C] is temperature and $\theta$ is the dimensionless temperature coefficient. This formulation is used to estimate the growth rates of all microbial species at 37℃. The results obtained at 20℃ and 37℃ under the condition that maximizes the total nitrogen removal ($N_{1\infty} = 100\; gm^{-3}$, $O_{\infty} = 5\; gm^{-3}$, and $A_{\infty} = 30\; gm^{-3}$) and the biofilm thickness $600\; \mu m$ are summarized in Table 5.

Results in Table 5 show that increasing temperature from 20℃ to 37℃ increases the total nitrogen removal by 31.5%, nitrification by 16.4%, and acetate consumption by 10.4% and reduces the oxygen consumption by 1.6%. The nitrogen removal fraction by PN/A, ND, and CND are 3.0%, 92.6%, and 4.8% at 20℃, respectively, while they change to 19.7%, 76.1%, and 4.5% at 37℃, respectively. The most modified pathway for nitrogen removal is degradation trough anammox due to more increase in the maximum growth rate of anammox. For further investigation of the effect of increasing temperature on nitrogen removal, we let biofilm grow to occupy 70% of the computational domain and record degradation of nitrogen via incorporated mechanisms in this study within the biofilm, i.e., with respect to biofilm thickness. Spatial spreading of biofilm occurs when total biomass density, $M$, approaches maximum cell density (1 after non-dimensionalization). Therefore, mass of the developed biofilm and its volume are equivalent and can be computed as the product of the biofilm thickness and area (i.e., biofilm thickness = mass/area). Using the total amount of biomass computed in Section 3.2, we can calculate the biofilm thickness and record the amount of biofilm species and substrate degradation within the biofilm. Results are depicted in Figures 8 and 9.

Increasing the temperature from 20℃ to 37℃ leads to an increase in the volume fractions of $B_a$, $B_{h2}$, $B_{h3}$, and $B_{an}$, while decreasing the volume fractions of $B_{h1}$ and $B_n$ within the biofilm. Notably, the volume fraction of anammox bacteria ($B_{an}$) grows exponentially with biofilm thickness, whereas $B_{h1}$ decreases sharply once the biofilm thickness reaches approximately 0.4 $mm$ at 20℃ and 0.2 $mm$ at 37℃. The reduction in $B_{h1}$ suggests the development of an effective anoxic zone within the biofilm, which promotes the growth of anoxic bacteria $B_{h2}$, $B_{h3}$, and anammox. Figure 9 illustrates the changes in nitrification and nitrogen removal rates with increasing biofilm thickness at both 20℃ and 37℃. At 20℃, nitrification, nitrogen removal via ND, and total nitrogen removal exhibit similar trends as biofilm grows. These trends are characterized by a rapid rise with the biofilm thickness increasing from 0.1 $mm$ to 0.42 $mm$, followed by a transition period between the biofilm thickness 0.42 $mm$ to 0.8 $mm$, and finally reaching a steady state for the biofilm thickness beyond about 0.8 $mm$. This alignment suggests that ND nitrogen removal plays a dominant role in total nitrogen removal, with a well-coordinated relationship between AOB activity and nitrite reduction bacteria ($B_{h2}$) at 20℃. At 37℃, both nitrification and total nitrogen removal are higher than those at 20℃. However, the rate of ND nitrogen removal decreases once biofilm thickness exceeds 0.4 $mm$. Despite this decline, total nitrogen removal continues to increase due to the sustained rise in PN/A nitrogen removal. CND contributes minimally to total nitrogen removal under both temperature conditions. Initially, denitrification rates increase but then quickly drop after reaching a peak, reflecting the variable growth rates of $B_{h3}$ and $B_n$ during the early stages of biofilm development. Thus, the increase in temperature promotes the dominance of nitrogen removal via ND and PN/A, with anammox becoming increasingly important as biofilm thickness increases.

5.   Discussion

This study developed a two-dimensional, multi-species biofilm model to investigate the roles of ND, PN/A, and CND in nitrogen removal within a membrane aerated biofilm reactor (MABR) under various operational conditions. At 20℃, the highest nitrogen removal rate of $0.54\; gN/m^2/d$ was achieved at bulk oxygen concentration ($O_{\infty}$) of $5 \; gm^{-3}$, ammonium nitrogen concentration ($N_{1\infty}$) of $100 \; gm^{-3}$, and acetate concentration ($A_{\infty}$) of $30 \; gm^{-3}$. Reducing $N_{1\infty}$ from $100 \; gm^{-3}$ to $15 \; gm^{-3}$ resulted in only a slight decrease in nitrogen removal (from $0.54 \; gN/m^2/d$ to $0.47\; gN/m^2/d$). However, decreasing the bulk oxygen concentrations to $2.5 \; gm^{-3}$, lowering $A_{\infty}$ to $6 \; gm^{-3}$, or increasing $A_{\infty}$ to $100 \; gm^{-3}$ led to reductions in nitrogen elimination by 72%, 70%, and 83.6%, respectively. Overall, nitrogen removal was limited by nitritation at low $O_{\infty}$ ($\leq 2.5 \; gm^{-3}$), denitritation at low $A_{\infty}$ ($\leq 6 \; gm^{-3}$), or the overgrowth of heterotrophic bacteria (HB) at high $A_{\infty}$ ($> 100 \; gm^{-3}$). These findings suggest that $O_{\infty}$ and $A_{\infty}$ exert a more significant influence on nitrogen removal rates and pathways than $N_{1\infty}$. This observation contrasts with previous studies, which emphasized the importance of the COD-to-nitrogen ratio as a key factor affecting nitrogen removal in MABRs, often based on variations in COD concentrations while maintaining constant nitrogen levels and surface loading rates ^[11,12]. The results of this study suggest that independent control of $A_{\infty}$ and $O_{\infty}$ is likely more effective for optimizing MABR nitrogen removal performance than controlling the ratio of $A_{\infty}$ to $N_{1\infty}$. From a practical perspective, $O_{\infty}$ can be effectively regulated through the supply of oxygen. Although $A_{\infty}$ is primarily governed by the intrinsic characteristics of the wastewater, it can be controlled to some extent by adjusting influent feed strategies, managing liquid recirculation between bioreactor chambers, and adding external carbon sources.

In this study, ND is identified as the dominant nitrogen removal pathway under most tested conditions. For moderate acetate concentrations ($A_{\infty}$ between $30 \; gm^{-3}$ and $100 \; gm^{-3}$), ND dominance is observed at oxygen concentrations ($O_{\infty}$) between $2.5 \; gm^{-3}$ and $5 \; gm^{-3}$. At low acetate levels ($A_{\infty} = 6\; gm^{-3}$), ND is the primary pathway for $O_{\infty}$ values between $1 \; gm^{-3}$ and $2.5 \; gm^{-3}$. CND becomes dominant only under specific conditions, namely $A_{\infty} = 6\; gm^{-3}$, $O_{\infty}\; gm^{-3}$, and ammonium nitrogen concentrations ($N_{1\infty}$) between $5 \; gm^{-3}$ and $15 \; gm^{-3}$. These results indicate that maintaining $O_{\infty} \leq 5 \; gm^{-3}$ promotes ND while effectively inhibiting the CND pathway in scenarios with moderate bulk acetate and ammonium concentrations. Under conditions yielding the highest nitrogen removal rate at 20℃ ($O_{\infty} = 5 \; gm^{-3}$, $N_{1\infty} = 100 \; gm^{-3}$, $A_{\infty} = 30 \; gm^{-3}$), approximately 92% of nitrogen removal occurred via the ND pathway, 3% through anammox, and 5% through conventional denitrification.

Nitrogen reduction processes are primarily determined by oxygen and acetate distribution within the biofilm. The simulation results obtained under the conditions of $O_{\infty} = 5 \; gm^{-3}$, $N_{1\infty} = 100 \; gm^{-3}$, $A_{\infty} = 30 \; gm^{-3}$, and biofilm thickness = 1 $mm$ show that oxygen is consumed within a 0.4 $mm$ depth from the membrane surface, while acetate is depleted within 0.6 $mm$ from the liquid-biofilm interface (c.f. Figure A1). The establishment of aerobic and anoxic/anaerobic zones in the biofilm facilitates ammonium oxidation to nitrite by AOB in the aerobic zone near the membrane. Controlling $O_{\infty}\leq 5 \; gm^{-3}$ create a balance between nitrite and oxygen concentrations in the aerobic zone for the ND pathway at a level that NOB may not be able to compete with AOB for oxygen and $B_{h2}$ for nitrite. This occurs because NOB have significantly lower affinities for both oxygen ($\kappa_{n_{o}} = 2.2 \; gm^{-3}$ vs. $\kappa_{a_{o}} = 0.6 \; gm^{-3}$) and nitrite ($\kappa_{n_{n2}} = 5.5 \; gm^{-3}$ vs. $\kappa_{h_{n2}} = 0.5 \; gm^{-3}$) compared to AOB and $B_{h2}$. As shown in Figure A1, with a balanced supply of oxygen and acetate, $NO_2^-$ is rapidly consumed by $B_{h2}$ in the anaerobic zone upon the availability of acetate as an electron donor. Thus, controlling oxygen and acetate supply is critical to synchronizing nitritation and denitritation, establishing ND as the dominant nitrogen removal pathway in MABR biofilms.

Under all conditions tested at 20℃ in this study, anammox biomass constitutes only a minor fraction of the total biofilm biomass. The simulation results demonstrate that oxygen concentration ($O_{\infty} = 0.25-5 \; gm^{-3}$) and nitrogen concentrations ($N_{1\infty} = 5-400\; gm^{-3}$) influence PN/A nitrogen removal under three acetate dosing conditions ($A_{\infty} = 6, \; 30, \; 100\; gm^{-3}$). At $A_{\infty} = 6\; gm^{-3}$, where carbon limitation may restrict conventional denitrification, the PN/A pathway is enhanced as $O_{\infty}$ increases from $0.25 \; gm^{-3}$ to $2.5 \; gm^{-3}$, particularly for $N_{1\infty}$ values in the range of $10 \; gm^{-3}$ and $400 \; gm^{-3}$ (c.f. Figure 3). At $A_{\infty} = 30$ and $100\; gm^{-3}$, PN/A nitrogen removal exhibits an increasing trend with $O_{\infty}$ ranging from $0.25 \; gm^{-3}$ to $5 \; gm^{-3}$ for nitrogen concentrations of $N_{1\infty} = 15$ and $100 \; gm^{-3}$ (c.f. Figure 6). PN/A becomes the dominant pathway only when nitrification is markedly limited by low oxygen availability, specifically at $O_{\infty} = 0.25 \; gm^{-3}$ across all conditions, and for oxygen levels between $0.25 \; gm^{-3}$ and $1\; gm^{-3}$ when $A_{\infty}$ is between $30 \; gm^{-3}$ and $100 \; gm^{-3}$ (c.f. Figures 3 and 6). As $O_{\infty}$ increases to $2.5 \; gm^{-3}$ or $5 \; gm^{-3}$, short-cut nitrogen removal pathways dominate in most scenarios, except under the conditions of $A_{\infty} = 6\; gm^{-3}$, $O_{\infty} = 5\; gm^{-3}$, and $N_{1\infty}\leq 15 \; gm^{-3}$, where conventional denitrification remains dominant. These findings underscore that precise control of oxygen supply at low concentrations is critical for promoting PN/A-dominated nitrogen removal. Nonetheless, even under conditions where PN/A performs optimally, its overall nitrogen removal rate remains significantly lower than that achieved by ND under higher oxygen levels. This limitation is primarily due to reduced nitrification efficiency at low $O_{\infty}$, which restricts the availability of nitrite for the anammox process. Increasing the temperature from 20℃ to 37℃ raises the contribution of PN/A to total nitrogen removal from 3% to 19.7%, driven by the increased growth rate of anammox bacteria at higher temperatures ^[47]. Nevertheless, this temperature increase does not substantially reduce oxygen or acetate consumption. This indicates that increasing temperature to promote the PN/A pathway offers limited benefits in terms of reducing energy and carbon requirements. Additionally, elevating the temperature of one cubic meter of water from 20℃ to 37℃ consumes approximately 20.0 $kWh$ of energy. Given the typical energy consumption of MABR systems around $6\; Kg-O_2/kWh$, the energy required for municipal wastewater treatment is only about $0.004\; kWh/m^3$, according to Peeters and Moran ^[48]. Therefore, raising the temperature of water to 37℃ is not a practical approach, except when treating already-heated water, such as return liquid from anaerobic digesters. While anammox is generally considered the most energy- and carbon-efficient nitrogen removal pathway, the nitrogen removal rate achievable through PN/A alone in MABR biofilms is limited. Therefore, operating MABRs to favor ND-dominant or ND and PN/A co-dominant nitrogen removal is likely to be more practical and beneficial, maximizing nitrogen removal while reducing energy and carbon inputs.

6.   Conclusions

In conclusion, this study successfully develops a comprehensive two-dimensional, multi-species biofilm model that elucidates the intricate dynamics of nitrogen removal processes in MABRs. By systematically varying environmental factors such as bulk concentrations of oxygen, ammonium, and acetate, as well as temperature and biofilm thickness, we have demonstrated their influence on the growth dynamics of various biomass species and the overall efficiency of nitrogen removal pathways, including shortcut denitrification (ND), anammox (PN/A), and conventional denitrification (CND). Our findings indicate that shortcut nitrogen removal is the predominant pathway under most operational conditions, with anammox thriving in oxygen-limited environments. Conversely, conventional denitrification becomes significant primarily at elevated oxygen levels and low ammonium concentrations. Notably, temperature emerged as a key factor enhancing overall nitrogen removal efficiency, primarily by promoting anammox activity. The simulation results further reveal that biofilm thickness plays a vital role in maintaining stable nitrogen elimination efficiencies, with thicker biofilms favoring the dominance of ammonium-oxidizing bacteria and shortcut denitrifiers.

The results obtained in this study demonstrate that the distribution and fate of oxygen, acetate, ammonium, $NO_3^-$, and $NO_2^-$ within the biofilm are influenced by the bulk concentrations of dissolved substrates and the biofilm thickness. While understanding the spatial distribution and transformation of these chemical species is essential for uncovering the nitrogen removal mechanisms in MABR biofilms, such insights are difficult to obtain through experimental studies alone due to the complexity and micro-scale structure of biofilms. This study demonstrates that mathematical modeling is a powerful tool for exploring these pathways. Moreover, the model developed here can be integrated with MABR process models to simulate the system's dynamic performance. The insights gained through the simulations provide a foundation for optimizing MABR design and operational strategies, ultimately contributing to more effective and sustainable wastewater treatment solutions.

Use of Generative-AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Conflict of interest

The authors declare no conflicts of interest.

Appendix

Spatial distribution of biomass volume fractions and substrates in biofilm