# Impact of Dual Substrate Limitation on Biodenitrification Modeling in Porous Media

Mostafa Abaali, Jérôme Harmand, Zoubida Mghazli

► To cite this version:

Mostafa Abaali, Jérôme Harmand, Zoubida Mghazli. Impact of Dual Substrate Limitation on Biodenitrification Modeling in Porous Media. Processes, 2020, 8 (8), pp.890. 10.3390/pr8080890. hal-02948952

HAL Id: hal-02948952

https://hal.inrae.fr/hal-02948952

Submitted on 25 Sep 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Impact of Dual Substrate Limitation on Biodenitrification Modeling in Porous Media

Mostafa Abaali<fnref n="1" />,* D, Jérôme Harmand<fnref n="2" /> and Zoubida Mghazli<fnref n="1" />

Received: 31 May 2020; Accepted: 21 July 2020; Published: 24 July 2020

<figure><img src="image_2.png" /></figure>

check for
updates

**Abstract:** In this work, we consider a model of the biodenitrification 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 finite element method that satisfies 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 find that the "minimum" model gives numerical results closer to the experimental results.

**Keywords:** modeling; biodenitrification; convection-diffusion-reaction equations; finite element method

## 1. 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₂ and/or NO₃) 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

<footnote n="1">¹ Lirne, Eima, Department of Mathematics, Faculty of Sciences, Ibn Tofail University, 14000 Kenitra, Morocco; zoubida.mghazli@uit.ac.ma</footnote>

<footnote n="2">² LBE, University of Montpellier, INRAE, F-11100 Narbonne, France; jerome.harmand@inrae.fr</footnote>

<footnote>* Correspondence: mostafa.abaali@uit.ac.ma; Tel.: +212-6665-15143</footnote>
interactions into account. We considered that the growth function $\mu(\cdot)$ depends only on the nutrient concentration $S$ following the Monod law:

$$ \mu(S) = \mu_{max}^{S} \frac{S}{S + K_{S}}, \qquad (1) $$

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 $\mu(\cdot)$ 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:

$$ \mu(S_1, S_2) = \left( \mu_{max}^{1} \frac{S_1}{S_1 + K_1} \right) \left( \mu_{max}^{2} \frac{S_2}{S_2 + K_2} \right) \qquad (2) $$

and the second, called the minimum formula, is given by:

$$ \mu(S_1, S_2) = \min \left( \mu_{max}^{1} \frac{S_1}{S_1 + K_1}; \mu_{max}^{2} \frac{S_2}{S_2 + K_2} \right), \qquad (3) $$

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 $\mu(\cdot)$ 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).

## 2. Mathematical Model

Let $\Omega$ be an open set of $\mathbb{R}^2$ with a regular enough boundary $\partial\Omega = \Gamma$, which is divided into three parts $\Gamma_1, \Gamma_2$ and $\Gamma_3$:

$$ \Gamma = \Gamma_1 \cup \Gamma_2 \cup \Gamma_3. $$

We suppose that $\Omega$ contains nitrified water, denitrifying bacteria, and a nutrient. The flow of the nitrified water comes from $\Gamma_1$ and goes through $\Gamma_2$; we assume that the flow is permanent with a velocity $u$. The impermeable part of the boundary is $\Gamma_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:

$$ Q_T := \Omega \times ]0, T], \quad \text{with the boundary} \quad \Sigma_T = \Gamma \times ]0, T]. $$
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_\infty$. The proportion of the occupation of the wall is a number between zero and one defined by the ratio:

$$ \bar{b}_f = \frac{b_f}{w_\infty}. $$

The function $\mu_i(\cdot)$ 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:

$$ \frac{\partial b_f}{\partial t} = (\mu_f(S)G(\bar{b}_f) - k_f - \beta)b_f + \alpha(1 - \bar{b}_f)\gamma b_m, \quad (4) $$

where:

* $k_f$ denotes the adherent bacteria mortality rate;

* $\beta$ 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 $\alpha$.

* $G(\bar{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]).

* $\gamma$ 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:

$$ \frac{\partial b_m}{\partial t} = (\mu_m(S) - k_m - \alpha)b_m + \gamma^{-1}b_f(\mu_f(S)(1 - G(\bar{b}_f)) + \beta), \quad (5) $$

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]):

$$ \frac{\partial S}{\partial t} = -\frac{1}{Y_m}b_m\mu_m(S) - \frac{1}{Y_f}\gamma^{-1}b_f\mu_f(S), \quad (6) $$

$$ \frac{\partial S_N}{\partial t} = -\frac{R}{Y_m}b_m\mu_m(S) - \frac{R}{Y_f}\gamma^{-1}b_f\mu_f(S), \quad (7) $$

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 $\mu_i$ for $i = f, m$. According to the form considered by Freter ([9,12]), let:

$$ G(X) = \frac{1 - X}{a + 1 - X}, $$

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 $\mu_i(\cdot)$, 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 $\mu_i(\cdot)$ in (1) and (2) satisfy the following hypothesis:

**Hypothesis 1.** The growth rate of bacteria $\mu_i(S, S_N)$, for $i = f, m$, satisfies:

$$ \mu_i \in C^1, \quad \mu_i(0, S_N) = \mu_i(S, 0) = 0. $$

**Hypothesis 2.** The function G satisfies:

$$ G \in C^1, \quad 0 < G(0) \le 1, \quad G(1) = 0. $$

The rate of change of concentrations due to diffusion is represented by terms of the form $D\nabla 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:

$$ \begin{aligned} \frac{\partial b_f}{\partial t} - \operatorname{div}(D_f \nabla b_f) &= \left( \mu_f(S, S_N)G(\bar{b}_f) - k_f - \beta \right) b_f + \alpha\gamma(1 - \bar{b}_f)b_m, \\ \frac{\partial b_m}{\partial t} - \operatorname{div}(D_m \nabla b_m - u \cdot b_m) &= \left( \mu_m(S, S_N) - k_m - \alpha(1 - \bar{b}_f) \right) b_m \\ &\quad + \left( \mu_f(S, S_N)(1 - G(\bar{b}_f)) + \beta \right) \gamma^{-1} b_f, \\ \frac{\partial S}{\partial t} - \operatorname{div}(D_1 \nabla S - u \cdot S) &= - \left( \frac{\mu_m(S, S_N)}{Y_m} \right) b_m - \left( \frac{\mu_f(S, S_N)}{Y_f} \right) \gamma^{-1} b_f, \\ \frac{\partial S_N}{\partial t} - \operatorname{div}(D_2 \nabla S_N - u \cdot S_N) &= - \left( \frac{R\mu_m(S, S_N)}{Y_m} \right) b_m \\ &\quad - \left( \frac{R\mu_f(S, S_N)}{Y_f} \right) \gamma^{-1} b_f. \end{aligned} $$

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 $\Gamma$. By the definition of fixed bacteria, there is no flux across the entire boundary of the domain:

$$ \frac{\partial b_f}{\partial n} = 0 \quad \text{on } \Gamma \times ]0, T[. \qquad (8) $$

In order to maintain the growth of the bacteria, we continuously inject from $\Gamma_1$ the substrate with the density $S^{in}$. The boundary conditions used on $\Gamma_1$, which models flow continuity, assuming that the concentrations are uniform outside, are Robin's type, also called Danckwerts [14]:

$$ - D_1 \nabla S \cdot n + (u \cdot n)S = (u \cdot n)S^{in} \quad \text{on } (\Gamma_1 \times ]0, T[). \qquad (9) $$

At the output, the flow of S is uniform, and the condition on $\Gamma_2$ is then given by:

$$ - D_1 \nabla S \cdot n + (u \cdot n)S = 0 \quad \text{on } (\Gamma_2 \times ]0, T[). \qquad (10) $$

On the impervious part $\Gamma_3$, it is natural to consider:

$$ \nabla S \cdot n = 0 \quad \text{on } (\Gamma_3 \times ]0, T[). \qquad (11) $$
For $b_m$, similar reasoning leads to the following boundary conditions:

$$ -D_m \nabla b_m \cdot \mathbf{n} + (\mathbf{u} \cdot \mathbf{n})b_m = 0 \quad \text{on} \; ((\Gamma_1 \cup \Gamma_2) \times ]0, T[), \qquad (12) $$

$$ \nabla b_m \cdot \mathbf{n} = 0 \quad \text{on} \; (\Gamma_3 \times ]0, T[). \qquad (13) $$

The nitrified water contained in $\Omega$ comes from $\Gamma_1$ with a velocity flow $\mathbf{u}$ and is withdrawn with the same flow from the part $\Gamma_2$ of the boundary. The concentration of nitrates thus verifies the following boundary conditions:

$$ -D_2 \nabla S_N \cdot \mathbf{n} + (\mathbf{u} \cdot \mathbf{n})S_N = (\mathbf{u} \cdot \mathbf{n})S_N^{in} \quad \text{on} \; (\Gamma_1 \times ]0, T[), \qquad (14) $$

$$ -D_2 \nabla S_N \cdot \mathbf{n} + (\mathbf{u} \cdot \mathbf{n})S_N = 0 \quad \text{on} \; (\Gamma_2 \times ]0, T[), \qquad (15) $$

$$ \nabla S_N \cdot \mathbf{n} = 0 \quad \text{on} \; (\Gamma_3 \times ]0, T[). \qquad (16) $$

To close this system, we define the initial conditions. At $t = 0$, the domain $\Omega$ contains nitrified water with an initial density $S_N^0$, planktonic bacteria with a density $b_m^0$, adherent bacteria with a density $b_f^0$, and a carbon source with an initial density $S^0$:

$$ S_N(x, 0) = S_N^0 \quad \text{in } \Omega, \qquad (17) $$

$$ b_m(x, 0) = b_m^0 \quad \text{in } \Omega \qquad (18) $$

$$ b_f(x, 0) = b_f^0 \quad \text{in } \Omega, \qquad (19) $$

$$ S(x, 0) = S^0 \quad \text{in } \Omega. \qquad (20) $$

We suppose that:

$$ S_N^0 \geq 0, \quad b_m^0 \geq 0, \quad b_f^0 \geq 0 \quad \text{and} \quad S^0 \geq 0. $$

### 3. Analysis and Approximation

The usual space $L^p(\Omega)$, for $1 \le p < \infty$, is defined by:

$$ L^p(\Omega) := \left\{ f : \|f\|_{p,\Omega} := \left(\int_{\Omega} |f|^p \right)^{1/p} < \infty \right\}. $$

The norm in $L^2(\Omega)$ is denoted by $\| f \|_{0,\Omega} = (\int_{\Omega} |f|^2)^{1/2}$. For $p = \infty$, the space $L^\infty(\Omega)$ is defined by $L^\infty(\Omega) = \{f : \sup_{\Omega} \text{ess} |f| < \infty\}$, equipped with the norm $\| f \|_\infty = \sup_{\Omega} \text{ess}|f|$. If $X$ is a Banach space, $L^p(0,T;X)$ is the set of the measurable functions in $X$ such that $\|f\|_{L^p(0,T;X)} := (\int_0^T \|f\|_X^p)^{1/p} < \infty$. For $p = \infty$ $\| f \|_{L^\infty(0,T;X)} := \sup_{X \in \Omega} \| f \|_X$. Let $H^1(\Omega)$ be the usual Sobolev space of the first order defined by:

$$ H^1(\Omega) := \{v \in L^2(\Omega) / \nabla v \in (L^2(\Omega))^n\}, $$

equipped with the standard norm: $\|v\|_{1,\Omega} := (\|v\|_{0,\Omega}^2 + \|\nabla v\|_{0,\Omega}^2)^{1/2}$, and $H^{-1}(\Omega)$ be its topological dual. From now on, we adopt the following notations:

$$ \mathbf{C} = (c_1, c_2, c_3, c_4) := (b_f, b_m, S, S_N) \quad \text{and} \quad \mathbf{C}^0 = (b_f^0, b_m^0, S^0, S_N^0). \qquad (21) $$

$$ \| \mathbf{C} \|_{0,\Omega} = \left( \sum_{i=1}^{4} \| c_i \|_{0,\Omega} \right) \quad \text{and} \quad \| \mathbf{C} \|_{1,\Omega} = \left( \sum_{i=1}^{4} \| c_i \|_{1,\Omega} \right). $$
We put:

$$F_1(\mathbf{C}) := \left( \mu_f(S, S_N)G(\bar{b}_f) - k_f - \beta \right) b_f + \alpha\gamma(1 - \bar{b}_f)b_m,$$

$$F_2(\mathbf{C}) := \left( \mu_m(S, S_N) - k_m - \alpha(1 - \bar{b}_f) \right) b_m + \left( \mu_f(S, S_N)(1 - G(\bar{b}_f)) + \beta \right) \gamma^{-1} b_f,$$

$$F_3(\mathbf{C}) := - \left( \frac{\mu_m(S, S_N)}{Y_m} \right) b_m - \left( \frac{\mu_f(S, S_N)}{Y_f} \gamma^{-1} \right) b_f,$$

$$F_4(\mathbf{C}) := - \left( \frac{R\mu_m(S, S_N)}{Y_m} \right) b_m - \left( \frac{R\mu_f(S, S_N)}{Y_f} \gamma^{-1} \right) b_f,$$

and:

$$\mathbf{F} := (F_1, F_2, F_3, F_4).$$

The global fluxes are defined by:

$$J_1(c_1) := D_f \nabla b_f,$$

$$J_2(c_2) := D_m \nabla b_m - \mathbf{u}b_m,$$

$$J_3(c_3) := D_1 \nabla S - \mathbf{u}S,$$

$$J_3(c_4) := D_2 \nabla S_N - \mathbf{u}S_N.$$

The boundary operator will be denoted by $\mathbf{B} := (B_1, B_2, B_3, B_4)$ with:

$$B_i(c_i) = \begin{cases} J_i(c_i) \cdot \mathbf{n} & \text{on } \Gamma_1 \cup \Gamma_2 \text{ for } i = 2, 3, 4 \\ \nabla c_i \cdot \mathbf{n} & \text{on } \Gamma_3 \text{ for } i = 2, 3, 4 \text{ and on } \Gamma \text{ for } i = 1. \end{cases}$$

Let $\mathbf{g} := (0, 0, g_3, g_4)$ with:

$$g_3 = \begin{cases} \mathbf{u} \cdot \mathbf{n} S^{in} & \text{on } \Gamma_1 \\ 0 & \text{on } \Gamma_2 \cup \Gamma_3 \end{cases} \quad \text{and} \quad g_4 = \begin{cases} \mathbf{u} \cdot \mathbf{n} S_N^{in} & \text{on } \Gamma_1 \\ 0 & \text{on } \Gamma_2 \cup \Gamma_3. \end{cases}$$

With these notations, the quasilinear diffusion-convection system (8)–(20) becomes:

$$\left\{ \begin{array}{l} \displaystyle \frac{\partial c_i}{\partial t} - \text{div}(J_i(c_i)) = F_i(\mathbf{C}) \quad \text{in } Q_T, \text{ for } i = 1, 2, 3, 4 \\ \mathbf{B}(\mathbf{C}) = \mathbf{g} \quad \text{on } \Sigma_T \\ \mathbf{C}(0, \cdot) = \mathbf{C}^0 \quad \text{in } \Omega. \end{array} \right. \qquad (22)$$

### 3.1. Existence Theorem

Since $\mathbf{C}$ is in $(L^\infty(\Omega))^4$ and $F_i$ is in $C^1(\mathbb{R}^4)$ (in the case when we use the Formulas (1) or (2)) for $1 \le i \le 4$, the classical local existence result is well known (see [15]). There exists $T > 0$ and a unique classical solution of (22) in $(C([0, T]; L^1(\Omega)))^4 \cap (L^\infty([0, T - \tau] \times \Omega))^4$, for all $\tau \in (0, T)$. To show the global existence and uniqueness of the weak solution of (22), the functions $F_i$ for $i = 1, \dots, 4$ have to satisfy some hypothesis. The first one is the quasi-positivity of $F_i$.

**Lemma 1.** For $1 \le i \le 4$, the functions $F_i$ are in $C^1(\mathbb{R}^4)$ and satisfy the quasi-positivity property:

$$\text{For } c_j \ge 0, \quad j = 1, \dots, 4, \quad F_i(\mathbf{C}) \ge 0 \quad \text{whenever } c_i = 0.$$

**Proof of Lemma 1.** Since $G$, $\mu_f$, and $\mu_m$ are in $C^1$, so is $F_i$, for $1 \le i \le 4$. The quasi-positivity is immediately obtained for non-negative densities $b_f$, $b_m$, $S$, and $S_N$. $\square$
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) \in Q_T$.

**Lemma 2.** There exists a lower triangular invertible matrix $\mathbb{M} = (m_{ij})_{1\le i,j\le 4}$ with non-negative diagonal entries and a vector $\mathbf{V}$ in $\mathbb{R}_+^4$ such that:

$$ \forall \mathbf{C} \in ([0, +\infty))^4, \quad \mathbb{M}\mathbf{F}(\mathbf{C}) \le \left(1 + \sum_{i=1}^{4} c_i\right) \mathbf{V} \qquad (23) $$

(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 \in \Omega$:

$$ \begin{aligned} F_1(\mathbf{C}) &\le \max\left((\mu_f G(\bar{b}_f) - k_f - \beta); (\alpha\gamma(1 - \bar{b}_f))(b_m + b_f)\right) \\ F_2(\mathbf{C}) &\le \max\left((\mu_m - k_m - \alpha(1 - \bar{b}_f); (\mu_f(1 - G(\bar{b}_f)) + \beta)\gamma^{-1}\right)(b_m + b_f). \end{aligned} $$

Let:

$$ \begin{aligned} v_1 &:= \max\left((\mu_f G(\bar{b}_f) - k_f - \beta); (\alpha\gamma(1 - \bar{b}_f))\right), \\ v_2 &:= \max\left((\mu_m - k_m - \alpha(1 - \bar{b}_f)); (\mu_f(1 - G(\bar{b}_f)) + \beta)\gamma^{-1}\right), \\ \mathbf{V} &:= (v_1, v_2, 1, 1). \end{aligned} \qquad (24) $$

Since $S$ and $S_N$ are positive and $F_3$ and $F_4$ are negative, we obtain for $i = 1, 2, 3, 4$:

$$ F_i(\mathbf{C}) \le v_i(b_f + b_m + S + S_N) \le v_i(1 + b_f + b_m + S + S_N). \qquad (25) $$

We take $\mathbb{M}$ as the identity matrix to obtain (23). $\square$

The third essential property considered for the reaction term F is the “polynomial growth” structure of F defined in the following lemma.

**Lemma 3.**

$$ \forall T > 0, \exists L > 0, p > 0 : \forall i, \forall \mathbf{y} \in ([0, +\infty]^4, |F_i(\mathbf{y})| \le L(1 + |\mathbf{y}|^p), \qquad (26) $$

where $|\cdot|$ is the euclidean norm in $\mathbb{R}^4$.

**Proof of Lemma 3.** This is a consequence of (25) by taking $p = 1$ and $M = \max_{1\le i\le 4} v_i$. $\square$

Lemmas 1–3 allow us to have the result of the global existence and uniqueness of the weak solution, in the following theorem (see Theorem 1 in [16]).

**Theorem 1.** The problem (22) has a unique global non-negative weak solution in the following: sense

$$ \left\{ \begin{array}{l} \forall T > 0, \mathbf{C} \in (C([0, T]; L^2(\Omega)) \cap L^\infty(Q_T) \cap L^2(0, T; H^1(\Omega)))^4, \\ \forall \Psi = (\psi_1, \psi_2, \psi_3, \psi_4) \in (C^\infty(\overline{Q_T}))^4 \text{ such that } \Psi(T) = 0, \text{ for } i = 1, 2, 3, 4 \\ - \int_{Q_T} c_i \frac{\partial \psi_i}{\partial t} + \int_{Q_T} J_i(c_i) \cdot \nabla \psi_i = \int_{Q_T} F_i(\mathbf{C}) \psi_i + \int_{\Omega} c_i^0 \psi_i(0) + \int_{\Sigma_T} g_i \psi_i. \end{array} \right. \qquad (27) $$
Moreover, for any $T > 0$, there exists: $M > 0$ depends on $T$, $c^0$, $\mathbf{u}$, and $\mathbf{V}$ defined by (24) such that:

$$ \|\mathbf{C}\|_{(L^\infty(Q_T))}^4 + \|\mathbf{C}\|_{(L^2(0,T;H^1(\Omega)))}^4 + \left\|\frac{\partial \mathbf{C}}{\partial t}\right\|_{(L^2(0,T;H^{-1}(\Omega)))}^4 \leq M \quad (28) $$

## 3.2. Approximation

By taking the test function $\Psi$ in $(D(Q_T))^4$, then in $(D(\overline{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:

$$ \left\{ \begin{array}{l} \text{Find } \mathbf{C} \in (L^2(0, T; L^2(\Omega)) \cap L^2(0, T; H^1(\Omega)))^4 \text{ such that for a.e. } t \in [0, T[ \\ \frac{d}{dt} \int_{\Omega} c_i(t) \psi_i + \int_{\Omega} J_i(c_i(t)) \cdot \nabla \psi_i = \int_{\Omega} F_i(\mathbf{C}(t)) \psi_i + \int_{\Gamma} g_i \psi_i, \quad \forall \Psi \in (H^1(\Omega))^4 \\ \mathbf{C}(0) = \mathbf{C}_0. \end{array} \right. \quad (29) $$

In order to have a $(H^1(\Omega))^4$-elliptic bilinear form, we transform the problem by using the following augmented bilinear form:

$$ \tilde{a}_i(u, v) = \int_{\Omega} J_i(u) \cdot \nabla v + \gamma \int_{\Omega} uv \quad (30) $$

where $\gamma = \inf_{1 \leq i \leq 4} D_i$, which is $H^1$-elliptic if $\|\mathbf{u}\|_{\infty} \leq 2\gamma$. Indeed, we have:

$$ \begin{array}{ll} \text{for } i = 1, & \tilde{a}_1(v, v) = D_1 \|\nabla v\|_{0,\Omega}^2 + \gamma \|v\|_{0,\Omega}^2 \geq \gamma \|v\|_{1,\Omega}^2, \\ \text{For } 2 \leq i \leq 4, & \tilde{a}_i(v, v) = D_i \|\nabla v\|_{0,\Omega}^2 - \int_{\Omega} v \mathbf{u} \cdot \nabla v + \gamma \|v\|_{0,\Omega}^2 \\ & \geq \left(D_i - \frac{\|\mathbf{u}\|_{\infty}}{2}\right) \|\nabla v\|_{0,\Omega}^2 + \left(\gamma - \frac{\|\mathbf{u}\|_{\infty}}{2}\right) \|v\|_{0,\Omega}^2 \\ & \geq \left(\gamma - \frac{\|\mathbf{u}\|_{\infty}}{2}\right) \|v\|_{1,\Omega}^2. \end{array} $$

The problem is then transformed as follows. Let $\Phi = (\phi_1, \phi_2, \phi_3, \phi_4)$, with $\phi_i = \exp(-\gamma t)c_i$, and let $\tilde{g}_i = \exp(-\gamma t)g_i$ and $G_i(\Phi) = \exp(-\gamma t)F_i(\exp(\gamma t)\Phi)$. The problem (29) is equivalent to the system: for a.e. $t \in [0, T]$,

$$ \left\{ \begin{array}{l} \text{find } \Phi(t) \in (H^1(\Omega))^4 \text{ such that for } i = 1, \dots, 4 \\ \frac{d}{dt} \int_{\Omega} \phi_i(t) \psi_i + \tilde{a}_i(\phi_i(t), \psi_i) = \int_{\Omega} G_i(\Phi) \psi_i + \int_{\Gamma} \tilde{g}_i \psi_i, \quad \forall \psi_i \in H^1(\Omega), \\ \Phi(0) = \mathbf{C}_0. \end{array} \right. \quad (31) $$

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.

### 3.2.1. Semi-Discretization

Let $\mathcal{T}_h = \bigcup_{\mathbf{T} \in \mathcal{T}_h} \mathbf{T}$ be a family of regular triangulations of $\Omega$, where $\mathbf{T}$ is a triangle and $h_{\mathbf{T}}$ its diameter. We denote $h = \max_{\mathbf{T} \in \mathcal{T}_h} \text{diam}(\mathbf{T})$. Let $\mathcal{P}_1(\mathcal{T}_h)$ and $\mathcal{V}_h$ be defined respectively by:

$$ \mathcal{P}_1(\mathcal{T}_h) = \{v \in C(\bar{\Omega}) \mid v|_{\mathbf{T} \in \mathcal{P}_1(\mathcal{T}_h)}, \forall \mathbf{T} \in \mathcal{T}_h\} \quad \text{and} \quad \mathcal{V}_h = (\mathcal{P}_1(\mathcal{T}_h))^4, $$
where, for each $T \in \mathcal{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 \in ]0, T]$, by:

$$ egin{cases} 
\text{Find } \Phi_h(t) \in \mathcal{V}_h \text{ such that for } i = 1, \dots, 4 \\ 
\frac{d}{dt} \int_{\Omega} \phi_{ih}(t) \psi_{ih} + \tilde{a}(\phi_{ih}, \psi_{ih}) = \int_{\Omega} G_i(\Phi_h(t)) \psi_{ih} + \int_{\Gamma} \tilde{\sigma}_i \psi_{ih}, \quad \forall \psi_{ih} \in P_1(\mathcal{T}_h), \\ 
\Phi_h(0) = \Pi_h \Phi_0, 
\end{cases} 
\quad (32) $$

where $\Pi_h$ is a projection operator on $\mathcal{V}_h$. Let $\bar{t} \in ]0, T]$. To give an a priori error estimate for a local solution of the semi-discretized problem (32), we consider a finite time interval $\bar{J} = (0, \bar{t}]$.

The error estimates in the $L^2$-norm are based on the elliptic projection: for all $t \in [0, \bar{t}]$ and for all $0 \le i \le 4$, let $\omega_{ih}(t)$ be the elliptic projection of the exact solution $\phi_i$ on $P_1(\mathcal{T}_h)$ defined by:

$$ \tilde{a}_i(\omega_{ih}(t) - \phi_i(t), \chi) = 0, \quad \forall \chi \in P_1(\mathcal{T}_h), \quad (33) $$

where $\tilde{a}_i(\cdot, \cdot)$ is defined in (30). Following the technique of Thomée [17], we can prove the following theorem.

**Theorem 2.** Let $\Phi$ and $\Phi_h$ be the solutions of (31) and (32), respectively, and let $\Pi_h \Phi_0$ be the elliptic operator defined by (33) for $t = 0$. We assume that $\|\mathbf{u}\|_{\infty} < 2\gamma$. Then, we have:

$$ \| \Phi(t) - \Phi_h(t) \|_{0,\Omega} \le L(\Phi)h^2, \quad \text{for } t \in \bar{J}, $$

where $L(\Phi)$ is a non-negative constant and depends on $t$ and the solution $\Phi$, but independent of $h$.

### 3.2.2. Full Discretization

We consider the discretization of $[0, T]$ given by $0 = t_0 < t_1 < \dots < t_N = T$ and put $\tau_n = t_n - t_{n-1}$ and $\tau = \max_{1 \le n \le N} \tau_n$. For all $k, 0 \le k \le N$, and all $i, 1 \le i \le 4$, we use the notation:

$$ egin{cases} 
\phi_i^k := \phi_i(t_k) \\ 
\Phi_h^k := \Phi_h(t_k) = (\phi_{h1}^k, \phi_{h2}^k, \phi_{h3}^k, \phi_{h4}^k). 
\end{cases} 
\quad (34) $$

The derivative with respect to time is approximated by the backward Euler scheme given by the following difference quotient:

$$ \frac{\partial \phi_i}{\partial t}(t_n) \approx \frac{\phi_i^n - \phi_i^{n-1}}{\tau_n}. $$

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:

$$ egin{cases} 
\text{for } 1 \le n \le N, \text{ find } \Phi_h^n \in \mathcal{V}_h \text{ such that for } i = 1, \dots, 4 \\ 
\int_{\Omega} \frac{\phi_{ih}^n - \phi_{ih}^{n-1}}{\tau_n} \psi_{ih} + \tilde{a}(\phi_{ih}^n, \psi_{ih}) = \int_{\Omega} G_i(\Phi_h^{n-1}) \psi_{ih} + \int_{\Gamma} \tilde{\sigma}_i \psi_{ih}, \quad \forall \Psi \in \mathcal{V}_h \\ 
\Phi_h^0 = \Pi_h(\Phi_0). 
\end{cases} 
\quad (35) $$

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 $\Phi_h^n$ and $\Phi$ be the solutions of (35) and (31), respectively. Then, there exists a constant $\bar{M}$ depending on $\Phi$, but independent of $h$ and $\tau$ such that:

$$ \|\Phi_h^n - \Phi(t_n)\|_{0,\Omega} \le M(\Phi)(h^2 + \tau). \quad (36) $$

**Remark 1.** The estimates of Theorem 3 remain valid for $C_h$ since we have $\Phi_h = \exp(-\gamma t)C_h$. Therefore, in the numerical tests, we use the following approximated problem:

$$ \left\{ \begin{array}{ll} \text{for } 1 \le n \le N, \text{ find } C_h^n \in \mathbb{V}_h \text{ such that for } i = 1, \dots, 4 \\ \displaystyle \int_{\Omega} \frac{c_{ih}^n - c_{ih}^{n-1}}{\tau_n} \psi_{ih} + \int_{\Omega} J_i(c_{ih}^n) \cdot \nabla \psi_{ih} = \int_{\Omega} F_i(C_h^{n-1}) \psi_{ih} \\ \displaystyle + \int_{\Gamma} g_i \psi_{ih}, \quad \forall \psi_{ih} \in \mathcal{P}_1(\mathcal{T}_h), \\ C_h^0 = \Pi_h(C_0), \end{array} \right. \quad (37) $$

## 4. Numerical Tests

In this section, we consider a domain $\Omega \subset \mathbb{R}^2$ defined by $\Omega = ]0, 3[ \times ]0, 1[, \Gamma_1 = \{0\} \times [0, 0.3]$ and $\Gamma_2 = \{3\} \times [0.7, 1]$ (see Figure 1).

<figure><img src="image_3.png" /><figcaption>Figure 1. Domain.</figcaption></figure>

We assume that $\Omega$ is a porous medium, hereafter denoted the “reactor” $\Omega$, into which we inject some water containing nutrients, sources of both nitrate and a carbon with concentrations $S_N^{in}$ and $S^{in}$, respectively. In this case, the flow in $\Omega$ is governed by Darcy’s equation, given in its mixed form by the system: find $(\mathbf{u}, p) \in H(\text{div}; \Omega) \times L^2(\Omega)$, with $\mathbf{u} \cdot \mathbf{n} = u_0$ on $\Gamma_1 \cup \Gamma_2$ such that:

$$ \left\{ \begin{array}{ll} \displaystyle \int_{\Omega} K^{-1} \mathbf{u} \cdot \mathbf{v} + \int_{\Omega} \text{div} \, \mathbf{v} \, p = \int_{\Gamma_2} \mathbf{v} \cdot \mathbf{n} p_D, & \forall \mathbf{v} \in \mathbf{V} \\ \\ \displaystyle \int_{\Omega} \text{div} \, \mathbf{u} \, q + \int_{\Omega} f q = 0, & \forall q \in L^2(\Omega), \end{array} \right. \quad (38) $$

where $H(\text{div}; \Omega) := \{\mathbf{v} \in (L^2(\Omega))^2 / \text{div} \, \mathbf{v} \in L^2(\Omega)\}$ and $\mathbf{V} := \{\mathbf{v} \in H(\text{div}; \Omega) / \mathbf{v} \cdot \mathbf{n} = 0 \text{ on } \Gamma_1 \cup \Gamma_2\}$, $K$ is the hydraulic conductivity of the medium, $u_0$ is the given flux on $\Gamma_1$, and $p_D$ is the given pressure on $\Gamma_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 $\mathbf{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.
<figure><img src="image_4.png" /><figcaption>Figure 2. Mesh.</figcaption></figure>

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.

<table><caption>Table 1. Model parameter values used for the simulations.</caption><thead><tr><th>Parameters</th><th>Values</th><th>Parameters</th><th>Values</th><th>Initial Conditions</th><th>Values</th></tr></thead><tbody><tr><td>D<sub>1</sub></td><td>0.3 cm<sup>2</sup>/h</td><td>μ<sup>f</sup><sub>max</sub></td><td>0.7 h<sup>-1</sup></td><td>b<sup>0</sup><sub>f</sub></td><td>10 g/cm<sup>2</sup></td></tr><tr><td>D<sub>2</sub></td><td>0.4 cm<sup>2</sup>/h</td><td>Y<sub>m</sub></td><td>0.5</td><td>b<sup>0</sup><sub>m</sub></td><td>10 g/mL</td></tr><tr><td>D<sub>b</sub></td><td>0.4 cm<sup>2</sup>/h</td><td>Y<sub>f</sub></td><td>0.5</td><td>S<sup>in</sup></td><td>104 mg/L</td></tr><tr><td>D<sub>f</sub></td><td>0.2 cm<sup>2</sup>/h</td><td>β</td><td>0.2 h<sup>-1</sup></td><td>S<sup>in</sup><sub>N</sub></td><td>100 mg/L</td></tr><tr><td>α</td><td>0.02 h<sup>-1</sup></td><td>T</td><td>100 days</td><td>K<sub>0</sub></td><td>0.865 cm/h</td></tr><tr><td>w<sub>∞</sub></td><td>50 g/cm<sup>2</sup></td><td>R</td><td>1.2</td><td>S<sup>0</sup></td><td>104 mg/L</td></tr><tr><td>—</td><td>—</td><td>K<sup>f</sup><sub>S</sub></td><td>54 mg/L</td><td>S<sup>0</sup><sub>N</sub></td><td>100 mg/L</td></tr><tr><td>k<sub>m</sub></td><td>0.005 h<sup>-1</sup></td><td>K<sup>m</sup><sub>S</sub></td><td>54 mg carbon/L</td><td>—</td><td>—</td></tr><tr><td>k<sub>b</sub></td><td>0.005 h<sup>-1</sup></td><td>K<sup>m</sup><sub>S<sub>N</sub></sub></td><td>50 mg NO<sub>3</sub><sup>-</sup>/L</td><td>—</td><td>—</td></tr><tr><td>μ<sup>m</sup><sub>max</sub></td><td>0.7 h<sup>-1</sup></td><td>K<sup>f</sup><sub>S<sub>N</sub></sub></td><td>50 mg carbon/L</td><td>—</td><td>—</td></tr></tbody></table>

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:

$$K(b_f) = K_0 \left( 1 - \frac{bf}{w_\infty} \right)^{nk},$$

where $K_0$ is the initial conductivity and $nk$ is a given parameter. At the first time step, we resolve the system (38) with $K_0$, which gives the first velocity $u$, and then, the system (37) to obtain $(b_f^{(1)}, b_m^{(1)}, S^{(1)}, S_N^{(1)})$. Then, at each subsequent time step $t_n$, for $n \ge 2$, the algorithm is as follows: given $(b_f^{(n-1)}, b_m^{(n-1)}, S^{(n-1)}, S_N^{(n-1)})$, we resolve the system (38) with $K(b_f^{(n-1)})$ to obtain the new velocity, and we resolve the system (37) to obtain $(b_f^{(n)}, b_m^{(n)}, S^{(n)}, S_N^{(n)})$, until $t_n = T$.

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_N^{in} = 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 $\Gamma_1$ and $\Gamma_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_N^{in} = 0$.
<figure><img src="image_5.png" /><figcaption>Figure 3. Concentrations for Models 1 and 2 at point (0.1, 0.1) for $S_N^{in} = 0$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 4. Concentrations for Models 1 and 2 at point (2.9, 0.8) for $S_N^{in} = 0$.</figcaption></figure>

In Figure 5, we plot the concentrations for Models 2 and 3 near $\Gamma_2$ for $S_N^{in} = 0$: for both models, the concentration of nitrate remains positive. These foregoing conclusions also hold for small values of $S_N^{in}$, as we can see in Figure 6, which represents the evolution of the concentrations of nitrate and carbon near $\Gamma_2$ with $S_N^{in} = 20$, for all three models.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Concentrations for Models 2 and 3 at point (2.9, 0.8) for $S_N^{in} = 0$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Concentrations of carbon and nitrates for the three models at point (2.9, 0.8) for $S_N^{in} = 20$.</figcaption></figure>

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.
<figure><img src="image_9.png" /><figcaption>Figure 7. $S_N$ with respect to S.</figcaption></figure>

We will now consider a larger $S_N^{in}$ 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.

<figure><img src="image_10.png" /><figcaption>Figure 8. Concentration for Models 1 and 2 at point (1.5, 0.5) for $S_N^{in} = 100$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 9. Concentrations of carbon and nitrates for the three models at point (2.9, 0.8) for $S_N^{in} = 100$.</figcaption></figure>

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_N^{in}$ is large enough, it prevents the nitrate from becoming negative, giving:

$$ \mu_{max}^{S} \frac{S}{S + K_S} = \min(\mu_{max}^{S} \frac{S}{S + K_S}; \mu_{max}^{S_N} \frac{S_N}{S_N + K_{S_N}}), $$

while small $S_N^{in}$ gives a negative concentration because:

$$ \min(\mu_{max}^{S} \frac{S}{S + K_S}; \mu_{max}^{S_N} \frac{S_N}{S_N + K_{S_N}}) = \mu_{max}^{S_N} \frac{S_N}{S_N + K_{S_N}}. $$

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.

<figure><img src="image_12.png" /><figcaption>Figure 10. Concentrations of the ODE model using the minimum and the multiplicative formulas (T = 50 day).</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Concentrations of the ODE model using the minimum and the multiplicative formulas (T = 300 day).</figcaption></figure>

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.

**5. 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_N^{in}$). 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_N^{in}$. 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.

**Author Contributions:** All authors were involved and fully aware of the work. Methodology, M.A., J.H., Z.M.; validation, M.A., J.H., Z.M.; formal analysis, M.A., J.H., Z.M.; writing, original draft preparation, M.A., J.H., Z.M.; writing, review and editing, M.A., J.H., Z.M. All authors read and agreed to the published version of the manuscript

**Funding:** This research received no external funding.

**Acknowledgments:** The authors thank PHC TOUBKAL TBK 17/47-36773WE project, the Euro-mediterranean TREASURE research network and the JPI project Control4reuse (cf. http://control4reuse.net) financed by the French Research National Agency under the contract ANR-18-IC4W-0002, for their financial support.

**Conflicts of Interest:** The authors declare no conflict of interest.
References

1. McCarthy, P.L. *Energetics of Organic Matter Degradation*; Mitchell, R., Ed.; Wiley Interscience: New York, NY, USA, 1972; pp. 91–118.

2. Payne, W.J. Reduction of nitrogenous oxides by microorganisms. *Bactiol. Rev.* **1973**, *37*, 409–452. [CrossRef]

3. Smith, H.L.; Waltman, P. *The Theory of the Chemostat: Dynamics of Microbial Competition*; Cambridge University Press: Cambridge, UK, 1995.

4. Bailey, J.E.; Ollis, D. Biochemical engineering fundamentals. *Chem. Eng. Educ.* **1976**, *10*, 162–165, 176.

5. Abaali, M.; Mghazli, Z. Mathematical modeling of biodegradation in situ application to biodenitrification. *Comput. Math. Appl.* **2020**, *79*, 1833–1844. [CrossRef]

6. Grover, J.P. *Resource Competition*; Springer Science & Business Media: Boston, NY, USA, 1997; Volume 18.

7. Chen-Charpentier, B.M.; Kojouharov, H.V. Mathematical modeling of bioremediation of trichloroethylene in aquifers. *Comput. Math. Appl.* **2008**, *56*, 645–656. [CrossRef]

8. Stewart, H.A.; Al-Omari, A.; Bott, C.; De Clippeleir, H.; Su, C.; Takacs, I.; Wett, B.; Massoudieh, A.; Murthy, S. Dual substrate limitation modeling and implications for mainstream deammonification. *Water Res.* **2017**, *116*, 95–105. [CrossRef] [PubMed]

9. Freter, R.; Brickner, H.; Fekete, J.; Vickerman, M.M.; Carey, K.E. Survival and implantation of Escherichia coli in the intestinal tract. *Infect. Immun.* **1983**, *39*, 686–703. [CrossRef] [PubMed]

10. Harmand, J.; Lobry, C.; Rapaport, A.; Sari, T. *Le chémostat: Théorie Mathématique de la Culture Continue de Micro-organismes*; ISTE Editions: Londres, UK, 2017.

11. Vanrolleghem, P.A.; Dochain, D. *Dynamical Modeling and Estimation in Wastewater Treatment Processes*; IWA Publishing: London, UK, 2001.

12. Freter, R.; Brickner, H.; Temme, S. An understanding of colonization resistance of the mammalian large intestine requires mathematical analysis. *Microecol. Ther.* **1986**, *16*, 147–155.

13. Michaelis, L.; Menten, M.L. Die Kinetic der Invertinwirkung. *Biochem. Z* **1913**, *49*, 333–369.

14. Aris, R. *Mathematical Modeling: A Chemical Engineer's Perspective*, 1st ed.; Elsevier: New York, NY, USA, 1999.

15. Pierre, M. Global Existence in Reaction-Diffusion Systems with Control of Mass: A Survey. *Milan J. Math.* **2010**, *78*, 417–455. [CrossRef]

16. Bothe, D.; Fischer, A.; Pierre, M.; Rolland, G. Global wellposedness for a class of reaction-advection-anisotropic-diffusion systems. *Evol. Equ.* **2017**, *17*, 101–130. [CrossRef]

17. Thomée, V. *Galerkin Finite Element Methods for Parabolic Problems*; Springer: Berlin, Germany, 1984.

18. Brezzi, F.; Fortin, M. *Mixed and Hybrid Finite Element Methods*; Springer Science & Business Media: Boston, NY, USA, 2012; Volume 15.

19. Hecht, F.; Le Hyaric, A.; Morice, J.; Ohtsuka, K.; Pironneau, O. FreeFEM++. Available online: http://www.freefem.org/ff++ (accessed on 14 July 2020).

20. Chevron, F. Dénitrification Biologique D'une Nappe Phréatique Polluée par des Composés Azotés D'origine Industrielle: Expérimentations en Laboratoire sur les Cinétiques, le Métabolisme et les Apports de Nutriments. Ph.D. Thesis, Université de Lille 1, Lille, France, 1996.

21. Clement, T.P.; Hooker, B.S.; Skeen, R.S. Microscopic models for the predicting changes in the saturated porous media properties caused by microbial growth. *Ground Water* **1996**, *34*, 934–942. [CrossRef]

22. Hambsch, B.; Werner, P. Die messung der wachstumsrate von bakterien zur optimierung, kontrolle und überwachung von biologischen denitrifikationsanlagen. *Vom Wasser* **1989**, *72*, 235–247.

<figure><img src="image_14.png" /></figure>

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).