Modelling the impacts of operational conditions on the performance of a full-scale membrane aerated biofilm reactor

(cid:129) A full-scale WWTP+MABR model is developed. (cid:129) The Model predicts the plant and MABR performance under two operational conditions. (cid:129) Additional simulations provide insights about bio ﬁ lm structure and gradients. (cid:129) Scenario analysis shows how predictions change when modifying modelling settings.


H I G H L I G H T S G R A P H I C A L A B S T R A C T
• A full-scale WWTP+MABR model is developed. • The Model predicts the plant and MABR performance under two operational conditions. • Additional simulations provide insights about biofilm structure and gradients. • Scenario analysis shows how predictions change when modifying modelling settings.

A B S T R A C T A R T I C L E I N F O Editor: Paola Verlicchi
Keywords: Mathematical modelling Nitrogen removal Process intensification Simulation Scenario analysis Membrane Aerated Biofilm Reactors (MABR) are gaining more and more acceptance in the plethora of wastewater process intensification technologies. Mathematical modelling has contributed to show their feasibility in terms of reduced energy consumption and footprint. Nevertheless, most simulation studies published until now are still focused on analyzing MABR as single units and not fully integrated within the flow diagram of the water treatment plant (WWTP). In this paper, the prediction capabilities of an integrated modelling approach is tested using full-scale data from Ejby Mølle WWTP+MABR site (Odense, Denmark). Mass balances, data reconciliation methods, process simulation and the different evaluation criteria were used to adjust influent, effluent and process indicators. Results show 10 % mismatch between flow, COD, N and P predictions and measurements in different plant locations. Using the adopted hydraulic retention time (HRT), nitrogen load (NL), membrane surface area (MA) and oxygen transfer rate (OTR), it was possible to predict nitrification rates (NR) within the interquartile range. This has been done under two different MABR operational conditions: with (#S2) and without (#S1) external aeration (EA) in the bulk liquid. The model provides additional process insights about biofilm structure, substrate gradients, weak acid base chemistry and precipitation potential. More specifically, simulations suggest the potential undesirable effects of sulfate (SRB) and iron reducing bacteria (IRB) on both microbial activity and composition of the biofilm. The latter may have a strong impact on ammonium (NH x ), sulfate (SO x ) and ferrous ion (Fe +2 ) conversion processes. The change of operational strategy in the scenario analysis highlights that the denitrifying activity of phosphorus accumulating organisms (PAOs) can enhance nutrient removal in MABR tanks. In addition, it was possible to assess the chance of success (in terms of energetic cost of nitrogen removal) of adding several MABR units in one tank of the WWTP under study before full-scale implementation.
Science of the Total Environment 856 (2023) 158980

Introduction
Membrane-aerated biofilm reactor (MABR) technology has gained considerable interest recently, emerging as an alternative to conventional biological nutrient removal with the potential to combine the benefits of biofilm-based technologies and low-energy aeration supply (Martin and Nerenberg, 2012). The MABR is based on the use of counter-diffusional biofilms, in which air is supplied through hollow-fiber or flat-sheet membranes, which at the same time function as biofilm support. Since microbial activity on the surface of the membrane acts as a catalyzer for oxygen transfer , the oxygen transfer efficiency (OTE) in the field can be higher than in clean water, with the potential of achieving significant reductions in energy requirements compared to diffused aeration (Castrillo et al., 2019). MABR aeration is based on diffusion, and (in terms of energy requirements) it should be compared to aeration based on bubbles. MABR technology is a bubbleless system, and that is the main reason of the reduction in energy requirements, compared to conventional systems such as activated sludge which is based on bubbles.
Besides the well-known advantages of biofilms for wastewater treatment, such as the possibility of working at high biomass concentrations, the ease of separation between biofilm / product and the coexistence of aerobic and an anoxic/anaerobic (Van Loosdrecht and Heijnen, 1993), the counter-diffusional nature of the MABR offers additional benefits provides more suitable conditions simultaneous nitrification and denitrification (Timberlake et al., 1988, Chen andZhou, 2022).Furthermore, this allows the MABR to be retrofitted into existing anoxic activated sludge zones to provide additional nitrification capacity. This type of MABR configuration is referred to as "hybrid", in which anoxic activated sludge surrounds the MABR membranes and can result in complete total nitrogen removal, i.e., achieving removals approaching 100 % at high organic loadings (Downing and Nerenberg, 2008).
Mathematical models have proven to be useful tools to virtually study MABRs. The scientific literature offers numerous works dealing with different aspects related to the MABR process and its operation. In this way, it has been possible to investigate the competition between microorganisms (Li et al., 2018;Chen et al., 2016;Sun et al., 2017), the effect of mass transfer limitations (Perez-Calleja et al., 2017;Liu et al., 2022) and how substrate/biomass change for different geometries (Acevedo-Alonso and Lackner, 2019) and axes (radial / longitudinal) (Chen et al., 2020). The main findings of these studies allowed improving nitrogen and/or sulfur removal, reduce N 2 O emissions or improve aeration efficiency of the systems (pilot / full-scale) under study.
Even though the previous works made a substantial progress in the field of water process systems engineering, there are important aspects that remained untouched. First and foremost, in all these studies, MABRs are studied as separate units, but never integrated within the context of a full wastewater treatment plant (WWTP). This is mainly attributed to software limitations for all the studies that have been carried out (Reichert, 1994). The latter made it for example impossible to study how changes in the operational conditions /control strategies may affect MABR performance (Gernaey et al., 2014). Secondly, earlier modelling studies did not consider how water chemistry may change through the biofilm. For example, pH may affect the P, S and Fe interactions in water treatment systems (Solon et al., 2017). The latter comes with the well-known stiffness and simulation, problems and special solving routines must be developed (Batstone et al., 2012, Feldman et al., 2017. Last but not the least, and strongly related to the first point, there are very few references where MABR has been fully modelled as intensification technology given a process layout. Only very few theoretical studies have been published up to date (Carlson et al., 2021). Hence, it is quite difficult to assess the extra value of adding MABR cassettes from a holistic perspective when it comes to nitrogen removal and/or energy savings.
The study presented herein goes beyond the state of the art in the field of mathematical modelling of water treatment technologies by presenting a simulation study predicting the performance of an MABR pilot plant included within a full-scale WWTP (Uri-Carreño et al., 2021. The main contribution of this work relies on establishing an integrated modelling approach where MABR are a truly process intensification technological option. Special emphasis is placed to account for P, S and Fe interactions both biologically and physico-chemically. These reactions are especially important in systems with chemically enhanced primary treatment (CEPT) and wastewater with high contents of sulfates. Biochemical reactions are based on the well-established ASM2d model (Henze et al., 2000) but with additional upgrades described in Solon et al. (2017) andKazadi Mbamba et al. (2019). As a result is possible to describe the effect Fe/S oxidation/reduction/precipitation processes may have in the overall process performance. Data reconciliation methods are used to adjust mass balances through the WWTP (Monje et al., 2022 (Feldman et al., 2017). Special numerical schemes have been included in order to make computer simulations feasible . The paper details the development of the proposed modelling approach. Prediction capabilities are verified with full-scale data. Additional scenario analyses further explore the effect of selected operational conditions (Section 3). Results are framed in the context of existing literature (Section 4). Lastly, opportunities and limitations that arise from utilization of these models are discussed in detail.

Case study: Ejby Mølle WWTP
The Ejby Mølle (EM) Wastewater treatment plant (WWTP) in Odense, Denmark, has a 410,000 population equivalents treatment capacity. The main liquid treatment train is comprised of grit removal and screening, chemically enhanced primary treatment with the addition of ferric sulfate, and Bio-Denipho™ nutrient removal including anaerobic zones for biological phosphorus removal. Nevertheless the activated sludge section is not operating in phases (the output of the anaerobic selector does not change). Hence, the plant is run as A 2 O/Bardenpho reactor. There is also a secondary line treating wet weather wastewater including storm tanks and first and second step trickling filters (only nitrification) with intermediate secondary clarifiers. Primary effluent and return activated sludge (RAS) stream are introduced before the anaerobic section. Wet weather flow (treated) and reject water are combined with the anaerobic reactor outflow just before the oxidation ditches (see Fig. 1, top left).
The MABR is adjacent to anaerobic zone at the EM WWTP and consisted of a 25 m 3 circular reactor. A full-scale hollow-fiber MABR cassette with a total volume of 11.3 m 3 and a total membrane surface area of 1920 m 2 was installed inside the pilot reactor in 2019 (see Fig. 1). The reactor was set up as a continuous stirred-tank reactor (CSTR) fed with mixed liquor from the full-scale anaerobic zones i.e. hybrid MABR. The feed (wastewater + return activated sludge) was pumped from the anaerobic zone and introduced to the tank near the bottom using a distribution grid. A recirculation pump was utilized to ensure proper mixing of the tank. The tank is controlled with state of the art sensing equipment. The effluent was located at the top of the tank and it is sent to the oxidation ditch. Hence, the whole sequence is in anaerobic -> hybrid MABR -> aerobic -> secondary settling. More details can be found in Uri-Carreño et al. (2021) (see Fig. 1, top right).
Low-pressure air was supplied to the MABR unit for process oxygen supply and mixing. Mixing airflow was maintained at 12 m 3 h −1 throughout the study period due to concerns regarding the feed's high suspended solids concentrations. Process and mixing airflows to the MABR were adjusted manually.
External aeration (EA) was provided with fine bubble diffusers on one side of the tank to increase the redox conditions inside the reactors. The EA was operated intermittently and controlled using two different strategies: time-based control during the first three months and using an oxidation reduction potential (ORP)-based feedback control during the study's remainder. The ORP-based feedback control also included managing the feed flow to the tank: Decreasing the flow when the ORP reached a predefined set-point.

Historical data, statistical reconciliation and calculation of (process) performance indicators
Plant data used in this study comprises 1 years of measurements between January and December 2020 (1-2 measurements per week). Samples were taken from 6 locations: 1 point at the raw wastewater inlet, 2 points after the primary clarifier (over and underflow), 2 points at the secondary clarifier (over and underflow), 1 point at the activated sludge section. There were some occasional measurements of the storm water and reject water line (see Fig. 1, top left). Measurements involved the determination of: pH, total suspended solids (TSS), volatile suspended solids (VSS), chemical oxygen demand (COD), total nitrogen (TN), total phosphorus (TP) and total iron (TFe) These analyses were complemented with filtered (0.45 μm) samples for quantifying ammonium / ammonia (NH x ), phosphates (H x PO 4 ), nitrates (NO x ), sulfates (SO x ), dissolved iron (Fe +2 ) and soluble calcium (Ca +2 ). Standard methods were used for the experimental determination of the different process parameters.
Flow-rate, COD, N and P data (totals) was reconciled using a five step approach suggested by Monje et al. (2022) where: 1) identity matrix was defined, 2) data were cleansed and curated, 3) missing fluxes were estimated, 4) optimal flows were calculated using Lagrange multipliers and 5) the new data set was verified.
In addition, it was possible to obtain 3 years (June 2018 to March 2021) of on-line measurements from the MABR unit, namely: input and output NH x and oxidation-reduction potential (ORP) signals. Gaseous carbon dioxide (CO 2 ) and dissolved oxygen (O 2 ) were also measured in the exhaust (see

Biochemical and physico-chemical reactions
Reaction rates in both bulk and biofilm are evaluated using the stoichiometry and kinetics as described in the Activated Sludge Model No 2 (ASM2d) (Henze et al., 2000). The operational temperature is constant at 15°C. Biological sulfate and iron reduction processes are considered by assuming multiple different organic substrates (acetate, fermentables) (Solon et al., 2017;Hauduc et al., 2019). In addition, conversion of sulfide (H 2 S / HS − ), to sulfur mineral (S 0 ) and then back to sulfates (SO x ) is carried out by sulfur oxidizing bacteria (SOB) involving different electron acceptors (O 2 , NO x ) (Chen et al., 2016). Chemical S oxidation is also accounted for (Gutierrez et al., 2010).
A common physico-chemical framework is implemented to simulate the acid-base system and is thereby able to predict pH (Solon et al., 2017). Multiple mineral precipitation potential is based on saturation Index (SI) calculations as described in Kazadi-Mbamba et al. (2015). Fe transformations are based on the principles described in Kazadi Mbamba et al. (2019) (oxidation, reduction, precipitation). Formation of hydrous ferric oxides (HFO) plus absorption and co-precipitation is modelled as described in Hauduc et al., 2019. Mass transfer between the liquid (bulk) and the gas phase is accounted for selected compounds (i = O 2 , NH 3 , H 2 S, N 2 and CO 2 ). The transport rates are formulated as a function of the difference between the saturation concentration and the actual concentration of the gas dissolved in the liquid (Batstone et al., 2012). Transport of O 2 , N 2 and CO 2 through the membrane is implemented as proposed in Casey et al. (2000).
Primary clarification is reactive to account for COD hydrolytic/ fermentation processes and Fe-P removal. TSS are separated using a splitter (Jeppsson et al., 2007). Different settling velocity speeds are assumed for both organic and inorganic material (Ekama and Wentzel, 2004). Secondary clarification is implemented following the principles stated by Takács et al. (1991) i.e. double exponential velocity function using a ten layer pattern. RAS and WAS are taken from the secondary settler underflow.

Multi-scale reactor model
The multi-scale approach described in Feldman et al. (2017) and Flores-Alsina et al. (2019) is used to model the EM WWTP. The Bio-Denipho™ reactor hydrodynamics was approximated by sequencing multiple CSTRs. Influent and return activated sludge enters reactor 1. The effluent of the trickling filter + reject water stream are introduced after reactor 4 (see Fig. 1, bottom left). The MABR tank includes both bulk and biofilm section. More specifically, the biofilm model is based on the one-dimensional approach reported in Wanner et al. (2005). The model contains both soluble and particulate state variables. The mass balance assumes that the transport of soluble compounds is governed solely by (homogenous) diffusion, whereas movement of particulate compounds takes place by convection (Saravanan and Sreekrishnan, 2006). Biofilm thickness (L) is given as the radial distance (z) and varies due to two phenomena: 1) the net growth of the particulate species, and 2) detachment from the biofilm surface (Lackner et al., 2008). The resulting system of partial differential equations (PDEs) is solved using the method of lines (Press et al., 2007). Details are depicted in Fig. 1 (bottom right).

Implementation and process simulation details
Both data reconciliation methods and process models are implemented in Matlab/Simulink. More specifically, biological and physico-chemical models are coded in C and run in Simulink as S-functions. Data reconciliation, model parameters, initial guesses and pre & post processing routines are coded as .m files. All the tools developed for this study are available for free upon request (see software availability section).
Two operational modes are considered: without (#S1) and with (#S2) external aeration (EA) in the bulk. In order to account for this, a PI controller is implemented in order to maintain the DO concentration at 0.01 and 0.15 g/m 3 respectively (see controller implementation details in Flores-Alsina et al., 2008). An average constant biofilm thickness (L) 0.75 mm and density (50,000 g/m 3 ) is assumed for both periods.
A two-step simulation procedure is adopted in order to avoid having convergence problems. First, the EM WWTP is run until reaching steady state conditions including the MABR tank but without accounting for the biofilm activity. COD, N and P are fractionated using the principles reported in Gernaey et al. (2011). Prediction capabilities are assessed calculating the deviation between the plant measurements (historical data) and computer simulations using the reconciled data. In the second step, previous steady state values are used as starting point for all the main units, i.e. primary clarifier, biological reactors, secondary settlers. With respect to the biofilm, MABR bulk value are used as initials. All particulates are assumed to be equally distributed. No changes in the density are considered. This ensures faster convergence and fewer numerical problems (mostly due the stiffness caused by the PCM). Simulations are run typically for 365 days to reach steady state conditions since the starting point is considered to be good already. The model could be solved all at once (both WWTP and MABR), but the initialization would take longer. Hydraulic retention time (HRT), membrane area (MA) and oxygen transfer rates (OTR) are obtained from historical data. Predicted nitrification rates (NR) are used to evaluate the performance of the MABR.

Scenario analysis
Two additional scenarios are used to test the capabilities of the proposed approach. The first set of simulations (SC1) shows the effect of changing the influent characteristics on the MABR process performance. A range of biofilm thickness values (L) are assumed which would correspond to different scouring strategies in the MABR tank. Next, different hydraulic retention time (HRT) values resulting from increasing/decreasing influent flow rate were assumed. The entire plant is re-simulated keeping the same OTR. EA in the bulk is maintained as defined in #S2. NR and phosphate removal efficiency (PRE) are used as (process) performance indicators. In the second set of (extra) simulations (SC2) the potential of using MABR as process intensification tool at the EM WWTP is assessed. The third anaerobic tank at EM WWTP is assumed to be re-converted into an MABR. Area/Volume (1920 m 2 /11.3 m 3 ) in the pilot plant is up-scaled to the full system. Again, OTR is assumed to be same as in # S2. The same applies here for EA. The plant in its current mode and the "intensified" plant are simulated assuming up to 30 % increase in loading rate. An ammonia controller is installed in both systems in order to ensure that total nitrification rates are comparable (see Flores-Alsina et al., 2008). RAS (external recycle) and WAS (secondary waste flow) are scaled up as flow increases. The ratio between nitrified nitrogen and energy consumption is used as (process) performance indicator. Aeration energy (AE) in the main aerobic reactor and EA in the MABR is quantified according to Gernaey et al. (2014). In the MABR tank, AE is estimated assuming the kgO 2 /kWh values reported in Uri-Carreño et al. (2021). Aeration for biofilm scouring is not accounted for. Fig. 2 a b, c & d shows both reconciled measurements and process simulations for flow, COD, N and P in the main points of the water line of the EM WWTP. Calculated average deviations are 2, 8, 15 and 7 % respectively. The latter demonstrates the capabilities of the proposed approach to describe in general terms organic carbon, nitrogen and phosphorus removal at EM WWTP. Mass balances indicate that 24 and 56 % of the incoming COD and N is removed within the activated sludge section. The remaining fraction ends up in the effluent (4 and 8 %) and waste (72 and 35 %) streams respectively (both primary and secondary). When it comes to P, a major part is accumulated in the sludge in form of polyphosphates and/or iron phosphates (97 %). A very minor fraction leaves with the secondary settler overflow (3 %). Default parameters reported in Solon et al., 2017 are used to run the simulations. Only Yields for OHO, AOB, NOB and PAO were increased 10 % to have a better match of the COD profiles within the reactor. Details about the model structure can be found in Supplemental information (with verified continuity checks).

Baseline simulation conditions
Using the adopted hydraulic retention time (HRT), nitrogen loading (NL), membrane area (MA) and oxygen transfer rate (OTR) it was possible to predict nitrification rates (NR) and phosphate removal efficiency (PRE) within the interquartile range. This has been done under the two considered MABR operational conditions: with (#S2) and without (#S1) external aeration (EA) in the bulk liquid (Fig. 3). In order to understand the differences in NR and PRE values it is necessary to examine what kind of processes are taking place within the biofilm and how those are related to the bulk conditions. This is further studied in the following sections.

Biofilm structure, microbial composition, distribution of inorganics & substrate profiles
The steady state simulations of the general biofilm structure and substrate profiles are illustrated in Fig. 4 (a, d). In both cases (#S1 and #S2), the higher content of biomass is found in the inner part of the biofilm (in close contact with the membrane). Inert material resulting from biological activity is generated within the biofilm and pushed towards the bulk by convective movement. Organics are rapidly consumed. Results also show how the different dissolved oxygen (DO) levels in the bulk (#S1 and #S2) affect the distribution and the quantity of inorganics. Indeed, the model suggests that the quantity of precipitates will be higher for #S1.
With respect to microbial composition (Fig. 4 b, e), ammonia oxidizing bacteria (AOB) are dominant inside the biofilm. This is attributed to the high affinity of these microbial groups for DO and the availability of substrate given the mass transfer limitations. Ordinary heterotrophic bacteria (OHO) are clearly favored in the outer parts of the biofilm. Phosphorus accumulating organisms (PAO) + sulfate (SRB) and iron (IRB) reducing bacteria are only present in the bulk. The latter two under low DO conditions (#S1) favor ferrous iron (Fe 2+ ) and sulfides (H 2 S/HS − ) formation, which will promote the growth of sulfide oxidizing bacteria (SOB) in the middle layers of the biofilm. Fig. 4 c and f shows the change of inorganics through the horizontal axis of the biofilm for #S1 and #S2. For example, hydrous ferric oxides (HFO) and sulfur mineral (S0) will be preferentially formed close to the membrane as result of the re-oxidation of Fe +2 and sulfates (SO X ) . The middle and outer parts of the biofilm will be occupied by iron sulfide precipitates (FeS) as product of SRB and IRB activity (particularly for #S1) Higher DO conditions in the bulk (#S2) will minimize H 2 S and Fe +2 presence i.e. IRB/SRB will be inhibited, hydrogen sulfide (H 2 S) stripped and Fe +2 H 2 S/HS − re-oxidized.
Substantial gradients for COD, N, S and Fe are predicted by the model (for each evaluated strategy). The ammonium is oxidized to nitrates by the AOB located close to the membrane (where DO concentration is higher) (Fig. 4 h, k). The latter is reduced to nitrogen gas (N 2 ) by the OHO present in the outer parts of the biofilm using organic substrates (acetate, fermentable) as electron donors (Fig. 4 g, j). Part of the NO X removal is also carried out by PAO in the bulk phase (this aspect is further explored in the scenario analysis section). Similar profiles can be observed for Fe 2+ and H 2 S/HS- (Fig. 4 m, p) i.e. rapidly oxidized to HFO/S0 when entering the biofilm. The model indicates a flat profile for phosphate (Fig. 4 i, l). It is important to highlight that both H 2 S and CO 2 can be stripped in the bulk.

Weak-acid base chemistry conditions and precipitation potential
More steady state model simulations show how the model predicts weak-acid base conditions as well as the risk of precipitation (Fig. 4 n trend towards the center. This is attributed to two main factors. First, the CO 2 (with acidic nature) is either consumed by AOB or transferred from the biofilm to the membrane. Second, there are organic acid gradients (as previously described, from high to low) (Fig. 4 n, q). Besides the already accounted FeS the model indicates potential saturation conditions for calcium phosphate (not included within the model). Calcium carbonate (calcite) and iron phosphates (vivianite) is under-saturated for both reactors ( Fig. 4 o, r). In general terms, biofilm structure, microbial composition, distribution of inorganics and substrate gradients are consistent to each other for the two evaluated operational strategies (#S1 and #S2). These elements helps to explain the NR and PRE results reported in Fig. 3. Low DO levels in the bulk (#S1) promotes the growth of SRB and IRB (in the bulk) that will reduce HFO to Fe 2+ and SO x to H 2 S/HS − respectively. These two, will diffuse within the biofilm and will be susceptible to: 1) be reoxidized to S 0 /SO x and HFO increasing the internal DO demand, 2) be precipitated as FeS (outer biofilm) and/or HFO-P (inner biofilm) competing with biomass for space, 3) H 2 S can inhibit AOB growth. It is important to highlight that FeS also coats the membrane and hinders the mass transfer processes (Uri-Carreño et al., 2022). Nevertheless, this is not accounted for by the model. All these aspects are reflected in the relatively lower NR (see Fig. 3c). Higher DO levels in the bulk (#S2) will promote H 2 S/HS − and Fe +2 oxidative reactions (also stripping for H 2 S). Consequently, there will be no growth of both SRB and IRB, which will be causing the above mentioned problems resulting from Fe and S reduction processes. Hence, AOB will be working under optimal conditions (no competition for DO, no inhibition, no intra-biofilm precipitation). This explains the relatively higher NR values during #S2 (see Fig. 3c). The latter have strong impact on the overall phosphate removal (see Fig. 3f). Indeed, higher NR will increase the quantity of NOx, which will be used by PAO to accumulate P as poly-phosphates (PP). The ratio between NR/PRE seems to indicate that NOx is also used by other groups of microorganisms (in our case we assume only OHO) (Henze et al., 2000). The reader should know that the authors did some testing (results not shown in the paper) and the model provided much high NR rates when S and Fe reactions were not included (or EA activated under S1 loading conditions).

Scenario analysis
The additional simulations presented in this section demonstrate the potential of using the proposed approach when performing model-based studies. Therefore, for exemplary purposes, the adjusted model will be used for different engineering purposes. More specifically the following will be investigated: 1) the effect of changing the operational conditions (influent load, biofilm length) in the pilot MABR performance and 2) the potential use of MABR as process intensification at EM WWTP.

Effects of influent characteristics on MABR process performance
The generated response surfaces for Fig. 5 (a, b) show the effect of biofilm thickness (L) on NR and phosphate removal efficiency (PRE). Thin biofilms (high scouring) have a detrimental effect on the NR. This is mainly because the high COD values give a competitive advantage to the OHO (instead of AOB). When biofilm thickness is increased, COD transfer is limited and NH X can easily reach the inner parts of the biofilm (due to its high diffusivity). As a result, AOB are clearly favored and the NR is much higher. When the biofilm is too thick NR decreases again. This is attributed to: 1) excessive accumulation of inert material (lower fraction of biofilm is alive), 2) NH x diffusivity problems, which does not allow AOB to grow. The same pattern can be observed for the range of evaluated loads. The latter also reveals that higher NR do not increase, which indicates that the system is running at the maximum capacity. Fig. 5a also reveals that when the nitrogen loading rate is decreased to similar levels as applied in #S1, the study reveals higher NR (≈1 gN/m 2 .d). The higher DO levels in the bulk does not allow S and Fe reduction processes, which trigger the previously described undesirable effects on the overall process performance (FeS precipitation, H 2 S inhibition and DO competition). With respect to PRE, results show that a thinner biofilm substantially increases PAO activity. This is mainly due to a lower anoxic region within the biofilm, which allows for a higher flux of NOx to the bulk (this is the place where PAO are present). At higher biofilm thickness, NOx is rapidly consumed by biofilm OHO and PAO cannot grow (with the consequent P uptake and storage in the form of PP). Last but not the least, the reader should be aware that the range of values was very high for MABR biofilm thickness. Nevertheless, the analysis shows very well how mass transfer conditions may affect the overall process performance. Fig. 6a show the energetic cost of nitrogen conversion obtained from 1) plant data, and 2) model predictions (Baseline). Both data and model reveal that the approximate cost of N removal is 10 kWh / Kg N. Model simulations suggest that this quantity will increase up to 17 kWh / Kg N in case the current load is increased by 50 %. The Figure also includes the computer simulations resulting from running the model when assuming that the third anaerobic reactor is reconverted into an MABR (V reactor / A membrane is the same as in the pilot case) The study reveals that the latter option (WWTP+MABR) could convert the same quantity of nitrogen at lower energetic expenditure (= 35 % reduction in kWh / kg N). Compared to the standard plant, this difference when using an MABR instead of the third anaerobic reactor increases at higher nitrogen loads (up to 50 % difference kWh / kg N). In other words, in case the nitrogen load is increased the cost of its treatment would be reduced almost by half by installing MABR capacity.

MABR as process intensification technology at EM WWTP (SC2)
A large part of N is handled in the MABR (25 %), which has much lower energy requirements (kgO 2 /kWh) than traditional activated sludge plants. In terms of global nitrogen balances, the controller ensures that both scenarios have the same nitrification efficiency. The system does not seem to suffer in terms of phosphorus removal. Most of the P release takes place in the first three anaerobic selectors. In the fourth reactor, the NOx produced is used by dPAO for P uptake. The lower AOB DO requirements in the aerobic zone also favors PAOs growth. Last but not the least, Fe (over) addition, also contributes to keep P low at all instances.
Assuming an electricity cost of 0.062 Euro/kWh, it is possible to calculate the yearly savings (in terms of electrical costs) when comparing the two treatment solutions. Hence, at the default loading conditions, the model reveals that the WWTP + MABR implementation would save around 25,000 Euro/year. This difference would be increased up to 70,000 Euro/ year in the event of having an increase of the incoming N load. This could be projected over several years in order to calculate the return on investment (ROI). The results of the study suggest that WWTP + MABR may be more energetically efficient at the current load but it would be even better at higher loads.
It is important to highlight that for this particular case, both alternatives were capable to nitrify all the incoming N (at the different loading rates). The reader should be aware that EM WWTP is operating at an HRT of 1 day. Therefore, potential benefits of conducting N removal at reduced footprint, which could be very interesting in other system, do not apply here. For this reason, when comparing the yearly saving in terms of OPEX, only derived from energy saving may be taking into account. These estimation should be properly assessed using appropriate electricity billing information (Aymerich et al., 2015).

Simulation results in the context of existing literature
Results reported in this study are in good agreement with previous modelling studies published in specialized literature. The same AOB and OHO distribution within the biofilm could be found in Lackner et al. (2008), Wang et al. (2009), Ni et al. (2013, Chen et al., 2016) and Li et al. (2018). The effect of organic substrates and the competition between AOB and OHO is also studied reaching similar conclusions as in this manuscript. When dealing with S transformations, similar findings were reported in Sun et al. (2017) with simultaneous reduction and re-oxidation of sulfur compounds within the same biofilm. Indeed, the role of S had to be explicitly modelled since it has profound effect on the overall process performance. The effect of biofilm thickness has been investigated by multiple authors (Acevedo-Alonso and Lackner, 2019;Chen et al., 2021) also pointing out the accumulation of inert material and loss of effectiveness. With respect to P removal in a hybrid MABR, Carlson et al. (2021) found higher P removal efficiencies when compared to a conventional A2O configuration due to the effect of denitrifying PAOs. The same study reports lower oxygen requirements and potential higher biogas production.

Experimental validation of the multi-scale modelling results
The simulation results reported in this study are in good agreement with the experimental observations reported in Uri-Carreño et al. (2022). The latter study revealed that all biofilm samples (for the very same WWTP and evaluated period) were mainly a mixture of organic compounds (e.g., microbes). There was presence of inorganic precipitates such as hydrous ferric oxides (HFO) and calcium (Ca) with absorbed and/or coprecipitated P (see Fig. 2S a, b). The period before the implementation of the external aeration (EA) there was an increase of the quantity of Fe-S substances, which could be pyrite and its precursors like mackinawite (see Fig. 2S a). These results are aligned with the experimental findings of Prot et al. (2022), where it is reported that FeS is preferentially formed over vivianite and soluble iron. Solid S compounds were drastically reduced after including EA. These results are in good agreement with the main findings reported in Sections 3.2 and 3.3, where the effect of DO was simulated at the EM MABR. The model suggested: i) higher quantity of inorganics (#S1), ii) potential precipitation of FeS within the biofilm (#S1), iii) no accumulation of S0 (#S1, #S2) (see Fig. 3 a, j, c, l). The risk of having Ca-P precipitates and HFO (with absorbed and precipitated P) is the same independently of the redox conditions ( Fig. 3 i, s). These compounds disappeared when the membranes were cleaned and/or replaced (see Fig. 2 a, b). In the same study (Uri-Carreño et al., 2022), the qPCR analysis and post statistical processing confirmed the presence of OHO, PAO, SRB and IRB in both bulk and biofilm (see Fig. 3 S). The study also stated an AOB increase from 0.2 to 6.7 % when the EA was included. The latter would perfectly explain the higher nitrification rates for strategy#2 and the rest of the conversion rates reported in Fig. 3. In this study, the authors report the mechanisms (included within the biofilm) the justify the differences in the NR values: 1) lower inhibition values. 2) less competition between AOB and DOB and 3) lower accumulation of precipitation within the biofilm.

Chances, limitations and future research
Last but not least, the reader should be aware that the WWTP+MABR model presented in this case study shows a good compromise between model complexity and prediction capabilities. Nevertheless, the list of relevant processes is far from being complete. It is important to highlight that no hydrodynamic test was carried out. The anaerobic selector is divided in 4 section and so is the model. The longitudinal discretization of the upper and lower part of the oxidation was made in order to account 2 aerobic and 2 anoxic zones. No information about biofilm thickness was available. In fact, the access to the biofilm part in the pilot was very difficult and a b Load increase it was necessary a crane a several operators. For this reason, it was done in counted occasions. The assumed values (L = 0.75 mm) is a pretty large thickness and seems not realistic in some typologies of MABR (for instance, in hollow fiber bundles). The lack of VSS/TSS data within the plant did not facilitate the validation of the precipitation dynamics. We have an indication that these processes are happing but we do not know how fast/slow they really are. Intra biofilm precipitation is a complex issue (Feldman et al., 2017) and more thorough measurements are necessary. The model did not include attachment either. This is the main reason why we could not find PAOs in the biofilm as it was reported in Uri-Carreño et al. (2022). The way the model is constructed does not allow the PAO growth within the biofilm since there are no alternating anaerobic/aerobic conditions. Biofilm density is assumed to be constant, even though this is not completely true since FeS accumulation is considered. The latter would require a complete reformulation of the particulate mass balance and the way convective movement is described. Since the biofilm model is a 1D model, the MABR model considers homogeneous spatial distribution of the biofilm in the whole reactor (in the second and third dimensions). This assumption is another limitation of the model, since the real biofilm will be differently exposed to the bulk liquid in different zones of the membranes and reactor. Uncertainty/risk analysis, a central concept for model-based evaluation of alternatives (Belia et al., 2009), is an important issue that should be considered as well. The consideration of uncertainty/risk factors within the evaluation procedure can substantially change the evaluation results, and preliminary selected technologies might come up as less desirable while other options retain a higher chance of success under uncertainty (Flores-Alsina et al., 2008). Dynamic profiles could substantially change the results reported in Section 4.2. Loading increase is only simulated under steady state conditions and the effect of peaks is not accounted for. All these elements will be explored and included in subsequent versions of the model as we have more data available.

Conclusions
The main findings of this study can be summarized in the following points: 1) A WWTP+ MABR model approach is developed. The WWTP includes primary + secondary clarification and activated sludge processes. The MABR describes diffusion, convection and mass transfer (from membrane to biofilm and from biofilm to the liquid bulk). Besides transformations of COD, N and P, the model also accounted for the interactions with the S and Fe cycle. 2) Both mass and model predictions describe flow, organic carbon, nitrogen and phosphorus removal at different plant locations. The MABR module predicts nitrification rates and phosphate removal efficiency based on hydraulic retention time, nitrogen load, membrane area and oxygen transfer rate. 3) Multi-scale simulations provided additional insights about the processes explaining these differences: Biofilm structure, microbial composition, distribution of inorganics, substrate profiles, weak acid base chemistry and precipitation potential 4) The model suggests the harmful effect that both SRB and IRB activity may have on the overall MABR performance. Both plant data and model confirm the beneficial effects of including an external aeration unit and thus increase DO levels to counteract the harmful effects of H 2 S and Fe 2+ 5) The scenario analysis reveals how NR and the denitrifying potential of PAO are affected by the biofilm thickness. A very thin biofilm favors OHO growth (instead of AOB). A very thick biofilm increases the quantity of inert material (=less alive part). It also reduces PRE due to the lower flux of NOx arriving to the bulk since it is consumed by OHO in the biofilm anoxic section (and cannot be used by PAOs) 6) The potential reconversion of the third anaerobic reactor into an MABR system improves the kWh / kg N significantly. Hence the cost of treating N can be reduced by 35 % at the current load and may even reach 50 % when the load is higher.

Software availability
The Matlab/Simulink implementation of the data reconciliation files + WWTP + MABR are available on request. To express an interest for obtaining the code, please contact Dr. Xavier Flores-Alsina (xfa@kt.dtu. dk) or Prof. Krist V. Gernaey (kvg@kt.dtu.dk) at the Department of Chemical and Biochemical Engineering at the Technical University of Denmark. The authors will also make the code available through the GitHub platform (www.github.com/wwtmodels).

Data availability
Data will be made available on request.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.