Coupling a Simple and Generic Membrane Fouling Model with Biological Dynamics: Application to the Modeling of an Anaerobic Membrane BioReactor (AnMBR)

1. Introduction

The anaerobic Membrane BioReactor (AnMBR) is an interesting wastewater treatment technology, which couples anaerobic digestion treatment of organic pollutants with a physical separation between sludge and liquid allowing to obtain a highly purified effluent. To reach an optimal treatment efficiency, it is crucial to control both the biological and the separation processes. Thus, it is important to model biological dynamics and to couple them to a membrane filtration model, to predict membrane fouling, which remains, by far, the main drawback of MBR technology [1,2].

Recent work established that models developed to describe the conventional activated sludge process (Activated Sludge Models (ASM), [3]) and the anaerobic digestion (Anaerobic Digestion Model N.1, [4]), can be used to describe MBR dynamics when slightly modifying the model parameters ([5]). However, such models were developed above all for continuous and homogeneous reactors and are not able to account for specific components as Soluble Microbial Product (SMP) dynamics that are known to play an important role in membrane fouling ([6,7,8,9]). In a conventional bioreactor, the matter recycling (due for instance to biomass mortality) is not necessarily taken into account, because it can usually be neglected with respect to the dilution rate. In MBRs, this hypothesis no longer holds and variables describing the dynamics of certain classes of molecules, such as the SMP produced during biomass growth and mortality, must be added to the model [10]. For anaerobic systems, simpler models than the ADM1 have been coupled to SMP dynamics. For instance, we proposed an extension of the two steps-anaerobic digestion model (AM2) ([11]), in order to include their dynamics, ([12]).

Regarding membrane fouling, many models have been proposed such as the resistance-in-series model ([13]), models based on a sectional approach ([14]) and on the fractal geometry to describe the fouling cake permeability ([15]), models based on the local pressure and flux variation leading to the uneven fouling cake up on the membrane surface ([16,17]), and finally, physical models which have been proposed to study fundamental membrane properties [18,19,20]. Simpler models are proposed to describe fouling as the result of only the cake formation mechanism [21] or adding pores blocking phenomenon due to soluble matter [22]. These models are purely physical ones; they describe the dynamics of abiotic membrane parameters and completely neglect biological dynamics even if some authors, as [23], proposed to combine them to describe fouling in MBRs.

In order to properly describe the entire MBR dynamics, models describing the dynamics of biological processes in the reactor must be coupled to membrane fouling models. If a number of such integrated models have been proposed for aerobic MBRs (cf. for instance [13,24,25,26]), very few have been proposed for AnMBRs. In [27], authors developed a mathematical model for MBR, by considering together reversible and irreversible fouling. Mixed liquor suspended solids were assumed to be the major components of the reversible fouling layer and, dissolved organic matter is thought to be responsible for the long-term irreversible fouling. [28] proposed a filtration model based on the resistance-in-series model and able to reproduce the filtration process of a Submerged AnMBR (cake layer build-up and consolidation during filtration; membrane scouring by biogas sparging; removal of cake layer by back-flushing and irreversible fouling consolidation). This model was validated in the long-term under different operating conditions, using data obtained from a SAnMBR demonstration plant [29]. While such models have a high prediction capabilities, they are usually too complicated for being used for control purposes. [30] proposed a model for a Submerged MBR, with slow-fast dynamics and, they used this structure for the parameter estimation procedure. Thereafter, the model is used to develop a nonlinear predictive control.

In [12], the AM2b model was specifically proposed for control purposes. It is a simple model which describes only the two main biological processes of the anaerobic digestion while including the SMP dynamics. In the first step (acidogenesis) the acidogenic bacteria $X_1$ consume the organic substrate $S_1$ to produce Volatile Fatty Acids $S_2$ (VFA), SMP and $C{O}_2$, while in the second step (methanogenesis), the methanogenic population $X_2$ consumes VFA and produces SMP, methane and $C{O}_2$. However this model was not coupled with a membrane fouling model to completely describe the dynamics of an entire AnMBR in the simplest way we can think of for control synthesis purposes.

To summarize the state on the art about AnMBR models, one may say that available fouling models are either not coupled to the biological models, or they are too complicated to be used in process control [31]. It is precisely the aim of the present paper, where the novelty is to propose a simple and generic membrane fouling model of which the usefulness is illustrated in coupling it with biological model to completely describe an AnMBR and, to develop optimization tools and strategy. In addition, we consider only two fouling mechanisms for the membrane model, which depend on the biological model outputs (measurements of the biological model are inputs for the fouling model). Contrary to several literature studies on fouling modeling, we consider that the total filtering membrane surface area, is not constant during a filtration period nor after several filtration/stop cycles. It is described in a very general way as a decreasing function of deposited matters onto the membrane surface and retained into membrane pores.

The paper is organized as follows. First, the assumptions used for building the membrane model are given and discussed. Then, the dynamic equations of the model are presented with respect to the specific functioning phases considered: filtration and cleaning phases. Then, simulation results are presented in order to study the qualitative properties of the model. To illustrate the easiness with which it may be coupled to a biokinetic model, the fouling model is then coupled with the AM2b model and confronted to experimental data. Thereafter, some techniques for fouling control are discussed and simulated. Finally, conclusions and perspectives are formulated.

2. Mathematical Equations of the Proposed Membrane-Fouling Model

2.1. Model Development

[32] proposed a membrane fouling model including an explicit relationship between the mass of solid matter attached onto the membrane and the flux going through this membrane. The dynamics of the solid attachment to and detachment from the membrane were related to the filtration flux. Using a resistance-in-series model, they considered the total fouling resistance to be caused both by the pores clogging resistance due to the solute deposition inside the membrane pores, the dynamic sludge film resistance and the stable sludge cake resistance. However, the proposed model is not suitable for automatic control purposes since it is too complicated. The aim of the present section is to simplify this model to come up with a very simple fouling model while keeping realistic hypotheses. The idea is to include a feedback of the decreasing flux due to membrane fouling into the actual output flow rate $Q_{out}$ leaving the MBR. In other words, we propose to consider $Q_{out}$ as a decreasing function of the total mass solids attached onto the membrane surface and of the solute (as SMP) deposited inside the pores. As recalled in the introduction, as a matter of fact, many studies in the literature agreed that SMP have a crucial role in the membrane fouling, especially in pore clogging [6,7,8,9].

• Fouling mechanisms
  It is well known that the fouling dynamics is different depending on the fouling mode considered, namely pore constriction, cake formation, complete blocking and intermediate blocking [19]. In our simple model, we consider only the two main membrane fouling mechanisms, as defined in [18] (see Figure 1):

  –

  The first one is caused by the mass $m_c\left(t\right)$ of solids which attach onto the membrane surface also called cake formation or cake fouling. According to the particles concentration and solids attachment rate, particles are retained leading to a decrease of the filtering area of the membrane.

  –

  The second is due to the mass of particles $m_p\left(t\right)$ retained inside the membrane pores as SMP, called hereafter pore constriction. Their size may be smaller than the pore sizes and they are known to progressively clog the membrane pores. This phenomenon typically reduces the porous area of the membrane.

The proposed modeling approach allows us to decouple dynamically the different fouling mechanisms (solids attached onto the membrane vs SMP clogging the pores). Using a resistance-in-series model based on Darcy’s law, the membrane fouling model for plants operating at constant transmembrane pressure (TMP) is thus given by (1):

$$J\left(t\right)=\frac{Q_{out}\left(t\right)}{A\left(t\right)}=\frac{\Delta P}{\mu \left(R_0+R\left(t\right)\right)},$$

(1)

where $Q_{out}\left(t\right)$ the output flow rate, $A\left(t\right)$ the membrane surface area, $\Delta P$ the transmembrane pressure (assuming constant), $\mu$ the permeate viscosity, $R_0$ the intrinsic membrane resistance and $R\left(t\right)$ the fouling resistance given by (2):

$$R\left(t\right)=R_c\left(t\right)+R_p\left(t\right),$$

(2)

where $R_c\left(t\right)$ and $R_p\left(t\right)$ are the cake and the pore clogging resistances, respectively.

• Models of membrane resistance and membrane area
  $R_c\left(t\right)$ and $R_p\left(t\right)$ are typically dependent on masses $m_c\left(t\right)$ and $m_p\left(t\right)$, respectively, and are modeled by (3) (adapted from [13]):

  $$R_c\left(t\right)=\alpha \frac{m_c\left(t\right)}{A\left(t\right)},R_p\left(t\right)={\alpha }^{\prime }\frac{m_p\left(t\right)}{ϵA\left(t\right)},$$

  (3)

  with $\alpha$ and ${\alpha }^{\prime }$ the specific resistances and $ϵA$ ($0<ϵ<1$) the porous area which is a fraction of the total useful surface area A, (see Figure 2 for a flat sheet membrane for instance).
  Contrary to several literature studies on fouling modeling, we consider that the total filtering membrane surface area $A\left(t\right)$, is not constant during a filtration period nor after several filtration/cleaning cycles: it is described by (4) in a very general way as a decreasing function of $m_c\left(t\right)$ and $m_p\left(t\right)$:

  $$A\left(t\right)=\frac{A_0}{1+\frac{m_c\left(t\right)}{\sigma }+\frac{m_p\left(t\right)}{{\sigma }^{\prime }}},$$

  (4)

  with $A_0$ the initial membrane surface and $\sigma$ and ${\sigma }^{\prime }$ parameters used to model the contribution of $m_c$ and $m_p$ to the surface reduction. Such a function is well adapted if we assume that the total useful filtration area is composed of two parts: a filtering surface and a porous surface. If $m_p\left(t\right)$ increases, then the porous surface decreases leading to the total loss of $A\left(t\right)$ even if the cake fouling is not yet significant. Likewise, if $m_c\left(t\right)$ increases then the filtering surface decreases because attached particles may prevent the flux to circulate freely, even if the pores clogging fouling reaches its equilibrium or if it is not yet significant. In short, the membrane surface $A\left(t\right)$ decreases when $m_p\left(t\right)$ and/or $m_c\left(t\right)$ increase (mathematically, $A\left(t\right)$ tends to zero as $m_c\left(t\right)$ and/or $m_p\left(t\right)$ tend to infinity).
  The function (4) is also able to model the fact that the initial filtering surface $A_0$ is not totally recovered after a backwash or a chemical cleaning. Theoretically, in equation (4), if $m_c\left(t\right)=0$ and $m_p\left(t\right)=0$ when we operate the MBR plant for the first time, or after each perfect backwash of the membrane, then the area $A\left(t\right)$ is equal to its initial value $A_0$. However, in practice, after each membrane backwash or cleaning, there is small remaining quantities of $m_c\left(t\right)$ and $m_p\left(t\right)$ which are not detached, causing progressively an irreversible fouling effect. On the long term, the surface $A\left(t\right)$ continuously decreases, leading to the membrane degeneration.

• Models of attached solids on the membrane surface and blocked SMP into pores
  Both compounds $m_c\left(t\right)$ and $m_p\left(t\right)$ have their own dynamics: they increase during the filtration phase and decrease during the relaxation (or backwash) phase. Since it is assumed that the mixed liquor is homogeneous, we assume that all soluble components ($S_T=\sum S_i$, $i=1,2,...$ and SMP) and particulate components ($X_T=\sum X_i$, $i=1,2,...$) may contribute, at different degrees, to the membrane fouling by cake formation (solids attachment). Thus, the dynamic of the mass $m_c\left(t\right)$ can be described by (5):

  $${\dot{m}}_c=Q_{out}(C_s{S}_T+C_x{X}_T+C_{smp}SMP),$$

  (5)

  where $C_s$, $C_x$ and $C_{smp}$ are weighting parameters used to model the contribution and the rate of each variable to the cake formation. In practice, they must be adjusted using calibration data (see the experimental results Section 4).
  The membrane has a selective rejection: particulate components (biomass) and large solute compounds (as macro-molecules of SMP) are totally retained by the membrane (their size being supposed to be greater than the pores diameter), while part of the solute components (substrates and a fraction of SMP) go through the membrane without retention (their size is assumed to be smaller than the pores diameter). We propose the following dynamic model (6) for the pores clogging by $m_p\left(t\right)$:

  $${\dot{m}}_p=Q_{out}\left({\beta }_{1}SMP+{\beta }_2{S}_T\right),$$

  (6)

  where ${\beta }_1$ is a parameter used to calibrate the rate of pores clogging, by the fraction of SMP leaving the bioreactor, while ${\beta }_2$ is used to model the contribution of others solute substrates to the pores clogging.
  On the other hand, no back-diffusion of $m_c\left(t\right)$ and $m_p\left(t\right)$ to the bulk solution is considered: we assume it is negligible with respect to the remaining attached and blocked matter.

• Additional hypothesis: There is no biomass growth on the membrane surface and detached solids do not affect matter concentration in the bulk liquid
  For simplicity, we assume that the biological growth of the attached biomass on the membrane (as well as in the pores) is neglected. This hypothesis is justified by the fact that backwash or relaxation periods arise quite often. In addition, we assume that if there are detached quantities of $m_c\left(t\right)$ and $m_p\left(t\right)$ during relaxation, which return into the bioreactor, they can be neglected with respect to their corresponding concentrations in the bulk liquid (see Figure 3). In many operated MBR, detached matter by backwash is not returned into the reaction medium and is rejected elsewhere. Finally, both fouling mechanisms are considered to be partially irreversible, but at different degrees i.e, fouling by pores clogging is more irreversible than fouling by cake formation.

2.2. Fouling Model for the Filtration Phase

Starting from the previous equations and hypotheses, the complete fouling model for the filtration phase ($\Delta P>0$) is given by equations (7)-(11). The output flow rate $Q_{out}$ is a decreasing function: after a certain period of functioning when the permeate flux has dramatically decreased, the process must be stopped and membrane cleaning must be performed.

$$\begin{array}{ccc}\hfill {\dot{m}}_c& =& Q_{out}(C_s{S}_T+C_x{X}_T+C_{smp}SMP),\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill {\dot{m}}_p& =& Q_{out}\left({\beta }_{1}SMP+{\beta }_2{S}_T\right),\hfill \end{array}$$

(8)

$$\begin{array}{ccc}\hfill R& =& \alpha \frac{m_c}{A}+{\alpha }^{\prime }\frac{m_p}{ϵA},\hfill \end{array}$$

(9)

$$\begin{array}{ccc}\hfill A& =& \frac{A_0}{1+\frac{m_c}{\sigma }+\frac{m_p}{{\sigma }^{\prime }}},\hfill \end{array}$$

(10)

$$\begin{array}{ccc}\hfill Q_{out}& =& J.A=\frac{\Delta P.A}{\mu \left(R_0+R\right)}.\hfill \end{array}$$

(11)

In equation (7) the dynamics of $m_c\left(t\right)$ are proportional to the total soluble ($S_T$ and $SMP$) and particulate matter ($X_T$) deposited onto the membrane surface. If we can measure separately the components of $S_T$ ($S_1$, $S_2$, ...) and/or those of $X_T$ ($X_1$, $X_2$, ...), where each one contributes differently to the cake formation then, equation (7) can be written by (12):

$${\dot{m}}_c=Q_{out}(\sum C_{si}{S}_i+\sum C_{xj}{X}_j+C_{smp}SMP),$$

(12)

where $C_{si}$ ($i=1,2,...$) and $C_{xj}$ ($j=1,2,...$) are weighting coefficients.

In equation (), the dynamics of $m_p\left(t\right)$ are essentially proportional to the fraction ${\beta }_1.SMP$ going through the membrane and blocked into the pores (macro-molecules of SMP are retained by the membrane, see [12]). However, after a long enough filtration time, porosity decreases, and pores clogging is also assumed to be dependent on small deposited quantities of $S_T$, but less significantly than SMP. This is modeled by ${\beta }_2{S}_T$ , with ${\beta }_2$ chosen lower than ${\beta }_1$.

The dynamics of $m_c\left(t\right)$ and $m_p\left(t\right)$ depend on the values of variables $S_T$, $X_T$, $SMP$ and $Q_{out}$. When the system operates at steady state, $S_T$, $X_T$, $SMP$ and $Q_{out}$ reach theirs equilibria and thus, $m_c$ and $m_p$ converge to their steady state values where, their accumulation values are constant. Consequently, the fouling resistance $R\left(t\right)$ reaches its equilibrium corresponding to the maximum fouling. This interpretation corresponds to reality since membrane fouling increases with time. For a long enough filtration time, there is no longer attachment of matter onto the membrane, because of the cake layer already formed: this equilibrium corresponds to a functioning mode in which detachment by shear forces compensates the attachment of matter by permeation forces.

2.3. Fouling Model for the Relaxation (or Backwash) Phase

To clean the membrane by backwash, the feed of the MBR is stopped and the transmembrane pressure is inverted ($\Delta P<0$), so that permeate flow backs into the feed, lifting the fouling layer from the pores and the surface of the membrane. In certain cases, MBR are operated with a relaxation period instead of a backwash. In others terms, the flux is simply stopped ($\Delta P=0$) for a given short time allowing the natural detachment of matters and particles, which can be modeled by equations (13)-(14), with $\omega$ and ${\omega }^{\prime }$ positive constants to be adjusted with respect to experimental data, $m_{c_{irr}}$ and $m_{p_{irr}}$ positive quantities to model irreversible fouling after cleaning operation (quantities of irreversible $m_c$ and $m_p$). For instance, we assume that we have $x\%$ of irreversible fouling after relaxation or backwash, which means that after membrane cleaning, $x\%$ of $m_c$ and $m_p$ reached at the end of the previous filtration period is irreversible.

$$\begin{array}{ccc}\hfill {\dot{m}}_c& =& -\omega m_c+m_{c_{irr}},\hfill \end{array}$$

(13)

$$\begin{array}{ccc}\hfill {\dot{m}}_p& =& -{\omega }^{\prime }{m}_p+m_{p_{irr}}.\hfill \end{array}$$

(14)

The relaxation (or backwash) time is neglected compared to the filtration time and it is expected that after this period, $m_c\left(t\right)$ and $m_p\left(t\right)$ are approximately equal to their initial values. However, there is always a certain quantity of attached matter which may remain onto the membrane surface and/or blocked inside the pores, yielding irreversible fouling. Using hypotheses discussed hereabove, one can also neglect concentrations of detached matter returning to the reactional medium. It is important to emphasize that all (or part of) model parameters values must be adjusted using experimental data. In the next section, we choose arbitrary values of parameters for properties numerical investigation of the system (7)-(11) and (13)-(14).

3. Simulation Results

3.1. Coupling the Membrane Fouling Model with the AM2b Model

The proposed integrated model combines a biological anaerobic model and the fouling model for an homogeneous bioreactor as illustrated in Figure 4. For the biological compartment, the anaerobic digestion model is not specified here. Unless the effects of SMP on fouling are neglected as in [33], neither the ADM1 model [4] nor the AM2 model [11] are good candidates since they do not include SMP dynamics. Instead, we suggest to use the AM2b model [12] which includes SMP dynamics and that has been precisely developed for control purposes. In any case, to couple it with the proposed fouling model, soluble and particulate matters must be related to state variables of the AM2b model. In addition, two functioning phases are considered: filtration and relaxation or backwash. By convention, $\Delta P>0$ holds for filtration; $\Delta P=0$ holds for relaxation and $\Delta P<0$ holds for backwash.

3.2. Investigating the Qualitative Behavior of the Model

In this section, we investigate numerically the model dynamics consisting in the coupling of equations (7)-(11) and (13)-(14) with the anaerobic digestion model AM2b ([12]). The AM2b model was developed to describe anaerobic digestion as a two step process, including dynamics of SMP. In the first step (acidogenesis) the acidogenic bacteria $X_1$ consume the organic substrate $S_1$ to produce Volatile Fatty Acids $S_2$ (VFA), $SMP$ and $C{O}_2$, while in the second step (methanogenesis), the methanogenic population $X_2$ consumes VFA and produces $SMP$, methane and $C{O}_2$.

The integrated model is given by equations (15)-(16). The set of parameter values used in the simulations is given in Table A1. Numerical simulations are performed for two cycles of filtration/relaxation and for three distinct combinations of parameters values $C_s$, $C_x$ and $C_{smp}$ (for simulation, their values are equal to: 0.1, 0.4 or 0.7). Quantities of soluble and particulate matters $C_s{S}_1$, $C_s{S}_2$, $C_{smp}SMP$, $C_x{X}_1$ and $C_x{X}_2$ attracted by permeation forces $Q_{out}/V$ are assumed fully attached onto the membrane surface. We consider that the relaxation period (5 min) is negligible compared to the filtration period (2h). During relaxation or backwash period, particulate and soluble matters partially re-injected into the bulk are taken into account in model (16) through the parameters $C_s^{\prime }$, $C_x^{\prime }$ and $C_{smp}^{\prime }$. Between two cleaning cycles, the bulk volume $V\left(t\right)$ is constant thanks to the volume dynamics equation and the flux balance given in model (15). When $Q_{out}$ decreases, then, $Q_{in}$ should decrease or $Q_w$ should increase or both should vary. Experimentally, both $Q_{in}$ and $Q_w$ are judiciously fixed in such a way we obtain an optimal ratio of organic matters (COD)/mixed liquor volatile suspended solids concentration (MLVSS)/day.

It should be noticed here that the precise adjustment of the filtration/relaxation (or backwash) periods in MBRs is an open problem of control where an optimal solution must be searched for. This problem is usually solved in applying short filtration sequences followed by relaxation/backwash periods in such a way the clogging is very limited. However, it must be underlined that such sequences are not optimized and are probably quite far from an optimal.

Coupled model for the filtration phase:

$$\begin{array}{ccc}{\dot{X}}_1\hfill & =& \left({\mu }_1\left(S_1\right)+{\mu }_{smp}\left(SMP\right)-k_{d1}-\frac{Q_w}{V}-\frac{Q_{out}}{V}{C}_x\right)X_1,\hfill \\ {\dot{X}}_2\hfill & =& \left({\mu }_2\left(S_2\right)-k_{d2}-\frac{Q_w}{V}-\frac{Q_{out}}{V}{C}_x\right)X_2,\hfill \\ {\dot{S}}_1\hfill & =& \frac{Q_{in}}{V}{S}_{1in}-\left(\frac{Q_{out}}{V}+\frac{Q_w}{V}\right)S_1-k_1{\mu }_1\left(S_1\right)X_1-\frac{Q_{out}}{V}{C}_s{S}_1,\hfill \\ {\dot{S}}_2\hfill & =& \frac{Q_{in}}{V}{S}_{2in}-\left(\frac{Q_{out}}{V}+\frac{Q_w}{V}\right)S_2-k_3{\mu }_2\left(S_2\right)X_2+\left(k_2{\mu }_1\left(S_1\right)+b_2{\mu }_{smp}\left(SMP\right)\right)X_1\hfill \\ \hfill & & -\frac{Q_{out}}{V}{C}_s{S}_2,\hfill \\ \dot{SMP}\hfill & =& -\left(\beta \frac{Q_{in}}{V}+(1-\beta )\frac{Q_w}{V}\right)SMP+\left(b_3{\mu }_1\left(S_1\right)+k_{d1}-b_1{\mu }_{smp}\left(SMP\right)\right)X_1\hfill \\ \hfill & & +\left(b_4{\mu }_2\left(S_2\right)+k_{d2}\right)X_2-\frac{Q_{out}}{V}{C}_{smp}SMP,\hfill \\ {\dot{m}}_c\hfill & =& Q_{out}(C_s{S}_T+C_x{X}_T+C_{smp}SMP),\hfill \\ {\dot{m}}_p\hfill & =& Q_{out}\left({\beta }_1.SMP+{\beta }_2{S}_T\right),\hfill \\ R\hfill & =& \alpha \frac{m_c}{A}+{\alpha }^{\prime }\frac{m_p}{ϵA},\hfill \\ A\hfill & =& \frac{A_0}{1+\frac{m_c}{\sigma }+\frac{m_p}{{\sigma }^{\prime }}},\hfill \\ Q_{out}\hfill & =& J.A=\frac{\Delta P.A}{\mu \left(R_0+R\right)},\hfill \\ Q_{in}\hfill & =& Q_{out}+Q_w,\hfill \\ \dot{V}\hfill & =& Q_{in}-Q_{out}-Q_w,\hfill \\ X_T\hfill & =& X_1+X_2,\hfill \\ S_T\hfill & =& S_1+S_2+SMP.\hfill \end{array}$$

(15)

Coupled model for the relaxation (or backwash) phase:

$$\begin{array}{ccc}{\dot{X}}_1\hfill & =& \left({\mu }_1\left(S_1\right)+{\mu }_{smp}\left(SMP\right)-k_{d1}+C_x^{\prime }\right)X_1,\hfill \\ {\dot{X}}_2\hfill & =& \left({\mu }_2\left(S_2\right)-k_{d2}+C_x^{\prime }\right)X_2,\hfill \\ {\dot{S}}_1\hfill & =& -k_1{\mu }_1\left(S_1\right)X_1+C_s^{\prime }{S}_1,\hfill \\ {\dot{S}}_2\hfill & =& -k_3{\mu }_2\left(S_2\right)X_2+\left(k_2{\mu }_1\left(S_1\right)+b_2{\mu }_{smp}\left(SMP\right)\right)X_1+C_s^{\prime }{S}_2,\hfill \\ \dot{SMP}\hfill & =& \left(b_3{\mu }_1\left(S_1\right)+k_{d1}-b_1{\mu }_{smp}\left(SMP\right)\right)X_1+\left(b_4{\mu }_2\left(S_2\right)+k_{d2}\right)X_2+C_{smp}^{\prime }SMP,\hfill \\ {\dot{m}}_c\hfill & =& -\omega m_c+m_{c_{irr}},\hfill \\ {\dot{m}}_p\hfill & =& -{\omega }^{\prime }{m}_p+m_{p_{irr}}.\hfill \end{array}$$

(16)

Simulation results are reported in Figure 5, where we have plotted the dynamic evolution of the attached mass $m_c\left(t\right)$ on the membrane surface, the blocked soluble matter $m_p\left(t\right)$ (SMP in the majority) inside the pores, the fouling resistances $R_c\left(t\right)$, $R_p\left(t\right)$ and $R\left(t\right)$, the output flow rate $Q_{out}\left(t\right)$, the permeate flux $J\left(t\right)$ and the membrane surface $A\left(t\right)$. Dynamic responses are simulated for three different combinations of parameters values $C_s$, $C_x$, $C_{smp}$ and ${\beta }_1$, to emphasize effects of deposited and blocked matter rates on the fouling dynamic. These rates depend on many parameters as concentrations of soluble and particulate matters, characteristics of mixed liquor and its viscosity or still temperature and matters specific capability to contribute to fouling.

The trajectories of the main variables are plotted in the case of rapid and strong fouling due, for example, to a high concentration of solid matter. In such a case, if we define a threshold flux $J_s$, over which the process can operate, then the process will be stopped very often and be switched in relaxation or backwash phases (at $t_{f1}$ for the first operating cycle). Dashed plots correspond to a slower and softer fouling: the slower the fouling, the longer the time period ($t_{f3}$) during which the process may operate without switching in a relaxation or backwash mode. Such simulations show that $C_s$, $C_x$, $C_{smp}$ and ${\beta }_1$ may be adjusted to match a large range of experimental data.

During the first minutes of the filtration process, the fouling is fast and significant. All variables have fast dynamics (increasing or decreasing) at the beginning and then attain gradually their equilibrium. At steady state, quantities of the attached mass $m_c\left(t\right)$ on the membrane surface (around $15g$ per $1m^2$ of area), and the SMP mass blocked into the pores $m_p\left(t\right)$ (around $8g$) are sufficient to cause membrane fouling. Resistances of pore blocking $R_p\left(t\right)$ and cake formation $R_c\left(t\right)$ are within the order of ${10}^{11}$. The useful filtering surface area $A\left(t\right)$, the output flow $Q_{out}\left(t\right)$ and the permeate flux $J\left(t\right)$ decrease significantly, especially during the first minutes of filtration as it is often the case in practice. During the relaxation (or backwash) phase, we have an exponential decreasing of $m_c\left(t\right)$ and $m_p\left(t\right)$, with detachment from the membrane surface and pores (shear-induced diffusion). We notice that even after the relaxation (or the backwash), the permeate flux $J\left(t\right)$ and the membrane area $A\left(t\right)$ are not equal to their initial values (see Figure 5), because of the irreversible fouling taken into account in the model. From cycle to cycle, the total resistance $R\left(t\right)$ has increasing equilibrium values. After several cycles, it will be necessary to clean chemically the membrane or to change it. Regarding these qualitative results, we claim that if a good and accurate prediction of fouling can be achieved using the model, new fouling control strategy can be proposed and filtration-relaxation (or backwash) sequences can be optimized.

4. Experimental Results

In this section, we are interested in the calibration of the fouling model (7)-(11), by using experimental data collected at the Center of Biotechnology of Sfax, Tunisia.

4.1. Pilot Plant (AnMBR) and Data Used for the Model Validation

The schematic representation of the pilot plan is shown in Figure 6. It was installed in the “Center of Biotechnology of Sfax”, Tunisia and, it was used for the treatment of municipal wastewaters. The system is composed of an anaerobic bioreactor coupled with an ultrafiltration membrane module, with a filtering area of 1 $m^2$ and 100 $kDa$ cut-off. The working volume of the bioreactor is 50 L and its temperature is maintained constant at 37 ^oC. The cross-flow velocity is fixed at 3 $m/s$ and the transmembrane pressure may vary until reaching 2 $bars$. The membrane was frequently chemically cleaned, approximately every 40 days, in order to maintain an acceptable flux (see Figure 8). The cleaning step was performed for 1 h at 35 ^oC and was followed by water cycle. For more details on the experimental process and analytical methods, the reader is referred to [33] and [34].

Two cycles of filtration/backwash are considered over 84 days in total: the first cycle for the period [t=1..39 days] and the second one for the period [t=41..84 days], with a backwash at the 40^th day (Figure 8). Collected data of the total COD ($S_T$ and $SMP$) and the total biomass ($X_T$) are represented in Figure 7 (red markers), while data for the permeate flux are shown in Figure 8 (black dots).

4.2. Experimental Identification and Validation of the Fouling Model

4.2.1. Parameters Estimation Procedure

In the present case, we apply the least-squares method for the parameters estimation of the model (7)-(11). The objective function (17) is minimized in adjusting model parameters such that the simulated flux $\tilde{J}\left(t\right)$ best fits the experimental flux $J\left(t\right)$. A nonlinear optimization algorithm was used (functions “fmincon” and “ode45” of matlab).

$$F=\sum _{i=1}^N{\left(J_i\left(t\right)-{\tilde{J}}_i\left(t\right)\right)}^2.$$

(17)

where N is the number of measurements.

Before running the optimization algorithm, we must fix the values of some constants of the model, as the intrinsic resistance membrane $R_0$ of equation (). It can be estimated from the initial value of the measured flux $J\left(0\right)$, as long as $R\left(t\right)$ is still negligible by (18):

$$R_0=\frac{\Delta P}{\mu J\left(0\right)}.$$

(18)

For simplicity, we fixed arbitrarily some parameters values of the model (7)-(11), notably $ϵ$ ($0<ϵ<1$), $\sigma$ and ${\sigma }^{\prime }$ as given in Table 1. Parameters defined in the literature as $\mu$, $\alpha$ and ${\alpha }^{\prime }$ have defaults values reported in Table A1.

Model parameters to be estimated by the least-square method are: $C_s$, $C_x$, $C_{smp}$ and ${\beta }_1$, (${\beta }_2$ is chosen smaller than ${\beta }_1$, for instance ${\beta }_2={\beta }_1/15$). They are considered as key parameters used to describe the rate of membrane fouling.

Our identification procedure consists in using a first part of data of the flux (16 measures for $t=1..39$, see Figure 8 on the left) to estimate parameters $C_s$, $C_x$, $C_{smp}$ and ${\beta }_1$ and, then the last part (18 measures for $t=41..84$, see Figure 8 on the right) to validate the model. Data of total $COD$ ($S_T$ and $SMP$) and total biomass $X_T$ (Figure 7) are inserted as inputs for the identification algorithm of the model (7)-(11), such that $S_T$ and $SMP$ are assumed to be as follows:

$$S_T=85\%CODandSMP=15\%COD.$$

Let us emphasize here that this is a simplifying assumption based on what is proposed in the literature, and that in practice these proportions can change with time and according environmental conditions ([35,36,37]). Here it is just essential that the sum of all soluble matters equals the $COD$ and that the sum of all particulate matters equals the $X_T$. Also, we notice that times and frequencies of $COD$ and $X_T$ measurements are different (29 measures of $COD$ and only 13 measures $X_T$). To solve this problem, an interpolation of $COD$ and $X_T$ is performed by adding more intermediate points as illustrated in Figure 7 (blue markers).

4.2.2. Results and Discussion

Simulation and experimental results are plotted in Figure 8. On the left, the simulated flux $\tilde{J}\left(t\right)$ (solid line) is compared with the real measured flux $J\left(t\right)$ (red dots) used for calibration. There is a good matching of the model simulations and real data (measures for $t=1..39$) and an accurate correlation coefficient $r^2=0.968$ is obtained. On the right, the identified model is validated on the second data set (measures for $t=41..84$), with a satisfied correlation coefficient $r^2=0.938$. Table 2 shows dimensionless parameter values that have been estimated by minimizing the criterion (17).

Figure 9 shows the dynamic evolution of the estimated solids mass $m_c\left(t\right)$ attached on the membrane surface and the estimated soluble matter $m_p\left(t\right)$ deposited inside the pores. $m_c\left(t\right)$ and $m_p\left(t\right)$ increase significantly during the first days of filtration process before their slopes decline when the membrane becomes more and more heavily fouled. At the end of the first filtration period (day 39), the estimated mass $m_c\left(t\right)$ is around of 43 $g/m^2$ and it is about 45 $g/m^2$ at the end of the second filtration period (day 84). Estimated quantities of $m_p\left(t\right)$ are about $7.4$$g/m^2$ and $7.75$$g/m^2$ at the end of the first and the second filtration periods respectively. We emphasize that these values are just estimations of $m_c\left(t\right)$ and $m_p\left(t\right)$ by numerical simulation during the calibration of the model (7)-(11) and that, their experimental values may be probably different from those simulated. Of course, if one has experimental measures of $m_c\left(t\right)$ and $m_p\left(t\right)$, then one can use them to calibrate accurately more parameters of the model.

Figure 10 represents the experimental and the simulated total fouling resistance $R\left(t\right)$ due to both $m_c\left(t\right)$ and $m_p\left(t\right)$. Experimental resistance values are deduced from flux data using equation (1). Resistance $R\left(t\right)$ during the two filtration cycles has the tendency to increase rapidly in the first period, before attaining values of about $12e13$ and $14e13$ at the end of the calibration and the validation period respectively. If the system is functioning with a long-term filtration, then $R\left(t\right)$ should converge towards a quasi-constant values (its equilibrium).

The proposed model is simple from the mathematical point of view and it reproduces quite well the fouling behavior of the AnBRM: it can then be used for control purposes. However, the estimated parameters values and/or those arbitrarily selected could change with time. Thus, it is necessary to readjust them regularly. For example, one could re-identify on a regular basis parameters in order to best fit the model for the last available measurements (for instance using the last data for $t=41..84$ days, see Figure 8, on the left). Also, if we decide to estimate the parameters of the biological model AM2b, then we will need more informative measurements.

Figure 8. Left: Data used to calibrate model parameters and simulated response of flux $J\left(t\right)$. Right: Validation of the model using a second set data.

Figure 8. Left: Data used to calibrate model parameters and simulated response of flux $J\left(t\right)$. Right: Validation of the model using a second set data.

Figure 9. Left: Simulated solids mass $m_c\left(t\right)$ attached onto the membrane surface. Right: Simulated solute $m_p\left(t\right)$ deposited inside the membranes pores.

Figure 9. Left: Simulated solids mass $m_c\left(t\right)$ attached onto the membrane surface. Right: Simulated solute $m_p\left(t\right)$ deposited inside the membranes pores.

Figure 10. Experimental and simulated resistance $R\left(t\right)$ during, Left: calibration period, Right: validation period.

Figure 10. Experimental and simulated resistance $R\left(t\right)$ during, Left: calibration period, Right: validation period.

5. Discussions, Open Questions and Perspectives on the Process Control Using the Proposed Model

Membrane fouling is the major drawback of MBRs and one important challenge is to propose new control strategies to minimize fouling and improve treatment efficiency. In particular two important questions must be addressed:

• What parameters mainly influence membrane fouling? This is basically a modeling question, and,
• How minimizing fouling (filtering conditions)? which can be seen as a control problem as soon as a model describing the fouling dynamics is available.

If several fouling models have already been proposed in the literature, very few were used for minimizing fouling and optimizing treatment. One of them is recently proposed in [30], where one used a model of SMBR with slow/fast dynamics to develop a nonlinear predictive control. However very often, the control strategies are tuned heuristically and use available process actuators:

• Gas sparging: It consists in injecting bubbles (air for aerobic process or biogas for anaerobic systems) for membrane scouring in order to limit attachment and promote detachment of matter by shear forces. This control parameter is however costly because it consumes energy.
• Intermittent filtration: MBR is operated in alternating filtration/relaxation cycles. This functioning mode allows the detachment of matter responsible for the reversible fouling.
• Backwash: It must be used for short time compared with filtration time to detach matter involved in irreversible fouling and it is costly because it also needs energy.

One common question for all these techniques is to find a good control sequence ensuring good process performances while minimizing membrane fouling. In practice, it can be seen as an optimal control problem, since we need to optimize the filtration flux, the filtration time or still the energy consumption. A practical problem could be: what is the optimal sequence for intermittent filtration or for filtration/backwash cycles? What is the optimal operating time and mode for bubbles injection?

A study was performed in [38] with the final aim to reduce by different strategies the energy costs in MBR. In particular, authors investigated the influence of the aeration intensity, the duration of filtration/backwashing cycles, and the number of membrane cleanings on the MBR energy demand. However, the used model is integrated and complicated, which it divided into a biological sub-model (19 biological state variables and 79 parameters) and, a physical sub-model (membrane model). In the following, we investigate in simulation the influence of the filtering parameters mentioned above on the flux production and process performances, by using the simple model proposed in this paper.

5.1. Influence of Gas Sparging

In this section, we investigate how gas sparging can be used for limiting membrane fouling. To do so, we need to modify the proposed model (7)-(11) in adding negative terms on the right sides of equations (7) and (8). This way, the reversible and irreversible fouling rates are reduced by gas sparging as illustrated by equations (19) and (20). Functions $f\left(m_c\right)$ and $g\left(m_p\right)$ are positive and represent the controller effect on the detachment of matter (gas sparging, membrane scouring). In some fouling models, these terms are simply constants but modulating their magnitude, our idea is to add them here as control parameters (as mentioned above in the modeling section).

$$\begin{array}{ccc}\hfill {\dot{m}}_c& =& Q_{out}(C_s{S}_T+C_x{X}_T+C_{smp}SMP)-f\left(m_c\right),\hfill \end{array}$$

(19)

$$\begin{array}{ccc}\hfill {\dot{m}}_p& =& Q_{out}\left({\beta }_{1}SMP+{\beta }_2{S}_T\right)-g\left(m_p\right).\hfill \end{array}$$

(20)

In other terms, efficient control consists to propose functions $f\left(m_c\right)$ and $g\left(m_p\right)$ depending on the intensity of gas sparging (parameter control). A first simple form of $f\left(m_c\right)$ and $g\left(m_p\right)$ which is already used in the literature is $k_m{m}_c$ and $k_p{m}_p$, which represent quantities of $m_c$ and $m_p$ detached by shear forces caused by membrane scouring, where $k_m$ and $k_p$ depend on the intensity of injected bubbles used to detach fouling [27,39,40]. A higher aeration intensity can have a positive effect on cake layer removal by shear force and thus improves the membrane permeability [41]. Figure 11, illustrates the time evolution of the flux $J\left(t\right)$ with respect to different values of $k_m$ (here $k_p=0$, it is assumed that the irreversible fouling detachment is neglected, since it is not significantly affected by gas sparging). It can be seen that $m_c\left(t\right)$ and $R_c\left(t\right)$ are inversely proportional to the control parameter $k_m$, for higher values of this later, accumulated matter on membrane surface and its corresponding resistance take small values. The output flow $Q_{out}\left(t\right)$ and permeate flux $J\left(t\right)$ are increasing proportionally to $k_m$. This is a classical result, but the given question is, how to optimally calibrate $k_m$ (and $k_p$) in order to best control the fouling with minimum energy consumption?

On the other hand, one sees on Figure 11, that deposited matters $m_p\left(t\right)$ inside the pores and its relative resistance $R_p\left(t\right)$ are proportional to $k_m$ and inversely proportional to $m_c\left(t\right)$. If the value of this parameter increases, then the quantity of $m_p\left(t\right)$ and the value of $R_p\left(t\right)$ increase likewise leading to a flux loss especially at the end of the filtration time (around steady-state). One can explain this result as follows: if the formed cake layer ($m_c\left(t\right)$) represents a second biological membrane (which prevents pores from fouling ($m_p\left(t\right)$) as it is known in the literature [25,42,43]) then when this layer is detached by gas sparging, more particles of different sizes go through pores and cause further fouling. So, a second question of such a control strategy may be asked: how can we control and favour the cake formation until acceptable level, to protect pores from fouling, but at the same time, without influencing permeate flux? This question, actually, remains open.

5.2. Influence of the Number of Filtration/Relaxation (Backwash) Cycles per Time Unit

Given a sufficiently large time horizon, what is the optimal number of filtration/relaxation or backwash cycles allowing a higher mean value for the MBR output flux (Figure 12). To illustrate the importance of this operational functioning mode, we performed numerical simulations by changing the number of filtration/relaxation cycles over a given functioning period with a constant ratio between filtration time and relaxation time ${\alpha }_t=T_{filtr}/T_{Relax}$ for all cycles. On Figure 13, results are given for: one cycle, two cycles, five cycles and ten cycles of filtration/relaxation with a period of 2 hours and ${\alpha }_t=7$. We are particularly interested by the mean value $J_{mean}$ of the produced flux on the given period. Using the simulations, we computed:

• for 1 cycle of filtration/relaxation: $T_{filtr}=105mn$, $T_{Relax}=15mn$ and $J_{mean}=17.9$$L/(h.m^2)$,
• for 2 cycles of filtration/relaxation: $T_{filtr}=52.2mn$, $T_{Relax}=7.5mn$ and $J_{mean}=22.9$$L/(h.m^2)$,
• for 5 cycles of filtration/relaxation: $T_{filtr}=21mn$, $T_{Relax}=3mn$ and $J_{mean}=29$$L/(h.m^2)$,
• for 10 cycles of filtration/relaxation: $T_{filtr}=10.5mn$, $T_{Relax}=1.5mn$ and $J_{mean}=31.5$$L/(h.m^2)$.

If the objective is to produce a maximum flux over the given period, then 10 filtration/relaxation cycles appears to be the best strategy. As a matter of fact, for a reasonable number of filtration/relaxation cycles and a constant ratio ${\alpha }_t$ between filtration time and relaxation time for all cycles, we have an increasing mean production of flux. But if the number of intermittent filtration cycles is too large on the considered functioning period, then it can damage the process by forcing it to operate in an On/Off mode with a high frequency. Simulations show that multiplying the number of cycles is beneficial but that this benefit does not increase anymore beyond 10 cycles. Obviously, increasing the number of filtration/relaxation cycles is useful but the higher the frequency, the lower the benefit.

It is thus suggested not to wait too long before proceeding to the membrane cleaning by relaxation or backwash to find the best ratio operated time by benefit in terms of flux produced.

5.3. Coupling Sparging Gas and Intermittent Filtration Controls

Figure 14 illustrates a control strategy based on gas sparging which is used at the beginning of the considered period when the flux is still higher than a threshold flux together with intermittent filtration as soon as the flux has reached the threshold flux.

Our idea here is to minimize the energy consumption when using gas sparging and the flux loss (resp. the permeate loss) when the process is in relaxation mode (resp. backwash). In others words, instead of using gas sparging and intermittent filtration simultaneously, we propose to use them sequentially for the following reasons:

• Gas sparging is used to detach the matter deposited on the membrane: this phenomenon occurs at the beginning of the filtration (fouling is soft and not yet dense). Here, one should control the gas sparging intensity, which may depend on different parameters as the mixed liquor characteristics, the concentration of soluble and particulate matters, ...
• Intermittent relaxation is used to detach a denser fouling (strong), which can occur after an enough long functioning time. These control parameters (typically the number and frequency of filtration/relaxation (backwash) cycles) may depend on the fouling characteristics as its irreversibility, its thickness...

To illustrate this idea, we performed numerical simulations plotted in Figure 14. The system is first simulated without any control (black line). Then this reference scenario is compared with the proposed control strategy. It means that gas sparging is first applied until the flux reaches the threshold flux (here $J_s=18$$L/(h.m^2)$). At this instant ($t=0.64h$), we apply intermittent control with $k_m=5$ in the equation (19), where $f\left(m\right)=k_{m}m$ with 4 cycles.

Simulations show that this control strategy allows one to increase the mean production flux to 33 $L/(h.m^2)$, whereas the mean flux without control was 18 $L/(h.m^2)$. As it is noticed in Figure 14, when applying the gas sparging control, it has increased favorably the permeate flux on the control period (until 0.64 hour). It should be noticed that even if we applied only the gas sparging all along the functioning period (without using intermittent filtration cycles), the mean flux is $28.76$$L/(h.m^2)$, lower than the produced flux when the two techniques are used together. Thus, intermittent filtration was an appropriate control strategy to obtain over the whole functioning period a maximum of flux.

Our study on control strategy is obviously inline with other studies as the work presented in [38]. Their main purpose was to investigate and select the best operating conditions in terms of aeration intensity, duration of filtration/backwashing cycles and number of membrane cleaning to optimize energy demand and operational costs. Thus, interesting further studies may focus on the optimization of the mean production flux, while saving energy consumption and costs to produce a good effluent. The first results we established here show that this could be achieved by using the simple control model, which integrates biological and membrane compartments of the MBR and evaluating on the considered time unit an appropriate criteria, as a function of the mean production flux, the energy demand and the operational cost. This work is under investigation.

We emphasize that in [44], parameters of the developed model (7)-(11) were identified using data generated by models proposed in [25] and [30] considered as virtual processes. It was shown that our generic model can capture important properties of these two models, as the mean value of the transmembrane pressure and the attached mass on the membrane and their dynamics.

6. Conclusion

In this paper we proposed a simple fouling model of Anaerobic Membrane BioReactor (AnMBR). The model was developed under certain classical hypotheses on the membrane fouling phenomena and was coupled with a reduced order anaerobic digestion model. Two mechanisms of fouling were considered, cake formation on the membrane surface and pores blocking. Contrary to many models of the literature which consider constant membrane surfaces A, we assumed that this later is a decreasing function of both the attached matter $m_c\left(t\right)$ and the deposited matter $m_p\left(t\right)$ (notably SMP). Our main idea was to introduce in the mass-balance model AM2b [12] a feedback of the decreasing flux due to membrane fouling into the actual output flow rate of the process. We performed simulations to investigate the qualitative behavior of the model and we validated it on experimental data. It was shown that the proposed model can predict quite well the fouling behavior for the considered AnMBR. It fitted accurately the real flux in both phases of identification: calibration phase and validation phase (see Figure 8). In a second part of the paper, we discussed the fouling control problem in focusing on the optimal control of gas sparging and intermittent filtration. Preliminary results were obtained about the results of different control strategies over a given time period: at the beginning stage of the process functioning, it appeared useful to use the gas sparging and the intermittent filtration at the end of the considered time period. Based on these results, we proposed to couple control benefits in order to produce the maximum mean flux over the total considered functioning period.

Perspectives of this work include (i) model extension to constant flux and variable TMP filtration and, its validation using other experimental informative data, where parameters identifiability should be studied, (ii) the design of an optimal control to minimize the fouling effects on the system performances while minimizing the energy requirements, (iii) the development of an optimization strategy to control and favour the cake formation until acceptable level, to protect pores from fouling, but at the same time, without influencing permeate flux and, (iv) the coupling of the fouling model with ADM1 model in integrating SMP and possibly EPS (Extracellular Polymeric Substances) dynamics.