Available online sensors can be used to create fingerprints for MABRs that characterize biofilm limiting conditions and serve as soft sensors

Membrane aerated biofilm reactors (MABRs) are a promising biological wastewater treatment technology that offers process intensification combined with energy efficiency advantages and unique opportunities for enhancing nutrient removal (Houweling & Daigger 2019; Lu et al. 2020; He et al. 2021). MABR uses gas-permeable membranes, immersed in the wastewater to be treated, to supply oxygen to the inner portion of a biofilm that grows naturally on the outer surface of the membrane. Complementary substrates (e.g. ammonia and soluble organics) diffuse in the opposite direction from the bulk of the liquid into the biofilm. This is known as ‘counter-diffusion’. The benefits of MABRs include, but are not limited to, high oxygen transfer efficiency and oxygen transfer rate, excellent biological nutrient removal, little stripping of volatile organic compounds, resilience to loading shocks, reduction of greenhouse gas emissions, and compact footprints (Côté et al. 2015; Nerenberg 2016; Kinh et al. 2017; Houweling & Daigger 2019; Lu et al. 2020; Carlson et al. 2021; He et al. 2021; Pérez-Calleja et al. 2022). These benefits have fueled, beginning in the early 2010s, the rapid development and adaptation of MABRs in industrial applications to enhance carbon oxidation and nitrogen removal (Sathyamoorthy et al. 2018; Guglielmi et al. 2020; Shechter et al. 2020; Corsino & Torregrossa 2022; Pérez-Calleja et al. 2022). A sharp increase in pilot- and full-scale applications was seen in 2018, and at the end of April 2022 there were 335 contracted installations globally (BlueTech Research 2022).

Oxygen limitation refers to the situation when the oxygen flux into the nitrifying biofilm cannot satisfy the demand for full nitrification and results in residual ammonia in the biofilm. This often happens when (1) the partial pressure of oxygen in the lumen is lower than needed, or (2) soluble organics penetrate deeply into the biofilm so that heterotrophs outcompete nitrifiers for the oxygen supply (Casey et al. 2000; Hwang et al. 2009; Pérez-Calleja et al. 2022). A typical solution is increasing the lumen gas pressure to boost the oxygen diffusional driving force (Houweling et al. 2017; Houweling & Daigger 2019).

Ammonia limitation refers to a situation where the driving force for ammonia diffusion is insufficient to overcome the mass transfer resistance of the liquid–biofilm boundary layer plus diffusional resistance through the biofilm (gray rectangles). Little ammonia can reach the deep zone where nitrifiers are most active. This often happens when excessive biofilm introduces flow maldistribution, channeling, and clogging (Syron & Casey 2008; Wagner et al. 2022), and therefore, biofilm thickness control is crucial. The last limitation is co-limitation, where both substrate fluxes into the active nitrifier region are limited. Co-limitation is rarely seen in continuously operated MABRs. One possible reason for such limitation could be turning off both air supply and influent feed for an extended period.

Identifying which limitation the system faces is crucial for MABR operation and control, and to optimize performance. Several approaches can be used to directly or indirectly monitor the substrate concentrations in the biofilm and biofilm thickness, thereby quantifying the mass transfer rate or resistance. These include microsensors (Downing & Nerenberg 2008; Tan et al. 2014; Liu et al. 2020) and pressure decay tests using inert gas (Casey et al. 2014). Although these methods are promising, they require extra equipment installation and additional maintenance. It would be appealing if alternative approaches can be established using sensor measurements regularly provided in MABR systems. In addition to conventional water quality sensor monitoring (e.g. ammonia, nitrate), MABRs provide additional signals from the gas phase of the membrane lumen and an opportunity to better characterize system performance by direct assessment of oxygen transferred to the biofilm (Kunetz et al. 2016; Peeters et al. 2017; Rieger et al. 2020). This is accomplished by measurement of the gas flow to the MABR units, pressure, and oxygen purity in the feed (if needed) and exhaust gas. Experimental studies have reported strong correlations between gas-phase and liquid-phase signals in MABR applications (Houweling & Daigger 2019; Uri-Carreño et al. 2021). Moreover, observations of commercial scale MABR installations have shown a distinctive and repeatable pattern relating these gas-phase signals to reactor loading patterns (Houweling & Daigger 2019). This provides an obvious opportunity for process monitoring and control, yet little study that has systematically examined useful information could be extracted from these patterns and correlations.

This paper uses model simulations to: (1) quantitatively analyze the relationships between gas-phase and liquid-phase signals and proposes fingerprint plots for MABR system characterization, (2) identify the limiting conditions for MABR biofilms based on the fingerprint plots, (3) relate the fingerprint plots with mass transfer resistances, and (4) propose a potential opportunity for soft sensor development and their use along with the fingerprint plots to control MABR biofilm thickness.

Rather than being set as a fixed value, the oxygen concentration at the base of the biofilm is calculated based on a mass balance integrating oxygen partial pressure, volumetric air supply, and gas diffusion rates into the biofilm. The result is a model with the capability to simulate exhaust gas composition and the transition from oxygen-abundant to oxygen-limiting conditions as a function of biofilm activity, reactor loading rates, and volumetric air supply.

Key model assumptions include: (1) a one-dimensional stratified biofilm (zF in Figure 2) with three layers; (2) the total biofilm biomass is constant throughout the study, yet the species (e.g. heterotrophs and nitrifiers) within the biofilm are viable as a result of operational conditions (see Supplementary Materials Figure S9);(3) transport from bulk liquid into biofilm regulated by a mass transfer boundary layer thickness (zBL in Figure 2) governed by Fick's Law (Wanner & Reichert 1996; Takács et al. 2007); (4) transport of gases between the membrane lumen and base layer of the biofilm governed by Fick's Law; (5) the gas partial pressures is assumed to be linear along the longitudinal axis in the lumen space (Chen et al. 2021), therefore, the average oxygen partial pressure of the inlet and outlet of membranes was used for Fick's Law calculation; and (6) the outlet oxygen partial pressure is calculated based on the airflow and the influent loads condition. The assumption (4) about the diffusion rate of gases between the membrane lumen and base biofilm layer requires relating the partial pressure of gases in the membrane lumen to an equivalent saturation concentration through Henry's Law. The diffusion gradient is calculated based on this equivalent saturation concentration and the actual concentration as simulated in the base layer of the biofilm. Côté et al. (1989) proposed a more detailed ‘resistors-in-series’ model to account for (1) diffusion from the lumen across the gas boundary layer, (2) diffusion through the membrane material, and (3) diffusion from the membrane across the liquid boundary layer into the base layer of the biofilm. In the present implementation, these resistances are combined into the effective media boundary layer thickness (zMedia in Figure 2) which can be understood as a lumped parameter to account for their combined effect.

The biokinetic model used to simulate growth of autotrophic (nitrifying) and heterotrophic organisms in the biofilm was SUMO1, which is based on the well-known activated sludge model framework (Henze et al. 2000).

The bioprocess model was implemented in SUMO 2021 and configured as a single MABR unit, as shown in Figure 2. A state-variable-based influent unit and a standard effluent unit were connected with the MABR unit to build a SUMO model. The influent flow rate was 10,000 m³/d, and the MABR unit had a volume of 1,000 m³ with a total media surface area of 125,000 m². For simplicity, the influent total chemical oxygen demand (COD) concentration was set to 200 mg/L (150 mg/L slowly biodegradable COD + 50 mg/L readily biodegradable COD), and ammonia, whose concentration varied in various simulations, was the only nitrogen species presented in the influent. Operational settings for the MABR follow the recommended values provided in the SUEZ Zeelung 2.0 manual. Parameters not mentioned above were left at their default values unless explicitly listed in Table S1. Kinetic parameters follow the default values in SUMO1.

The primary MABR gas-phase signals usually consist of air flow rate, inlet and outlet gas pressure, and oxygen purity in the exhaust gas. Secondary signals calculated from these primary signals include oxygen transfer rate and oxygen transfer efficiency, which are key indicators of MABR performance (Houweling et al. 2017). Primary signals are more favorable for process control and analysis than secondary signals because they are directly measured and less prone to error propagation.

In this study, signals for bulk ammonia concentration and oxygen purity in the membrane exhaust gas are chosen, and their relationship plots are defined here as the MABR fingerprint plots. This choice was derived based on field observations, engineering intuition, and elegant mathematical features. Specifically, bulk ammonia concentration reflects mass transfer of the electron donors and oxygen purity reflects mass transfer of the acceptor. As electrons should always be balanced, these two signals obey certain patterns. Therefore, it is hypothesized that the changes in signal patterns could reveal information about MABR biofilms, just as fingerprints identify individual people.

It is the nature of bioprocesses that the organisms in them will eventually reach an equilibrium, subject to process configurations, loadings, and operations, which is known as a steady-state. The performance of bioprocesses often fluctuates around steady-state points and, therefore, steady-state points can be viewed as a foundation for bioprocess performance. Numerous steady-state simulations were implemented to (1) investigate the steady-state shape of fingerprint plots; (2) explore factors that influence the shape of the fingerprint plots and describe them in a valid mathematical form; and (3) identify the relationships between the parameters in the mathematical form and the diffusional layer thicknesses (zF, zBL and zMedia, defined in Section 2.1).

A total of 1,200 steady-state simulations were implemented using the SUMO model (Section 2.2) with combinations of different influent ammonia concentrations (Inf_SNHx), biofilm thickness (zF), liquid–biofilm boundary layer thickness (zBL), and effective media layer thickness (zMedia). The parameter values used are listed in Table 1. Each simulation generated a point on the fingerprint plot. Points with the same zF, zBL, and zMedia were connected to generate a steady-state curve, indicating system equilibria under different loading conditions.

Analysis of the fingerprint plots leads to formulation of a series of hypotheses concerning the underlying mechanisms and functional relationships that can be inferred from the plots. A dynamic simulation with a composite influent ammonia profile consisting of multiple steps and impulses was conducted to test these hypotheses. The ammonia and oxygen concentrations in different biofilm layers plus total nitrifier biomass were investigated and used to assess the behaviors of the dynamic trajectories in the fingerprint plots.

Taking any curve in Figure 3 as an example, as the membrane gas pressure and air flow rate are fixed, the oxygen flux into the biofilm is inversely proportional to the exhaust gas oxygen purity. The lower the oxygen purity, the higher the transferred oxygen flux. When the influent ammonia concentration is low (<5 mg-N/L) and full nitrification is achieved, changes in the bulk ammonia concentration are subtle while substantial oxygen mass transfer occurs, resulting in the observed decrease in oxygen purity and the vertical segment. It should be noted that, theoretically, the vertical segment should be close to but not vertical; however, for simplicity, it is assumed to be vertical with an infinite slope. At this stage, the system is ammonia-limited. When influent ammonia concentration further increases, the equilibrium bulk ammonia concentration also increases. The increased bulk ammonia concentration creates an additional driving force for ammonia flux into the biofilm and, therefore, increased demand for oxygen flux. As a result, the oxygen purity decreases as the bulk ammonia concentration increases, forming the inclined segment. The controlling factor on this segment remains ammonia, and the system is still ammonia-limited. Oxygen limitation occurs when the steady-state point enters the horizontal segment, where the residual oxygen in the lumen is not enough to overcome the mass transfer resistance to reach the biofilm, and the unconsumed oxygen ends up in the exhaust gas, so the increase of bulk ammonia no longer causes a decrease in exhaust gas oxygen purity.

The effects of biofilm thickness on the shape of the steady-state curves are shown in Figure 3. The end values of the vertical and horizontal segments are almost the same, indicating biofilm thickness is not a significant contributor. Biofilm thickness mainly affects the slope of the inclined segments. By visual observation, the thicker the biofilm, the milder the slope. This occurs because (1) nitrification activities mainly occur in the deepest layer of biofilm, which is closest to the gas lumen. (2) A thicker biofilm has two conflicting effects. One is increasing nitrification by housing more nitrifier biomass (Acevedo Alonso & Lackner 2019; Elsayed et al. 2021), which usually happens in the MABR start-up phase. The other is decreasing nitrification due to the greater mass transfer resistances caused by thicker biofilm (Janczewski & Trusek-Holownia 2016; Li & Liu 2019). The latter is found to be more influential in this study. (3) The trans-biofilm diffusional resistance is proportional to the biofilm thickness. The thicker the biofilm, the greater the diffusional resistance. To overcome greater resistance, higher bulk ammonia concentrations are needed. Assuming the amount of oxygen transferred is fixed (an arbitrary horizontal line in Figure 3), the intersections of this horizontal lines with the sloping line represent the bulk ammonia concentrations needed to deliver the ammonia flux to consume this given oxygen flux. A thicker biofilm pushes the intersection to the right, indicating that higher bulk ammonia concentrations are needed to match the oxygen flux.

This negative flux–thickness relationship has been observed by others. Pérez-Calleja et al. (2022) reported that the ammonia flux first increases with biofilm thickness because of the accumulated biomass, and then decreases with increased biofilm thickness due to the greater mass transfer resistance. Their critical thickness was around 100 μm, which is smaller than the minimum biofilm thickness considered in this study. The selected biofilm thicknesses in this study are already beyond the critical thickness and only capture the monotonic behavior after exceeding the critical thickness. So, the priority effect of the increased biofilm thickness in this study is the greater diffusional resistance rather than greater ammonia uptake. Common agreement is reached about the effect of biofilm thickness (Chen et al. 2021; Elsayed et al. 2021; Sanchez-Huerta et al. 2022), and biofilm thickness control remains a vital and ongoing topic in MABR practical applications (Karna & Visvanathan 2019; Lu et al. 2020; He et al. 2021; Silveira et al. 2022).

The effects of liquid–biofilm boundary layer thickness on the shape of the steady-state curves are shown in Figure 4. Similar to the effect of biofilm thickness, liquid–biofilm boundary layer thickness is not a significant contributor to the end values of the vertical and horizontal segments. Moreover, the thicker the layer, the milder the slope and the greater the diffusional resistance. Higher bulk ammonia concentrations are needed to transfer the same ammonia flux that matches the oxygen supply. Similar results were reported in other studies (Shoji et al. 2020; Chen et al. 2021; Elsayed et al. 2021; Pérez-Calleja et al. 2022).

The effects of media layer thickness on the shape of the steady-state curves are shown in Figure 5. Unlike zF and zBL that mainly impact the inclined segments, the influence of zMedia covers all three segments. Decreasing zMedia increases the lower bound of the vertical segment, steepens the inclined segment, and lowers the horizontal segment. Lower values of the vertical segment indicates that a greater gradient is needed to maintain the same oxygen transfer. The lower horizontal segments indicate that a greater maximum oxygen flux can be drawn from the gas inside the membrane. The gentler slope reflects greater diffusional resistances from the bulk liquid to the biofilm.

In conclusion, among the three thickness parameters investigated in this study, i.e., biofilm thickness, liquid boundary layer thickness, and media thickness, a thinner biofilm and liquid boundary layer thickness would solely decrease diffusional resistant for ammonia from the bulk liquid to the biofilm, while a thinner media layer thickness would result in both decreased oxygen and ammonia transfer resistance.

In the earlier sections, it was observed that thicknesses influence the shape of steady-state curves, which can be described with transformed ReLU functions. Determining a transformed ReLU function only requires three parameters (k, m, c). A follow-on research question is whether quantitative relationships can be established between the ReLU function parameters and the thicknesses.

All steady-state curves were fitted with the transformed ReLU function, whose fitting parameters are provided in Table S2. All fits have high R² values (>0.975), indicating the transformed ReLU shape is a commonly shared shape. Statistics for the fitted parameters are provided in Table 2. Parameter k (⁠⁠) is the absolute value of the slope of the inclined segment, with greater values indicating thinner diffusional layers and lower diffusional resistance. Parameter m, calculated from ⁠, is the x coordinate of the intersection point between the inclined segment and the horizontal segment. This point marks the transition from ammonia limitation to oxygen limitation. The higher this value, the higher the bulk ammonia concentration needed to create an ammonia flux that matches the same oxygen flux. Parameter c is the y coordinate of the horizontal segment, which reveals the lowest steady-state oxygen purity in the exhaust gas in a continuous air supply mode. It represents the maximum capacity for oxygen flux given certain air pressure and flow rates.

Quantitative relationships are found when regressing these parameters with different thicknesses. Both training data (steady-state curves generated with Table 1) and testing data (additional 20 steady-state curves that are randomly generated) displays good fits between predicted and actual values (Figures S3–S5). The final regression models are shown in Table 3. Feature selection was performed by removing terms from the regression model when the coefficients of their corresponding terms are not statistically different from zero (P value > 0.5, Figures S6–S8).

As for the parameter k, the slope of the inclined segment, all three thicknesses are statistically significant, and higher second-order terms are essential for the model, especially the interactive terms. Analysis on the derivatives of the regression model helps characterize influences of each thickness on the ReLU parameters. For instance, to characterize how thicker boundary layers influence k, partial derivatives are needed, ⁠. In this study, given zBL = 100–500 μm, zF = 200–600 μm, and zMedia = 80–240 μm, is always negative. Therefore, increasing zBL results in milder slopes, indicating greater resistance, which aligns with observations in Section 3.1.3. Another important observation is that also depends on how thick zF and zMedia are, as these two variables are on the right-hand side of the derivative equation. The parameter m, the x coordinate of the intersection point between the inclined segment and the horizontal segment, has fewer terms (rows in Table 3) than k. Thicker zMedia shifts m towards the left (smaller value). On the contrary, the other two thicknesses shift m towards the right when they increase. Again, these match observations in Sections 3.1.2 and 3.1.3. The parameter c, the y coordinate of the horizontal segment, is only a function of zMedia, matching results in Section 3.1.4.

When the biofilm has acclimated to a high loading condition with increased nitrifier biomass, the dynamic trajectories behave differently as shown in Figure 7. Instead of deviating, trajectories move along the steady-state curve. Another interesting finding is that, as a response to high peak concentration (blue), the dynamic trajectory directly enters the oxygen-limited segment without deviations. Moving along steady-state curves indicates systems are acting like at steady-state, where performance is at equilibrium with the limiting substrate. Therefore, the limiting factor for Figure 7 is not biomass but substrate fluxes into the biofilm.

Figures 6 and 7 provide several insights into MABR systems, for instance, resilience. Resilience is defined as the ability of a system to quickly recover from stressed loadings back to normal conditions (Veleva et al. 2022). Figures 6 and 7 illustrate this for dynamic trajectories that move along the steady-state curve where the MABR is performing at its current best. Even when the influent concentration is below its normal concentration, the MABR can quickly recover its original performance. It is hypothesized that resilience capacity is proportional to the total nitrifier biomass in the biofilm.

As shown in Figure 9(a), deviations from the steady-state curve occurred during phase I. It was hypothesized that this occurred because of biomass limitation, as indicated in Figure 8(c) and 8(d), where the inner layer has both substrates substantially present. The nitrifier biomass also continued to grow, as observed in Figure 8(b). The new steady state was reached as the residual ammonia concentration in the inner layer dropped to a low level. Phases II and III tested the resilience of the biofilm. These impulses introduced ammonia concentration disturbances in the outer two layers, while the inner ones remained relatively static. Changes in the DO profiles in the inner layer were also observed. Meanwhile, the dynamic trajectories moved on the steady-state curve, demonstrating the ability of the system to handle the impulse peaks. Phases III and I had the same peak elevation, yet only Phase III was able to handle the peak because Phase III had more nitrifier biomass in the biofilm and thus more resilience. Lastly, in Phase IV, where the biofilm was approaching oxygen limitation (the horizontal segment), the DO concentration in the inner layer biofilm was almost zero and the nitrifier biomass entered a plateau.

Based on the results of the study, the authors propose the following hypothesis concerning use of the bulk ammonia concentrations and the oxygen purity in the exhaust gas to monitor the condition of a MABR biofilm:

• •

  Steady-state curves lay a foundation to assess MABR performance, which consists of three segments, namely, vertical, inclined, and horizontal. The shapes formed by these three segments can be mathematically characterized by transformed ReLU functions with three fitting parameters, k, m, and c.

• •

  Dynamic trajectories, together with steady-state curves, can help identify limiting conditions of MABR biofilms.

• • o If dynamic trajectories deviate from the steady-state curve, the system is currently biomass-limited.

  • o If dynamic trajectories move on the vertical and the inclined segments, the system is ammonia-limited. Ammonia concentrations are depleted in the inner layer of the biofilm.

  • o If dynamic trajectories move on the horizontal segments, the system is oxygen-limited. Oxygen concentrations are depleted in the inner layer of biofilm.

• •

  Values of these fitting parameters relate to three different diffusional thickness layers and the relationship can be described with second-order polynomial models.

The proposed hypothesis offers the following practical applications based on the online signals of bulk ammonia concentration and oxygen purity in exhaust gas.

The dynamic trajectories can be used to characterize steady-state curves, and therefore system performance. For a particular MABR product, as its membrane material, permeability, and thickness is determined, the value of the horizontal segment should be a constant, which can be determined experimentally and consistently used for all applications. The vertical and inclined linear segments can be determined based on dynamic trajectories. When the system is only ammonia-limited, dynamic trajectories are almost linear, and the inclined segment can be directly characterized. The intercept of the inclined segment should be less than 21%. Although non-linear profiles exist when the system is biomass-limited, the initial portion of the dynamic trajectory is linear and lies on the inclined segment, which can be used to determine the slope and intercept. Once the system steady-state curve is characterized, limiting conditions can be easily identified and inform decision-making.

As demonstrated in Section 3.1.5, the parameters of the steady-state curves, k, m, and c, can be predicted using three thicknesses. Although the specific thicknesses cannot be uniquely determined from these parameters mathematically, the determined parameters can be used for MABR process operation and biofilm thickness control. For instance, since zMedia is fixed by the physical design of the system, the inclined segment is a function of zBL and zF. Inferences on thicknesses can be made if the inclined segment is characterized. This is significant because, in practice, there are limited tools available to measure them. In addition, for both of these thicknesses, their influence on the slope is the same – the thicker they are, the milder the slope, which makes control design based on slope feasible. To be specific, available control actions reduce both thicknesses at the same time, that is through gas purging to the external membranes (aeration) to slough off excess biofilm and reducing the boundary layer thickness (Heffernan et al. 2017; Shechter et al. 2020). Therefore, it is advantageous to combine them and use the slope value as a soft-sensing signal to describe diffusional resistances and to design automatic control scheme.

Experimental investigations are needed to further develop and evaluate these potential uses of the identified fingerprints. However, this work indicates that such investigations can potentially be beneficial.

The relationship between the residual bulk ammonia concentration and exhaust gas oxygen purity exhibits distinctive and repetitive patterns that are frequently observed in practical MABR applications. This relationship can be visualized by assembling multiple steady-state points to form a steady-state curve and by plotting the dynamic trajectories describing how this relationship point moves as time elapses. The latter is defined as the dynamic trajectory, which can be captured by online sensors. The steady-state curves and the dynamic trajectory are informative in characterizing MABRs and are proposed as MABR fingerprint plots. With an established SUMO model, behaviors of fingerprint plots were systematically examined, quantified, and explained.

The steady-state curves share a common shape consisting of three linear segments characterizing steady-state performance of the MABR system. The shape can be generally characterized by a ReLU function that is uniquely determined by three parameters. Results revealed these three parameters describing the ReLU function closely correlate with liquid boundary layer thickness, biofilm thickness, and membrane media thickness. Thus, these three parameters can be utilized to infer MABR mass transfer resistances. The dynamic trajectory and the steady-state curve are closely related. On the one hand, overlaying the dynamic trajectory on the steady-state curve can reveal the factors limiting MABR performance. On the other hand, the initial linear portion of the dynamic trajectory can be used to derive the steady-state curve, making it possible to plot steady-state curves based on online signals. It is hypothesized that the defined fingerprint plots can form the basis for MABR process monitoring and control. To be specific, parameters that characterize the shape of the fingerprint plots can be used to develop soft sensors for monitoring MABR mass transfer resistances, to identify limiting conditions for improving process operations, and to develop process control design, for example biofilm thickness control based on the slope of the fingerprint plots.

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

The authors declare there is no conflict.