Impact of Dual Substrate Limitation on Biodenitriﬁcation Modeling in Porous Media

: In this work, we consider a model of the biodenitriﬁcation process taking place in a spatially-distributed bioreactor, and we take into account the limitation of the kinetics by both the carbon source and the oxidized nitrogen. This model concerns a single type of bacteria growing on nitrate, which splits into adherent bacteria or free bacteria in the liquid, taking all interactions into account. The system obtained consists of four diffusion-convection-reaction equations for which we show the existence and uniqueness of a global solution. The system is approximated by a standard ﬁnite element method that satisﬁes an optimal a priori error estimate. We compare the results obtained for three forms of the growth function: single substrate limiting, “multiplicative” form, and “minimum” form. We highlight the limitation of the ‘ single substrate limiting model”, where the dependency of the bacterial growth on the nitrate is neglected, and ﬁnd that the “minimum” model gives numerical results closer to the experimental results.


Introduction
The biodenitrification process (degradation of nitrite and nitrate into gaseous nitrogen) is realized by heterotrophic microbial ecosystems. In the absence of oxygen, such ecosystems use oxidized nitrogen (NO 2 and/or NO 3 ) instead of oxygen as an electron acceptor while they need an organic carbon source for their growth. Mccarthy [1] and Payn [2] described the process of biodenitrification as a respiratory process in which certain bacteria (so-called denitrifying bacteria) use nitrates instead of oxygen as an acceptor of electrons, intended to provide energy for cell activity and the synthesis of new cells. The process generally takes place under conditions called anoxic, i.e., when the dissolved oxygen is replaced by another electron acceptor. From a thermodynamical viewpoint, nitrates are the best acceptor of electrons that can replace oxygen so that they should be considered in the modeling as a limiting compound. Most standard models of microbial growth in laboratory bioreactors, such as the chemostat or the piston flow reactor, take into account the tendency of bacteria to adhere to surfaces and thus form biofilm (cf. [3,4]); however, such models neglect the possible diffusion of attached biomass. In reality, the bacterial population consists of cells suspended in the fluid termed planktonic or free cells, and cells adhering to the surface, termed adherent cells. At any time, planktonic cells can adhere to the walls forming biofilms, while adherent cells can detach from the biofilm (due to erosion and sloughing) and move into the planktonic cell compartment. Earlier (see [5]), we considered a model of the biodenitrification process taking place in a spatially-distributed bioreactor with a single type of bacteria growing on nitrate and that splits into adherent and free bacteria in the liquid, taking all interactions into account. We considered that the growth function µ(·) depends only on the nutrient concentration S following the Monod law: so the nitrates' concentration, S N , was not considered as a limiting substrate. Generally, in the case of the existence of two limiting substrates, the growth function µ(·) can take various forms (see [6]) depending on the relationship among its nutrients. Among the most frequently used forms, we cite two formulas that are used when the two limiting resources are both essential (i.e., are needed). The first, called the multiplicative formula, is given by: and the second, called the minimum formula, is given by: where S 1 , S 2 are the two limiting substrates and K 1 , K 2 are respectively the associated half-saturation indices. Charpentier, Ch. et al. [7] (2008) used a slightly different multiplicative formula. Stewart, H.A. et al. [8] (2017) compared three dual limitation models (multiplicative, minimum, and Bertolazzi) based on experiments considering two bacteria types, where the growth of the first one is limited by dissolved oxygen and nitrite, whereas the growth of the second by ammonium and nitrite.
In this work, we provide a comparison between the three formulas. We will first study the limit of the model with µ(·) given by (1), by comparing it to the case where nitrates are considered as a limiting substrate; it emerges that this first model remains valid up to a threshold beyond which the results are no longer valid. Currently, the multiplicative model is the most commonly used as it is continuous, smooth, and easy to handle in numerical simulations. We will compare the two Formulas (3) and (2), especially in the presence of the diffusion, which pertains to our case, and show, via numerical tests, that (3) gives better results than others reported in the literature.
In the second section, we recall the mathematical model represented by a non-linear coupled system of four equations, before introducing the hypotheses. In the third section, we give the analysis of the existence of solutions and the approximated problem by a standard finite element method. In the last section, some numerical tests are presented where the advantage of the growth function (3) is highlighted by comparisons with previous simulations obtained with (1) and (2).

Mathematical Model
Let Ω be an open set of R 2 with a regular enough boundary ∂Ω = Γ, which is divided into three parts Γ 1 , Γ 2 and Γ 3 : We suppose that Ω contains nitrified water, denitrifying bacteria, and a nutrient. The flow of the nitrified water comes from Γ 1 and goes through Γ 2 ; we assume that the flow is permanent with a velocity u. The impermeable part of the boundary is Γ 3 . In the reactor, the bacteria will be divided into two categories: those that adhere to the walls of the reactor to form a biofilm, which is assumed to be a monolayer, and those that remain mobile and free in the medium, called planktonic cells. For a given T, let the space-time domain defined by: We denote by b m the density of the mobile bacteria and by b f the surface density of bacteria adhering to the walls, with a maximum value denoted by w ∞ . The proportion of the occupation of the wall is a number between zero and one defined by the ratio: The function µ i (·) is the growth rate of b i , for i = f , m, and S is the concentration of the limiting substrate. According to [5], the equation modeling the evolution of adherent bacteria is given by: where: • k f denotes the adherent bacteria mortality rate; • β denotes a term corresponding to the detachment rate from the wall.

•
A portion of the mobile category can attach to the walls with a certain rate that is denoted by α. • G(b f ) denotes the proportion of daughter cells of the fixed bacteria able to find a place to attach onto the wall, the remainders being washed out by the liquid flow (cf. [9]). • γ is the coefficient of conversion of the volume density to the surface density.
For the mobile bacteria, the equation modeling their evolution according to [5] is given by: where k f denotes the free bacteria mortality rate.
The degradation of the nutrient S and of the contaminant S N is governed by the equations (cf. [10,11]): where: • Y i , for i = m, f , is respectively the coefficient rate of yield of mobile and fixed bacteria, defined as the ratio of the bacterial mass produced (in g or mol) by the mass of the substrate consumed (in g or mol), • R is the rate of the degradation of nitrates.
We will give now some comments about the coefficients G and µ i for i = f , m. According to the form considered by Freter ( [9,12]), let: where a is a small positive number. The most used expression for the growth rate is the Monod law [13] given by (1). In the present work, we take into account the fact that the bacteria need both carbon and nitrate to grow. We consider the growth rates depending on two limiting nutrients S and S N according to the Formulas (3) and (2). Equations (4)-(7) are now considered with the growth rates µ i (·), for i = m, f , given in (2) for the analysis part. In the numerical simulations, we consider a comparison between the Formulas (3) and (2).
The functions G and µ i (·) in (1) and (2) satisfy the following hypothesis: Hypothesis 1. The growth rate of bacteria µ i (S, S N ), for i = f , m, satisfies: Hypothesis 2. The function G satisfies: The rate of change of concentrations due to diffusion is represented by terms of the form D∇c, where c is the concentration or density and D is the diffusion coefficient. The transport is represented by terms of the form cu, where u is the velocity. We obtain a system of four equations, defined in Q T . The first one is a reaction-diffusion equation, and the others are reaction-convection-diffusion equations in which the transport velocity is u. Equations (4)-(7) become in the time-space variables: In order to describe the time evolution of the variables b f , b m , S, and S N completely, we have to specify the behavior of these variables on the boundary of the domain. Let n be the outward normal vector to the boundary Γ. By the definition of fixed bacteria, there is no flux across the entire boundary of the domain: In order to maintain the growth of the bacteria, we continuously inject from Γ 1 the substrate with the density S in . The boundary conditions used on Γ 1 , which models flow continuity, assuming that the concentrations are uniform outside, are Robin's type, also called Danckwerts [14]: At the output, the flow of S is uniform, and the condition on Γ 2 is then given by: On the impervious part Γ 3 , it is natural to consider: For b m , similar reasoning leads to the following boundary conditions: The nitrified water contained in Ω comes from Γ 1 with a velocity flow u and is withdrawn with the same flow from the part Γ 2 of the boundary. The concentration of nitrates thus verifies the following boundary conditions: To close this system, we define the initial conditions. At t = 0, the domain Ω contains nitrified water with an initial density S 0 N , planktonic bacteria with a density b 0 m , adherent bacteria with a density b 0 f , and a carbon source with an initial density S 0 : We suppose that:

Analysis and Approximation
The usual space L p (Ω), for 1 ≤ p < ∞, is defined by: Let H 1 (Ω) be the usual Sobolev space of the first order defined by: , and H −1 (Ω) be its topological dual. From now on, we adopt the following notations: We put: The global fluxes are defined by: The boundary operator will be denoted by B : Let g := (0, 0, g 3 , g 4 ) with: With these notations, the quasilinear diffusion-convection system (8)- (20) becomes:

Existence Theorem
Since C is in (L ∞ (Ω)) 4 and F i is in C 1 (R 4 ) (in the case when we use the Formulas (1) or (2)) for 1 ≤ i ≤ 4, the classical local existence result is well known (see [15]). There exists T > 0 and a unique classical solution of (22)  The second essential property of the reaction term F is the "triangular" structure of F defined in the following lemma. We recall that in our problem, F does not depend explicitly on (t, x) ∈ Q T .

Lemma 2.
There exists a lower triangular invertible matrix M = (m ij ) 1≤i,j≤4 with non-negative diagonal entries and a vector V in R 4 + such that: (where the inequality stands component by component).

Proof of Lemma 2.
In the definition of F 1 and F 2 , the coefficients of the variables b m and b f are both positive, and we can consequently write, for all x ∈ Ω: Since S and S N are positive and F 3 and F 4 are negative, we obtain for i = 1, 2, 3, 4: We take M as the identity matrix to obtain (23).
The third essential property considered for the reaction term F is the "polynomial growth" structure of F defined in the following lemma.
where | · | is the euclidean norm in R 4 .

Proof of Lemma 3.
This is a consequence of (25) by taking p = 1 and M = max Moreover, for any T > 0, there exists: M > 0 depends on T, c 0 , u, and V defined by (24) such that:

Approximation
By taking the test function Ψ in (D(Q T )) 4 , then in (D(Q T )) 4 , and using some integration by parts, it is standard to see that the weak solution of the last theorem is a solution of the system (22) in (D (Q T ) 4 ) (the distribution space). Now, to give an approximation of this solution by a finite element method, we consider the following formulation obtained from (22) by integration by parts. It reads: In order to have a (H 1 (Ω)) 4 -elliptic bilinear form, we transform the problem by using the following augmented bilinear form: The problem is then transformed as follows. Let Φ = (φ 1 , φ 2 , φ 3 , φ 4 ), with φ i = exp (−γt)c i , and letg i = exp (−γt)g i and G i (Φ) = exp (−γt)F i (exp (γt)Φ). The problem (29) is equivalent to the system: for a.e. t ∈]0, T], The discretization of the problem (31) is based on two steps: the space discretization (or semi-discretization), which is made by a finite element method of degree one, and the time discretization, which uses the backward Euler scheme.

Semi-Discretization
Let T h = T T be a family of regular triangulations of Ω, where T is a triangle and h T its diameter.
We denote h = max T∈T h diam(T). let P 1 (T h ) and V h be defined respectively by: where, for each T ∈ T h , P 1 (T) stands for the space of restriction to T of polynomial functions of degree one. The semi-discrete problem associated with (31) is given, for a.e. t ∈]0, T], by: where Π h is a projection operator on V h . Lett ∈]0, T]. To give an a priori error estimate for a local solution of the semi-discretized problem (32), we consider a finite time intervalJ = (0,t]. The error estimates in the L 2 -norm are based on the elliptic projection: for all t ∈ [0,t] and for all 0 ≤ i ≤ 4, let ω ih (t) be the elliptic projection of the exact solution φ i on P 1 (T h ) defined by: whereã i (·, ·) is defined in (30). Following the technique of Thomée [17], we can prove the following theorem.

Theorem 2.
Let Φ and Φ h be the solutions of (31) and (32), respectively, and let Π h Φ 0 be the elliptic operator defined by (33) for t = 0. We assume that u ∞ < 2γ. Then, we have: where L(Φ) is a non-negative constant and depends on t and the solution Φ, but independent of h.

Full Discretization
We consider the discretization of [0, T] given by 0 = t 0 < t 1 < ... < t N = T and put τ n = t n − t n−1 and τ = max 1≤n≤N τ n . For all k, 0 ≤ k ≤ N, and all i, 1 ≤ i ≤ 4, we use the notation: The derivative with respect to time is approximated by the backward Euler scheme given by the following difference quotient: To linearize the reaction term and have decoupled equations, we consider at time t = t n the reaction term at time t = t n−1 . The system (32) is then fully approximated by the following implicit Euler scheme: The error estimate for the fully approximated problem is given in the following theorem. The proof follows also the technique of Thomée. Theorem 3. Let Φ n h and Φ be the solutions of (35) and (31), respectively. Then, there exists a constantM depending on Φ, but independent of h and τ such that: Remark 1. The estimates of Theorem 3 remain valid for C h since we have Φ h = exp(−γt)C h . Therefore, in the numerical tests, we use the following approximated problem:

Numerical Tests
In this section, we consider a domain Ω ⊂ R 2 defined by Figure 1). We assume that Ω is a porous medium, hereafter denoted the "reactor" Ω, into which we inject some water containing nutrients, sources of both nitrate and a carbon with concentrations S in N and S in , respectively. In this case, the flow in Ω is governed by Darcy's equation, given in its mixed form by the system: find (u, p) ∈ H(div; Ω) × L 2 (Ω), with u · n = u0 on Γ 1 ∪ Γ 2 such that: where H(div; Ω) := v ∈ (L 2 (Ω)) 2 /div v ∈ L 2 (Ω) and V := {v ∈ H(div; Ω)/v · n = 0 on Γ 1 ∪ Γ 2 }, K is the hydraulic conductivity of the medium, u0 is the given flux on Γ 1 , and p D is the given pressure on Γ 2 . This problem will be approximated by the mixed finite element method of Raviart-Thomas with the lowest order (see [18]). The resolution of this system gives the velocity u, which is needed to resolve our system (37). All the resolutions will be done with the FreeFem ++ software [19] with the mesh of Figure 2 and the curves drawn using MATLAB. We consider the initial data given in [20] and the parameters of the system (22) given in [7]. These parameters are summarized in Table 1.
The hydraulic conductivity changes due to the evolution of the fixed bacteria. We use the following relation, given in [21], to describe the evolution of the hydraulic conductivity with respect to this evolution: where K0 is the initial conductivity and nk is a given parameter. At the first time step, we resolve the system (38) with K0, which gives the first velocity u, and then, the system (37) to obtain (b (1) N ). Then, at each subsequent time step t n , for n ≥ 2, the algorithm is as follows: given (b ) to obtain the new velocity, and we resolve the system (37) to obtain (b We call Models 1, 2, and 3 the models defined with the functions (1), (2), and (3), respectively, and we give some comparison between these models. Using growth function (1), which is dependent only on carbon as in [5], can give rise to problems for some situations. For example, in the case where S in N = 0, we come across negative values of the nitrate concentration, which is nonsense. Figures 3 and 4 give the evolution of the concentrations obtained with Models 1 and 2, near Γ 1 and Γ 2 , respectively. In both sides of the reactor, the concentration of nitrate becomes negative for Model 1, while it remains positive for Model 2. Model 3 also gives a non-negative concentration for S in N = 0.  In Figure 5, we plot the concentrations for Models 2 and 3 near Γ 2 for S in N = 0: for both models, the concentration of nitrate remains positive. These foregoing conclusions also hold for small values of S in N , as we can see in Figure 6, which represents the evolution of the concentrations of nitrate and carbon near Γ 2 with S in N = 20, for all three models.  In addition, it was shown in [22] that the ratio between the carbon used and the nitrate degraded remains constant throughout the experiment. In the curve representing S N with respect to S given in Figure 7, this conclusion is satisfied for Model 2, but not for Model 1, at least after about 80 days. We conclude that the domain of validity of Model 1 is more reduced than that of Model 2 or Model 3. We will now consider a larger S in N than in the previous cases. The simulations of the evolution of the concentrations of the four components-the nitrate, the nutrient, the free, and the adherent-at the point (1.5,0.5) with T = 300 days, are represented in Figure 8 for Models 1 and 2. We observe that the results achieved using Model 1 better agree qualitatively with the experimental results given by Chevron [20] (see Figure 25 and page 104) compared to the other two models. Specifically, the percentage of nitrate removal stabilizes by 95 percent in 80 days of reactor operation. Indeed, Figure 8 shows that, in Model 1, ninety percent of nitrates are removed in about 80 days, whereas Model 2 seems to be less consistent with these experimental results and would need additional parameter calibration. The comparison between the concentrations of nitrate and carbon for the three models is given in Figure 9.  We observe that Model 1 and Model 3 concur with the experimental results. From all these observations, we can posit that the growth function (3) (the minimum formula) seems to be more suitable for our problem than (1) and (2). We explain the fact that Model 1 also gives good results by the fact that when S in N is large enough, it prevents the nitrate from becoming negative, giving: while small S in N gives a negative concentration because: In [8], for another problem, the authors gave a comparison between three formulas, including the minimum and multiplicative formulas, and found that the minimum formula concurs better with the experimental results. In the literature, the Formulas (3) and (2) are generally presented as equivalent. If we consider the system (22) with only the reaction terms, i.e., the non-spatialized model, which is a system of Ordinary Differential Equations (ODE), we notice that these formulas do indeed give equivalent results, as can be seen in Figures 10 and 11.  These last figures are obtained using a Runge-Kutta scheme for ODE. The difference between the two cases is found mainly at the start of the calculations (first five days). Furthermore, it can be seen that the nitrate elimination rate stabilizes at 75%, which is a long way from the experimental results. We conclude that: firstly, the introduction into the equation of the diffusion and transport gives a more realistic model, and secondly, Formula (3) is better than (2) (as mentioned above). This can be explained by the fact that (2) can create a numerical (not physical) diffusion, as shown in Figures 8 and 9.

Conclusions
In this work, we present the impact of different forms of growth rate functions (1)-(3) on the numerical results of a spatialized model of the denitrification process in a porous media. In particular, we take into account both free and adherent bacteria. This model is composed of four coupled non-linear partial differential equations (PDEs) related to biological phenomena. We establish the existence and uniqueness of the weak solution of this system with Formulas (1) and (2). We show that the first one, which depends on one substrate only, gives satisfactory results as long as the nitrate concentration remains positive (for large values of S in N ). For the second growth rate expression, which takes both substrate concentrations into account and which is the most commonly used in the literature, the numerical results are valid whatever the values of S in N . However, this second expression seems to be less consistent with the experimental results, even in the case where ODE models using Formulas (2) and (3) give equivalent results. We found that the growth function (3) is the most suitable for our PDE model since numerical tests have shown that the minimum form gives results closer to those of the experiment.