# Comparison of two iterative schemes to solve the Chemotaxis - Biodegradation System

Mostafa Abaali, Salih Ouchtout

► To cite this version:

Mostafa Abaali, Salih Ouchtout. Comparison of two iterative schemes to solve the Chemotaxis - Biodegradation System. Boletim da Sociedade Paranaense de Matematica, 2025, 43, pp.77995. 10.5269/bspm.77995. hal-05472046

HAL Id: hal-05472046

https://ifp.hal.science/hal-05472046v1

Submitted on 22 Jan 2026

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 CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

# Comparison of two iterative schemes to solve the Chemotaxis - Biodegradation System

Mostafa Abaali* and Salih Ouchtout

**ABSTRACT:** The mathematical model describing the biodegradation process in porous media by bacteria—divided into planktonic and adherent types—incorporates the chemotaxis effect, wherein the flow velocity depends on the concentration of adherent bacteria. This model consists of a system of five strongly coupled nonlinear parabolic equations, including a nonlinear advection term. Using the fixed-point theorem, we establish the existence, uniqueness, and non-negativity of the solution.

To approximate the solution, we employ the finite element (FE) method. To linearize the system at each time step, we compare two iterative schemes. The first, called the Coupled Prediction Scheme (CPS), is shown to converge, while the second is the conventional Fixed-Point Method (FPM). Theoretical and numerical results show that CPS offers a faster and more efficient alternative to the conventional FPM.

**Key Words:** Chemotaxis, finite element, nonlinear advection term, fixed point methods.

## Contents

<table><tr><td>1</td><td>Introduction</td><td>1</td></tr><tr><td>2</td><td>Problem setting</td><td>2</td></tr><tr><td>3</td><td>Existence and regularity results</td><td>5</td></tr><tr><td>4</td><td>Approximation</td><td>6</td></tr><tr><td>5</td><td>The iterative scheme analysis</td><td>7</td></tr><tr><td>5.1</td><td>Coupled Prediction Scheme.</td><td>9</td></tr><tr><td>5.2</td><td>Fixed-point method.</td><td>15</td></tr><tr><td>6</td><td>Numerical algorithm for solving the coupled problem</td><td>17</td></tr><tr><td>7</td><td>Numerical results</td><td>19</td></tr><tr><td>8</td><td>Conclusion and perspectives</td><td>22</td></tr></table>

## 1. Introduction

Biological denitrification process, orchestrated by a diversity of anaerobic microorganisms, catalyzes the conversion of nitrates into harmless gaseous nitrogen, thereby helping to reduce nitrate concentrations in groundwater. However, the effective implementation of biological denitrification in aquifer systems requires a thorough understanding of its mechanisms. The purpose of this paper is to provide the numerical analysis of the mathematical model developed in [4], formed by partial differential equations (PDEs), modeling the denitrification phenomenon.

The mathematical framework we have constructed comprises a system of five coupled parabolic equations. These equations incorporate a non-linear advection term. To establish the existence, uniqueness, and non-negativity of the solution, we have employed the fixed point theorem, in the same way as in [4], which serves as a foundational aspect of our analysis.

In order to approximate the solution, we have opted for the finite elements (FE) method, a widely utilized numerical technique. Furthermore, to enhance computational efficiency, we have introduced an iterative scheme at each time step. This scheme, known as the Coupled Prediction Scheme (CPS) (see [10]), leverages a projection method to linearize the system and facilitate convergence. Through rigorous analysis,

<footnote>* Corresponding author.</footnote>

<footnote n="2010">2010 Mathematics Subject Classification: 34K28, 35J66, 35J50, 74H15, 74H20, 74H25.</footnote>

<footnote>Submitted July 22, 2025. Published September 18, 2025</footnote>
we have demonstrated the efficacy and reliability of this iterative approach in solving our mathematical model.

This article examines the application of the (CPS) to enhance the resolution of the system modeling the biodegradation process. The CPS was initially introduced in [10] within the context of the thermally coupled Navier-Stokes equation. In this study, we have adapted this approach to coupled nonlinear convection-diffusion equations. In [10], the authors demonstrated the efficiency and precision of the CPS compared to the basic projection scheme (BPS).

The modeling of natural phenomena most of the time gives strongly coupled systems that pose difficulties during the resolution and often leads to very complex time dependent dynamics requiring efficient solvers. To remedy this problem, several strategies have been developed. Among the strategies developed in this context is the use of explicit terms allowing a “numerical” decoupling of equations (thus eliminating the need of a fixed point) [1,2,14,17] or the explicit treatment of the nonlinear terms yielding a numerically less expansive coupled system. Explicit methods impose restrictions on the time step (or the mesh size) which can be severely detrimental to the overall performance; even if the implicit/semi-implicit approach is suboptimal in certain cases, it is generally more robust and gives good performance in most cases. In [10] authors introduce a new iterative scheme based on a projection method called the coupled prediction scheme(CPS) they also presented the basic projection scheme (BPS). In [15] authors proposed a BPS based on the finite volume method.

## 2. Problem setting

Let $\Omega$ be an open set of $\mathbb{R}^d$, $d = 2$ or 3 with an enough regular boundary $\partial\Omega = \Gamma$, which is divided into three parts $\Gamma_1$, $\Gamma_2$ and $\Gamma_3$ (see Fig.1):

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

<figure><img src="image_3.png" /><figcaption>Figure 1: $\Omega$ and $\Gamma_1$, $\Gamma_2$ et $\Gamma_3$.</figcaption></figure>

The domain $\Omega$ contains polluted water, denitrifying bacteria, and a nutrient. The flow of polluted water enters through $\Gamma_1$ and exits through $\Gamma_2$; it is steady with a constant velocity $\mathbf{u}$. The boundary $\Gamma_3$ represents the impermeable part of the domain. Within the reactor, the bacteria are divided into two categories: adherent bacteria, which attach to the reactor walls and form a monolayer biofilm, and planktonic bacteria, which remain mobile and suspended in the fluid. We denote by $x_2$ the density of planktonic bacteria, and by $x_1$ the surface density of adherent bacteria, with a maximum surface density $\omega_\infty$. The concentration of the limiting substrate is denoted by $S_1$, and the concentration of the contaminant by $S_2$. The biodenitrification process is modeled by the system given in [1,2]. 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]. $$

In [5] D. Bothe et al proved existence of global weak solutions a problem which describes the evolution of concentrations of $p$ charged species subject to Fickian diffusion and chemical reactions in the presence of an electrical field called, the Nernst-Planck-Poisson problem, including in particular the Boltzmann
statistics case. A typical model describing chemotaxis is the Keller-Segel equations derived by Keller and Segel [13] which have been studied extensively. We give the following model proposed by Abaali M. and Ouchtout S. [4] to describe the biodenitrification phenomena:

$$ F: \left\{ \begin{array}{l} \frac{\partial x_1}{\partial t} - \nabla \cdot (D_1 \nabla x_1) = (\mu_a(S_1, S_2)G(\bar{x}_1) - k_a - \beta)x_1 + \alpha\gamma(1 - \bar{x}_1)x_2 \quad \text{in } Q_T, \\ \frac{\partial x_2}{\partial t} - \nabla \cdot (D_2 \nabla x_2) + \mathbf{u} \cdot \nabla x_2 + \nabla \cdot (x_2\chi(S_2)\nabla S_2) = \\ \quad \left(\mu_p(S_1, S_2) - k_p - \alpha(1 - \bar{x}_1)\right)x_2 + \left(\gamma^{-1}\mu_a(S_1, S_2)(1 - G(\bar{x}_1)) + \gamma^{-1}\beta\right)x_1 \quad \text{in } Q_T, \\ \frac{\partial S_1}{\partial t} + \mathbf{u} \cdot \nabla S_1 - \nabla \cdot (D_3 \nabla S_1) = -\left(\frac{\mu_p(S_1, S_2)}{Y_p}\right)x_2 - \left(\frac{\mu_a(S_1, S_2)}{Y_a}\gamma^{-1}\right)x_1 \quad \text{in } Q_T, \\ \frac{\partial S_2}{\partial t} + \mathbf{u} \cdot \nabla S_2 - \nabla \cdot (D_4 \nabla S_2) = -\left(\frac{R\mu_p(S_1, S_2)}{Y_p}\right)x_2 - \left(\frac{R\mu_a(S_1, S_2)}{Y_a}\gamma^{-1}\right)x_1 \quad \text{in } Q_T, \end{array} \right. \quad (2.1) $$

where

* $D_i, i = 1, 2, 3, 4$ is the diffusion coefficient.

* $k_a, k_p$ denotes respectively the adherent bacteria mortality rate and the planktonic bacteria mortality rate.

* $\beta$ denotes a term corresponding to the detachment rate from the wall.

* A portion of the planktonic category can attach to the walls with a certain rate that is denoted by $\alpha$.

* $G(\bar{x}_1)$ denotes the proportion of daughter cells of the adherent bacteria able to find a place to attach onto the wall, the remainders being washed out by the liquid flow (cf. [11]), where $\bar{x}_1 = \frac{x_1}{\omega_\infty}$.

* $\gamma$ is the coefficient of conversion of the volume density to the surface density.

* $Y_i$, for $i = a, p$, is respectively the coefficient rate of yield of adherent and planktonic 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 degradation of contaminant.

* $\mu_i(S_1, S_2) = \mu_{max} \frac{S_1}{K_{S_1}^i + S_1} \frac{S_2}{K_{S_2}^i + S_2}$, $i = a, p$ represent the growth function(see [1]).

* The nonnegative function $\chi(\cdot)$ denote the chemotaxis sensitivity.

With the following boundary conditions

$$ \left\{ \begin{array}{l} \nabla x_1 \cdot \mathbf{n} = 0 \quad \text{on } \Gamma, \\ -D_2 \nabla x_2 \cdot \mathbf{n} + (\mathbf{u} \cdot \mathbf{n})x_2 = 0 \quad \text{on } \Gamma_1 \cup \Gamma_2, \\ \nabla x_2 \cdot \mathbf{n} = 0 \quad \text{on } \Gamma_3, \\ -D_3 \nabla S_1 \cdot \mathbf{n} + (\mathbf{u} \cdot \mathbf{n})S_1 = (\mathbf{u} \cdot \mathbf{n})S_1^{in}, \quad \text{on } \Gamma_1, \\ -D_3 \nabla S_1 \cdot \mathbf{n} + (\mathbf{u} \cdot \mathbf{n})S_1 = 0 \quad \text{on } \Gamma_2, \\ \nabla S_1 \cdot \mathbf{n} = 0 \quad \text{on } \Gamma_3, \\ -D_4 \nabla S_2 \cdot \mathbf{n} + (\mathbf{u} \cdot \mathbf{n})S_2 = (\mathbf{u} \cdot \mathbf{n})S_2^{in} \quad \text{on } \Gamma_1, \\ -D_4 \nabla S_2 \cdot \mathbf{n} + (\mathbf{u} \cdot \mathbf{n})S_2 = 0 \quad \text{on } \Gamma_2, \\ \nabla S_2 \cdot \mathbf{n} = 0 \quad \text{on } \Gamma_3. \end{array} \right. \quad (2.2) $$

The obtained system is a nonlinear coupled system. For a stationary system in [4] author used the fixed point strategy to show the existence of the steady states. To linearize the advection term, let
$Z = (z_1, z_2, z_3, z_4) \in (H^1(\Omega))^4$ and the system

$$
F_Z: \left\{
\begin{array}{l}
\displaystyle \frac{\partial x_1}{\partial t} - \nabla \cdot (D_1 \nabla x_1) = (\mu_a(S_1, S_2)G(\bar{x}_1) - k_a - \beta)x_1 + \alpha\gamma(1 - \bar{x}_1)x_2 \quad \text{in } Q_T, \\[2ex]
\displaystyle \frac{\partial x_2}{\partial t} - \nabla \cdot (D_2 \nabla x_2) + \mathbf{u} \cdot \nabla x_2 + \nabla \cdot (x_2\chi(z_4)\nabla z_4) = \\
\displaystyle \quad \left(\mu_p(S_1, S_2) - k_p - \alpha(1 - \bar{x}_1)\right)x_2 + \left(\gamma^{-1}\mu_a(S_1, S_2)(1 - G(\bar{x}_1)) + \gamma^{-1}\beta\right)x_1 \quad \text{in } Q_T, \\[2ex]
\displaystyle \frac{\partial S_1}{\partial t} + \mathbf{u} \cdot \nabla S_1 - \nabla \cdot (D_3 \nabla S_1) = -\left(\frac{\mu_p(S_1, S_2)}{Y_p}\right)x_2 - \left(\frac{\mu_a(S_1, S_2)}{Y_a}\gamma^{-1}\right)x_1 \quad \text{in } Q_T, \\[2ex]
\displaystyle \frac{\partial S_2}{\partial t} + \mathbf{u} \cdot \nabla S_2 - \nabla \cdot (D_4 \nabla S_2) = -\left(\frac{R\mu_p(S_1, S_2)}{Y_p}\right)x_2 - \left(\frac{R\mu_a(S_1, S_2)}{Y_a}\gamma^{-1}\right)x_1 \quad \text{in } Q_T.
\end{array}
\right.
\tag{2.3}
$$

which implies

$$
F_Z: \left\{
\begin{aligned}
\frac{\partial x_1}{\partial t} - \nabla \cdot (D_1 \nabla x_1) &= \left( \mu_a(S_1, S_2)G(\bar{x}_1) - k_a - \beta \right) x_1 + \alpha \gamma (1 - \bar{x}_1) x_2 \quad \text{in } Q_T, \\
\frac{\partial x_2}{\partial t} - \nabla \cdot \left( D_2 \nabla x_2 - (\mathbf{u} + \chi(z_4) \nabla z_4) x_2 \right) &= \\
&\quad \left( \mu_p(S_1, S_2) - k_p - \alpha (1 - \bar{x}_1) \right) x_2 + \left( \gamma^{-1} \mu_a(S_1, S_2) (1 - G(\bar{x}_1)) + \gamma^{-1} \beta \right) x_1 \quad \text{in } Q_T, \quad (2.4) \\
\frac{\partial S_1}{\partial t} - \nabla \cdot (D_3 \nabla S_1 - \mathbf{u} S_1) &= - \left( \frac{\mu_p(S_1, S_2)}{Y_p} \right) x_2 - \left( \frac{\mu_a(S_1, S_2)}{Y_a} \gamma^{-1} \right) x_1 \quad \text{in } Q_T, \\
\frac{\partial S_2}{\partial t} - \nabla \cdot (D_4 \nabla S_2 - \mathbf{u} S_2) &= - \left( \frac{R \mu_p(S_1, S_2)}{Y_p} \right) x_2 - \left( \frac{R \mu_a(S_1, S_2)}{Y_a} \gamma^{-1} \right) x_1 \quad \text{in } Q_T.
\end{aligned}
\right.
$$

From now on we adopt the following notations:

$$
\mathbf{C}_Z = (c_1, c_2, c_3, c_4) := (x_1, x_2, S_1, S_2) \quad \text{and} \quad \mathbf{C}_Z^0 = (x_1^0, x_2^0, S_1^0, S_2^0). \tag{2.5}
$$

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

We put:

$$
f_1(x, \mathbf{C}_Z) := \left( \mu_a(c_3(x), c_4(x))G\left(\frac{c_1}{\omega_\infty}(x)\right) - k_a - \beta \right) c_1(x) + \alpha\gamma\left(1 - \frac{c_1}{\omega_\infty}(x, t)\right) c_2(x),
$$

$$
f_2(x, \mathbf{C}_Z) := \left\{
\begin{aligned}
& \left( \mu_p(c_3(x), c_4(x)) - k_p - \alpha \left(1 - \frac{c_1}{\omega_\infty}(x)\right) \right) c_2(x) \\
& + \left( \mu_a(c_3(x), c_4(x)) \left(1 - G\left(\frac{c_1}{\omega_\infty}(x)\right)\right) + \beta \right) \gamma^{-1} c_1(x),
\end{aligned}
\right.
$$

$$
f_3(x, \mathbf{C}_Z) := - \left( \frac{\mu_p(c_3, c_4)}{Y_p} \right) c_2 - \left( \frac{\mu_a(c_3, c_4)}{Y_a} \gamma^{-1} \right) c_1,
$$

$$
f_4(x, \mathbf{C}_Z) := - \left( \frac{R\mu_p(c_3, c_4)}{Y_p} \right) c_2 - \left( \frac{R\mu_a(c_3, c_4)}{Y_a} \gamma^{-1} \right) c_1,
$$

$$
\mathbf{F} := (f_1, f_2, f_3, f_4).
$$

The global fluxes are defined by:

$$
J_1(c_1) := D_1 \nabla c_1,
$$

$$
J_2(c_2) := D_2 \nabla c_2 + u_2 c_2,
$$

$$
J_3(c_3) := D_3 \nabla c_3 - u_3 c_3,
$$

$$
J_4(c_4) := D_4 \nabla c_4 - u_4 c_4,
$$

where $u_2 = \mathbf{u} + \chi(z_4)\nabla z_4$; $u_3 = u_4 = \mathbf{u}$. 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} c_3^{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} c_4^{in} & \text{on } \Gamma_1 \\ 0 & \text{on } \Gamma_2 \cup \Gamma_3. \end{cases}$$

With these notations the quasilinear diffusion-convection system (2.1)-(2.2) becomes:

$$\left\{ \begin{array}{l} \displaystyle \frac{\partial c_i}{\partial t} - \operatorname{div}(J_i(c_i)) = f_i(x, \mathbf{C}) \quad \text{in } Q_T, \text{ for } 1 \le i \le 4, \\ \mathbf{B}(\mathbf{C}) = \mathbf{g} \quad \text{in } \Sigma_T, \\ \mathbf{C}(0, \cdot) = \mathbf{C}^0 \quad \text{in } Q_T. \end{array} \right. \qquad (2.6)$$

and the quasilinear diffusion-convection system(2.3)-(2.2) becomes:

$$\left\{ \begin{array}{l} \displaystyle \frac{\partial c_i}{\partial t} - \operatorname{div}(J_i(c_i)) = f_i(x, \mathbf{C}_Z) \quad \text{in } Q_T, \text{ for } 1 \le i \le 4, \\ \mathbf{B}(\mathbf{C}_Z) = \mathbf{g} \quad \text{in } \Sigma_T, \\ \mathbf{C}_Z(0, \cdot) = \mathbf{C}_Z^0 \quad \text{in } Q_T. \end{array} \right. \qquad (2.7)$$

### 3. Existence and regularity results

**Assumption 3.1** The growth rate of bacteria $\mu_i(x, y)$, for $i = a, p$, satisfies

$$\mu_i \in C^1, \quad \mu(0, y) = \mu(x, 0) = 0, \text{ and } \exists L_{\mu_i} > 0 \\ \text{such that } \forall (x, x', y, y') \in [0, +\infty]^4; \|\mu_i(x, y) - \mu_i(x', y')\|_{\infty} \le L_{\mu_i} (\|x - x'\|_{0,\Omega} + \|y - y'\|_{0,\Omega}).$$

**Remark 1** If we consider the multiplicative formula $\mu_i(S_1, S_2) = \mu_{\max}^i \frac{S_1}{K_{S_1}^i + S_1} \frac{S_2}{K_{S_2}^i + S_2}$, $i = a, p$, then $L_{\mu_i} \le \mu_{\max}^i, i = a, p$.

**Assumption 3.2** The function $G$ satisfies

$$G \in C^1, \quad 0 < G(0) \le 1, \quad G(1) = 0, \text{ and } \exists L_G > 0 \\ \text{such that } \forall (x, y) \in [0, +\infty]^2; \|G(x) - G(y)\|_{\infty} \le L_G \|x - y\|_{0,\Omega}.$$

**Remark 2** Since $G(X) = \frac{1+X}{1+a+X}$ for $X \in [0, \omega_\infty]$ (see [12]) we have

$$G'(X) = \frac{a}{(1 + a + X)^2} \le \frac{a}{(1 + a)^2} < 1,$$

then we can take $L_G = \frac{a}{(1 + a)^2}$.

**Assumption 3.3** The function $\chi$ is a continuous differentiable function satisfying $\chi' \ge 0$.

The existence and uniqueness of the solution of the problem (2.7) is given by the following theorem (Theorem 3.1) given in [4].

**Theorem 3.4** Problem (2.7) has a unique global nonnegative weak solution in the following sense

$$\forall T > 0, \quad \mathbf{C}_Z \text{ belongs to } \left( C([0, T]; L^2(\Omega)) \cap L^\infty(Q_T) \cap L^2(0, T; H^1(\Omega)) \right)^4 \\ \text{and for all } \Psi = (\psi_1, \psi_2, \psi_3, \psi_4) \in \left( C^\infty(\overline{Q}_T) \right)^4 \text{ such that } \Psi(T) = 0$$

$$- \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}_Z) \psi_i + \int_{\Omega} c_i^0 \psi_i(0) + \int_{\Sigma_T} g_i \psi_i. \qquad (3.1)$$

Moreover there exists $M > 0$ depending on system parameters (see [2]) such that

$$\|\mathbf{C}_Z\|_{(L^\infty(Q_T))^4} + \|\mathbf{C}_Z\|_{(L^2(0,T;H^1(\Omega)))^4} + \left\|\frac{\partial \mathbf{C}}{\partial t}\right\|_{(L^2(0,T;H^{-1}(\Omega)))^4} \le M. \qquad (3.2)$$
The problem (3.1) is equivalent to the following system:

Find $\mathbf{C}(t) \in (H^1(\Omega))^4$ for a.e. $t \in [0, T]$ such that and for a.e. $t \in [0, T]$ and for all $\psi_i \in H^1(\Omega)$, $1 \le i \le 4$

$$ \left\{ \begin{array}{l} \displaystyle \int_{\Omega} \frac{dc_i(t)}{dt} \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 (3.3) \\ \mathbf{C}_0 = \mathbf{C}(0). \end{array} \right. $$

## 4. Approximation

Several strategies are available for the treatment of nonlinear terms and for solving the problem (3.3). The choice of a scheme for both time derivatives as well as the treatment of the nonlinear terms will lead to different resolution schemes. The discretization of the problem (3.3) is based on two steps: the space discretization (or semi-discretization) which is made by a finite element method of degree 1, and the time discretization which uses the backward Euler scheme.

### Semi-discretization

Let $\Omega_h = \bigcup_{\mathbf{T} \in \mathbf{T}} \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 \mathbf{T}} \text{diam}(\mathbf{T})$ and let $\mathcal{P}_1(\Omega_h)$ and $\mathcal{V}_h$ defined respectively by:

$$ \mathcal{P}_1(\Omega_h) = \{v \in C(\bar{\Omega}) \mid v|_{\mathbf{T} \in \mathcal{P}_1(\Omega_h)}, \forall \mathbf{T} \in \mathbf{\Omega}_h\} \quad \text{and} \quad \mathcal{V}_h = (\mathcal{P}_1(\Omega_h))^4, $$

where, for each $\mathbf{T} \in \mathbf{\Omega}_h$, $P_1(\mathbf{T})$ stands for the space of restriction to $\mathbf{T}$ of polynomial functions of degree 1. The semi-discrete problem associated to (3.3) is given, for a.e. $t \in [0, T]$, by:

$$ \left\{ \begin{array}{l} \text{find } \mathbf{C}_h(t) \in \mathcal{V}_h \text{ such that for } 1 \le i \le 4 \\ \displaystyle \int_{\Omega} \frac{dc_{ih}(t)}{dt} v_{ih} + \int_{\Omega} J_i(c_{ih}(t)) \cdot \nabla v_{ih} = \int_{\Omega} F_i(\mathbf{C}_h(t)) v_{ih} + \int_{\Gamma} g_i v_{ih}, \quad \forall v_{ih} \in \mathcal{P}_1(\Omega_h), \\ \mathbf{C}_h(0) = \Pi_h \mathbf{C}_0, \end{array} \right. \quad (4.1) $$

where $\Pi_h$ is a projection operator on $\mathcal{V}_h$.

### Full-discretization

Let 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:

$$ \left\{ \begin{array}{l} c_{(ih,k)} := c_{ih}(t_k) \\ \mathbf{C}_h^k := \mathbf{C}_h(t_k) = (c_{(1h,k)}, c_{(2h,k)}, c_{(3h,k)}, c_{(4h,k)}). \end{array} \right. \quad (4.2) $$

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

$$ \frac{\partial c_{ih}}{\partial t}(t_n) \approx \frac{c_{(ih,n)} - c_{(ih,n-1)}}{\tau_n}. $$

With this notation at each time step ($t = t_n$) the problem (3.3), is then fully approximated using the following implicit Euler scheme:

$$ \left\{ \begin{array}{l} \text{for } 1 \le n \le N, \text{ find } \mathbf{C}_h^n \in \mathcal{V} \text{ such that for } 1 \le i \le 4 \\ \displaystyle \int_{\Omega} \frac{c_{(ih,n)} - c_{(ih,n-1)}}{\tau_n} v_{ih} + \int_{\Omega} J_i(c_{(ih,n)}) \cdot \nabla v_{ih} = \int_{\Omega} F_i(\mathbf{C}_h^{n-1}) v_{ih} + \int_{\Gamma} g_i \psi_{ih}, \quad (4.3) \\ \forall v_{ih} \in \mathcal{P}_1(\Omega_h), \\ \mathbf{C}_h^0 = \Pi_h(\mathbf{C}_0). \end{array} \right. $$
5. The iterative scheme analysis

To approximate the solution of problem (4.3), we propose an iterative procedure based on a decoupled computation of the equations. In this section, we focus only on the study of the convergence of the iterative procedure for one time step of the coupled problem (4.3). In the following and for simplification, we omit the spatial step index. The technique used for the proof can be found in [10,16,18]. At each time step ($t = t_n$): the system (4.3) consists to find $C^n \in (H^1(\Omega))^4$ such that $\forall (v_1, v_2, v_3, v_4) \in (H^1(\Omega))^4$,

$$ 
\left\{
\begin{aligned}
& \int_{\Omega} \tau_n D_1 \nabla c_{(1,n)} \cdot \nabla v_1 + \int_{\Omega} \left(1 + \tau_n(k_a + \beta)\right) c_{(1,n)} v_1 \\
& - \int_{\Omega} \tau_n(\mu_a(c_{(3,n)}, c_{(4,n)}) G(\frac{c_{(1,n)}}{\omega_{\infty}}) c_{(1,n)} - \alpha \gamma(1 - \frac{c_{(1,n)}}{\omega_{\infty}}) c_{(2,n)}) v_1 \\
& - \int_{\Omega} c_{(1,n-1)} v_1 = 0, \\
& \int_{\Omega} \tau_n (D_2 \nabla c_{(2,n)} - \mathbf{u} c_{(2,n)}) \cdot \nabla v_2 + \int_{\Omega} \left(1 + \tau_n(k_p + \alpha)\right) c_{(2,n)} v_2 \\
& - \int_{\Omega} \tau_n \left(\mu_p(c_{(3,n)}, c_{(4,n)}) + \alpha \frac{c_{(1,n)}}{\omega_{\infty}}\right) c_{(2,n)} v_2 \\
& - \int_{\Omega} \tau_n \left(\gamma^{-1} \mu_a(c_{(3,n)}, c_{(4,n)}) (1 - G(\frac{c_{(1,n)}}{\omega_{\infty}})) + \gamma^{-1} \beta\right) c_{(1,n)} v_2 \\
& + \int_{\Omega} \tau_n(c_{(2,n)} \chi(c_{(4,n)}) \nabla c_{(4,n)}) \nabla v_2 - \int_{\Omega} c_{(2,n-1)} v_2 = 0, \\
& \int_{\Omega} \tau_n (D_3 \nabla c_{(3,n)} - \mathbf{u} c_{(3,n)}) \cdot \nabla v_3 + \int_{\Omega} c_{(3,n)} v_3 + \int_{\Omega} \left(\frac{\tau_n \mu_p(c_{(3,n)}, c_{(4,n)})}{Y_p}\right) c_{(2,n)} v_3 \\
& + \int_{\Omega} \left(\frac{\tau_n \mu_a(c_{(3,n)}, c_{(4,n)})}{Y_a} \gamma^{-1}\right) c_{(1,n)} v_3 + \int_{\Gamma_1} ((\tau_n \mathbf{u} \cdot n)c_3^{in} - c_{(3,n-1)}) v_3 = 0, \\
& \int_{\Omega} \tau_n (D_4 \nabla c_{(4,n)} - \mathbf{u} c_{(4,n)}) \cdot \nabla v_4 + \int_{\Omega} c_{(4,n)} v_4 + \int_{\Omega} \left(\frac{\tau_n R \mu_p(c_{(3,n)}, c_{(4,n)})}{Y_p}\right) c_{(2,n)} v_4 \\
& + \int_{\Omega} \left(\frac{\tau_n R \mu_a(c_{(3,n)}, c_{(4,n)})}{Y_a} \gamma^{-1}\right) c_{(1,n)} v_4 + \int_{\Gamma_1} (\tau_n (\mathbf{u} \cdot n)c_4^{in} - c_{(4,n-1)}) v_4 = 0.
\end{aligned}
\right.
\quad (5.1) $$

To linearize the reaction term and have decoupled equations, we introduce the iterative algorithm based on the fixed-point method at each time step in order to solve the coupled problem efficiently. For a given initial conditions $C^{(n,0)} = (c_{(1,n)}^0, c_{(2,n)}^0, c_{(3,n)}^0, c_{(4,n)}^0)$, at each time step, we compute $C^{(n,j+1)} = (c_{(1,n)}^{j+1}, c_{(2,n)}^{j+1}, c_{(3,n)}^{j+1}, c_{(4,n)}^{j+1})$, with the following steps:

* Initialisation: $(c_{(1,n)}^0, c_{(2,n)}^0, c_{(3,n)}^0, c_{(4,n)}^0) = (c_{(1,n-1)}, c_{(2,n-1)}, c_{(3,n-1)}, c_{(4,n-1)})$,

* $c_{(1,n)}^{j+1}$ solution of the following linear equation

$$ 
\begin{aligned}
& \int_{\Omega} \tau_n D_1 \nabla c_{(1,n)}^{j+1} \cdot \nabla v_1 + \int_{\Omega} \left(1 + \tau_n(k_a + \beta)\right) c_{(1,n)}^{j+1} v_1 \\
& - \int_{\Omega} \tau_n(\mu_a(c_{(3,n)}^j, c_{(4,n)}^j) G(\frac{c_{(1,n)}^j}{\omega_{\infty}}) c_{(1,n)}^{j+1} - \alpha \gamma(1 - \frac{c_{(1,n)}^j}{\omega_{\infty}}) c_{(2,n)}^j v_1 \\
& - \int_{\Omega} c_{(1,n-1)} v_1 = 0,
\end{aligned}
\quad (5.2) $$
• $c_{(2,n)}^{j+1}$ solution of the following linear equation

$$ 
\begin{aligned}
& \int_{\Omega} \tau_n (D_2 \nabla c_{(2,n)}^{j+1} - \mathbf{u} c_{(2,n)}^{j+1}) \cdot \nabla v_2 + \int_{\Omega} \left(1 + \tau_n (k_p + \alpha)\right) c_{(2,n)}^{j+1} v_2 \\
& - \int_{\Omega} \tau_n \left(\mu_p (c_{(3,n)}^j, c_{(4,n)}^j) + \alpha \left(1 - \frac{c_{(1,n)}^j}{\omega_{\infty}}\right)\right) c_{(2,n)}^{j+1} v_2 \\
& - \int_{\Omega} \tau_n \left(\gamma^{-1} \mu_a (c_{(3,n)}^j, c_{(4,n)}^j) \left(1 - G\left(\frac{c_{(1,n)}^{j+1}}{\omega_{\infty}}\right)\right) + \gamma^{-1} \beta\right) c_{(1,n)}^{j+1} v_2 \\
& + \int_{\Omega} \tau_n c_{(2,n)}^{j+1} \chi(c_{(4,n)}^j) \nabla c_{(4,n)}^j \cdot \nabla v_2 - \int_{\Omega} c_{(2,n-1)} v_2 = 0,
\end{aligned}
\quad (5.3) $$

• $c_{(3,n)}^{j+1}$ solution of the following linear equation

$$ 
\begin{aligned}
& \int_{\Omega} \tau_n (D_3 \nabla c_{(3,n)}^{j+1} - \mathbf{u} c_{(3,n)}^{j+1}) \cdot \nabla v_3 + \int_{\Omega} c_{(3,n)}^{j+1} v_3 + \int_{\Omega} \left( \frac{\tau_n \mu_p (c_{(3,n)}^j, c_{(4,n)}^j)}{Y_p} \right) c_{(2,n)}^{j+1} v_3 \\
& + \int_{\Omega} \left( \frac{\tau_n \mu_a (c_{(3,n)}^j, c_{(4,n)}^j)}{Y_a} \gamma^{-1} \right) c_{(1,n)}^{j+1} v_3 + \int_{\Gamma_1} \left( (\tau_n \mathbf{u} \cdot \mathbf{n}) c_3^{in} - c_{(3,n-1)} \right) v_3 = 0,
\end{aligned}
\quad (5.4) $$

• $c_{(4,n)}^{j+1}$ solution of the following linear equation

$$ 
\begin{aligned}
& \int_{\Omega} \tau_n (D_4 \nabla c_{(4,n)}^{j+1} - \mathbf{u} c_{(4,n)}^{j+1}) \cdot \nabla v_4 + \int_{\Omega} c_{(4,n)}^{j+1} v_4 + \int_{\Omega} \left( \frac{\tau_n R \mu_p (c_{(3,n)}^j, c_{(4,n)}^j)}{Y_p} \right) c_{(2,n)}^{j+1} v_4 \\
& + \int_{\Omega} \left( \frac{\tau_n R \mu_a (c_{(3,n)}^j, c_{(4,n)}^j)}{Y_a} \gamma^{-1} \right) c_{(1,n)}^{j+1} v_4 + \int_{\Gamma_1} (\tau_n (\mathbf{u} \cdot \mathbf{n}) c_4^{in} - c_{(4,n-1)}) v_4 = 0.
\end{aligned}
\quad (5.5) $$

**Assumption 5.1** Assuming that:

• $\mu_{max}^a > k_a + \beta$

• $\mathbf{u}^0 + 2\mu_{max}^p > 2k_p$,

• $\min(D_2, D_3, D_4) > \frac{\mathbf{u}^0}{2}$,

these asymptotics are natural when the bacteria successfully colonize the reactor (see [3]).

Let $\tau_n < \min(\tau^1; \tau^2; \tau^3)$ where

$$ *) \quad \tau^1 = \frac{1}{\mu_{max}^a - (k_a + \beta)}, $$

$$ *) \quad \tau^2 = \frac{2}{\mathbf{u}^0 + 2\mu_{max}^p - 2k_p}, $$

$$ *) \quad \tau^3 = \frac{2}{\mathbf{u}^0}. $$

To show the scheme convergence we proof that under some hypothesis the sequence $C^{(n,j)}$ is contraction and then its convergence.
## 5.1. Coupled Prediction Scheme.

**Theorem 5.2** Assuming that the assumptions (3.1), (3.2), (3.3), and (5.1) hold and using (3.2), the following inequality holds:

$$ \sum_{i=1}^{4} \|c_{(i,n)}^{k+1} - c_{(i,n)}^{m+1}\|_{1,\Omega} \le (\kappa)^m \sum_{i=1}^{4} \alpha_i \|c_{(i,n)}^{k-m} - c_{(i,n)}^{0}\|_{0,\Omega}, \quad (5.6) $$

where $A_2$ is given in (5.13), $X_2$ given in (5.15), $\kappa$ given in (5.24), $\alpha_1 = 1 + \frac{A_2}{X_2}$, $\alpha_2 = 0$ and $\alpha_i = 1$ for $i = 3, 4$. Moreover the choice of $\tau_n$ such that $\kappa < 1$ grantees the strong convergence of the sequence $(c_{(1,n)}^k; c_{(2,n)}^k; c_{(3,n)}^k; c_{(4,n)}^k)$ in $(H^1(\Omega))^4$.

**Proof:** For this we follow the next steps:

* **Step 1:** In (5.2) by the linearity and for $v_1 = c_{(1,n)}^{k+1} - c_{(1,n)}^{m+1}$ and computing the difference between the equation taken at iteration $m+1$ and $k+1$ we get:

$$ \begin{aligned} T_1 &= \int_{\Omega} \tau_n D_1 (\nabla (c_{(1,n)}^{k+1} - c_{(1,n)}^{m+1}))^2 + \int_{\Omega} \left(1 + \tau_n (k_a + \beta)\right) \left(c_{(1,n)}^{k+1} - c_{(1,n)}^{m+1}\right)^2 \\ &= \int_{\Omega} \tau_n \left( \mu_a(c_{(3,n)}^k, c_{(4,n)}^k) G\left(\frac{c_{(1,n)}^k}{\omega_\infty}\right) c_{(1,n)}^{k+1} - \mu_a(c_{(3,n)}^m, c_{(4,n)}^m) G\left(\frac{c_{(1,n)}^m}{\omega_\infty}\right) c_{(1,n)}^{m+1} \right) \left(c_{(1,n)}^{k+1} - c_{(1,n)}^{m+1}\right) \\ &+ \int_{\Omega} \tau_n \alpha \gamma \left( \left(1 - \frac{c_{(1,n)}^k}{\omega_\infty}\right) c_{(2,n)}^k - \left(1 - \frac{c_{(1,n)}^m}{\omega_\infty}\right) c_{(2,n)}^m \right) \left(c_{(1,n)}^{k+1} - c_{(1,n)}^{m+1}\right) \\ &= \int_{\Omega} \tau_n \mu_a(c_{(3,n)}^k, c_{(4,n)}^k) G\left(\frac{c_{(1,n)}^k}{\omega_\infty}\right) \left(c_{(1,n)}^{k+1} - c_{(1,n)}^{m+1}\right)^2 \\ &+ \int_{\Omega} \tau_n \mu_a(c_{(3,n)}^k, c_{(4,n)}^k) c_{(1,n)}^{m+1} \left( G\left(\frac{c_{(1,n)}^k}{\omega_\infty}\right) - G\left(\frac{c_{(1,n)}^m}{\omega_\infty}\right) \right) \left(c_{(1,n)}^{k+1} - c_{(1,n)}^{m+1}\right) \\ &+ \int_{\Omega} \tau_n G\left(\frac{c_{(1,n)}^m}{\omega_\infty}\right) c_{(1,n)}^{m+1} \left( \mu_a(c_{(3,n)}^k, c_{(4,n)}^k) - \mu_a(c_{(3,n)}^m, c_{(4,n)}^m) \right) \left(c_{(1,n)}^{k+1} - c_{(1,n)}^{m+1}\right) \\ &+ \int_{\Omega} \tau_n \alpha \gamma \left( c_{(2,n)}^k - c_{(2,n)}^m \right) \left( c_{(1,n)}^{k+1} - c_{(1,n)}^{m+1} \right) \\ &+ \int_{\Omega} \frac{\tau_n \alpha \gamma}{\omega_\infty} c_{(1,n)}^m \left( c_{(2,n)}^m - c_{(2,n)}^k \right) \left( c_{(1,n)}^{k+1} - c_{(1,n)}^{m+1} \right) \\ &+ \int_{\Omega} \frac{\tau_n \alpha \gamma}{\omega_\infty} c_{(2,n)}^k \left( c_{(1,n)}^m - c_{(1,n)}^k \right) \left( c_{(1,n)}^{k+1} - c_{(1,n)}^{m+1} \right) \end{aligned} $$

which implies that (note that $c_{(1,n)}^k \le \omega_\infty$)

$$ \begin{aligned} T_1 &\le \tau_n \|\mu_a(\cdot, \cdot)\|_{\infty} \|G(\cdot)\|_{\infty} \|c_{(1,n)}^{k+1} - c_{(1,n)}^{m+1}\|_{0,\Omega}^2 \\ &+ \omega_\infty \tau_n \|\mu_a(\cdot, \cdot)\|_{\infty} \int_{\Omega} \left( G\left(\frac{c_{(1,n)}^k}{\omega_\infty}\right) - G\left(\frac{c_{(1,n)}^m}{\omega_\infty}\right) \right) \left(c_{(1,n)}^{k+1} - c_{(1,n)}^{m+1}\right) \\ &+ \omega_\infty \tau_n \|G(\cdot)\|_{\infty} \int_{\Omega} \left( \mu_a(c_{(3,n)}^k, c_{(4,n)}^k) - \mu_a(c_{(3,n)}^m, c_{(4,n)}^m) \right) \left(c_{(1,n)}^{k+1} - c_{(1,n)}^{m+1}\right) \\ &+ \tau_n \alpha \gamma \int_{\Omega} \left(c_{(2,n)}^k - c_{(2,n)}^m\right) \left(c_{(1,n)}^{k+1} - c_{(1,n)}^{m+1}\right) \\ &+ \tau_n \alpha \gamma \int_{\Omega} \left(c_{(2,n)}^m - c_{(2,n)}^k\right) \left(c_{(1,n)}^{k+1} - c_{(1,n)}^{m+1}\right) \\ &+ \frac{M^1 \tau_n \alpha \gamma}{\omega_\infty} \int_{\Omega} \left(c_{(1,n)}^m - c_{(1,n)}^k\right) \left(c_{(1,n)}^{k+1} - c_{(1,n)}^{m+1}\right), \end{aligned} $$
where $M^1 = \|c_{(2,n)}^k\|_{\infty}$ using Assumptions (3.1), (3.2) and (3.3) we have

$$
\begin{aligned}
T_1 \le {}& \tau_n \|\mu_a(\cdot, \cdot)\|_{\infty} \|G(\cdot)\|_{\infty} \|c_{(1,n)}^{k+1} - c_{(1,n)}^{m+1}\|_{0,\Omega}^2 \\
& + \tau_n L_G \|\mu_a(\cdot, \cdot)\|_{\infty} \|c_{(1,n)}^k - c_{(1,n)}^m\|_{0,\Omega} \|c_{(1,n)}^{k+1} - c_{(1,n)}^{m+1}\|_{0,\Omega} \\
& + \omega_{\infty} \tau_n L_{\mu_a} \|G(\cdot)\|_{\infty} \left( \|c_{(3,n)}^k - c_{(3,n)}^m\|_{0,\Omega} + \|c_{(4,n)}^k - c_{(4,n)}^m\|_{0,\Omega} \right) \|c_{(1,n)}^{k+1} - c_{(1,n)}^{m+1}\|_{0,\Omega} \\
& + \tau_n \alpha \gamma \|c_{(2,n)}^k - c_{(2,n)}^m\|_{0,\Omega} \|c_{(1,n)}^{k+1} - c_{(1,n)}^{m+1}\|_{0,\Omega} \\
& + \tau_n \alpha \gamma \|c_{(2,n)}^m - c_{(2,n)}^k\|_{0,\Omega} \|c_{(1,n)}^{k+1} - c_{(1,n)}^{m+1}\|_{0,\Omega} \\
& + \frac{M^1 \tau_n \alpha \gamma}{\omega_{\infty}} \|c_{(1,n)}^m - c_{(1,n)}^k\|_{0,\Omega} \|c_{(1,n)}^{k+1} - c_{(1,n)}^{m+1}\|_{0,\Omega},
\end{aligned}
$$

since $\|\mu_i(\cdot, \cdot)\|_{\infty} < \mu_{max}^i$, $i = a, p$ and $\|G(\cdot)\|_{\infty} < 1$ that implies

$$
\begin{aligned}
T_1 \le {}& \tau_n \mu_{max}^a \|c_{(1,n)}^{k+1} - c_{(1,n)}^{m+1}\|_{1,\Omega}^2 \\
& + (\tau_n L_G \mu_{max}^a + \tau_n \alpha \gamma) \|c_{(1,n)}^k - c_{(1,n)}^m\|_{0,\Omega} \|c_{(1,n)}^{k+1} - c_{(1,n)}^{m+1}\|_{0,\Omega} \\
& + \omega_{\infty} \tau_n L_{\mu_a} \|c_{(3,n)}^k - c_{(3,n)}^m\|_{0,\Omega} \|c_{(1,n)}^{k+1} - c_{(1,n)}^{m+1}\|_{0,\Omega} \\
& + \omega_{\infty} \tau_n L_{\mu_a} \|c_{(4,n)}^k - c_{(4,n)}^m\|_{0,\Omega} \|c_{(1,n)}^{k+1} - c_{(1,n)}^{m+1}\|_{0,\Omega} \\
& + \left( \tau_n \alpha \gamma + \frac{M^1 \tau_n \alpha \gamma}{\omega_{\infty}} \right) \|c_{(2,n)}^k - c_{(2,n)}^m\|_{0,\Omega} \|c_{(1,n)}^{k+1} - c_{(1,n)}^{m+1}\|_{0,\Omega},
\end{aligned}
$$

by Assumption (5.1) we have $\tau_n \mu_{max}^a < 1 + \tau_n(k_a + \beta)$, let

$$A_1 = L_G \mu_{max}^a + \alpha\gamma; \quad B_1 = \alpha\gamma + \frac{M^1 \alpha\gamma}{\omega_{\infty}}; \quad Q_1 = P_1 = \omega_{\infty} L_{\mu_a}, \qquad (5.7)$$

which implies

$$
\begin{aligned}
X_1 \|c_{(1,n)}^{k+1} - c_{(1,n)}^{m+1}\|_{1,\Omega} \le {}& A_1 \|c_{(1,n)}^k - c_{(1,n)}^m\|_{0,\Omega} + B_1 \|c_{(2,n)}^k - c_{(2,n)}^m\|_{0,\Omega} + Q_1 \|c_{(3,n)}^k - c_{(3,n)}^m\|_{0,\Omega} \\
& + P_1 \|c_{(4,n)}^k - c_{(4,n)}^m\|_{0,\Omega},
\end{aligned}
\qquad (5.8)
$$

which implies

$$
\begin{aligned}
\|c_{(1,n)}^{k+1} - c_{(1,n)}^{m+1}\|_{1,\Omega} \le {}& \frac{A_1}{X_1} \|c_{(1,n)}^k - c_{(1,n)}^m\|_{0,\Omega} + \frac{B_1}{X_1} \|c_{(2,n)}^k - c_{(2,n)}^m\|_{0,\Omega} + \frac{Q_1}{X_1} \|c_{(3,n)}^k - c_{(3,n)}^m\|_{0,\Omega} \\
& + \frac{P_1}{X_1} \|c_{(4,n)}^k - c_{(4,n)}^m\|_{0,\Omega},
\end{aligned}
\qquad (5.9)
$$

where

$$X_1 = \min\left(D_1; (\tau_n)^{-1} + k_a + \beta - \mu_{max}^a\right). \qquad (5.10)$$

For the partial strategy for $i = 2, 3, 4$ we take $k = m$ then

$$X_1 \|c_{(1,n)}^{k+1} - c_{(1,n)}^{m+1}\|_{1,\Omega} \le A_1 \|c_{(1,n)}^k - c_{(1,n)}^m\|_{0,\Omega}, \qquad (5.11)$$

which implies

$$\|c_{(1,n)}^{k+1} - c_{(1,n)}^{m+1}\|_{1,\Omega} \le \frac{A_1}{X_1} \|c_{(1,n)}^k - c_{(1,n)}^m\|_{0,\Omega},$$

which implies

$$\|c_{(1,n)}^{k+1} - c_{(1,n)}^{m+1}\|_{1,\Omega} \le (\frac{A_1}{X_1})^m \|c_{(1,n)}^{k-m} - c_{(1,n)}^0\|_{0,\Omega}, \qquad (5.12)$$

then the algorithm converges if $\frac{A_1}{X_1} < 1$.
* **Step 2:** In (5.3) by the linearity and the fact that $1 - \frac{c_1}{\omega_\infty} \le 1$ and with the fact that the conductivity decreases over time (see [1,2,4]) then the velocity of the flow also decreases over time then $\|\mathbf{u}\|_\infty \le \mathbf{u}^0$, where $\mathbf{u}^0$ is the initial velocity. For $v_2 = c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1}$ and computing the difference between the equation taken at iteration $m+1$ and $k+1$ and since $\frac{c_1}{\omega_\infty} < 1$ we get:

$$
\begin{aligned}
T_2 &= \int_{\Omega} \tau_n D_2 \left( \nabla (c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1}) \right)^2 + \int_{\Omega} \left( 1 + \tau_n (k_p + \alpha) \right) \left( c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1} \right)^2 \\
&\le \int_{\Omega} \tau_n \mathbf{u} \left( c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1} \right) \nabla \left( c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1} \right) \\
&\quad + \int_{\Omega} \tau_n \left( \mu_p (c_{(3,n)}^k, c_{(4,n)}^k) c_{(2,n)}^{k+1} - \mu_p (c_{(3,n)}^m, c_{(4,n)}^m) c_{(2,n)}^{m+1} \right) \left( c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1} \right) \\
&\quad + \int_{\Omega} \tau_n \alpha \left( c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1} \right)^2 \\
&\quad + \int_{\Omega} \tau_n \gamma^{-1} \left( \mu_a (c_{(3,n)}^k, c_{(4,n)}^k) c_{(1,n)}^{k+1} - \mu_a (c_{(3,n)}^m, c_{(4,n)}^m) c_{(1,n)}^{m+1} \right) \left( c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1} \right) \\
&\quad + \int_{\Omega} \tau_n \gamma^{-1} \left( \mu_a (c_{(3,n)}^m, c_{(4,n)}^m) G\left(\frac{c_{(1,n)}^{m+1}}{\omega_\infty}\right) c_{(1,n)}^{m+1} - \mu_a (c_{(3,n)}^k, c_{(4,n)}^k) G\left(\frac{c_{(1,n)}^{k+1}}{\omega_\infty}\right) c_{(1,n)}^{k+1} \right) \left( c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1} \right) \\
&\quad + \int_{\Omega} \tau_n \gamma^{-1} \beta \left( c_{(1,n)}^{k+1} - c_{(1,n)}^{m+1} \right) \left( c_{2,n}^{k+1} - c_{2,n}^{m+1} \right) \\
&\quad + \int_{\Omega} \tau_n \left( c_{(2,n)}^{k+1} \chi(c_{(4,n)}^k) \nabla c_{(4,n)}^k - c_{(2,n)}^{m+1} \chi(c_{(4,n)}^m) \nabla c_{(4,n)}^m \right) \nabla \left( c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1} \right),
\end{aligned}
$$

which implies

$$
\begin{aligned}
T_2 &= \int_{\Omega} \tau_n D_2 \left( \nabla (c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1}) \right)^2 + \int_{\Omega} \left( 1 + \tau_n (k_p + \alpha) \right) \left( c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1} \right)^2 \\
&\le \int_{\Omega} \tau_n \mathbf{u} \left( c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1} \right) \nabla \left( c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1} \right) \\
&\quad + \int_{\Omega} \tau_n \left( \mu_p (c_{(3,n)}^k, c_{(4,n)}^k) c_{(2,n)}^{k+1} - \mu_p (c_{(3,n)}^m, c_{(4,n)}^m) c_{(2,n)}^{m+1} \right) \left( c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1} \right) \\
&\quad + \int_{\Omega} \tau_n \alpha \left( c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1} \right)^2 \\
&\quad + \int_{\Omega} \tau_n \gamma^{-1} \left( \mu_a (c_{(3,n)}^k, c_{(4,n)}^k) c_{(1,n)}^{k+1} - \mu_a (c_{(3,n)}^m, c_{(4,n)}^m) c_{(1,n)}^{m+1} \right) \left( c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1} \right) \\
&\quad + \int_{\Omega} \tau_n \gamma^{-1} \left( \mu_a (c_{(3,n)}^m, c_{(4,n)}^m) G\left(\frac{c_{(1,n)}^{m+1}}{\omega_\infty}\right) c_{(1,n)}^{m+1} - \mu_a (c_{(3,n)}^k, c_{(4,n)}^k) G\left(\frac{c_{(1,n)}^{k+1}}{\omega_\infty}\right) c_{(1,n)}^{k+1} \right) \left( c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1} \right) \\
&\quad + \int_{\Omega} \tau_n \gamma^{-1} \beta \left( c_{(1,n)}^{k+1} - c_{(1,n)}^{m+1} \right) \left( c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1} \right) \\
&\quad + \int_{\Omega} \tau_n \left( c_{(2,n)}^{k+1} \chi(c_{(4,n)}^k) \nabla c_{(4,n)}^k - c_{(2,n)}^{m+1} \chi(c_{(4,n)}^m) \nabla c_{(4,n)}^m \right) \nabla \left( c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1} \right),
\end{aligned}
$$
which implies

$$
\begin{aligned}
T_2 &= \int_{\Omega} \tau_n D_2 \left( \nabla (c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1}) \right)^2 + \int_{\Omega} \left( 1 + \tau_n (k_p + \alpha) \right) \left( c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1} \right)^2 \\
&\leq \frac{\tau_n \mathbf{u}^0}{2} \left( \|c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1}\|_{0,\Omega}^2 + \|\nabla(c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1})\|_{0,\Omega}^2 \right) + \tau_n \mu_{max}^p \|c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1}\|_{0,\Omega}^2 \\
&\quad + M_2 \tau_n L_{\mu_p} \left( \|c_{(3,n)}^k - c_{(3,n)}^m\|_{0,\Omega} + \|c_{(4,n)}^k - c_{(4,n)}^m\|_{0,\Omega} \right) \|c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1}\|_{0,\Omega} \\
&\quad + \tau_n \alpha \|c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1}\|_{0,\Omega}^2 + \tau_n \gamma^{-1} \mu_{max}^a \|c_{1,n}^{k+1} - c_{1,n}^{m+1}\|_{0,\Omega} \|c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1}\|_{0,\Omega} \\
&\quad + \omega_{\infty} \tau_n \gamma^{-1} L_{\mu_a} \left( \|c_{(3,n)}^k - c_{(3,n)}^m\|_{0,\Omega} + \|c_{(4,n)}^k - c_{(4,n)}^m\|_{0,\Omega} \right) \|c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1}\|_{0,\Omega} \\
&\quad + \tau_n \gamma^{-1} \mu_{max}^a \|c_{1,n}^{k+1} - c_{1,n}^{m+1}\|_{0,\Omega} \|c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1}\|_{0,\Omega} \\
&\quad + \gamma^{-1} \tau_n L_G \mu_{max}^a \|c_{(1,n)}^k - c_{(1,n)}^m\|_{0,\Omega} \|c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1}\|_{0,\Omega} \\
&\quad + \omega_{\infty} \gamma^{-1} \tau_n L_{\mu_a} \left( \|c_{(3,n)}^k - c_{(3,n)}^m\|_{0,\Omega} + \|c_{(4,n)}^k - c_{(4,n)}^m\|_{0,\Omega} \right) \|c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1}\|_{0,\Omega} \\
&\quad + \tau_n \gamma^{-1} \beta \|c_{(1,n)}^{k+1} - c_{(1,n)}^{m+1}\|_{0,\Omega} \|c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1}\|_{0,\Omega} \\
&\quad + M_2 \tau_n \|\chi(\cdot)\|_{\infty} \|\nabla c_{(4,n)}^k - \nabla c_{(4,n)}^m\|_{0,\Omega} \|\nabla(c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1})\|_{0,\Omega} \\
&\quad + (M_2)^2 \tau_n L_{\chi} \|c_{(4,n)}^k - c_{(4,n)}^m\|_{0,\Omega} \|\nabla(c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1})\|_{0,\Omega} \\
&\quad + M_2 \tau_n \|\chi(\cdot)\|_{\infty} \|c_{(1,n)}^{k+1} - c_{(1,n)}^{m+1}\|_{0,\Omega} \|\nabla(c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1})\|_{0,\Omega},
\end{aligned}
$$

where $M_2 = \|c_{(2,n)}^{k+1}\|_{\infty}$. For the partial strategy for $i=3,4$ we take $k=m$ then which implies

$$
\begin{aligned}
T_2 &\le \frac{\tau_n \mathbf{u}^0}{2} \left( \|c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1}\|_{0,\Omega}^2 + \|\nabla(c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1})\|_{0,\Omega}^2 \right) + \tau_n \mu_{max}^p \|c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1}\|_{0,\Omega}^2 \\
&\quad + \tau_n \alpha \|c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1}\|_{0,\Omega}^2 + \frac{2A_1 \tau_n \gamma^{-1} \mu_{max}^a}{X_1} \|c_{1,n}^k - c_{1,n}^m\|_{0,\Omega} \|c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1}\|_{0,\Omega} \\
&\quad + \gamma^{-1} \tau_n L_G \mu_{max}^a \|c_{(1,n)}^k - c_{(1,n)}^m\|_{0,\Omega} \|c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1}\|_{0,\Omega} \\
&\quad + \frac{A_1 \tau_n \gamma^{-1} \beta}{X_1} \|c_{1,n}^k - c_{1,n}^m\|_{0,\Omega} \|c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1}\|_{0,\Omega} \\
&\quad + \frac{A_1 M_2 \tau_n \|\chi(\cdot)\|_{\infty}}{X_1} \|c_{1,n}^k - c_{1,n}^m\|_{0,\Omega} \|\nabla(c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1})\|_{0,\Omega}.
\end{aligned}
$$

Using Assumption (5.1), let

$$
\begin{gather}
\zeta_1 = 1 + \tau_n k_p - \left( \frac{\tau_n \mathbf{u}^0}{2} + \tau_n \mu_{max}^p \right) > 0 \nonumber \\
A_2 = \frac{2A_1 \gamma^{-1} \mu_{max}^a + A_1 \gamma^{-1} \beta + A_1 M_2 \|\chi(\cdot)\|_{\infty}}{X_1} + \gamma^{-1} L_G \mu_{max}^a \tag{5.13}
\end{gather}
$$

which implies

$$
\|c_{(2,n)}^{k+1} - c_{(2,n)}^{m+1}\|_{1,\Omega} \le \frac{A_2}{X_2} \left(\frac{A_1}{X_1}\right)^m \|c_{(1,n)}^{k-m} - c_{(1,n)}^0\|_{0,\Omega}, \quad (5.14)
$$

where

$$
X_2 = \min \left( D_2 - \frac{\mathbf{u}^0}{2}; \tau_n^{-1} \zeta_1 \right), \qquad (5.15)
$$

then the algorithm converges if $\frac{A_1}{X_1} < 1$.
* **Step 3:** In (5.4), by the linearity, then for $v_3 = c_{(3,n)}^{k+1} - c_{(3,n)}^{m+1}$ and computing the difference between the equation taken at iteration $m+1$ and $k+1$ we get:

$$ 
\begin{aligned}
T_3 &= \int_{\Omega} \tau_n D_3 \left( \nabla(c_{(3,n)}^{k+1} - c_{(3,n)}^{m+1}) \right)^2 + \int_{\Omega} \left( c_{(3,n)}^{k+1} - c_{(3,n)}^{m+1} \right)^2 \\
&= \int_{\Omega} \tau_n \mathbf{u} \left( c_{(3,n)}^{k+1} - c_{(3,n)}^{m+1} \right) \nabla \left( c_{(3,n)}^{k+1} - c_{(3,n)}^{m+1} \right) \\
&\quad + \frac{\tau_n}{Y_p} \int_{\Omega} \left( \mu_p(c_{(3,n)}^m, c_{(4,n)}^m) c_{(2,n)}^{m+1} - \mu_p(c_{(3,n)}^k, c_{(4,n)}^k) c_{(2,n)}^{k+1} \right) \left( c_{(3,n)}^{k+1} - c_{(3,n)}^{m+1} \right) \\
&\quad + \frac{\tau_n \gamma^{-1}}{Y_a} \int_{\Omega} \left( \mu_p(c_{(3,n)}^k, c_{(4,n)}^k) c_{(1,n)}^{k+1} - \mu_p(c_{(3,n)}^m, c_{(4,n)}^m) c_{(1,n)}^{m+1} \right) \left( c_{(3,n)}^{k+1} - c_{(3,n)}^{m+1} \right),
\end{aligned}
$$

For the partial strategy for $i=4$ we take $k=m$ then which implies

$$ 
\begin{aligned}
T_3 &= \int_{\Omega} \tau_n D_3 \left( \nabla(c_{(3,n)}^{k+1} - c_{(3,n)}^{m+1}) \right)^2 + \int_{\Omega} \left( c_{(3,n)}^{k+1} - c_{(3,n)}^{m+1} \right)^2 \\
&\leq \frac{\tau_n \mathbf{u}^0}{2} \left( \|c_{(3,n)}^{k+1} - c_{(3,n)}^{m+1}\|_{0,\Omega}^2 + \|\nabla(c_{(3,n)}^{k+1} - c_{(3,n)}^{m+1})\|_{0,\Omega}^2 \right) \\
&\quad + \frac{\tau_n}{Y_p} \mu_{max}^p \|c_{(2,n)}^{m+1} - c_{(2,n)}^{k+1}\|_{0,\Omega} \|c_{(3,n)}^{k+1} - c_{(3,n)}^{m+1}\|_{0,\Omega} \\
&\quad + \frac{M_2 \tau_n L_{\mu_p}}{Y_p} \|c_{(3,n)}^k - c_{(3,n)}^m\|_{0,\Omega} \|c_{(3,n)}^{k+1} - c_{(3,n)}^{m+1}\|_{0,\Omega} \\
&\quad + \frac{\tau_n \gamma^{-1}}{Y_a} \mu_{max}^a \|c_{(1,n)}^{m+1} - c_{(1,n)}^{k+1}\|_{0,\Omega} \|c_{(3,n)}^{k+1} - c_{(3,n)}^{m+1}\|_{0,\Omega} \\
&\quad + \frac{\omega_{\infty} \gamma^{-1} \tau_n L_{\mu_a}}{Y_a} \|c_{(3,n)}^k - c_{(3,n)}^m\|_{0,\Omega} \|c_{(3,n)}^{k+1} - c_{(3,n)}^{m+1}\|_{0,\Omega} \\
&\leq \frac{\tau_n \mathbf{u}^0}{2} \left( \|c_{(3,n)}^{k+1} - c_{(3,n)}^{m+1}\|_{0,\Omega}^2 + \|\nabla(c_{(3,n)}^{k+1} - c_{(3,n)}^{m+1})\|_{0,\Omega}^2 \right) \\
&\quad + \frac{A_2 \tau_n}{X_2 Y_p} \mu_{max}^p \|c_{(1,n)}^m - c_{(1,n)}^k\|_{0,\Omega} \|c_{(3,n)}^{k+1} - c_{(3,n)}^{m+1}\|_{0,\Omega} \\
&\quad + \frac{M_2 \tau_n L_{\mu_p}}{Y_p} \|c_{(3,n)}^k - c_{(3,n)}^m\|_{0,\Omega} \|c_{(3,n)}^{k+1} - c_{(3,n)}^{m+1}\|_{0,\Omega} \\
&\quad + \frac{A_1 \tau_n \gamma^{-1}}{X_1 Y_a} \mu_{max}^a \|c_{(1,n)}^m - c_{(1,n)}^k\|_{0,\Omega} \|c_{(3,n)}^{k+1} - c_{(3,n)}^{m+1}\|_{0,\Omega} \\
&\quad + \frac{\omega_{\infty} \gamma^{-1} \tau_n L_{\mu_a}}{Y_a} \|c_{(3,n)}^k - c_{(3,n)}^m\|_{0,\Omega} \|c_{(3,n)}^{k+1} - c_{(3,n)}^{m+1}\|_{0,\Omega}.
\end{aligned}
$$

Using Assumption (5.1), let

$$ A_3 = \frac{A_2}{X_2 Y_p} \mu_{max}^p + \frac{A_1 \gamma^{-1}}{X_1 Y_a} \mu_{max}^p; \quad Q_3 = \frac{M_2 L_{\mu_p}}{Y_p} + \frac{\omega_{\infty} \gamma^{-1} L_{\mu_a}}{Y_a}, \qquad (5.16) $$

which implies

$$ \|c_{(3,n)}^{k+1} - c_{(3,n)}^{m+1}\|_{1,\Omega} \leq \frac{A_3}{X_3} \|c_{(1,n)}^k - c_{(1,n)}^m\|_{0,\Omega} + \frac{Q_3}{X_3} \|c_{(3,n)}^k - c_{(3,n)}^m\|_{0,\Omega}. \qquad (5.17) $$

With the asymptotic $\frac{A_1}{X_1} < 1$, $\|c_{(1,n)}^k - c_{(1,n)}^m\|_{0,\Omega} \to 0$ this implies

$$ \|c_{(3,n)}^{k+1} - c_{(3,n)}^{m+1}\|_{1,\Omega} \leq \left(\frac{Q_3}{X_3}\right)^m \|c_{(3,n)}^{k-m} - c_{(3,n)}^0\|_{0,\Omega}, \qquad (5.18) $$
where

$$X_3 = \min \left( D_3 - \frac{\mathbf{u}^0}{2}; \tau_n^{-1} - \frac{\mathbf{u}^0}{2} \right) = \min \left( D_3; \tau_n^{-1} \right) - \frac{\mathbf{u}^0}{2}, \quad (5.19)$$

the algorithm converges if $\frac{Q_3}{X_3} < 1$ then a sufficient condition for convergence is given by $\frac{Q_3}{X_3} < 1$.

* **Step 4:** In (5.5), by the linearity, then for $v_4 = c_{(4,n)}^{k+1} - c_{(4,n)}^{m+1}$ and computing the difference between the equation taken at iteration $m+1$ and $k+1$ we get:

$$ 
\begin{aligned}
T_4 &= \int_{\Omega} \tau_n D_4 \left( \nabla(c_{(4,n)}^{k+1} - c_{(4,n)}^{m+1}) \right)^2 + \int_{\Omega} \left( c_{(4,n)}^{k+1} - c_{(4,n)}^{m+1} \right)^2 \\
&= \int_{\Omega} \tau_n \mathbf{u} \left( c_{(4,n)}^{k+1} - c_{(4,n)}^{m+1} \right) \nabla \left( c_{(4,n)}^{k+1} - c_{(4,n)}^{m+1} \right) \\
&\quad + \frac{\tau_n R}{Y_p} \int_{\Omega} \left( \mu_p(c_{(3,n)}^m, c_{(4,n)}^m) c_{(2,n)}^{m+1} - \mu_p(c_{(3,n)}^k, c_{(4,n)}^k) c_{(2,n)}^{k+1} \right) \left( c_{(3,n)}^{k+1} - c_{(3,n)}^{m+1} \right) \\
&\quad + \frac{\tau_n R \gamma^{-1}}{Y_a} \int_{\Omega} \left( \mu_p(c_{(3,n)}^k, c_{(4,n)}^k) c_{(1,n)}^{k+1} - \mu_p(c_{(3,n)}^m, c_{(4,n)}^m) c_{(1,n)}^{m+1} \right) \left( c_{(3,n)}^{k+1} - c_{(3,n)}^{m+1} \right),
\end{aligned}
$$

which implies

$$ 
\begin{aligned}
T_4 &= \int_{\Omega} \tau_n D_4 \left( \nabla(c_{(4,n)}^{k+1} - c_{(4,n)}^{m+1}) \right)^2 + \int_{\Omega} \left( c_{(4,n)}^{k+1} - c_{(4,n)}^{m+1} \right)^2 \\
&\le \frac{\tau_n \mathbf{u}^0}{2} \left( \|c_{(4,n)}^{k+1} - c_{(4,n)}^{m+1}\|_{0,\Omega}^2 + \|\nabla(c_{(4,n)}^{k+1} - c_{(4,n)}^{m+1})\|_{0,\Omega}^2 \right) \\
&\quad + \frac{\tau_n R}{Y_p} \mu_{max}^p \|c_{(2,n)}^{m+1} - c_{(2,n)}^{k+1}\|_{0,\Omega} \|c_{(4,n)}^{k+1} - c_{(4,n)}^{m+1}\|_{0,\Omega} \\
&\quad + \frac{M_2 R \tau_n L_{\mu_p}}{Y_p} \left( \|c_{(3,n)}^k - c_{(3,n)}^m\|_{0,\Omega} + \|c_{(4,n)}^k - c_{(4,n)}^m\|_{0,\Omega} \right) \|c_{(4,n)}^{k+1} - c_{(4,n)}^{m+1}\|_{0,\Omega} \\
&\quad + \frac{R \tau_n \gamma^{-1}}{Y_a} \mu_{max}^a \|c_{(1,n)}^{m+1} - c_{(1,n)}^{k+1}\|_{0,\Omega} \|c_{(4,n)}^{k+1} - c_{(4,n)}^{m+1}\|_{0,\Omega} \\
&\quad + \frac{\omega_{\infty} R \gamma^{-1} \tau_n L_{\mu_a}}{Y_a} \left( \|c_{(3,n)}^k - c_{(3,n)}^m\|_{0,\Omega} + \|c_{(4,n)}^k - c_{(4,n)}^m\|_{0,\Omega} \right) \|c_{(4,n)}^{k+1} - c_{(4,n)}^{m+1}\|_{0,\Omega} \\
&\le \frac{\tau_n \mathbf{u}^0}{2} \left( \|c_{(4,n)}^{k+1} - c_{(4,n)}^{m+1}\|_{0,\Omega}^2 + \|\nabla(c_{(4,n)}^{k+1} - c_{(4,n)}^{m+1})\|_{0,\Omega}^2 \right) \\
&\quad + \frac{A_2 \tau_n R}{X_2 Y_p} \mu_{max}^p \|c_{(1,n)}^m - c_{(1,n)}^k\|_{0,\Omega} \|c_{(4,n)}^{k+1} - c_{(4,n)}^{m+1}\|_{0,\Omega} \\
&\quad + \frac{M_2 R \tau_n L_{\mu_p}}{Y_p} \left( \|c_{(3,n)}^k - c_{(3,n)}^m\|_{0,\Omega} + \|c_{(4,n)}^k - c_{(4,n)}^m\|_{0,\Omega} \right) \|c_{(4,n)}^{k+1} - c_{(4,n)}^{m+1}\|_{0,\Omega} \\
&\quad + \frac{A_1 R \tau_n \gamma^{-1}}{X_1 Y_a} \mu_{max}^a \|c_{(1,n)}^m - c_{(1,n)}^k\|_{0,\Omega} \|c_{(4,n)}^{k+1} - c_{(4,n)}^{m+1}\|_{0,\Omega} \\
&\quad + \frac{\omega_{\infty} R \gamma^{-1} \tau_n L_{\mu_a}}{Y_a} \left( \|c_{(3,n)}^k - c_{(3,n)}^m\|_{0,\Omega} + \|c_{(4,n)}^k - c_{(4,n)}^m\|_{0,\Omega} \right) \|c_{(4,n)}^{k+1} - c_{(4,n)}^{m+1}\|_{0,\Omega}.
\end{aligned}
$$

Using Assumption (5.1), let

$$A_4 = R A_3, \quad Q_4 = R Q_3, \quad P_4 = Q_4, \qquad (5.20)$$

which implies

$$ \|c_{(4,n)}^{k+1} - c_{(4,n)}^{m+1}\|_{1,\Omega} \le \frac{A_4}{X_4} \|c_{(1,n)}^k - c_{(1,n)}^m\|_{0,\Omega} + \frac{Q_4}{X_4} \|c_{(3,n)}^k - c_{(3,n)}^m\|_{0,\Omega} + \frac{P_4}{X_4} \|c_{(4,n)}^k - c_{(4,n)}^m\|_{0,\Omega}. \quad (5.21) $$
With asymptotics $\frac{A_1}{X_1} < 1$ and $\frac{Q_3}{X_3} < 1$ we have $\|c_{(1,n)}^k - c_{(1,n)}^m\|_{0,\Omega} \to 0$ and $\|c_{(3,n)}^k - c_{(3,n)}^m\|_{0,\Omega} \to 0$ then

$$ \|c_{(4,n)}^{k+1} - c_{(4,n)}^{m+1}\|_{1,\Omega} \le \left(\frac{P_4}{X_4}\right)^m \|c_{(4,n)}^{k_m} - c_{(4,n)}^{0}\|_{0,\Omega} \le \left(R\frac{Q_3}{X_4}\right)^m \|c_{(4,n)}^{k_m} - c_{(4,n)}^{0}\|_{0,\Omega}, \quad (5.22) $$

where

$$ X_4 = \min \left( D_4 - \frac{\mathbf{u}^0}{2}; \tau_n^{-1} - \frac{\mathbf{u}^0}{2} \right) = \min (D_4; \tau_n^{-1}) - \frac{\mathbf{u}^0}{2}. \quad (5.23) $$

By (5.12), (5.14), (5.18) and (5.22), we obtain

$$ \begin{aligned} \sum_{i=1}^{4} \|c_{(i,n)}^{k+1} - c_{(i,n)}^{m+1}\|_{1,\Omega} &\le \left(\frac{A_1}{X_1}\right)^m \left(1 + \frac{A_2}{X_2}\right) \|c_{(1,n)}^{k-m} - c_{(1,n)}^{0}\|_{0,\Omega} + \\ &\quad + \left(\frac{Q_3}{X_3}\right)^m \|c_{(3,n)}^{k_m} - c_{(3,n)}^{0}\|_{0,\Omega} + \left(R\frac{Q_3}{X_4}\right)^m \|c_{(4,n)}^{k-m} - c_{(4,n)}^{0}\|_{0,\Omega}. \end{aligned} $$

Let

$$ \kappa = \max \left( \frac{A_1}{X_1}; \frac{Q_3}{X_3}; R \frac{Q_3}{X_4} \right) \quad (5.24) $$

**Theorem 5.3** Assuming that the assumptions (3.1), (3.2), (3.3), and (5.1) hold and using (3.2), the following inequality holds:

$$ \sum_{i=1}^{4} \|c_{(i,n)}^{k+1} - c_{(i,n)}^{m+1}\|_{1,\Omega} \le (\kappa)^m \sum_{i=1}^{4} \alpha_i \|c_{(i,n)}^{k-m} - c_{(i,n)}^{0}\|_{0,\Omega}, \quad (5.25) $$

where $\kappa$ is given in (5.24) and $\alpha_1 = 1 + \frac{A_2}{X_2}$, $\alpha_2 = 0$ and $\alpha_i = 1$ for $i = 3, 4$. Moreover the choice of $\tau_n$ such that $\kappa < 1$ guarantees the strong convergence of the sequence $(c_{(1,n)}^k; c_{(2,n)}^k; c_{(3,n)}^k; c_{(4,n)}^k)$ in $(H^1(\Omega))^4$.

If $\kappa < 1$ the sequence $(C^{(n,k)})_{k \in \mathbb{N}}$ is a Cauchy sequence in $(H^1(\Omega))^4$ and it converges to $C^n$ in $(H^1(\Omega))^4$ strongly. Since the functions $\mu_a, \mu_p, G$ and $\chi$ are continuous and bounded we have a.e in $\Omega$

* $\lim_{k \to \infty} \mu_a(c_{(3,n)}^k, c_{(4,n)}^k) = \mu_a(c_{(3,n)}, c_{(4,n)})$,

* $\lim_{k \to \infty} \mu_p(c_{(3,n)}^k, c_{(4,n)}^k) = \mu_p(c_{(3,n)}, c_{(4,n)})$,

* $\lim_{k \to \infty} G(\frac{c_{(1,n)}^k}{\omega_\infty}) = G(\frac{c_{(1,n)}^k}{\omega_\infty})$,

* $\lim_{k \to \infty} \chi(c_{(4,n)}^k) = \chi(c_{(4,n)})$,

then $C^n$ is solution of (5.1).

□

## 5.2. Fixed-point method.

**Theorem 5.4** Assuming that the assumptions (3.1), (3.2), (3.3), and (5.1) hold and using (3.2), the following inequality holds:

$$ \|C^{n,k+1} - C^{n,m+1}\|_{1,\Omega} \le (\kappa^1)^m \|C^{n,k-m} - C^{n,0}\|_{0,\Omega}, \quad (5.26) $$

where $\kappa^1$ is given in (5.30). Moreover if $\kappa^1 < 1$ then the sequence $(C^{n,k})_{k \in \mathbb{N}}$ convergence strongly in $(H^1(\Omega))^4$.
**Proof:** Using (5.9), (5.10), (5.17) and (5.21) we have

$$ \sum_{i=1}^{4} \|c_{(i,n)}^{k+1} - c_{(i,n)}^{m+1}\|_{1,\Omega} \leq \sum_{i=1}^{4} \beta_i \|c_{(1,n)}^{k} - c_{(1,n)}^{m}\|_{0,\Omega}, \quad (5.27) $$

where

$$ \beta_1 = \sum_{i=1}^{4} \frac{A_i}{X_i}; \quad \beta_2 = 0; \quad \beta_3 = \frac{Q_3}{X_3} + \frac{Q_4}{X_4}; \quad \beta_4 = \frac{RQ_3}{X_4}, $$

which implies

$$ \sum_{i=1}^{4} \|c_{(i,n)}^{k+1} - c_{(i,n)}^{m+1}\|_{1,\Omega} \leq \max_{1 \leq i \leq 4} \beta_i \sum_{i=1}^{4} \|c_{(1,n)}^{k} - c_{(1,n)}^{m}\|_{0,\Omega}, \quad (5.28) $$

which implies

$$ \|\mathbf{C}^{k+1} - \mathbf{C}^{m+1}\|_{1,\Omega} \leq \kappa^1 \|\mathbf{C}^k - \mathbf{C}^m\|_{0,\Omega}, \quad (5.29) $$

where

$$ \kappa^1 = \max_{1 \leq i \leq 4} \beta_i, \quad (5.30) $$

which implies

$$ \|\mathbf{C}^{n,k+1} - \mathbf{C}^{n,m+1}\|_{1,\Omega} \leq (\kappa^1)^m \|\mathbf{C}^{n,k-m} - \mathbf{C}^{n,0}\|_{0,\Omega}. \quad (5.31) $$

The sequence $(\mathbf{C}^{(n,k)})_{k \in \mathbb{N}}$ is a Cauchy sequence in $(H^1(\Omega))^4$ and it converges to $\mathbf{C}^n$ in $(H^1(\Omega))^4$ strongly. Since the functions $\mu_a, \mu_p, G$ and $\chi$ are continuous and bounded we have a.e in $\Omega$

$$ \lim_{k \to \infty} \mu_a(c_{(3,n)}^k, c_{(4,n)}^k) = \mu_a(c_{(3,n)}, c_{(4,n)}), $$

$$ \lim_{k \to \infty} \mu_p(c_{(3,n)}^k, c_{(4,n)}^k) = \mu_p(c_{(3,n)}, c_{(4,n)}), $$

$$ \lim_{k \to \infty} G\left(\frac{c_{(1,n)}^k}{\omega_\infty}\right) = G\left(\frac{c_{(1,n)}^k}{\omega_\infty}\right), $$

$$ \lim_{k \to \infty} \chi(c_{(4,n)}^k) = \chi(c_{(4,n)}), $$

Then $\mathbf{C}^n$ solution of (5.1).

□

**Corollary 5.5** By (5.24) and (5.30) we have $\kappa \leq \kappa^1$ then Coupled Prediction Scheme (CPS) converges quickly than the Fixed-Point Methods (FPM).

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} = \mathbf{u}^0$ on $\Gamma_1 \cup \Gamma_2$ such that

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

where $H(\text{div}; \Omega) := \{\mathbf{v} \in (L^2(\Omega))^2 / \operatorname{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 Mixed Finite Element Method of Raviart-Thomas with the lowest order (see [6]). The resolution of this system gives the velocity $\mathbf{u}$ which is needed to resolve our system (4.3).

The hydraulic conductivity depends on the adherent bacteria evolution. The evolution of the hydraulic conductivity is given in [9] by the following relation:

$$ K^n = K(c_{(1,n)}) = K^{n-1} \left(1 - \frac{c_{(1,n)}}{w_\infty}\right)^{nk}, \quad (5.33) $$
where $K^{n-1}$ is the conductivity at time $t_{n-1}$ and $nk$ is a given parameter. At the first time step we resolve the system (5.32) with $K0$ which gives the first velocity $\mathbf{u}^1(t = t_1)$ and then the system (4.3) to obtain $(c_1^{(1)}, c_2^{(1)}, c_3^{(1)}, c_4^{(1)})$. Then, at each time step $t_n$, for $n \ge 2$, we will use a strategy for the stopping test.

**Remark 3** When the biodegradation process is carried out, the concentration of adherent bacteria increases, leading to a decrease in conductivity over time. This results in a reduced flow rate, so that the final flow rate $\mathbf{u}$ is less than the initial flow rate $\mathbf{u}_0$ (i.e $\|\mathbf{u}\|_{\infty} < \mathbf{u}^0$.)

## 6. Numerical algorithm for solving the coupled problem

The biodegradation system (2.1) is a non-linear coupled system of four equations. Its coupling with the flow system (5.32) makes the biodegradation and flow problem a non-linear coupled problem. In this section, we present the two different algorithms that will be applied and compared in the following numerical section. Let us refer to $c_1^h, c_2^h, c_3^h$ and $c_4^h$ as the four unknown components of the Biodegradation system (2.1), and $(\mathbf{u}^h, \mathbf{p}^h)$ as the unknown of the flow equations (5.32). The first strategy, called Coupled Prediction Scheme (CPS) (see Algorithm 1), involves solving the biodegradation system iteratively, one component after another. In contrast, the second strategy, called Fixed-point method (FPM) (see Algorithm 2), solves the biodegradation system with strong coupling, resulting in a linear system with a large matrix.

**Algorithm 1** Iterative algorithm for solving the coupled problem using (CPS).

1: **input** Values of the parameters and properties
2: **load** mesh {$\Omega^h$}
3: **define** finite element variables: $c_1^h, c_2^h, c_3^h, c_4^h, q_1^h, q_2^h, q_3^h, q_4^h$; $\mathbf{u}^h, \mathbf{v}^h, \mathbf{p}^h, \phi^h$
4: $k \leftarrow 0$
5: $c_{i,in} \leftarrow c_i(0)$; $i = 1, \dots, 4$
6: $c_{i,in0} \leftarrow c_i(0)$; $i = 1, \dots, 4$
7: **for** $t \leftarrow 1$ **to** $t_{max}$ **do**
8: &nbsp;&nbsp;&nbsp;&nbsp;**for** $k \leftarrow 1$ **to** $k_{max}$ **do**
9: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$l_1 \leftarrow 0$; $l_2 \leftarrow 0$; $l_3 \leftarrow 0$; $l_4 \leftarrow 0$
10: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$c_1^{(k,l1)} \leftarrow c_{1,in}$; $c_2^{(k,l2)} \leftarrow c_{2,in}$; $c_3^{(k,l3)} \leftarrow c_{3,in}$; $c_4^{(k,l4)} \leftarrow c_{4,in}$
11: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;calculate $K^k$ using $c_1^{(k,l1)}$
12: **solve** (5.32) **return** $\mathbf{u}^h$ and $\mathbf{p}^h$
13: **for** $l_1 \leftarrow 1$ **to** $l_1_{max}$ **do**
14: &nbsp;&nbsp;&nbsp;&nbsp;**solve** ((5.2), for $i=1$) using $c_{1,in0}, c_1^{(k,l1-1)}, c_2^{(k,l2-1)}, c_3^{(k,l3-3)}, c_4^{(k,l4-1)}$ **return** $c_1^{(k,l1)}$
15: &nbsp;&nbsp;&nbsp;&nbsp;$e_{rr}^{(l1)} \leftarrow \|c_1^{(k,l1)} - c_1^{(k,l1-1)}\|_{\infty}$
16: **if** $e_{rr}^{(l1)} < \varepsilon_{l1}$ **then**
17: &nbsp;&nbsp;&nbsp;&nbsp;**break**
18: **end if**
19: $c_1^{(k,l1-1)} \leftarrow c_1^{(k,l1)}$
20: **end for**
21: **for** $l_2 \leftarrow 1$ **to** $l_2_{max}$ **do**
22: &nbsp;&nbsp;&nbsp;&nbsp;**solve** ((5.3), for $i=2$) using $\mathbf{u}^h, c_{2,in0}, c_1^{(k,l1-1)}, c_2^{(k,l2-1)}, c_3^{(k,l3-3)}, c_4^{(k,l4-1)}$ **return** $c_2^{(k,l2)}$
23: &nbsp;&nbsp;&nbsp;&nbsp;$e_{rr}^{(l2)} \leftarrow \|c_2^{(k,l2)} - c_2^{(k,l2-1)}\|_{\infty}$
24: **if** $e_{rr}^{(l2)} < \varepsilon_{l2}$ **then**
25: &nbsp;&nbsp;&nbsp;&nbsp;**break**
26: **end if**
27: $c_2^{(k,l2-1)} \leftarrow c_2^{(k,l2)}$
28: **end for**
29: **for** $l_3 \leftarrow 1$ **to** $l_3_{max}$ **do**
30: &nbsp;&nbsp;&nbsp;&nbsp;**solve** ((5.4), for $i=3$) using $\mathbf{u}^h, c_{3,in0}, c_1^{(k,l1-1)}, c_2^{(k,l2-1)}, c_3^{(k,l3-3)}, c_4^{(k,l4-1)}$ **return** $c_3^{(k,l3)}$
31: &nbsp;&nbsp;&nbsp;&nbsp;$e_{rr}^{(l3)} \leftarrow \|c_3^{(k,l3)} - c_3^{(k,l3-1)}\|_{\infty}$
32: **if** $e_{rr}^{(l3)} < \varepsilon_{l3}$ **then**

33: **break**

34: **end if**

35: $c_3^{(k,l3-1)} \leftarrow c_3^{(k,l3)}$

36: **end for**

37: **for** $l4 \leftarrow 1$ **to** $l4_{max}$ **do**

38: **solve** ((5.5), for i=4) using $\mathbf{u}^h$, $c_{4,in0}$, $c_1^{(k,l1-1)}$, $c_2^{(k,l2-1)}$, $c_3^{(k,l3-3)}$, $c_4^{(k,l4-1)}$ **return** $c_4^{(k,l4)}$

39: $e_{rr}^{(l4)} \leftarrow \|c_4^{(k,l4)} - c_4^{(k,l4-1)}\|_{\infty}$

40: **if** $e_{rr}^{(l4)} < \varepsilon_{l4}$ **then**

41: **break**

42: **end if**

43: $c_4^{(k,l4-1)} \leftarrow c_4^{(k,l4)}$

44: **end for**

45: $\hat{l1}(k) \leftarrow l1$; $\hat{l2}(k) \leftarrow l2$; $\hat{l3}(k) \leftarrow l3$; $\hat{l4}(k) \leftarrow l4$

46: $c_1^k \leftarrow c_1^{(k,\hat{l1})}$; $c_2^k \leftarrow c_2^{(k,\hat{l2})}$; $c_3^k \leftarrow c_3^{(k,\hat{l3})}$; $c_4^k \leftarrow c_4^{(k,\hat{l4})}$

47: $E^{(k)} \leftarrow \|c_1^{(k)} - c_1^{(k-1)}\|_{\infty} + \|c_2^{(k)} - c_2^{(k-1)}\|_{\infty} + \|c_3^{(k)} - c_3^{(k-1)}\|_{\infty} + \|c_4^{(k)} - c_4^{(k-1)}\|_{\infty}$

48: **if** $E^{(k)} < \varepsilon_k$ **then**

49: **break**

50: **end if**

51: $c_{i,in} \leftarrow c_i^{(k)}$; $i = 1, ..., 4$

52: **end for**

53: $c_{i,in0} \leftarrow c_{i,in}$; $i = 1, ..., 4$

54: **end for**

**Algorithm 2** Iterative algorithm for solving the coupled problem using (FPM).

1: **input** Values of the parameters and properties

2: **load mesh** {$\Omega^h$}

3: **define** finite element variables: $c_1^h, c_2^h, c_3^h, c_4^h, q_1^h, q_2^h, q_3^h, q_4^h, \mathbf{u}^h, \mathbf{v}^h, p^h, \phi^h$

4: $\mathbf{C} \leftarrow (c_1^h, c_2^h, c_3^h, c_4^h)$

5: $k \leftarrow 0$

6: $\mathbf{C}_{in} \leftarrow \mathbf{C}(0)$

7: $\mathbf{C}_0 \leftarrow \mathbf{C}(0)$

8: **for** $t \leftarrow 1$ **to** $t_{max}$ **do**

9: **for** $k \leftarrow 1$ **to** $k_{max}$ **do**

10: $l \leftarrow 0$;

11: $\mathbf{C}^{(k,l)} \leftarrow \mathbf{C}_0$

12: calculate $K^k$ using $\mathbf{C}^{(k,l)}$

13: **solve** (5.32) **return** $\mathbf{u}^h$ and $p^h$

14: **for** $l \leftarrow 1$ **to** $l_{max}$ **do**

15: **solve** (5.1) using $\mathbf{u}^h$, $\mathbf{C}_{in}$, $\mathbf{C}^{(k,l)}$ **return** $\mathbf{C}^{(k,l+1)}$

16: $e_{rr}^{(l)} \leftarrow \|\mathbf{C}^{(k,l+1)} - \mathbf{C}^{(k,l)}\|_{\infty}$

17: **if** $e_{rr}^{(l)} < \varepsilon_l$ **then**

18: **break**

19: **end if**

20: $\mathbf{C}^{(k,l)} \leftarrow \mathbf{C}^{(k,l+1)}$

21: **end for**

22: $\hat{l}(k) \leftarrow l$

23: $\mathbf{C}^k \leftarrow \mathbf{C}^{(k+1,\hat{l})}$

24: $E^{(k)} \leftarrow \|\mathbf{C}^{(k+1)} - \mathbf{C}^{(k)}\|_{\infty}$

25: **if** $E^{(k)} < \varepsilon_k$ **then**

26: **break**

27: **end if**
28: $C_0 \leftarrow C^{(k)}$

29: end for

30: $C_{in0} \leftarrow C_{in}$

31: end for

## 7. Numerical results

In this section, we present and discuss the numerical results of the presented nonlinear coupled model of the biodenitrification process. In all the following experiments, we take $\Omega \subset \mathbb{R}^2$, $\partial\Omega = \Gamma_1 \cup \Gamma_2 \cup \Gamma_3$ such that: $\Omega =]0, L_1[\times]0, L_2[$, $L_1 = 3$, $L_2 = 1$, where $\Gamma_1 = \{(x, y) \in \mathbb{R}^2$, with $x = 0$; $0 \le y \le 0, 3L_2\}$, $\Gamma_2 = \{(x, y) \in \mathbb{R}^2$, with $x = L_1$; $0, 7L_2 \le y \le L_2\}$ and $\Gamma_3 = \partial\Omega - \{\Gamma_1 \cup \Gamma_2\}$ (see Figure 1). The initial conditions are provided in [8] and the parameter values are taken from [1,2,7]. All these data are summarized in Table 1 and Table 2.

<table><caption>Table 1: Parameters values used in the simulations.</caption><thead><tr><th>Parameters</th><th>Values</th><th>Parameters</th><th>Values</th></tr></thead><tbody><tr><td>$D_1$</td><td>0.8cm<sup>2</sup>/h</td><td>$\mu_{max}$</td><td>0.4h<sup>-1</sup></td></tr><tr><td>$D_2$</td><td>0.8cm<sup>2</sup>/h</td><td>$Y_p$</td><td>0.8</td></tr><tr><td>$D_3$</td><td>0.7cm<sup>2</sup>/h</td><td>$Y_a$</td><td>0.8</td></tr><tr><td>$D_4$</td><td>0.7cm<sup>2</sup>/h</td><td>$\beta$</td><td>0.2h<sup>-1</sup></td></tr><tr><td>$\alpha$</td><td>0.02h<sup>-1</sup></td><td>$T$</td><td>20 days</td></tr><tr><td>$w_{\infty}$</td><td>0.05g/cm<sup>2</sup></td><td>$R$</td><td>0.9</td></tr><tr><td>$\mathbf{u}^0$</td><td>0, 001</td><td>$p_D$</td><td>0</td></tr><tr><td>$a$</td><td>0, 1</td><td>$K_{S_1}^a$</td><td>54 mg/l</td></tr><tr><td>$k_p$</td><td>0.005h<sup>-1</sup></td><td>$K_{S_1}^p$</td><td>54 mg carbon/l</td></tr><tr><td>$k_a$</td><td>0.005h<sup>-1</sup></td><td>$K_{S_2}^p$</td><td>50mgNO<sub>3</sub><sup>-</sup>/l</td></tr><tr><td>$\mu_{max}^p$</td><td>0.4h<sup>-1</sup></td><td>$K_{S_2}^a$</td><td>50 mg carbon/l</td></tr><tr><td>$\gamma$</td><td>10cm<sup>3</sup>/cm<sup>2</sup></td><td></td><td></td></tr></tbody></table>

<table><caption>Table 2: Initial conditions.</caption><thead><tr><th>$X_1^0$</th><th>$X_2^0$</th><th>$S_1^{in}$</th><th>$S_2^{in}$</th><th>$K_0$</th><th>$S_1^0$</th><th>$S_2^0$</th></tr></thead><tbody><tr><td>10mg/cm<sup>2</sup></td><td>10 mg/l</td><td>104 mg/l</td><td>100 mg/l</td><td>0.865 cm/h</td><td>104 mg/l</td><td>100 mg/l</td></tr></tbody></table>

<table><caption>Table 3: Values.</caption><tbody><tr><td>$A_1$</td><td>0.232</td><td>$X_1$</td><td>0.8</td></tr><tr><td>$A_2$</td><td>0.1192</td><td>$X_2$</td><td>0.795</td></tr><tr><td>$A_3$</td><td>0.08945</td><td>$X_3$</td><td>0.695</td></tr><tr><td>$A_4$</td><td>0.0805</td><td>$X_4$</td><td>0.695</td></tr><tr><td>$\tau_n$</td><td>0.67</td><td>$Q_3$</td><td>0.11525</td></tr><tr><td>$Q_4$</td><td>0.13725</td><td></td><td></td></tr></tbody></table>

Using these data and equations (5.9), (5.10), (5.17), (5.15), (5.16), (5.19), (5.20) and (5.23), we have Table 3.

Since $(\tau_n)^{-1} > 1.492537$, then $(\tau_n)^{-1} + k_a + \beta - \mu_{max}^a > 0.995 > D_1$, which implies that $X_1 = D_1$, leading to $A_1 < 0.38$. With the data given in Table 3, we have $\kappa = 0.29$ and $\kappa^1 = 0.6844$, which implies that $\kappa < \kappa^1$. Therefore, the (CPS) scheme converges faster than the (FPM) scheme.

Our main aim in this section is to compare the two algorithms we presented and analyzed earlier. To do this, we decided to deal with different situations, depending on the step size of the spatial discretization $h$ and the step size of the time discretization $\tau$ (see Table 4).
<table><caption>Table 4: Values of the mesh step $h_i$ and the time step $\tau_i$.</caption><thead><tr><th>h<sub>1</sub></th><th>h<sub>2</sub></th><th>h<sub>3</sub></th><th>τ<sub>1</sub></th><th>τ<sub>1</sub></th><th>τ<sub>3</sub></th></tr></thead><tbody><tr><td>1/15</td><td>1/25</td><td>1/35</td><td>0.05</td><td>0.075</td><td>0.1</td></tr></tbody></table>

## COMPUTATION TIMES

Figure 2 shows three different quasi-uniform meshes of the domain. Sub-figure 2-(a) with mesh step $h_1 = 1/15$ generating 1444 triangles and 783 vertices, Sub-figure 2-(b) with mesh step $h_2 = 1/25$ generating 4036 triangles and 2119 vertices, and Sub-figure 2-(c) with mesh step $h_3 = 1/35$ generating 8096 triangles and 4189 vertices. In addition, we consider three time discretization steps, $\tau_1 = 0.05$, $\tau_2 = 0.075$, and $\tau_1 = 0.1$ (see Table 4).

<figure><img src="image_4.png" /><figcaption>Figure 2: Three quasi-uniform meshes of the domain.</figcaption></figure>

First, by setting the time step to, for example, the smallest value $\tau = \tau_1$, Figure 3 shows the evolution of the computational cost of each algorithm as a function of mesh size $h$, over the interval $[t_0, T_f]$. We can see that for the initial mesh ($h = h_1$), Algorithm 1 executes 333 s versus 536 s for Algorithm 2. For the highly refined mesh ($h = h_3$), Algorithm 1 executes 1243 s versus 4385 s for Algorithm 2. As can be seen in the figure, the red computation cost curve for Algorithm 2 is superior to the blue computation cost curve for Algorithm 1. Consequently, in all these situations, we deduce that Algorithm 1 is faster than Algorithm 2.

In the second experimental test, by fixing the mesh step, taking, for example, the highly refined mesh $h = h_3$, Figure 4 shows the evolution of the computational cost of each algorithm, during the interval $[t_0, T_f]$, as a function of time discretization steps. We observe that for the case of time step $\tau = \tau_1$, Algorithm 1 executes 1243 s versus 4385 s for Algorithm 2. For the case of time step $\tau = \tau_3$, Algorithm 1 executes 728 s versus 1658 s for Algorithm 2. Similar to what we obtained in the first test, we can see in Figure 4 that the red curve of Algorithm 2 is superior to the blue curve of Algorithm 1. Consequently, in all these situations we also deduce that Algorithm 1 is faster than Algorithm 2.

## SPATIAL DISTRIBUTION
<figure><img src="image_5.png" /><figcaption>Figure 3: Evolution of the computational cost of the two algorithms as a function of mesh size, over the interval $[t_0, T_f]$ with a fixed time step $\tau=\tau_1$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 4: Evolution of the computational cost of the two algorithms, during the interval $[t_0, T_f]$, as a function of time steps, using a highly refined mesh $h = h_3$.</figcaption></figure>

<figure><img src="image_7.png" /></figure>
<figure><img src="image_8.png" /><figcaption>Figure 6: Spatial distribution of $x_1, x_2, S_1$ and $S_2$ involved in the coupled biodenitrification problem at time $T_f$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7: Spatial distribution of Pressure and Velocity in the coupled biodenitrification problem at time $T_f$.</figcaption></figure>

Applying Algorithm 1 using a uniform descritization step in time $\tau = \tau_1$ together with a quasi-uniform mesh ($h = h_1$), Figures 6 and 7 shows the numerical simulations of the spatial distribution of all the quantities of nonlinear coupled model of the biodenitrification process at $T_f$. The first two subfigures represent the concentrations of the two types of bacteria: adherent bacteria ($x_1$) and planktonic (mobile) bacteria ($x_2$), respectively. The next two subfigures show the concentrations of the nutrient ($S_1$) and the contaminant (nitrate, $S_2$), respectively. The final two subfigures display the components of the flow system: pressure and velocity field. The values of the bacterial and substrate quantities at the initial time $t_0$ are uniform at all points of the spatial domain. Whereas at time $T_f$, we observe that the concentrations of bacteria and substrate quantities are low at the flow outlet and high at the flow inlet. Such behavior is already obtained in the reference [4].

## 8. Conclusion and perspectives

This article presents a comparison between two algorithms for linearizing a mathematical system: the first is the CPS method, and the second is the FPM method. The system describes the biodegradation process in porous media by bacteria, which are divided into planktonic and adherent types, with the chemotaxis effect taken into account. The mathematical model consists of a system of five strongly coupled parabolic equations, incorporating a nonlinear advection term. To establish the existence and uniqueness of the solution, we employed the fixed-point theorem.

For the solution approximation, we utilized the finite element (FE) method. To linearize the system at each time step, we introduced an iterative scheme referred to as the Coupled Prediction Scheme (CPS) and demonstrated its convergence. This approach provides a faster and more efficient numerical algorithm. The backward Euler scheme was used for time discretization. Several numerical tests in the 2D case, commonly used in the literature for simplified models, validate the effectiveness of the Coupled Prediction Scheme. The numerical results indicate that the CPS scheme converges faster than the Fixed-Point Method (FPM).

In future work, we plan to develop an algorithm that accounts for the dependence of the last three variables on the fluid velocity, enabling the simultaneous solution of the two governing equations.

## References

[1] Abaali M., and Mghazli Z., Mathematical modelling of biodegradation in situ application to biodenitrification. Computers & Mathematics with Applications 79.6 : 1833-1844(2020)
[2] Abaali M., Harmand J., and Mghazli Z., Impact of Dual Substrate Limitation on Biodenitrification Modeling in Porous Media. Processes 8.8, 890(2020)

[3] Abaali M., Existence and stability of a steady state for a coupled system of elliptic equation modelling the biodenitri-fication process. Math Meth Appl Sci., 1-18(2022) https://doi:10.1002/mma.8235

[4] Abaali M., and Ouchtout S., Mathematical modeling of the chemotaxis-biodenitrification process. Computers & Math-ematics with Applications 169, 56-67(2024)

[5] Bothe D., et al., Global existence for diffusion-electromigration systems in space dimension three and higher. Nonlinear Analysis: Theory, Methods & Applications 99, 152-166(2014)

[6] Brezzi F., and Fortin M., Mixed and hybrid finite element methods. Springer Science & Business Media, Vol. 15. (2012)

[7] Chen-Charpentier B.M., Kojouharov H.V., Mathematical modeling of bioremediation of trichloroethylene in aquifers. Comput. Math. Appl. 56 (3), 645-656(2008)

[8] 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. Diss. Université de Lille 1, (1996)

[9] 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 34 (5), 934-942 (1996)

[10] Deteix J., Jendoubi A., and Yakoubi D., A Coupled Prediction Scheme for Solving the Navier-Stokes and Convection-Diffusion Equations. SIAM Journal on Numerical Analysis 52.5, 2415-2439 (2014)

[11] Freter, R., Brickner, H., Fekete, J., Vickerman, M.M., Carey, K.E, Survival and implantation of Escherichia coli in the intestinal tract. Infect. Immun. 39, 686-703(1983)

[12] Freter R., Brickner H. and Temme S., An understanding of colonization resistance of the mammalian large intestine requires mathematical analysis. Microecology and Therapy 16, 147-155(1986)

[13] Keller, Evelyn F., and Lee A. Segel., Initiation of slime mold aggregation viewed as an instability. Journal of theoretical biology 26.3, 399-415(1970)

[14] Nithiarasu P., A unified fractional step method for compressible and incompressible flows, heat transfer and incom-pressible solid mechanics. Internat. J. Numer. Methods Heat Fluid Flow, 18, 111-130(2008)

[15] Ouertatani, N., Cheikh, N. B., Beya, B. B., & Lili, T., Numerical simulation of two-dimensional Rayleigh-Bénard convection in an enclosure. Comptes Rendus Mécanique, 336(5), 464-470(2008)

[16] Rebollo, T. C., Del Pino, S., & Yakoubi, D., An iterative procedure to solve a coupled two-fluids turbulence model. ESAIM: Mathematical Modelling and Numerical Analysis, 44(4), 693-713(2010)

[17] Shishkina O., Shishkin A., and Wagner C., Simulation of turbulent thermal convection in complicated domains. J. Comput. Appl. Math., 226, 336-344(2009)

[18] Yakoubi, D., Analyse et mise en œuvre de nouveaux algorithmes en méthodes spectrales. (Doctoral dissertation, Université Pierre et Marie Curie-Paris VI) (2007)

Mostafa Abaali,
Engineering Sciences Laboratory, ESL ENSA,
Ibn Tofail University, Kenitra
Morocco.
ORCID: https://orcid.org/0000-0001-5225-465X.
E-mail address: mostafaabaali@gmail.com

and

Salih Ouchtout,
IFPEN, 1 et 4 avenue de Bois-Preau, 92852 Rueil-Malmaison,
France
ORCID: https://orcid.org/0000-0003-0792-4081.
E-mail address: salih.ouchtout@gmail.com