Directional migration and competition in fluid media: Global existence, uniqueness, and robust simulation of chemotaxis and tumor growth dynamics

Georges Chamoun

doi:10.2478/ijmce-2026-0018

Introduction

This paper is dedicated to the examination of the interaction between two biological species, whether they are in competition or collaboration, within a fluid environment. This study lies in analyzing their dynamics in the presence of anisotropic heterogeneous nonlinear degenerate diffusion and convection tensors. These species can ensure food, nutrients, and oxygen levels, or they can collaborate to avoid predators and threats. The evolution of these densities over time is influenced by various scenarios: Lotka-Volterra competition, anisotropic degenerate diffusion, nonlinear chemotactic movement, and fluid dynamics. We denote the densities of species 1 (resp. species 2) by M(x, t) (resp. W(x, t)). Since the main biological applications of this paper relate to interspecies competition and fluid interactions in real biological regions, we consider W as an open bounded spatial domain in ℝ^d; d ≤ 3 corresponding to 1D fibers as nerve channels, 2D tissues as membranes, or 3D organs as tumor microenvironments.

Over the past decade, there has been a surge in the exploration of biological mechanisms. However, in this field, mathematical and numerical analysis have advanced to the point where they often surpass real experiments in predictive power. Indeed, it helps to explore all the consequences of manipulating the various parameters associated with any particular scenario of Figure 1. Before delving into our full model, which gathers all possible constraints, we now navigate the relevant literature related to each scenario.

To take competition and spatial linear diffusion of species, mathematicians and ecologists studied first the following system of reaction-diffusion equations: 1 ${\begin{array}{r} \partial_{t} M - {\tilde{D}}_{1} Δ M = μ_{1} M (1 - M - α_{1} W), in Ω \times (0, T] \\ \partial_{t} W - {\tilde{D}}_{2} Δ M = μ_{2} W (1 - W - α_{2} M), in Ω \times (0, T] \\ \nabla M \cdot η = 0, \nabla W \cdot η = 0, on \partial Ω \times (0, T], \end{array}$ where ${\tilde{D}}_{1}$ , ${\tilde{D}}_{2}$ , μ₁, μ₂, α₁, α₂ are strictly positive constants and η is the outward-pointing unit normal to the smooth boundary ∂,Ω. If ${\tilde{D}}_{1} = {\tilde{D}}_{2} = 0$ , then the system (1) is reduced to the well-known predator-prey systems, initially introduced by Lotka in [1] and further developed by Volterra in [2], as a set of two ordinary differential equations: 2 ${\begin{array}{l} \dot{M} & = & μ_{1} M (1 - M - α_{1} W) \\ \dot{W} & = & μ_{2} W (1 - W - α_{2} M) . \end{array}$

System (2) has four critical points O(0, 0), N₁(1, 0), N₂(0, 1) and $N_{3} (\frac{1 - α_{1}}{1 - α_{1} α_{2}}, \frac{1 - α_{2}}{1 - α_{1} α_{2}})$ . According to the Jacobian matrix at each critical point, one may determine whether it is a stable node with real negative eigenvalues or an unstable saddle point with real eigenvalues of opposite signs. Solutions to (2) are convergent to the stable node N₃ in the weak competition regime when α₁, α₂ ∈ (0, 1) and towards the stable nodes N₁ or N₂ in the strong competition regime when 0 < α₁ < 1 < α₂ or 0 < α₂ < 1 < α₁. Next, for the system (1), the authors in [3] show that N₃ remains the global attractor in the weak competition regime. Moreover, [4,5] prove that, besides unstable N₃ and stable N₁ or N₂ in the strong regime, system can have stable inhomogeneous steady states if ${\tilde{D}}_{1}$ and ${\tilde{D}}_{2}$ are not too large.

Various population dynamics can arise based on the species’ capacity to orient their movement in response to chemical gradients. This attraction or repulsion is called chemotactic motion. For a single species, Keller and Segel introduced in [6] the most classical model, and there was a vast literature in [7, 8]. In the context of the Lotka-Volterra system (2), chemotactic effects were integrated in [9, 10] to simulate behaviors such as prey evasion or predator pursuit. Extensive research has also been conducted on coupling chemotaxis effects with SIR-type equations to investigate the spatial spread of pandemics [11, 12]. Recent advances in mathematical oncology have also emphasized the role of memory effects and anomalous diffusion in tumor-immune interactions. In this context, fractional-order models have been proposed to better capture nonlocal biological responses and treatment dynamics. For instance, authors in [13, 14] investigated an immune system-cancer interaction model involving Mittag-Leffler kernels, demonstrating the relevance of non-classical diffusion operators in immunotherapy modeling. Furthermore, uniqueness properties for a class of fractal-fractional differential systems were analyzed in [15], providing new theoretical insights into the well-posedness of nonlocal biological models. Chemotaxis scenarios were also introduced to model host-parasite interaction in [16, 17]. Further inquiry into the persistence of competitive exclusion, incorporating chemotaxis and competitive interactions between two species, can be found in the relevant literature [18–21].

Traditionally, numerous studies mentioned above focused on isotropic homogeneous species dynamics within a stationary fluid environment. However, recent works highlighted in [22] underscore the significant influence of the surrounding fluid dynamics on the species’ behavior. Indeed, transportation effects change organism behavior by affecting capture rates and prey searching strategies. In the context of a one-species case, a wide range of literature exists on chemotaxis within the Navier-Stokes equation, both without or with logistic terms [23,24]. On the other hand, the introduction of full anisotropic diffusion tensors is biologically motivated since tumor tissues and extracellular matrices often facilitate migration along specific directions while limiting movement in others. Unfortunately, classical numerical methods cannot help to construct a convergent numerical scheme to discretize models including anisotropic diffusive tensors without any restriction on meshes. In [25], the realistic representation of spatial anisotropy and heterogeneous models was studied without logistic sources and for one biological organism. In this paper, a system comprising two swimming species with degenerate diffusive properties, competitive kinetics, and chemotaxis responses react to significant directional variations and specific inhomogeneous spatial orientations through the following system is studied.3 ${\begin{matrix} \partial_{t} M - \nabla \cdot (Q (x) a (M) \nabla M) + \nabla \cdot (Q (x) χ_{1} (M) \nabla C) + U \cdot \nabla M = μ_{1} M (1 - M - Φ_{1} (W)), \\ \partial_{t} W - \nabla \cdot (R (x) b (W) \nabla W) + \nabla \cdot (R (x) χ_{2} (W) \nabla C) + U \cdot \nabla W = μ_{2} W (1 - W - Φ_{2} (M)), \\ \partial_{t} C - \nabla \cdot (S (x) \nabla C) + U \cdot \nabla C = - (α M + β W) C, \\ \partial_{t} U - v Δ U + κ (U \cdot \nabla) U + \nabla P = - (γ M + λ W) \nabla ϕ, \\ \nabla \cdot U = 0, in Ω \times (0, T], \\ Q (x) a (M) \nabla M \cdot η = 0, R (x) b (W) \nabla W \cdot η = 0, S (x) \nabla C \cdot η = 0, U = 0, on \partial Ω \times (0, T], \\ M (x, 0) = M_{0} (x), W (x, 0) = W_{0} (x), C (x, 0) = C_{0} (x), U (x, 0) = U_{0} (x), in Ω . \end{matrix}$

The details, regarding all the variables of the system (3), are provided in the Table 1.

Table 1

Model variables and parameters associated with System (3).

Variables	Description
M, W	The density of species 1 and 2
C	The concentration of the chemical
U, P	The velocity and the pressure of the fluid
Q, R and S	The anisotropic diffusive tensors
a(M), b(W)	The density-dependent diffusion coefficients
χ₁(M), χ₂(W)	The chemotactic sensitivity functions
μ₁, μ₂ > 0	The growth rate of populations 1 and 2
α₁, α₂ > 0	The strength of populations 1 and 2 in competition
α, β > 0	The consumption rate of chemicals by populations 1 and 2
$\vec{g} = - (γ M + λ W) \nabla ϕ$	The external gravitational force exerted by the species 1 and 2 on the fluid

The functional responses Φ₁ and Φ₂ were originally proposed by Holling [26, 27] based on the time a predator spends between two activities, "prey searching" and "prey handling". Although there are three types of these functions, this paper focuses on the classical predator-prey dynamical system, which exhibits linear dependence on the density of other species [28]. Consequently, the relationship is expressed as: Φ₁(W) = α₁W and Φ₂( M) = α₂M.

From a biological perspective, the coefficients μ₁ and μ₂ represent intrinsic proliferation rates of the competing populations, while the parameters α₁ and α₂ quantify the strength of interspecies competition through resource limitation effects. The anisotropic tensors Q(x) and R(x) encode directional migration patterns induced by heterogeneous tissue structures such as aligned fibers in extracellular matrices. The chemotactic sensitivities χ₁ and χ₂ determine the ability of each species to respond to chemical gradients, whereas the consumption parameters α and β model the depletion of nutrients or oxygen by the interacting populations. Finally, the coupling coefficients γ and λ characterizes the feedback influence of the population densities on the surrounding fluid motion.

Moreover, the homogeneous Neumann boundary conditions reflect the assumption of an isolated biological environment, where no external influx or efflux of cells or signaling molecules occurs across the boundary. This setting is consistent with confined tumor regions or ecological habitats surrounded by impermeable membranes. The no-slip Dirichlet condition imposed on the fluid velocity models, the interaction between the surrounding medium and a rigid extracellular matrix or tissue boundary, ensures that the fluid remains stationary at the domain interface.

The diffusion equations exhibit a "volume-filling" effect, which ensures biological admissibility. This effect is achieved through the diffusion coefficients a(M) and b(W), which approach zero as the densities M and W are close to a threshold or when there is no cell population present. Upon normalization, the threshold density is set to 1, and a common example of a diffusion-density coefficient is a(M) = M(1 – M) for M ∈ [0, 1]. The sign of the chemotactic sensitivity, denoted as χ₁ or χ₂, signifies whether the species is attracted to or repelled from the chemical signals.

Furthermore, the dynamics of the fluid are governed by the Navier-Stokes equations, or by the Stokes equations when the parameter κ equals zero. The gravitational force of the species on the fluid is determined by the product of the species volume V_b, the gravitational acceleration g (equal to 9,8m/s²), and the density difference between the living cells and the fluid. This force acts in the upward direction through a fixed vertical gradient field ∇ϕ such that $\nabla ϕ = ” V_{b} (ρ_{b} - ρ) g ” \vec{z}$ along the upwards unit vector $\vec{z}$ . The total force $\vec{g}$ remains dependent on the evolving densities M and W.

For Chemotaxis-Navier-Stokes systems with linear non-degenerate isotropic diffusion, various criteria for global well-posedness and blow-up were derived in [29]. Additionally, the authors in [30] found the same asymptotic behavior of system (2) in the weak and strong competition regimes, where the chemical concentration and the fluid velocity converge exponentially to zero in all cases. Recently, the global existence of two-species chemotaxis-Stokes equations has been demonstrated under sufficient conditions [31], with subsequent studies establishing the asymptotic dynamics and methods for preventing blow-up in both two-dimensional 2D [32] and three-dimensional 3D [33] cases. Secondly, this work lies in demonstrating the global well-posedness and the uniqueness of weak solutions of system (3) over time. Indeed, the presence of singular degenerate diffusive terms may exhibit limited regularity properties and difficulties in establishing the uniform bounds of solutions While global existence, boundedness and asymptotic behavior results have been obtained for isotropic homogeneous models, the biological realism of these models remain limited in heterogeneous tissues where cell migration is strongly direction-dependent. Moreover, the model accounts for competitive interactions between two chemotactic species evolving in a moving fluid environment. This coupling introduces additional analytical challenges due to the loss of uniform parabolicity and the strong nonlinear feedback between species densities and fluid transport, which are not captured by classical isotropic chemotaxis-Navier-Stokes formulations.

To our knowledge, there are no large numerical studies for the competitive two-species chemotaxis-Navier-Stokes model in fluids. Classical methods, including finite volume techniques, impose restrictions on meshes while finite element methods often encounter instabilities, particularly in cases dominated by convection. Although a finite volume method has been proposed for species in a fluid at rest [34], it does not apply to anisotropic tensors on general meshes. Similarly, approaches using the classical finite element method discretization lead to many instabilities in some dominated situations. Indeed, such standard methods may exhibit unbiological negative densities and numerical blow-up solutions in general meshes due to the lack of discrete maximum principles and to the absence of monotonicity-preserving properties. Hence, this article focus on constructing an efficient and robust convergent numerical method to discretize system (3) on general meshes. This method must satisfy the discrete maximum principle, ensuring the confinement of discrete solutions, which is a mathematical translation of the volume-filling effect. The development of the numerical algorithm to consider two species instead of one and to incorporate competitive kinetics instead of zero logistic source terms has been executed through a Fortran code.

In this paper, Section 2 outlines the main assumptions regarding parameters and presents the main results of this work. Sections 3 and 4 are devoted to the proofs of the main theoretical results of global existence and uniqueness of weak solutions. Moreover, a recent numerical scheme originally designed for single-species trajectories without logistic sources is expanded to provide an efficient computational framework and an accurate predictive code for the two-species competitive model in a fluid in Section 5. Then, the Section 6 conducts the numerical simulations to enhance the understanding of the swimming competing species dynamics, in the presence of anisotropic and heterogeneous diffusive tensors. A conclusion is finally given in Section 7.

Set of the problem and main results

This section is devoted to detail the formal problem, to define the key parameters and to present the main objectives and hypotheses. Initially, the main assumptions are expressed as,

The experimental prevention of overcrowding is represented through degenerate nonlinear diffusive and convective functions: 4 $a, b, χ_{1}, χ_{2} belong to the set {w \in C ([0, 1], ℝ^{+}) such that w (0) = w (1) = 0} .$
The potential: 5 $ϕ = ϕ (x); \nabla ϕ \in {(L^{\infty} (Ω))}^{d} .$
∀x ∈ Ω, the matrices Q, R and S in ℳ_d(ℝ) are symmetric and verify 6 $Q_{i, j} \in L^{\infty} (Ω), R_{i, j} \in L^{\infty} (Ω), S_{i, j} \in L^{\infty} (Ω), \forall i, j \in {1, \dots, d} .$
There exist $c_{Q} \in ℝ_{+}^{*}, c_{R} \in ℝ_{+}^{*}$ and $c_{S} \in ℝ_{+}^{*}$ such that a.e. x ∈ Ω, ∀ξ ∈ ℝ^d, 7 $Q (x) ξ \cdot ξ \geq c_{Q} | ξ |^{2}, R (x) ξ \cdot ξ \geq c_{R} | ξ |^{2}, S (x) ξ \cdot ξ \geq c_{S} | ξ |^{2} .$
$\exists \bar{D}$ and $\bar{D_{1}} \in ℝ_{+}^{*}$ such that a.e. x ∈ Ω, ∀M ∈ [0, 1], 8 $‖ D (x, M) ‖_{ℳ_{d} (ℝ)} = ‖ Q (x) a (M) ‖_{ℳ_{d} (ℝ)} \leq \bar{D} and {‖ D_{1} (x, M) ‖}_{ℳ_{d} (ℝ)} = {‖ Q (x) χ_{1} (M) ‖}_{ℳ_{d} (ℝ)} \leq \bar{D_{1}} .$
$\exists \bar{E}$ and ${\bar{E}}_{1} \in ℝ_{+}^{*}$ such that a.e. x ∈ Ω, ∀W ∈ [0, 1], 9 $‖ E (x, W) ‖_{ℳ_{d} (ℝ)} = ‖ R (x) b (W) ‖_{ℳ_{d} (ℝ)} \leq \bar{E} and {‖ E_{1} (x, W) ‖}_{ℳ_{d} (ℝ)} = {‖ R (x) χ_{2} (W) ‖}_{ℳ_{d} (ℝ)} \leq \bar{E_{1}} .$
10 $0 \leq M_{0} \leq 1, 0 \leq W_{0} \leq 1, C_{0} \geq 0 a .e . in Ω and C_{0} \in L^{\infty} (Ω) .$
11 $U_{0} \in H and g \in L^{2} (0, T; V^{'}) such that H = {\bar{{U \in D (Ω), \nabla \cdot U = 0}}}^{L^{2} (Ω)} and V = {\bar{{U \in D (Ω), \nabla \cdot U = 0}}}^{H_{0}^{1} (Ω)} .$
12 $ℬ (U, V_{1}, V_{2}) = \int_{Ω} (U \cdot \nabla V_{1}) V_{2} d x is continuous over H_{0}^{1} (Ω) \times H^{1} (Ω) \times H^{1} (Ω) and ℬ (U, V, V) = 0 if \nabla \cdot U = 0.$
13 $M_{0}, W_{0} \in L^{\infty} (Ω), U_{0} \in W^{2 - \frac{2}{p}, p} (Ω), C_{0} \in W^{2 - \frac{2}{p}, p} (Ω); p > d and d \leq 3.$ 14 $\nabla ϕ \in W^{1, \infty} (Ω), χ_{1}, χ_{2} of class C^{1} and the coefficients Q_{i, j}, R_{i, j}, S_{i, j} \in C^{1} (\bar{Ω}) .$
15 $\exists κ_{0} > 0; {(χ_{1} (M_{1}) - χ_{1} (M_{2}))}^{2} \leq κ_{0} (M_{1} - M_{2}) (A (M_{1}) - A (M_{2})), \forall M_{1}, M_{2} \in [0, 1] .$
16 $\exists κ_{1} > 0; {(χ_{2} (W_{1}) - χ_{2} (W_{2}))}^{2} \leq κ_{1} (W_{1} - W_{2}) (B (W_{1}) - B (W_{2})), \forall W_{1}, W_{2} \in [0, 1] .$

Remark 1.

Through these assumptions, we ensure the boundedness of the tensors in (6) that can be homogeneous isotropic or heterogeneous anisotropic. Then, the coercivity and the continuity of the diffusive and convective operators for the proof of global existence are guaranteed in (7), (8), and (9). Moreover, the boundedness of the initial data is given in (10), and the initial conditions for the Navier-Stokes equations are given in (11) and (12). Further assumptions (13) till (16) are introduced for the uniqueness proof section. Finally, note that the model (3) can be reduced to the standard Keller-Segel model with linear diffusion and chemotaxis sensitivity. It can also be generalized to porous medim where a(M) = M^α(1 – M)^γ for α, γ ≥ 0.

We recall that C_c(0, T; V) denotes the space of V-valued continuous functions with compact support in the time interval (0, T), endowed with the uniform topology. Moreover, C_w(0, T;L²(Ω)) denotes the space of functions that are weakly continuous from [0, T] into L²(Ω).

Definition 1.

A quadruple (M, W, C, U) is called as weak solution of (3) if $\begin{array}{r} 0 \leq M (x, t) \leq 1, 0 \leq W (x, t) \leq 1, C (x, t) \geq 0 a .e . in Q_{T} = Ω \times [0, T], \\ M, W \in C_{w} (0, T; L^{2} (Ω)), \partial_{t} M, \partial_{t} W \in L^{2} (0, T; {(H^{1} (Ω))}^{'}), \\ A (M) : = \int_{0}^{M} a (r) d r \in L^{2} (0, T; H^{1} (Ω)), B (W) : = \int_{0}^{W} b (r) d r \in L^{2} (0, T; H^{1} (Ω)), \\ C \in L^{\infty} (Q_{T}) \cap L^{2} (0, T; H^{1} (Ω)) \cap C (0, T; L^{2} (Ω)); \partial_{t} C \in L^{2} (0, T; {(H^{1} (Ω))}^{'}), \\ U \in L^{\infty} (0, T; H) \cap L^{2} (0, T; V) \cap C_{w} (0, T; H); \frac{d U}{d t} \in L^{1} (0, T; V^{'}), \end{array}$ and (M, W, C, U) satisfy 17 $\begin{matrix} \int_{0}^{T} < \partial_{t} M, ψ_{1} >_{{(H^{1})}^{'}, H^{1}} d t + \iint_{Q_{T}} [Q (x) a (M) \nabla M - Q (x) χ_{1} (M) \nabla C - M U] \cdot \nabla ψ_{1} d x d t \\ = \iint_{Q_{T}} μ_{1} M (1 - M - α_{1} W) ψ_{1} d x d t, \end{matrix}$ 18 $\begin{matrix} \int_{0}^{T} < \partial_{t} W, ψ_{2} >_{{(H^{1})}^{'}, H^{1}} d t + \iint_{Q_{T}} [R (x) b (W) \nabla W - R (x) χ_{2} (W) \nabla C - W U] \cdot \nabla ψ_{2} d x d t \\ = \iint_{Q_{T}} μ_{2} W (1 - W - α_{2} M) ψ_{2} d x d t, \end{matrix}$ 19 $\int_{0}^{T} < \partial_{t} C, ψ_{3} >_{{(H^{1})}^{'}, H^{1}} d t + \iint_{Q_{T}} [S (x) \nabla C - C U] \cdot \nabla ψ_{3} d x d t = \iint_{Q_{T}} - (α M + β W) C ψ_{3} d x d t,$ 20 $\begin{matrix} \int_{0}^{T} < \partial_{t} U, ψ >_{V^{'}, V} d t + v \iint_{Q_{T}} \nabla U \cdot \nabla ψ d x d t + κ \iint_{Q_{T}} (U \cdot \nabla) U ψ d x d t = \iint_{Q_{T}} - (γ M + λ W) \nabla Φ ψ d x d t \\ = \iint_{Q_{T}} g ψ d x d t, \end{matrix}$ for all ψ₁, ψ₂, ψ₃ ∈ L²(0, T; H¹(Ω)) and ψ ∈ C_c(0, T;V).

It is known that in dimension 1, nonlinear degenerate systems exhibit global strong solutions. In 2D, blowup is usually avoidable under milder conditions, and we have better regularity for Navier-Stokes due to the Ladyzhenskaya inequality. However, in 3D, the global well-posedness becomes more delicate. With degenerate nonlinear diffusive fluxes, we lose the H¹ regularity of species, but suitable operators are introduced to overcome the problem. The first main aim of this paper is to prove the global existence of (3) given in the following theorem:

Theorem 1.

Given assumptions (10)-(11), the system (3) admits at least one global weak solution in the sense of Definition 1.

This theorem establishes the existence of global weak solutions, but the uniqueness remains a big challenge even in the context of classical Navier-Stokes equations in three dimensions. Hence, under further regularity assumptions on the initial data and for κ = 0, we state the following important result.

Theorem 2.

Under further assumptions (13) till (16), the system (3) coupled with the linear Stokes equation (κ = 0) has a global unique weak solution in space dimensions d ≤ 3.

To validate all the theoretical results, we then turned to the numerical study of system (3). The third main result is a generalized convergent numerical scheme that will be studied and detailed in Section 5. Moreover, an efficient Fortran code will be constructed to study species dynamics under all possible constraints.

Proof of Theorem 1

The main difficulties of this proof are the nonlinearity and the degeneracy of the diffusive and convective tensors. Hence, to avoid the strong degeneracy, we fix ε > 0 and we replace a(M) by a_ε(M) = a(M) + ε and b(W) by b_ε(W) = b(W) + ε.

The construction of classical solutions would require higher regularity of both the unknown functions and the data, which is generally not available. Although we consider here a non-degenerate version of our model, the problem remains highly nonlinear and strongly coupled with the Navier-Stokes equations. In particular, convective terms involve products of functions with limited regularity. For that, we will use a weak formulation to avoid differentiating nonlinearities that may not be smooth. We prove the global existence of non degenerate solutions (M_ε, W_ε, C_ε, U_ε) in the weak sense of Definition 1, such that ∀ψ₁, ψ₂, ψ₃ ∈ L²(0, T ; H¹(Ω)) and ψ ∈ C_c (0, T ; V), 21 $\begin{matrix} \int_{0}^{T} < \partial_{t} M_{ε}, ψ_{1} > d t + \iint_{Q_{T}} [Q (x) a_{ε} (M_{ε}) \nabla M_{ε} - Q (x) χ_{1} (M_{ε}) \nabla C_{ε} - M_{ε} U_{ε}] \cdot \nabla ψ_{1} d x d t \\ = \iint_{Q_{T}} μ_{1} M_{ε} (1 - M_{ε} - α_{1} W_{ε}) ψ_{1} d x d t, \end{matrix}$ 22 $\begin{matrix} \int_{0}^{T} < \partial_{t} W_{ε}, ψ_{2} > d t + \iint_{Q_{T}} [R (x) b_{ε} (W_{ε}) \nabla W_{ε} - R (x) χ_{2} (W_{ε}) \nabla C_{ε} - W_{ε} U_{ε}] \cdot \nabla ψ_{2} d x d t = \\ \iint_{Q_{T}} μ_{2} W_{ε} (1 - W_{ε} - α_{2} M_{ε}) ψ_{2} d x d t, \end{matrix}$ 23 $\iint_{Q_{T}} [\partial_{t} C_{ε} ψ_{3} + [S (x) \nabla C_{ε} - C_{ε} U_{ε}] \cdot \nabla ψ_{3}] d x d t = \iint_{Q_{T}} - (α M_{ε} + β W_{ε}) C_{ε} ψ_{3} d x d t,$ 24 $\int_{0}^{T} < \partial_{t} U_{ε}, Ψ > V^{'}, V d t + v \iint_{Q_{T}} \nabla U_{ε} \cdot \nabla ψ d x d t + κ \iint_{Q_{T}} (U_{ε} \cdot \nabla) U_{ε} ψ d x d t = \iint_{Q_{T}} - (γ M_{ε} + λ W_{ε}) \nabla ϕ ψ d x d t,$ and from the definition of V, one has 25 $\nabla \cdot U_{ε} = 0.$

3.1

Weak solution of the non-degenerate problem

In this subsection, we discretize the time variable interval [0, T] to get approximate solutions for the non-degenerate problem in the sense of (21)-(25). For that, we define $h = \frac{T}{Ñ}$ as a constant time step such that $Ñ \in ℕ^{*}$ . Then, we define gradually from $n \in {0, \dots, Ñ - 1}$ the following problems: Given the initial data $(M_{0, ε}^{Ñ}, W_{0, ε}^{Ñ}, C_{0, ε}^{Ñ}, U_{0, ε}^{Ñ}) = (M_{0}, W_{0}, C_{0}, U_{0})$ , and when $(M_{n, ε}^{Ñ}, W_{n, ε}^{Ñ}, C_{n, ε}^{Ñ}, U_{n, ε}^{Ñ})$ are known, we can define $(M_{n + 1, ε}^{Ñ}, W_{n + 1, ε}^{Ñ}, C_{n + 1, ε}^{Ñ}, U_{n + 1, ε}^{Ñ})$ that satisfy: 26 $\begin{matrix} \frac{1}{h} \int_{Ω} (M_{n + 1, ε}^{Ñ} - M_{n, ε}^{Ñ}) ψ_{1} d x + \int_{Ω} (Q (x) a_{ε} (M_{n + 1, ε}^{Ñ}) \nabla M_{n + 1, ε}^{Ñ} - Q (x) χ_{1} (M_{n + 1, ε}^{Ñ}) \nabla C_{n + 1, ε}^{Ñ}) \cdot \nabla ψ_{1} d x \\ - \int_{Ω} M_{n + 1, ε}^{Ñ} U_{n, ε}^{Ñ} \cdot \nabla ψ_{1} d x = \int_{Ω} μ_{1} M_{n + 1, ε}^{Ñ} (1 - M_{n + 1, ε}^{Ñ} - α_{1} W_{n, ε}^{Ñ}) ψ_{1} d x, \end{matrix}$ 27 $\begin{matrix} \frac{1}{h} \int_{Ω} (W_{n + 1, ε}^{Ñ} - W_{n, ε}^{Ñ}) ψ_{2} d x + \int_{Ω} (R (x) b_{ε} (W_{n + 1, ε}^{Ñ}) \nabla W_{n + 1, ε}^{Ñ} - R (x) χ_{2} (W_{n + 1, ε}^{Ñ}) \nabla C_{n + 1, ε}^{Ñ}) \cdot \nabla ψ_{2} d x \\ - \int_{Ω} W_{n + 1, ε}^{Ñ} U_{n, ε}^{Ñ} \cdot \nabla ψ_{2} d x = \int_{Ω} μ_{2} W_{n + 1, ε}^{Ñ} (1 - W_{n + 1, ε}^{Ñ} - α_{2} M_{n, ε}^{Ñ}) ψ_{2} d x, \end{matrix}$ 28 $\begin{matrix} \frac{1}{h} \int_{Ω} (C_{n + 1, ε}^{Ñ} - C_{n, ε}^{Ñ}) ψ_{3} d x + \int_{Ω} S (x) \nabla C_{n + 1, ε}^{Ñ} \cdot \nabla ψ_{3} d x - \int_{Ω} C_{n + 1, ε}^{Ñ} U_{n, ε}^{Ñ} \cdot \nabla ψ_{3} d x \\ = \int_{Ω} - (α M_{n, ε}^{Ñ} + β W_{n, ε}^{Ñ}) C_{n + 1, ε}^{Ñ} ψ_{3} d x, \end{matrix}$ 29 $\begin{matrix} \frac{1}{h} \int_{Ω} (U_{n + 1, ε}^{Ñ} - U_{n, ε}^{Ñ}) ψ d x + v \int_{Ω} \nabla U_{n + 1, ε}^{Ñ} \cdot \nabla ψ d x + κ \int_{Ω} (U_{n + 1, ε}^{Ñ} \cdot \nabla) U_{n + 1, ε}^{Ñ} ψ d x . \\ = \int_{Ω} - (γ M_{n, ε}^{Ñ} + λ W_{n, ε}^{Ñ}) \nabla ϕ ψ d x, \end{matrix}$ ∀ψ₁, ψ₂, ψ₃ ∈ H¹(Ω) and ∀ψ ∈ V.

3.1.1

Confinement of approximate solutions

To ensure the admissibility of the approximate solutions and their biological "Volume-filling" interpretation, we prove the following result.

Proposition 3.

There exists ζ > 0, independent of ε, such that for all n = 0, ⋯, $Ñ - 1$ , $0 \leq M_{n + 1, ε}^{Ñ} \leq 1, 0 \leq W_{n + 1, ε}^{Ñ} \leq 1 and 0 \leq C_{n + 1, ε}^{Ñ} \leq ζ a .e . x \in Ω .$

Proof.

The initialization is well satisfied due to (10). Suppose that the assertion is true for the rank n and prove v it for the rank n + 1. Let us start by proving by induction that $M_{n + 1, ε}^{Ñ} \geq 0$ a.e. x ∈ Ω. Suppose that $M_{n, ε}^{Ñ} \geq 0$ , $M_{n + 1, ε}^{Ñ} < 0$ and take $Ψ_{1} = - {(M_{n + 1, ε}^{Ñ})}^{-}$ in (26). Thus, $\begin{matrix} - \frac{1}{h} \int_{Ω} (M_{n + 1, ε}^{Ñ} - M_{n, ε}^{Ñ}) {(M_{n + 1, ε}^{Ñ})}^{-} d x - \int_{Ω} (D_{ε} (x, M_{n + 1, ε}^{Ñ}) \nabla M_{n + 1, ε}^{Ñ} - Q (x) χ_{1} (M_{n + 1, ε}^{Ñ}) \nabla C_{n + 1, ε}^{Ñ}) \cdot \nabla {(M_{n + 1, ε}^{Ñ})}^{-} d x \\ + \int_{Ω} M_{n + 1, ε}^{Ñ} U_{n, ε}^{Ñ} \cdot \nabla {(M_{n + 1, ε}^{Ñ})}^{-} d x = - \int_{Ω} μ_{1} h (M_{n + 1, ε}^{Ñ}) {(M_{n + 1, ε}^{Ñ})}^{-} d x + \int_{Ω} μ_{1} α_{1} M_{n + 1, ε}^{Ñ} W_{n, ε}^{Ñ} {(M_{n + 1, ε}^{Ñ})}^{-} d x . \end{matrix}$

We remark that $\int_{Ω} M_{n + 1, ε}^{Ñ} U_{n, ε}^{Ñ} \cdot \nabla {(M_{n + 1, ε}^{Ñ})}^{-} d x = \frac{1}{2} \int_{Ω} (U_{n, ε}^{Ñ} \cdot \nabla {[{(M_{n + 1, ε}^{Ñ})}^{-}]}^{2}) d x = 0$ . Consequently, $- \frac{1}{h} \int_{Ω} M_{n + 1, ε}^{Ñ} {(M_{n + 1, ε}^{Ñ})}^{-} d x \leq 0$ due to the coercivity of the diffusive operator, the positivity of $W_{n, ε}^{Ñ}$ and the continuous extension of χ₁ and h(M) = M(1 – M) by zero outside [0, 1]. Therefore, the last inequality achieved leads to a contradiction. The proof of this Proposition can be easily achieved using similar arguments.

3.1.2

Existence of solutions for (26)-(29)

Due to Proposition 3, Hypothesis (5) and in space dimension d ≤ 3, the existence of a discrete solution $U_{n + 1, ε}^{Ñ} \in V$ of (29) is a very well-known result. Moreover, there exists a unique solution $C_{n + 1, ε}^{Ñ} \in H^{1} (Ω)$ solution of (28) due to Lax-Migram theorem.

The existence of $M_{n + 1}^{Ñ}$ (resp. $W_{n + 1}^{Ñ}$ ) solution of (26) (resp. (27)) is proved using the Schauder fixed point theorem. Let us introduce an application θ defined over and into a closed convex subset $D = {M_{n + 1, ε}^{Ñ} \in L^{2} (Ω); 0 \leq M_{n + 1, ε}^{Ñ} (x) \leq 1$ , for a.e. x ∈ Ω} as $θ (w) = M_{n + 1, ε}^{Ñ}$ solution in H¹(Ω) of 30 $\begin{matrix} \frac{1}{h} \int_{Ω} (M_{n + 1, ε}^{Ñ} - M_{n, ε}^{Ñ}) ψ_{1} d x + \int_{Ω} Q (x) a_{ε} (w) \nabla M_{n + 1, ε}^{Ñ} \cdot \nabla ψ_{1} d x - \int_{Ω} Q (x) χ_{1} (w) \nabla C_{n + 1, ε}^{Ñ} \cdot \nabla ψ_{1} d x . \\ + \int_{Ω} M_{n + 1, ε}^{Ñ} U_{n, ε}^{Ñ} \cdot \nabla ψ_{1} d x = \int_{Ω} μ_{1} w (1 - w - α_{1} W_{n, ε}^{Ñ}) ψ_{1} d x . \end{matrix}$

The application θ is well defined due to the Lax-Milgram theorem. Now, we prove that θ(𝒟) is relatively compact in L²(Ω). We take $ψ_{1} = M_{n + 1, ε}^{Ñ}$ as a test function in (30). Due to Proposition 3, (6)-(9) and Young inequality, there exist a₁, a₂, a₃ and a₄ strictly positive numbers such that $\begin{matrix} \frac{1}{h} \int_{Ω} {| M_{n + 1, ε}^{Ñ} |}^{2} d x + γ^{'} \int_{Ω} {| \nabla M_{n + 1, ε}^{Ñ} |}^{2} d x \\ \leq \bar{D} 1 a_{1} {‖ \nabla M_{n + 1, ε}^{Ñ} ‖}_{L^{2}}^{2} + \frac{D 1}{a_{1}} {‖ \nabla C_{n + 1, ε}^{Ñ} ‖}_{L^{2}}^{2} + \frac{a_{2}}{h} {‖ M_{n + 1, ε}^{Ñ} ‖}_{L^{2}}^{2} + \frac{1}{a_{2} h} {‖ M_{n, ε}^{Ñ} ‖}_{L^{2}}^{2} \\ + μ_{1} a_{3} {‖ M_{n + 1, ε}^{Ñ} ‖}_{L^{2}}^{2} + \frac{μ_{1}}{a_{3}} ‖ w (1 - w) ‖_{L^{2}}^{2} + μ_{1} α_{1} a_{4} {‖ M_{n + 1, ε}^{Ñ} ‖}_{L^{2}}^{2} + \frac{μ_{1} α_{1}}{a_{4}} {‖ W_{n, ε}^{Ñ} ‖}_{L^{2}}^{2} . \end{matrix}$

One can choose $a_{1} = \frac{γ^{'}}{2 \bar{D 1}}, a_{2} = \frac{1}{6}, a_{3} = \frac{1}{6 h μ_{1}}$ and $a_{4} = \frac{1}{6 h μ_{1} α_{1}}$ to show finally that 31 ${‖ M_{n + 1, ε}^{Ñ} ‖}_{H^{1} (Ω)} \leq \frac{1}{2 h} \int_{Ω} {| M_{n + 1, ε}^{Ñ} |}^{2} d x + \frac{γ^{'}}{2} \int_{Ω} {| \nabla M_{n + 1, ε}^{Ñ} |}^{2} d x \leq C,$ where C is a constant independent of w.

It remains now to prove that θ is a continuous mapping. For that, we consider any convergent sequence (w_n)_n to w in 𝒟 and prove that θ(w_n) converges to θ(w) in L²(Ω) as n → +∞. Remark that for any subsequence (W_{n_j})_j of (W_n)_n and due to (8), D_ε(x, w_{n_j}) = Q(x)a_ε(w_{n_j}) (resp. D₁(x, w_{n_j}) = Q(x)χ(w)) converges to D_ε(x, w) (resp. D₁(x, w)) in (L²(Ω))^d² using dominated convergence theorem. Next, the equation (30) can be written in terms of subsequences (M_{n_j})_q from the bounded sequence (M_n)_n in H¹(Ω) as: $\begin{matrix} \frac{1}{h} \int_{Ω} (M_{{(n_{j})}_{q}} - M_{n, ε}^{Ñ}) ψ_{1} d x + \int_{Ω} D_{ε} (x, w_{{(n_{j})}_{q}}) \nabla M_{{(n_{j})}_{q}} \cdot \nabla ψ_{1} d x - \int_{Ω} D_{1} (x, w_{{(n_{j})}_{q}}) \nabla C_{n + 1, ε}^{Ñ} \cdot \nabla ψ_{1} d x \\ + \int_{Ω} M_{{(n_{j})}_{q}} (U_{n, ε}^{Ñ} \cdot \nabla ψ_{1}) d x = μ_{1} \int_{Ω} w_{{(n_{j})}_{q}} (1 - w_{{(n_{j})}_{q}} - α_{1} W_{n, ε}^{Ñ}) ψ_{1} d x . \end{matrix}$

By passing to the limit as q → +∞ and due to the weak (resp. strong) convergence of M_{(n_j)}_q to M in the Hilbert space H¹(Ω) (resp. in L²(Ω) and a.e. in Ω), we can deduce that M = θ(w). Therefore, θ(w) is an accumulation point of the sequence (M_n)_n which belongs to a relatively compact subset. Then, the entire sequence (M_n)_n tends to θ(w) in L²(Ω). Finally, the existence of a fixed point $M_{n + 1, ε}^{N}$ for θ, solution of (26), is a straightforward consequence of the Schauder fixed point theorem.

3.1.3

Estimates

Let us recall the constant interpolation operator ℐ⁰ defined as $ℐ^{0} w (t) = \sum_{n = 0}^{Ñ - 1} w_{n + 1} ϑ_{]_{n h, (n + 1) h}]} (t)$ and the linear interpolation operator ℐ¹ defined as $ℐ^{1} w (t) = \sum_{n = 0}^{Ñ - 1} [(1 + n - \frac{t}{h}) w_{n} + (\frac{t}{h} - n) w_{n + 1}] ϑ_{] n h, (n + 1) h]} (t)$ , where ϑ_]nh,(n+1)h](t) being the characteristic equation. Consequently, $\frac{d}{d t} (ℐ^{1} w (t)) = \sum_{n = 0}^{Ñ - 1} [\frac{1}{h} (w_{n + 1} - w_{n})] ϑ_{]_{n h, (n + 1) h}]} (t)$ .

First of all, we start by recalling some estimates for the equation (29). Indeed, the regularity of $g = - (λ ℐ^{0} M_{ε}^{Ñ} + γ ℐ^{0} W_{ε}^{Ñ}) \nabla ϕ \in L^{\infty} (Q_{T})$ is sufficient to deduce the existence of U_ε ∈ L^∞(0, T;H) ⋂ L²(0, T;V) such that, modulo a subsequence, as $Ñ$ tends to +∞, 32 ${‖ \frac{\partial}{\partial t} (ℐ^{1} U_{ε}^{Ñ}) ‖}_{L^{2} (0, T; {(H^{1} (Ω))}^{'})} \leq d_{1},$ 33 $ℐ^{0} U_{ε}^{Ñ} ⇀ U_{ε} weak * in L^{\infty} (0, T; H) and weakly in L^{2} (0, T; V),$ 34 $\forall h > 0, {‖ τ_{- h} ℐ^{0} U_{ε}^{Ñ} - ℐ^{0} U_{ε}^{Ñ} ‖}_{L^{2} (0, T - h; (L^{2} (Ω)))} \leq d_{1} h; d_{1} = {‖ U_{0} ‖}_{H}^{2} + \int_{0}^{T} ‖ g (s) ‖_{V^{'}}^{2} d s .$

Consequently, 35 $ℐ^{0} U_{ε}^{Ñ} \to U_{ε} strongly in L^{2} (0, T; H), as Ñ tends to + \infty .$

Next, uniform a priori estimates on the interpolation of $M^{Ñ}, W^{Ñ}$ and $C^{Ñ}$ with respect to $Ñ$ are obtained. Proposition 3 implies the boundedness in L^∞(Q_T). Thus, there exist M_ε, W_ε and C_ε such that, modulo a subsequence, as $Ñ$ tends to +∞, 36 $ℐ^{0} C_{ε}^{Ñ} ⇀ C_{ε}, ℐ^{0} M_{ε}^{Ñ} ⇀ M_{ε} and ℐ^{0} W_{ε}^{Ñ} ⇀ W_{ε} weakly- * in L^{\infty} (Q_{T}) .$

Choose $ψ_{3} = C_{n + 1, ε}^{Ñ}$ as a function test in (28). Consider (8), Proposition 3 and the inequality $(a - b) a \geq \frac{1}{2} (a^{2} - b^{2})$ to deduce that, $\frac{1}{h} \int_{Ω} [{(C_{n + 1, ε}^{Ñ})}^{2} - {(C_{n, ε}^{Ñ})}^{2}] d x + c_{S} \int_{Ω} {(\nabla C_{n + 1, ε}^{Ñ})}^{2} d x \leq C_{1}; C_{1} is a constant independent of Ñ .$

Now, we multiply by h and we sum from n = 0 to $n = Ñ - 1$ , $\int_{Ω} {(C_{Ñ, ε}^{Ñ})}^{2} d x + c_{S} {‖ \nabla ℐ^{0} C_{ε}^{Ñ} ‖}_{L^{2} (0, T; L^{2} (Ω))}^{2} \leq C_{1} T + \int_{Ω} {(C_{0, ε}^{Ñ})}^{2} d x = d_{2}; d_{2} independent of Ñ .$

The same strategy can be followed in the equations (26) and (27) to conclude that: 37 $ℐ^{0} C_{ε}^{Ñ} ⇀ C_{ε}, ℐ^{0} M_{ε}^{Ñ} ⇀ M_{ε} and ℐ^{0} W_{ε}^{Ñ} ⇀ W_{ε} weakly in L^{2} (0, T; H^{1} (Ω)) as Ñ tends to + \infty .$

To obtain strong convergence in L²(Q_T), we need to ensure time translate estimates for approximate solutions. Let h > 0, one has, $\begin{matrix} I : = {‖ τ_{- h} ℐ^{0} M_{ε}^{Ñ} - ℐ^{0} M_{ε}^{Ñ} ‖}_{L^{2} (0, T - h; (L^{2} (Ω)))}^{2} = \int_{0}^{T - h} \int_{Ω} {(ℐ^{0} M_{ε}^{Ñ} (t + h, x) - ℐ^{0} M_{ε}^{Ñ} (t, x))}^{2} d x d t \\ = \int_{0}^{T - h} \int_{Ω}^{n_{1} - 1} (M_{n = n_{0}}^{Ñ} (x) - M_{n, ε}^{Ñ} (x)) (M_{n_{1}, ε}^{Ñ} (x) - M_{n_{0}, ε}^{Ñ} (x)) d x d t, \end{matrix}$ where $n_{0} (t) = [\frac{t}{h}] (resp . n_{1} (t) = [\frac{t + h}{h}])$ is the integer part of the real $\frac{t}{h} (resp . \frac{t}{h} + 1)$ . We consider the equation (26) with $ψ_{1} = (M_{n_{1}, ε}^{Ñ} - M_{n_{0}, ε}^{Ñ})$ and add from n₀ to n₁ – 1 in order to obtain that $\begin{matrix} I \leq \int_{0}^{T - h} (\sum_{n = n_{0}}^{n_{1} - 1} [- \int_{Ω} D_{ε} (x, M_{n + 1, ε}^{Ñ}) \nabla M_{n + 1, ε}^{Ñ} \cdot \nabla (M_{n_{1}, ε}^{Ñ} - M_{n_{0}, ε}^{Ñ}) d x \\ - \int_{Ω} S (x) χ (M_{n + 1, ε}^{Ñ}) \nabla C_{n + 1, ε}^{Ñ} \cdot \nabla (M_{n_{1}, ε}^{Ñ} - M_{n_{0}, ε}^{Ñ}) d x - \int_{Ω} M_{n + 1, ε}^{Ñ} U_{n, ε}^{Ñ} \cdot \nabla (M_{n_{1}, ε}^{Ñ} - M_{n_{0}, ε}^{Ñ}) d x \\ - \int_{Ω} μ_{1} M_{n + 1, ε}^{Ñ} (1 - M_{n + 1, ε}^{Ñ} - α_{1} W_{n, ε}^{Ñ}) (M_{n_{1}, ε}^{Ñ} - M_{n_{0}, ε}^{Ñ}) d x]) \leq d_{4} \int_{0}^{T - h} χ_{n} (t, t + h) d t, \end{matrix}$ where $\int_{0}^{T - h} χ_{n} (t, t + h) d t \leq h$ . Therefore, 38 $\forall h > 0, {‖ τ_{- h} ℐ^{0} M_{ε}^{Ñ} - ℐ^{0} M_{ε}^{Ñ} ‖}_{L^{2} (0, T - h; (L^{2} (Ω)))} \leq d_{3} h; d_{3} independent of Ñ .$

Due to Kolmogorov’s compactness criterion under translate estimates (38), we can conclude with similar steps for (26) and (27) that the sequences $(ℐ^{0} M_{ε}^{Ñ}), (ℐ^{0} W_{ε}^{Ñ})$ and $(ℐ^{0} C_{ε}^{Ñ})$ are relatively compact in L²(0, T;L²(Ω)). Hence, modulo subsequences, as $Ñ$ tends to +∞, 39 $ℐ^{0} C_{ε}^{Ñ} \to C_{ε}, ℐ^{0} M_{ε}^{Ñ} \to M_{ε} and ℐ^{0} W_{ε}^{Ñ} \to W_{ε} strongly in L^{2} (Q_{T}) and a .e . in Q_{T} .$

By choosing ψ₁ ∈ H¹(Ω) as a test function in (26) and using property (8) and Cauchy-Schwarz inequality, one can deduce that: 40 $\begin{matrix} | \frac{1}{h} \int_{Ω} (M_{n + 1, ε}^{Ñ} - M_{n, ε}^{Ñ}) ψ_{1} d x | \leq \bar{D} {‖ \nabla M_{n + 1, ε}^{Ñ} ‖}_{{(L^{2} (Ω))}^{d}} {‖ \nabla ψ_{1} ‖}_{{(L^{2} (Ω))}^{d}} + {\bar{D}}_{1} {‖ \nabla C_{n + 1, ε}^{Ñ} ‖}_{{(L^{2} (Ω))}^{d}} {‖ \nabla ψ_{1} ‖}_{{(L^{2} (Ω))}^{d}} \\ + ζ {‖ U_{n, ε}^{Ñ} ‖}_{{(L^{2} (Ω))}^{d}} {‖ \nabla ψ_{1} ‖}_{{(L^{2} (Ω))}^{d}} + μ_{1} {‖ ψ_{1} ‖}_{{(L^{2} (Ω))}^{d}} . \end{matrix}$

Then, we can simplify by ∥ψ₁∥_H¹(ω), raise to the square, and use Young’s inequality to show that: $\begin{array}{l} \forall n \in {0, \dots Ñ} \\ \frac{1}{h^{2}} {‖ M_{n + 1, ε}^{Ñ} - M_{n, ε}^{Ñ} ‖}_{{(H^{1} (Ω))}^{'}}^{2} \leq 2 ({\bar{D}}^{2} {‖ \nabla M_{n + 1, ε}^{Ñ} ‖}_{{(L^{2} (Ω))}^{d}}^{2} + {\bar{D}}_{1}^{2} {‖ \nabla C_{n + 1, ε}^{Ñ} ‖}_{{(L^{2} (Ω))}^{d}}^{2} + ζ^{2} {‖ U_{n, ε}^{Ñ} ‖}_{{(L^{2} (Ω))}^{d}}^{2} + μ_{1}^{2}) . \end{array}$

Hence, after adding from n = 0 to $n = Ñ - 1$ and multiplying by h, one has, 41 ${‖ \frac{\partial}{\partial t} (ℐ^{1} M_{ε}^{Ñ}) ‖}_{L^{2} (0, T; {(H^{1} (Ω))}^{'})}^{2} = \sum_{n = 0}^{Ñ - 1} \frac{1}{h} {‖ M_{n + 1, ε}^{Ñ} - M_{n, ε}^{Ñ} ‖}_{{(H^{1} (Ω))}^{'}}^{2} \leq d_{4}; d_{4} independent of Ñ .$

From the estimate (41) and $\int_{n h}^{(n + 1) h} {(1 + n - \frac{t}{h})}^{2} d t = \frac{h}{3}$ , one can deduce that: 42 ${‖ ℐ^{1} M_{ε}^{Ñ} - ℐ^{0} M_{ε}^{Ñ} ‖}_{L^{2} (0, T; {(H^{1} (Ω))}^{'})} = \sum_{n = 0}^{Ñ - 1} \int_{n h}^{(n + 1) h} {‖ (1 + n - \frac{t}{h}) [M_{n, ε}^{Ñ} - M_{n + 1, ε}^{Ñ}] ‖}_{{(H^{1} (Ω))}^{'}}^{2} d t \leq \frac{d_{4} h^{2}}{3} .$

The sequence $(ℐ^{1} M_{ε}^{Ñ})$ is bounded in L^∞(Q_T) and the sequence $\frac{\partial}{\partial t} (ℐ^{1} M_{ε}^{Ñ})$ is bounded in L²(0, T; (H¹ (Ω))′), then $ℐ^{1} M_{ε}^{Ñ} \to w$ strongly in C(0, T; (H¹(Ω))′), as $Ñ \to + \infty$ due to the compact embedding of the space L^∞(Ω) into (H¹(Ω))′.

Otherwise, estimates (39) and (42) imply that $ℐ^{1} M_{ε}^{Ñ} \to M_{ε}$ strongly in L²(0, T; (H¹(Ω))′). Due to the uniqueness of the limit, w = M_ε. Therefore and with similar reasoning, we have 43 $ℐ^{1} C_{ε}^{Ñ} ⇀ C_{ε}, ℐ^{1} M_{ε}^{Ñ} ⇀ M_{ε} and ℐ^{1} W_{ε}^{Ñ} ⇀ W_{ε} weakly-*
in L^{\infty} (Q_{T}) as Ñ tends
to + \infty .$

Similar arguments for (27) and (28) imply that, 44 $\frac{\partial}{\partial t} (ℐ^{1} C_{ε}^{Ñ}) ⇀ \frac{\partial C_{ε}}{\partial t}, \frac{\partial}{\partial t} (ℐ^{1} M_{ε}^{Ñ}) ⇀ \frac{\partial M_{ε}}{\partial t} and \frac{\partial}{\partial t} (ℐ^{1} W_{ε}^{Ñ}) ⇀ \frac{\partial W_{ε}}{\partial t} weakly in L^{2} (0, T; {(H^{1} (Ω))}^{'}) .$

Moreover, M_ε belongs to C(0, T; (H 1(Ω))′) ⊂ C_w(0, T; (H¹(Ω))′) and to L^∞(0, T; L²(Ω)), then M_ε ∈ C_w(0, T; L²(Ω)). We know also that $ℐ^{1} M_{ε}^{Ñ} (0)$ converges strongly to M_ε(0, x) in (H¹(Ω))′ and by definition $ℐ^{1} M_{ε}^{Ñ} (0) = M_{0, ε} (x)$ . Thus, M_ε(0, x) = M_{0, ε}(x).

Finally and for $n_{0} = [\frac{t}{h}] + 1$ , we have by definition $ℐ^{0} M_{ε}^{Ñ} (t, x) = M_{n_{0}, ε} (x)$ that converges to M_ε(t, x) as $Ñ$ tends to +∞. As we have M_n₀,ε(x) ∈ [0, 1] then M_ε(t, x) ∈ [0, 1]. Using similar arguments, we also obtain the confinement of W_ε and C_ε.

3.1.4

Passing to the limit

It remains to be shown that the limit functions (M_ε, W_ε, C_ε, U_ε) constitute a weak solution of the continuous problem (21)-(25).

Using interpolation operators, we can re-write the equation (29) as, $\begin{matrix} \int_{0}^{T} < \partial_{t} (ℐ^{1} U_{ε}^{Ñ}), ψ > d t + v \iint_{Q_{T}} \nabla (ℐ^{0} U_{ε}^{Ñ}) \cdot \nabla ψ d x d t + κ \iint_{Q_{T}} (ℐ^{0} U_{ε}^{Ñ} \cdot \nabla) ℐ^{0} U_{ε}^{Ñ} ψ d x d t \\ = - \iint_{Q_{T}} (γ τ_{h} ℐ^{0} M_{ε}^{Ñ} + λ τ_{h} ℐ^{0} W_{ε}^{Ñ}) \nabla ϕ \cdot ψ d x d t, \forall ψ \in L^{2} (0, T; V) . \end{matrix}$

This equation tends to (24), when $Ñ$ goes to +∞. Indeed, we use (32) for the time evolution term, (33) for the diffusive term, (33) and (35) for the convective term and (5), (38) and (39) for the reaction term.

For the equation (28), we can write the terms under the following form: $\begin{matrix} \int_{0}^{T} < \partial_{t} (ℐ^{1} C_{ε}^{Ñ}), ψ_{3} > d t + \iint_{Q_{T}} S (x) \nabla (ℐ^{0} C_{ε}^{Ñ}) \cdot \nabla ψ_{3} d x d t - \iint_{Q_{T}} ℐ^{0} C_{ε}^{Ñ} (τ_{h} ℐ^{0} U_{ε}^{Ñ} \cdot \nabla ψ_{3}) d x d t \\ = - \iint_{Q_{T}} (α τ_{h} ℐ^{0} M_{ε}^{Ñ} + β τ_{h} ℐ^{0} W_{ε}^{Ñ}) ℐ^{0} C_{ε}^{Ñ} ψ_{3} d x d t, \forall ψ_{3} \in L^{2} (0, T; H^{1} (Ω)) . \end{matrix}$

All the terms of this equation converge to their corresponding terms in (23), when $Ñ$ goes to +∞. Indeed, estimate (44) (resp. (37) and the boundedness of the tensor R) is needed for the convergence of the time evolution term (resp. diffusive flux). For the convective term and due to (33), (34) and (35), one has $\begin{matrix} | \iint_{Q_{T}} ℐ^{0} C_{ε}^{Ñ} (τ_{h} ℐ^{0} U_{ε}^{Ñ} \cdot \nabla ψ_{3}) d x d t - \iint_{Q_{T}} C_{ε} U_{ε} \cdot \nabla ψ_{3} d x d t | \leq {‖ (ℐ^{0} C_{ε}^{Ñ} - C_{ε}) \cdot \nabla ψ_{3} ‖}_{L^{2} (Q_{T})} {‖ τ_{h} ℐ^{0} U_{ε}^{Ñ} ‖}_{L^{2} (Q_{T})} \\ + {‖ C_{ε} ‖}_{L^{\infty} (Q_{T})} {‖ τ_{h} ℐ^{0} U_{ε}^{Ñ} - U_{ε} ‖}_{L^{2} (Q_{T})} {‖ \nabla ψ_{3} ‖}_{L^{2} (Q_{T})} \to 0. \end{matrix}$

For the reaction term, Proposition (3) and estimates (37), (38) and (39) imply that $\begin{matrix} | \iint_{Q_{T}} (α τ_{h} ℐ^{0} M_{ε}^{Ñ} + β τ_{h} ℐ^{0} W_{ε}^{Ñ}) ℐ^{0} C_{ε}^{Ñ} ψ_{3} d x d t - \iint_{Q_{T}} (α M_{ε} + β W_{ε}) C_{ε} ψ_{3} d x d t | \\ \leq (α {‖ τ_{h} ℐ^{0} M_{ε}^{Ñ} ‖}_{L^{\infty} (Q_{T})} + β {‖ τ_{h} ℐ^{0} W_{ε}^{Ñ} ‖}_{L^{\infty} (Q_{T})}) {‖ ℐ^{0} C_{ε}^{Ñ} - C_{ε} ‖}_{L^{2} (Q_{T})} {‖ ψ_{3} ‖}_{L^{2} (Q_{T})} \\ + {‖ C_{ε} ‖}_{L^{\infty} (Q_{T})} (α {‖ τ_{h} ℐ^{0} M_{ε}^{Ñ} - M_{ε} ‖}_{L^{2} (Q_{T})} + β {‖ τ_{h} ℐ^{0} W_{ε}^{Ñ} - W_{ε} ‖}_{L^{2} (Q_{T})} {‖ ψ_{3} ‖}_{L^{2} (Q_{T})} \to 0. \end{matrix}$

For the equation (26), we have $\begin{matrix} \int_{0}^{T} < \partial_{t} (ℐ^{1} M_{ε}^{Ñ}), ψ_{1} > d t + \iint_{Q_{T}} D_{ε} (x, ℐ^{0} M_{ε}^{Ñ}) \nabla (ℐ^{0} M_{ε}^{Ñ}) \cdot \nabla ψ_{1} d x d t \\ = \iint_{Q_{T}} D_{1} (x, ℐ^{0} M_{ε}^{Ñ}) \nabla (ℐ^{0} C_{ε}^{Ñ}) \cdot \nabla ψ_{1} d x d t + \iint_{Q_{T}} ℐ^{0} M_{ε}^{Ñ} (τ_{h} ℐ^{0} U_{ε}^{Ñ} \cdot \nabla ψ_{1}) d x d t \\ + μ_{1} \iint_{Q_{T}} ℐ^{0} M_{ε}^{Ñ} (1 - ℐ^{0} M_{ε}^{Ñ} - α_{1} τ_{h} ℐ^{0} W_{ε}^{Ñ}) ψ_{1} d x d t, \forall ψ_{1} \in L^{2} (0, T; H^{1} (Ω)) . \end{matrix}$

At the limit, we obtain the equation (21). Indeed, estimate (44) guides the time evolution term for convergence. For the diffusive flux, assumption (8), the a.e. convergence from estimate (39) and the dominated convergence theorem imply that $D_{ε} (x, ℐ^{0} M_{ε}^{Ñ}) \nabla ψ_{1}$ converges to D_ε(x, M_ε)∇ψ₁ in L²(Q_T). Consequently, $\begin{matrix} | \iint_{Q_{T}} ℐ^{0} M_{ε}^{Ñ} (τ_{h} ℐ^{0} U_{ε}^{Ñ} \cdot \nabla ψ_{1}) d x d t - \iint_{Q_{T}} D_{ε} (x, M_{ε}) \nabla M_{ε} \cdot \nabla ψ_{1} d x d t | \leq \\ {‖ (D_{ε} (x, ℐ^{0} M_{ε}^{Ñ}) - D_{ε} (x, M_{ε})) \nabla ψ_{1} ‖}_{L^{2} (Q_{T})} {‖ \nabla ℐ^{0} M_{ε}^{Ñ} ‖}_{L^{2} (Q_{T})} + \iint_{Q_{T}} D_{ε} (x, M_{ε}) (\nabla ℐ^{0} M_{ε}^{Ñ} - \nabla M_{ε}) \cdot \nabla ψ_{1} d x d t \to 0. \end{matrix}$

Similar reasons are followed for the convective terms. Finally, the strong convergences in estimate (39) imply that $μ_{1} \iint_{Q_{T}} ℐ^{0} M_{ε}^{Ñ} (1 - ℐ^{0} M_{ε}^{Ñ} - α_{1} τ_{h} ℐ^{0} W_{ε}^{Ñ}) ψ_{1} d x d t$ converges to $μ_{1} \iint_{Q_{T}} M_{ε} (1 - M_{ε} - α_{1} W_{ε}) ψ_{1} d x d t$ , as $Ñ$ tends to +∞.

3.2

Weak solution of the degenerate problem

To prove the existence of weak solutions for the original degenerate problem, we tend the regularization parameter ε to zero. In the sequel, a priori estimates will be proved to conclude.

3.2.1

Estimates

For the U equation: Take the initial condition U_{ε, 0} ∈ H and remark that the reaction term g = – (γM_ε + λW_ε)∇ϕ ∈ L^∞(Q_T) ⊂ L²(0, T, V′). Then, for ψ = U_ε as a test function in (24), one can conclude that, as ε tends to 0, there exists U ∈ L²(0, T; V) ⋂ L^∞(0, T; H); 45 $U_{ε} ⇀ U^{weakly ⋆} in L^{\infty} (0, T; H) and U_{ε} ⇀ U weakly in L^{2} (0, T; V) .$

Moreover, an application of the Aubin-Simon compactness theorem implies that, modulo a subsequence, 46 $U_{ε} \to U in L^{2} (0, T; H), as ε \to 0.$

Furthermore, we can easily prove that $\frac{d U}{d t} \in L^{1} (0, T; V^{'})$ and $P \in W^{- 1, \infty} (0, T; L_{0}^{2} (Ω))$ .

For the C equation: Remark that the reaction term (αM_ε + βW_ε)C_ε ∈ L^∞(Q_T) and take ψ₃ = C_ε as a test function in (23). Therefore, classical results for parabolic equations imply that there exists a solution C ∈ L^∞(Q_T) ⋂ L²(0, T; H¹(Ω)) such that, modulo a subsequence and when ε goes to 0, 47 $C_{ε} ⇀ C^{weakly✶ ⋆} in L^{\infty} (Q_{T}) and C_{ε} ⇀ C weakly in L^{2} (0, T; H^{1} (Ω)) .$

In addition to that, 48 $\frac{\partial C_{ε}}{\partial t} ⇀ \frac{\partial C}{\partial t} weakly in L^{2} (0, T; {(H^{1} (Ω))}^{'}) .$

For the M equation: Take ψ₁ = A_ε(M_ε) = A(M_ε) + εM_ε ∈ (21) to obtain $\begin{matrix} \int_{0}^{T} < \partial_{t} (M_{ε}), A_{ε} (M_{ε}) > d x + \iint_{Q_{T}} Q (x) \nabla A_{ε} (M_{ε}) \cdot \nabla A_{ε} (M_{ε}) d x d t = \iint_{Q_{T}} M_{ε} U_{ε} \cdot \nabla A_{ε} (M_{ε}) d x d t \\ + \iint_{Q_{T}} Q (x) χ_{1} (M_{ε}) \nabla C_{ε} \cdot \nabla A_{ε} (M_{ε}) d x d t - μ_{1} \iint_{Q_{T}} M_{ε} (1 - M_{ε} - α_{1} W_{ε}) A_{ε} (M_{ε}) d x d t . \end{matrix}$

Using Young’s inequality, Proposition (3), hypothesis (7) and the uniform continuity of χ, one has 49 $\begin{matrix} \sup_{0 \leq t \leq T} \int_{Ω} A (M_{ε}) (x, t) d x + ε \sup_{0 \leq t \leq T} \int_{Ω} \frac{{| M_{ε} (x, t) |}^{2}}{2} d x + \frac{c_{Q}}{2} \iint \iint_{Q_{T}} | \nabla A_{ε} (M_{ε})) |^{2} d x d t \\ + \frac{ε}{2} \iint_{Q_{T}} {| \nabla M_{ε} |}^{2} d x d t \leq ζ_{2}, \end{matrix}$ where ζ₂ is a constant independent of ε and $A (s) = \int_{0}^{s} A (r) d r$ . Therefore, the sequence (M_ε) (resp. A(M_ε)) is bounded in L^∞(Q_T) (resp. L²(0, T; H¹(Ω))) and the sequence $(\sqrt{ε} M_{ε})$ is bounded in L²(0, T; H¹(Ω)). Consequently, $\begin{matrix} | \int_{0}^{T} < \partial_{t} M_{ε}, ψ_{1} > d t | \leq {‖ \nabla A (M_{ε}) ‖}_{L^{2} (Q_{T})} {‖ \nabla ψ_{1} ‖}_{L^{2} (Q_{T})} + {‖ \sqrt{ε} \nabla M_{ε} ‖}_{L^{2} (Q_{T})} {‖ \nabla ψ_{1} ‖}_{L^{2} (Q_{T})} \\ + | | Q (x) χ (M_{ε}) | | L^{\infty} (Q_{T}) | | \nabla C_{ε} | | L^{2} (Q_{T}) | | \nabla ψ_{1} | | L^{2} (Q_{T}) + | | U_{ε} | | L^{2} (Q_{T}) | | \nabla ψ_{1} ‖_{L^{2} (Q_{T})} \\ + μ_{1} {‖ M_{ε} ‖}_{L^{2} (Q_{T})} {‖ ψ_{1} ‖}_{L^{2} (Q_{T})} \leq ζ_{3} {‖ ψ_{1} ‖}_{L^{2} (0, T; H^{1} (Ω))}, \end{matrix}$ where ζ₃ is a constant independent of ε. Thus, as ε → 0, 50 $\frac{\partial M_{ε}}{\partial t} ⇀ \frac{\partial M}{\partial t} in L^{2} (0, T; {(H^{1} (Ω))}^{'}) .$

We have proved till now that (A(M_ε)) is uniformly bounded in L²(0, T; H¹(Ω)) and ∂_tA(M_ε) in L²(0, T; (H¹(Ω))′). Therefore, by compactness (Theorem 2.5.2 in [35]), modulo a subsequence, A(M_ε) → L in L²(Q_T). Since the function A is strictly increasing, there exists M such that L = A(M). Using the continuous function A⁻¹ and applying dominated convergence theorem to M_ε, we can conclude that 51 $M_{ε} \to M in L^{2} (Q_{T}) and a .e . in Q_{T} .$

Similar arguments are used for the W-equation.

Remark 2.

Under all these estimates, one can easily deduce that C_ε → C in L²(0, T; H¹(Ω)) as ε tends to 0.

3.2.2

Passing to the limit

We will only check for the M-equation to show that the limits are solutions of the continuous problem. Let us now tend the regularization parameter ε to 0. Hence, $| \int_{0}^{T} < \partial_{t} (M_{ε}), ψ_{1} > d t - \int_{0}^{T} < \partial_{t} M, ψ_{1} > d t | \to 0, as ε \to 0.$

Since Q(x)∇A_ε(M_ε) ⇀ Q(x)∇A(M) in L²(Q_T) then $| \iint_{Q_{T}} Q (x) \nabla A_{ε} (M_{ε}) \cdot \nabla ψ_{1} d x d t - \iint_{Q_{T}} Q (x) \nabla A (M) \cdot \nabla ψ_{1} d x d t | \to 0, as ε \to 0.$

Due to the Remark 2 and the boundedness of the sensitivity χ₁ and M_ε, one has $| \iint_{Q_{T}} Q (x) χ (M_{ε}) \nabla C_{ε} \cdot \nabla ψ d x d t - \iint_{Q_{T}} Q (x) χ (M) \nabla C \cdot \nabla ψ d x d t | \to 0, as ε \to 0.$

Next, dominated convergence theorem and estimates (46) and (51) imply that $\begin{matrix} | \iint_{Q_{T}} M_{ε} U_{ε} \cdot \nabla ψ d x d t - \iint_{Q_{T}} M U \cdot \nabla ψ d x d t | \\ \leq {‖ (M_{ε} - M) \nabla ψ ‖}_{L^{2} (Q_{T})} {‖ U_{ε} ‖}_{L^{2} (Q_{T})} + {‖ U_{ε} - U ‖}_{L^{2} (Q_{T})} ‖ M \nabla ψ ‖_{L^{2} (Q_{T})} \to 0. \end{matrix}$

Finally, using estimate (51) and similar result for W leads to $\begin{matrix} \iint_{Q_{T}} | M_{ε} (1 - M_{ε} - α_{1} W_{ε}) - M (1 - M - α_{1} W) | ψ d x d t \\ \leq {‖ M_{ε} - M ‖}_{L^{2} (Q_{T})} ‖ ψ ‖_{L^{2} (Q_{T})} + α_{1} {‖ W_{ε} - W ‖}_{L^{2} (Q_{T})} ‖ ψ ‖_{L^{2} (Q_{T})} \to 0. \end{matrix}$

Proof of theorem 2

This section aims to prove the existence and the uniqueness of a weak solution of (3) with κ = 0. Coupling with the Stokes equation also leads to global existence with further regularities on the weak solutions. Indeed, we obtain that $\frac{d U}{d t} \in L^{2} (0, T; V^{'})$ and hence U ∈ C(0, T; H). Moreover, if $U_{0} \in W^{2 - \frac{2}{p}, p} (Ω)$ and ∇ · U₀ = 0 from assumption (13) then $U \in W_{p}^{2, 1} (Q_{T})$ for 1 < p < +∞. Consequently, there is p ≥ 2 large enough for which ∇U ∈ L²(0, T; L^∞(Ω)) [36]. On the other hand, if $C_{0} \in W^{2 - \frac{2}{p}, p} (Ω)$ and under the further regularity of U, we can conclude that there exists p ≥ 2 such that ∇C ∈ L²(0, T; L^∞(Ω)) (Theorem 6 in [37]). Next, a duality technique is used to prove the uniqueness of weak solutions. Denote 𝒩w₁ and $N w_{2} \in H^{2} (Ω) \cap L_{0}^{2} (Ω)$ as unique solutions of 52 ${\begin{matrix} - \nabla \cdot (Q (x) \nabla N w_{1}) & = w_{1} \\ Q (x) \nabla N w_{1} \cdot η & = 0 \end{matrix} and {\begin{matrix} - \nabla \cdot (R (x) \nabla N w_{2}) & = w_{2} \\ R (x) \nabla N w_{2} \cdot η & = 0 \end{matrix} .$

Consider two weak solutions (M₁, W₁, C₁, U₁) and (M₂, W₂, C₂, U₂) of the system (3) (with κ = 0). For T > 0 fixed and for (t, x) ∈ [0, T] × Ω, take $\begin{matrix} M (t, x) = M_{1} (t, x) - M_{2} (t, x), W (t, x) = W_{1} (t, x) - W_{2} (t, x), C (t, x) = C_{1} (t, x) - C_{2} (t, x), \\ U (t, x) = U_{1} (t, x) - U_{2} (t, x) . \end{matrix}$

For the U-equation: We subtract equations related to U₁ and U₂ and we choose ψ = U in the variational equations in order to obtain: $< \partial_{t} U, U > + v \int_{Ω} \nabla U \cdot \nabla U d x = - \int_{Ω} (γ M + λ W) \nabla ϕ \cdot U d x .$

The right-hand side can be written using (52) as: $\begin{matrix} - \int_{Ω} (γ M + λ W) \nabla ϕ \cdot U d x = - γ \int_{Ω} Q (x) \nabla N M \cdot \nabla (\nabla ϕ \cdot U) d x - λ \int_{Ω} R (x) \nabla N W \cdot \nabla (\nabla ϕ \cdot U) d x \\ = - γ \int_{Ω} Q (x) \nabla N M \cdot [(\nabla ϕ \cdot \nabla) U + (U \cdot \nabla) \nabla ϕ + \nabla ϕ \times curl (U) + U \times curl (\nabla ϕ)] d x \\ - λ \int_{Ω} R (x) \nabla N W \cdot [(\nabla ϕ \cdot \nabla) U + (U \cdot \nabla) \nabla ϕ + \nabla ϕ \times curl (U) + U \times curl (\nabla ϕ)] d x . \end{matrix}$

Recall that curl(∇ϕ) = 0, ∥q₁ × q₂∥_L²(Ω) ≤ ∥q₁∥_L²(Ω)∥q₂∥_L²(Ω) and ∃c′ ≥ 0; ∥curl(U)∥_L²(Ω) ≤ c′∥∇U∥_L²(Ω)). Then Poincaré and Young inequalities lead to 53 $\frac{d}{d t} ‖ U (t) ‖_{L^{2} (Ω)}^{2} + v ‖ \nabla U (t) ‖_{L^{2} (Ω)}^{2}$ $\begin{array}{l} \leq (γ ‖ Q (x) ‖_{L^{\infty} (Ω)} ‖ \nabla N M (t) ‖_{L^{2} (Ω)} + λ ‖ R (x) ‖_{L^{\infty} (Ω)} ‖ \nabla N W (t) ‖_{L^{2} (Ω)}) ((1 + c_{P} + c^{'}) ‖ \nabla ϕ ‖_{W^{1, \infty} (Ω)}) ‖ \nabla U (t) ‖_{L^{2} (Ω)} \\ \leq (\frac{1}{δ} c_{1}^{2} (‖ Q (x) ‖_{L^{\infty} (Ω)}^{2} ‖ \nabla N M (t) ‖_{L^{2} (Ω)}^{2} + ‖ R (x) ‖_{L^{\infty} (Ω)}^{2} ‖ \nabla N W (t) ‖_{L^{2} (Ω)}^{2}) ‖ \nabla ϕ ‖_{W^{1, \infty} (Ω)}^{2}) + (γ + λ) δ ‖ \nabla U (t) ‖_{L^{2} (Ω),}^{2} \end{array}$ where c_P, c′ and c₁ are positive constants.

For the C-equation: One has, $\begin{matrix} < \partial_{t} C, C > + \int_{Ω} S (x) \nabla C \cdot \nabla C d x + \int_{Ω} (U \cdot \nabla C_{1}) C d x + \int_{Ω} (U_{2} \cdot \nabla C) C d x \\ = \int_{Ω} (α M_{2} + β W_{2}) (C_{2} - C_{1}) C d x - \int_{Ω} (α M + β W) C_{1} C d x \\ = \int_{Ω} (α M_{2} + β W_{2}) C^{2} d x - \int_{Ω} (α Q (x) \nabla N M + β R (x) \nabla N W) \cdot \nabla (C_{1} C) d x . \end{matrix}$

Due to the assumptions (7), (12) and to Young’s inequality, we obtain 54 $\begin{matrix} \frac{d}{d t} ‖ C (t) ‖_{L^{2} (Ω)}^{2} + c_{S} ‖ \nabla C (t) ‖_{L^{2} (Ω)}^{2} \leq ζ ‖ U (t) ‖_{L^{2} (Ω)} ‖ \nabla C (t) ‖_{L^{2} (Ω)} + (α + β) ‖ C (t) ‖_{L^{2} (Ω)}^{2} \\ + α ζ ‖ Q ‖_{L^{\infty} (Ω)} | | \nabla N M (t) ‖_{L^{2} (Ω)} (‖ \nabla C (t) ‖_{L^{2} (Ω)} + {‖ \nabla C_{1} (t) ‖}_{L^{2} (Ω)} \\ + β ζ ‖ R ‖_{L^{\infty} (Ω)} ‖ \nabla N W (t) ‖_{L^{2} (Ω)} ‖ \nabla C_{1} (t) | | L^{2} (Ω) + | | \nabla C (t) ‖_{L^{2} (Ω)} \\ \leq ζ^{2} δ ‖ \nabla C (t) ‖_{L^{2} (Ω)}^{2} + \frac{1}{δ} ‖ U (t) ‖_{L^{2} (Ω)}^{2} + (α + β) ‖ C (t) ‖_{L^{2} (Ω)}^{2} + \frac{1}{δ} ‖ Q ‖_{L^{\infty} (Ω)}^{2} ‖ \nabla N M (t) ‖_{L^{2} (Ω)}^{2} \\ + \frac{1}{δ} ‖ R ‖_{L^{\infty} (Ω)}^{2} ‖ \nabla N W (t) ‖_{L^{2} (Ω}^{2} + \frac{2}{δ} ‖ C (t) ‖_{L^{2} (Ω)}^{2} + δ (α^{2} + β^{2}) ζ^{2} ‖ \nabla C (t) ‖_{L^{2} (Ω)}^{2} \\ + δ (α^{2} ‖ Q ‖_{L^{\infty} (Ω)}^{2} ‖ \nabla N M (t) ‖_{L^{2} (Ω)}^{2} + β^{2} ‖ R ‖_{L^{\infty} (Ω)}^{2} ‖ \nabla N W (t) ‖_{L^{2} (Ω)}^{2} {‖ \nabla C_{1} ‖}_{L^{\infty} (Ω)}^{2} . \end{matrix}$

For the M-equation: The subtraction of the two variational equalities related to M₁ and M₂ leads to: 55 $\begin{matrix} \int_{0}^{t} < \frac{\partial M}{\partial t}, N M > d s = - \int_{0}^{t} \int_{Ω} Q (x) \nabla (A (M_{1}) - A (M_{2})) \cdot \nabla N M d x d s \\ + \int_{0}^{t} \int_{Ω} Q (x) (χ_{1} (M_{1}) - χ_{1} (M_{2})) \nabla C_{1} \cdot \nabla N M d x d s + \int_{0}^{t} \int_{Ω} Q (x) χ (M_{2}) \nabla C \cdot \nabla N M d x d s \\ + \int_{0}^{t} \int_{Ω} M_{1} U \cdot \nabla N M d x d s + \int_{0}^{t} \int_{Ω} M U_{2} \cdot \nabla N M d x d s + \int_{0}^{t} \int_{Ω} (μ_{1} M (1 - M_{2} - α_{1} W_{2}) N M d x d s \\ + \int_{0}^{t} \int_{Ω} (μ_{1} M_{1} (- M - α_{1} W) N M d x d s . \end{matrix}$

Then, the boundedness of M₁ implies that 56 $\int_{Ω} M_{1} U \cdot \nabla N M d x \leq ‖ \nabla N M (t) ‖_{L^{2} (Ω)} ‖ U (t) ‖_{L^{2} (Ω)} .$

Now, we use the dual problem and the Leibniz rule to deduce that 57 $\begin{matrix} \int_{Ω} M (U_{2} \cdot \nabla N M) d x = \int_{Ω} Q (x) \nabla N M \cdot \nabla (U_{2} \cdot \nabla N M) d x \\ = \int_{Ω} Q (x) \nabla N M \cdot [(U_{2} \cdot \nabla) \nabla N M + (\nabla N M \cdot \nabla) U_{2} + \nabla N M \times curl (U_{2}) + U_{2} \times curl (\nabla N M)] d x . \end{matrix}$

The last integral is zero since curl(∇𝒩M) = 0. Let us bound the first integral: $\int_{Ω} Q (x) \nabla N M \cdot (U_{2} \cdot \nabla) \nabla N M d x = \sum_{i} \sum_{j, k} \int_{Ω} U_{i} Q_{j, k} \frac{\partial}{\partial x_{k}} (N M) \cdot \frac{\partial}{\partial x_{i}} (\frac{\partial}{\partial x_{j}} (N M)) d x .$

For U_i ∈ L^∞(Ω), the C¹-coefficients and the symmetry of the tensor Q, Green Formula and ∇ · U₂ = 0, one obtains: 58 $\begin{matrix} 2 \int_{Ω} Q (x) \nabla N M \cdot (U_{2} \cdot \nabla) \nabla N M d x = - \int_{Ω} (\nabla \cdot U_{2}) Q (x) \nabla N M \cdot \nabla N M d x \\ - \int_{Ω} (U_{2} \cdot \nabla (Q (x))) \nabla N M \cdot \nabla N M d x \leq {‖ U_{2} ‖}_{L^{\infty} (Ω)} ‖ \nabla Q ‖_{L^{\infty} (Ω)} ‖ \nabla N M (t) ‖_{L^{2} (Ω)}^{2} . \end{matrix}$

Again, dual problem (52) and uniform bounds imply that 59 $\int_{Ω} (μ_{1} M (1 - M_{2} - α_{1} W_{2}) N M d x d s \leq μ_{1} \int_{Ω} M N M d x d s \leq μ_{1} ‖ Q ‖_{L^{\infty} (Ω)} ‖ \nabla N M (t) ‖_{L^{2} (Ω)}^{2},$ and, 60 $\int_{Ω} (μ_{1} M_{1} (- M - α_{1} W) N M d x d s \leq - μ_{1} \int_{Ω} M N M d x d s \leq μ_{1} ‖ Q ‖_{L^{\infty} (Ω)} ‖ \nabla N M (t) ‖_{L^{2} (Ω)}^{2} .$

Moreover, we have 61 $\int_{Ω} Q (x) \nabla N M (t) \cdot \nabla N M (t) d x = \int_{Ω} Q (x) \nabla N M (0) \cdot \nabla N M (0) d x + 2 \int_{0}^{t} < \frac{\partial M}{\partial t}, N M > d s .$ because –∇ · (Q(x)∇∂_t(𝒩M)) = ∂_tM in (H¹(Ω))′ and S is symmetric. Recalling (56) till (61), we deduce from the equation (55), Cauchy Schwarz and Young’s inequality that 62 $\begin{matrix} \int_{Ω} Q (x) \nabla N M (t) \cdot \nabla N M (t) d x \leq - 2 \int_{0}^{t} \int_{Ω} (M_{1} - M_{2}) (A (M_{1}) - A (M_{2})) d x d s \\ + 2 δ \int_{0}^{t} \int_{Ω} {(χ_{1} (M_{1}) - χ_{1} (M_{2}))}^{2} d x d s + \frac{2}{δ} \int_{0}^{t} ‖ Q ‖_{L^{\infty} (Ω)}^{2} {‖ \nabla C_{1} ‖}_{L^{\infty} (Ω)}^{2} ‖ \nabla N M (t) ‖_{L^{2} (Ω)}^{2} d s \\ + 2 δ c_{χ_{1}}^{2} \int_{0}^{t} ‖ \nabla C (t) ‖_{L^{2} (Ω)}^{2} d s + \frac{2}{δ} \int_{0}^{t} ‖ Q ‖_{L^{\infty} (Ω)}^{2} ‖ \nabla N M (t) ‖_{L^{2} (Ω)}^{2} d s + 2 δ \int_{0}^{t} ‖ U (t) ‖_{L^{2} (Ω)}^{2} \\ + \frac{2}{δ} \int_{0}^{t} ‖ \nabla N M (t) ‖_{L^{2} (Ω)}^{2} d s + 2 \int_{0}^{t} {‖ U_{2} ‖}_{L^{\infty} (Ω)} ‖ \nabla Q ‖_{L^{\infty} (Ω)} ‖ \nabla N M (t) ‖_{L^{2} (Ω)}^{2} d s \\ + 2 \int_{0}^{t} ‖ Q ‖_{L^{\infty} (Ω)}^{2} ({‖ \nabla U_{2} (t) ‖}_{L^{\infty} (Ω)}^{2} + {‖ \nabla U_{2} (t) ‖}_{L^{2} (Ω)}^{2}) ‖ \nabla N M (t) ‖_{L^{2} (Ω)}^{2} d s + 4 μ_{1} \int_{0}^{t} ‖ Q ‖_{L^{\infty} (Ω)} ‖ \nabla N M (t) ‖_{L^{2} (Ω)}^{2} d s . \end{matrix}$

For the W-equation: Same arguments imply that 63 $\begin{matrix} \int_{Ω} R (x) \nabla N W (t) \cdot \nabla N W (t) d x \leq - 2 \int_{0}^{t} \int_{Ω} (W_{1} - W_{2}) (B (M_{1}) - B (M_{2})) d x d s \\ + 2 δ \int_{0}^{t} \int_{Ω} {(χ_{2} (W_{1}) - χ_{2} (W_{2}))}^{2} d x d s + \frac{2}{δ} \int_{0}^{t} ‖ Q ‖_{L^{\infty} (Ω)}^{2} {‖ \nabla C_{1} ‖}_{L^{\infty} (Ω)}^{2} ‖ \nabla N W (t) ‖_{L^{2} (Ω)}^{2} d s \\ + 2 δ c_{χ_{2}}^{2} \int_{0}^{t} ‖ \nabla C (t) ‖_{L^{2} (Ω)}^{2} d s + \frac{2}{δ} \int_{0}^{t} ‖ R ‖_{L^{\infty} (Ω)}^{2} ‖ \nabla N W (t) ‖_{L^{2} (Ω)}^{2} d s \\ + 2 δ \int_{0}^{t} ‖ U (t) ‖_{L^{2} (Ω)}^{2} + \frac{2}{δ} \int_{0}^{t} ‖ \nabla N W (t) ‖_{L^{2} (Ω)}^{2} d s + 2 \int_{0}^{t} {‖ U_{2} ‖}_{L^{\infty} (Ω)} ‖ \nabla R ‖_{L^{\infty} (Ω)} ‖ \nabla N W (t) ‖_{L^{2} (Ω)}^{2} d s \\ + 2 \int_{0}^{t} ‖ R ‖_{L^{\infty} (Ω)}^{2} ({‖ \nabla U_{2} (t) ‖}_{L^{\infty} (Ω)}^{2} + {‖ \nabla U_{2} (t) ‖}_{L^{2} (Ω)}^{2}) ‖ \nabla N W (t) ‖_{L^{2} (Ω)}^{2} d s + 4 μ_{2} \int_{0}^{t} ‖ Q ‖_{L^{\infty} (Ω)} ‖ \nabla N W (t) ‖_{L^{2} (Ω)}^{2} d s . \end{matrix}$

Then, we consider $0 < δ < \min (\frac{1}{κ_{0}}, \frac{1}{κ_{1}}, \frac{v}{λ + γ}, \frac{c_{M}}{2 c_{χ_{1}}^{2} + 2 c_{χ_{2}}^{2} + ζ^{2} (1 + (α^{2} + β^{2}))}, 1)$ and we add the integrated inequalities with respect to the time for all equations to deduce that: $\begin{matrix} c_{Q} ‖ \nabla N M (t) ‖_{L^{2} (Ω)}^{2} + c_{R} ‖ \nabla N W (t) ‖_{L^{2} (Ω)}^{2} + ‖ C (t) ‖_{L^{2} (Ω)}^{2} + ‖ U (t) ‖_{L^{2} (Ω)}^{2} \\ \leq \int_{0}^{t} (\frac{2}{δ} + α + β) ‖ C (t) ‖_{L^{2} (Ω)}^{2} d s + \int_{0}^{t} (\frac{1}{δ} + 4) ‖ U (t) ‖_{L^{2} (Ω)}^{2} d s \\ + \int_{0}^{t} (\frac{1}{δ} [2 + 3 ‖ Q ‖_{L^{\infty} (Ω)}^{2} + c_{1}^{2} ‖ Q ‖_{L^{\infty} (Ω)}^{2} ‖ \nabla ϕ ‖_{W^{1, \infty} (Ω)}^{2} + 2 ‖ Q ‖_{L^{\infty} (Ω)}^{2} {‖ \nabla C_{1} (t) ‖}_{L^{\infty} (Ω)}^{2}] + 4 μ_{1} ‖ Q ‖_{L^{\infty} (Ω)} \\ + 2 ‖ Q ‖_{L^{\infty}}^{2} ({‖ \nabla U_{2} (t) ‖}_{L^{\infty}}^{2} + {‖ \nabla U_{2} (t) ‖}_{L^{2}}^{2}) + α^{2} ‖ Q ‖_{L^{\infty}}^{2} {‖ \nabla C_{1} (t) ‖}_{L^{2}}^{2} + 2 {‖ U_{2} ‖}_{L^{\infty}} ‖ \nabla Q ‖_{L^{\infty} (Ω)}) ‖ \nabla N M (t) ‖_{L^{2} (Ω)}^{2} d s . \\ + \int_{0}^{t} (\frac{1}{δ} [2 + 3 ‖ R ‖_{L^{\infty} (Ω)}^{2} + c_{1}^{2} ‖ R ‖_{L^{\infty} (Ω)}^{2} ‖ \nabla ϕ ‖_{W^{1, \infty} (Ω)}^{2} + 2 ‖ R ‖_{L^{\infty} (Ω)}^{2} {‖ \nabla C_{1} (t) ‖}_{L^{\infty} (Ω)}^{2}] + 4 μ_{2} ‖ R ‖_{L^{\infty} (Ω)} \\ + 2 ‖ R ‖_{L^{\infty}}^{2} ({‖ \nabla U_{2} (t) ‖}_{L^{\infty}}^{2} + {‖ \nabla U_{2} (t) ‖}_{L^{2}}^{2}) + β^{2} ‖ R ‖_{L^{\infty}}^{2} {‖ \nabla C_{1} (t) ‖}_{L^{2})}^{2} + 2 {‖ U_{2} ‖}_{L^{\infty}} ‖ \nabla R ‖_{L^{\infty} (Ω)}) ‖ \nabla N W (t) ‖_{L^{2} (Ω)}^{2} d s . \end{matrix}$

Consequently, $\begin{matrix} ‖ \nabla N M (t) ‖_{L^{2} (Ω)}^{2} + ‖ \nabla N W (t) ‖_{L^{2} (Ω)}^{2} + ‖ C (t) ‖_{L^{2} (Ω)}^{2} + ‖ U (t) ‖_{L^{2} (Ω)}^{2} \\ \leq \int_{0}^{t} μ (s) [‖ \nabla N M ‖_{L^{2} (Ω)}^{2} + ‖ \nabla N W ‖_{L^{2} (Ω)}^{2} + ‖ C ‖_{L^{2} (Ω)}^{2} + ‖ U ‖_{L^{2} (Ω)}^{2}] d s, \end{matrix}$ where μ(s) is a positive integrable function. The proof of Theorem 2 is well achieved since U(t) = C(t) = ∇𝒩M(t) = ∇𝒩W(t) = 0 for every t ∈ [0, T] due to Gronwall Lemma.

Numerical scheme

5.1

Overview and motivation

Numerical modeling of interspecies competition with chemotactic behavior poses significant challenges when diffusion is anisotropic. In the isotropic setting, finite volume (FV) schemes have been successfully applied, as in [38], to simulate competition in quiescent fluids. However, it is well known that classical two-point flux approximations become inconsistent in the presence of anisotropic diffusion.

To accurately capture the competitive anisotropic dynamics of two interacting species M and W, we adopt a robust and convergent numerical scheme initially introduced for a one-species anisotropic Keller-Segel model and further developed for a chemotaxis-fluid coupling without logistic terms. In the present work, we extend this methodology to a two-species chemotaxis-fluid system incorporating Lotka-Volterra-type interspecies competition.

The numerical approach, implemented on Fortran, combines:

A finite volume method for discretizing convective fluxes for species densities and chemoattractant concentrations.
A finite element method for discretizing the diffusive fluxes for species densities and for solving the incompressible Navier-Stokes equations governing fluid flow.

The underlying computational domain is discretized, as seen in Figure 2, using a triangular primal mesh 𝒯_h, from which a dual mesh 𝒟_h is constructed, with diamond-shaped control volumes centered on each interior edge.

5.2

Discrete variables and notation

Let τ > 0 denote the time step, and define discrete time levels t_n = nτ for $n = 0, 1, \dots, Ñ$ . The species densities M and W, and the chemoattractant concentration C, are approximated as piecewise constant functions on the dual mesh. Initial conditions for each variable are defined by cell averages over each dual control volume at each time step. The fluid velocity U and pressure P are approximated using piecewise linear and piecewise constant functions, respectively, defined over the primal mesh 𝒯_h.

Initial conditions for each variable are defined by cell averages over each dual control volume D ∈ 𝒟_h as follows: 64 $M_{D}^{0} = \frac{1}{| D |} \int_{D} M_{0} (x) d x, W_{D}^{0} = \frac{1}{| D |} \int_{D} W_{0} (x) d x, C_{D}^{0} = \frac{1}{| D |} \int_{D} C_{0} (x) d x .$

5.3

Numerical algorithm

The time-stepping procedure proceeds in two main stages:

First step: Velocity-pressure coupling (Navier-Stokes solver)

Given the densities ${\tilde{M}}_{h}^{n} = {M_{D}^{n}}_{D \in D_{h}}$ and ${\tilde{W}}_{h}^{n} = {W_{D}^{n}}_{D \in D_{h}}$ , and the previous time-step velocity field $U_{h}^{n}$ we compute the updated fluid velocity $U_{h}^{n + 1} \in X_{h}$ and pressure $P_{h}^{n + 1} \in Y_{h}$ , where:

X_h is the space of piecewise linear functions on 𝒯_h, continuous at interior edge midpoints.
Y_h is the space of piecewise constant pressures.

The system is solved using a variational formulation with linear constraints, and the velocity-pressure pair is obtained via a classical Uzawa iterative algorithm.

Second step: Scalar field updates (Finite volume scheme)

With $U_{h}^{n + 1}$ computed, we update the discrete densities and chemoattractant concentrations by solving the following finite volume equations for each diamond control volume D ∈ 𝒟_h:

Species 1 (M): 65 $\begin{matrix} | D | \frac{M_{D}^{n + 1} - M_{D}^{n}}{Δ t} - \sum_{E \in D_{h}} Q_{D, E} A (M_{E}^{n + 1}) + \sum_{E \in N (D)} G (M_{D}^{n + 1}, M_{E}^{n + 1}; δ C_{D, E}^{n + 1}) \\ + \sum_{E \in N (D)} G_{1} (M_{D}^{n + 1}, M_{E}^{n + 1}; U_{D, E}^{n + 1}) = μ_{1} M_{D}^{n + 1} (1 - M_{D}^{n + 1} - α_{1} W_{D}^{n}), \end{matrix}$
Species 2 (W): 66 $\begin{matrix} | D | \frac{W_{D}^{n + 1} - W_{D}^{n}}{Δ t} - \sum_{E \in D_{h}} ℛ_{D, E} B (W_{E}^{n + 1}) + \sum_{E \in N (D)} G (W_{D}^{n + 1}, W_{E}^{n + 1}; δ C_{D, E}^{n + 1}) \\ + \sum_{E \in N (D)} G_{1} (W_{D}^{n + 1}, W_{E}^{n + 1}; U_{D, E}^{n + 1}) = μ_{2} W_{D}^{n + 1} (1 - W_{D}^{n + 1} - α_{2} M_{D}^{n}), \end{matrix}$
Chemoattractant (C): 67 $| D | \frac{C_{D}^{n + 1} - C_{D}^{n}}{Δ t} - \sum_{E \in D_{h}} S_{D, E} C_{E}^{n + 1} + \sum_{E \in N (D)} G_{1} (C_{D}^{n + 1}, C_{E}^{n + 1}; U_{D, E}^{n + 1}) = - (α M_{D}^{n + 1} + β W_{D}^{n + 1}) C_{D}^{n + 1} .$

Here:

𝒬_{D, E}, ℛ_{D, E} and 𝒮_D,E are stiffness matrix entries from the finite element method.
$Q_{D, E} = - \sum_{K \in T_{h}} {(Q (x) \nabla φ_{E}, \nabla φ_{D})}_{0, K}; | Q_{D, E} | \leq \frac{c_{Q}}{k_{T}} \frac{(d i a m {(K_{D, E})}^{d - 2}}{{(d - 1)}^{2}}, \forall D \in D_{h}, E \in N (D)$ .
$δ C_{D, E}^{n + 1} = Q_{D, E} (C_{E}^{n + 1} - C_{D}^{n + 1})$ represents chemoattractant gradients.
Convective terms G and G₁ are constructed to be conservative, consistent, and coercive, with G₁ defined using an upwind flux. It can be written as $G_{1} (M_{D}^{n + 1}, M_{E}^{n + 1}, U_{D, E}^{n + 1}) = {(U_{D, E}^{n + 1})}^{+} M_{D}^{n + 1} - {(U_{D, E}^{n + 1})}^{-} M_{E}^{n + 1} = - {(U_{D, E}^{n + 1})}^{-} (M_{E}^{n + 1} - M_{D}^{n + 1})$ such that $U_{D, E}^{n + 1} = \int_{σ_{D, E}} U_{h}^{n + 1} \cdot η_{D, E} d γ$ .

We assume that all transmissibilities are positive; if not, a correction of the diffusive fluxes will be needed to rejoin the positive case.

5.4

Convergence analysis

We consider two types of discrete approximations:

A finite volume solution $({\tilde{M}}_{h, Δ t} (x, t), {\tilde{W}}_{h, Δ t} (x, t), {\tilde{C}}_{h, Δ t} (x, t)) = (M_{D}^{n + 1}, W_{D}^{n + 1}, C_{D}^{n + 1})$ , piecewise constant in space and time.
A nonconforming finite element solution (M_h,Δt, W_h,Δt, C_h,Δt, U_h,Δt), piecewise linear in space and piecewise constant in time.

These approximations satisfy the following maximum principle and boundedness properties: $0 \leq M_{D}^{k}, W_{D}^{k} \leq 1, 0 \leq C_{D}^{k} \leq ζ, {‖ U_{h}^{k} ‖}_{L^{\infty} (Ω)} \leq C_{1}, \forall D \in D_{h}, \forall k \in {0, 1, \dots, Ñ} .$

5.5

Convergence of the full scheme

We summarize the main outcomes in the following convergence result inspired from [39] for the fluid part and [40] for the chemotaxis part. Let the model assumptions hold, including (11), and assume that the regularization parameter satisfies $0 < ρ < \frac{2 v}{d}$ . Then:

1)
Velocity-Pressure solver:
- If d = 2, the discrete Navier-Stokes problem admits a unique velocity field U_h,Δt.
- For d ≤ 4, at least one solution exists.
- As h, Δt → 0, the velocity converges: $U_{h, Δ t} \to U strongly in L^{2} (Q_{T}), U_{h, Δ t} \overset{*}{⇀} U in L^{\infty} (0, T; L^{2} (Ω)) .$
2)
Chemotaxis-Fluid system:
- There exists a solution $({\tilde{M}}_{h, Δ t}, {\tilde{W}}_{h, Δ t}, {\tilde{C}}_{h, Δ t})$ to the discrete system (65)-(67).
- For any sequence h_m → 0, there exists a subsequence such that (M_{h_m}, W_{h_m}, C_{h_m}, U_{h_m}) converges a.e. on Q_T to a solution (M, W, C, U), which is a weak solution of the continuous problem (3) as in Definition 1.

By establishing discrete gradient estimates on $A ({\tilde{M}}_{h, Δ t})$ and $B ({\tilde{W}}_{h, Δ t})$ , and using Kolmogorov’s compactness theorem, we prove strong compactness of the sequences in L²(Q_T). Due to the monotonicity of A and B, modulo subsequences, $(A ({\tilde{M}}_{h, Δ t}))$ and $(B ({\tilde{W}}_{h, Δ t}))$ converge to A(M) and B(W). Then, space translate estimates imply that the limits A(M) and B(W) are in L²(0, T, H¹(Ω)). We can now apply the dominated convergence theorem to ${\tilde{M}}_{h, Δ t} = A^{- 1} (A ({\tilde{M}}_{h, Δ t}))$ and to ${\tilde{W}}_{h, Δ t} = B^{- 1} (B ({\tilde{W}}_{h, Δ t}))$ since the operators A and B are invertible and continuous and the discrete solutions ${\tilde{M}}_{h, Δ t}$ and ${\tilde{W}}_{h, Δ t}$ are bounded in L^∞(Q_T). Consequently, there exist a subsequence of ${\tilde{M}}_{h, Δ t} (resp . {\tilde{W}}_{h, Δ t})$ that converges strongly in L²(Q_T) and a.e. in Q_T to the same function M (resp. W).

Numerical simulations

In this section, we present a series of numerical experiments designed to explore the anisotropic dynamics and competitive interactions of two biological species, M and W, within a fluid environment. These simulations aim to provide deeper insight into how flow-driven organisms interact, compete, and potentially collaborate depending on their environment and intrinsic properties. To accomplish this, we examine the following system: 68 ${\begin{matrix} \partial_{t} M - \nabla \cdot (Q (x) a (M) \nabla M) + \nabla \cdot (Q (x) χ_{1} (M) \nabla C) + c_{1} (U \cdot \nabla M) & = & μ_{1} M (1 - M - α_{1} W), \\ \partial_{t} W - \nabla \cdot (R (x) b (W) \nabla W) + \nabla \cdot (R (x) χ_{2} (W) \nabla C) + c_{2} (U \cdot \nabla W) & = & μ_{2} W (1 - W - α_{2} M), \\ \partial_{t} C - {\tilde{D}}_{3} \nabla \cdot (S (x) \nabla C) + c_{3} (U \cdot \nabla C) & = & - (α M + β W) C, \\ \partial_{t} U - ν ∆ U + κ (U \cdot \nabla) U + \nabla P & = & - (γ M + λ W) \nabla ϕ, \\ \nabla \cdot U & = & 0, \end{matrix}$ \[ \left\{ \begin{aligned} \partial_t M -\nabla\cdot\bigl(Q(x)a(M)\nabla M\bigr) +\nabla\cdot\bigl(Q(x)\chi_1(M)\nabla C\bigr) +c_1(U\cdot\nabla M) &=\mu_1M(1-M-\alpha_1W),\\[3pt] \partial_t W -\nabla\cdot\bigl(R(x)b(W)\nabla W\bigr) +\nabla\cdot\bigl(R(x)\chi_2(W)\nabla C\bigr) +c_2(U\cdot\nabla W) &=\mu_2W(1-W-\alpha_2M),\\[3pt] \partial_t C -\widetilde{D}_3\nabla\cdot\bigl(S(x)\nabla C\bigr) +c_3(U\cdot\nabla C) &=-(\alpha M+\beta W)C,\\[3pt] \partial_t U -\nu\Delta U +\kappa(U\cdot\nabla)U +\nabla P &=-(\gamma M+\lambda W)\nabla\phi,\\[3pt] \nabla\cdot U &= 0. \end{aligned} \right. \] with the following functions:

Degenerate diffusions: $a (M) = {\tilde{D}}_{1} M (1 - M)$ and $b (W) = {\tilde{D}}_{2} W (1 - W)$ ,
Sensitivities: χ₁(M) = c₄M(1 – M) and χ₂(W) = c₅W(1 – W).

This system is numerically approximated using the robust scheme detailed in Section 5. In addition to the qualitative validation of the model, the convergence behaviour of the proposed scheme is quantitatively assessed through a mesh refinement study. Numerical solutions are computed on two successively refined discretizations (h, Δt) and (h/2, Δt/2). Let $U_{h}^{n}$ and $U_{h / 2}^{n}$ denote the discrete solutions at time t_n. The relative L²-error between both approximations is defined as $E_{h} (t_{n}) = \frac{{‖ U_{h / 2}^{n} - U_{h}^{n} ‖}_{L^{2} (Ω)}}{{‖ U_{h / 2}^{n} ‖}_{L^{2} (Ω)}}$ . The experimental order of convergence (EOC) is then estimated by $p \approx \frac{\log (E_{h} / E_{h / 2})}{\log (2)}$ . The obtained results indicate a consistent first-order convergence rate of the scheme, i.e. 𝒪(h, Δt) ≈ 1.

6.1

Test 1: Classical Lotka-Volterra kinetics (No Fluid, No Diffusion)

This benchmark test isolates the intrinsic competitive behavior of two species by removing all spatial effects: no diffusion, chemotaxis, or fluid flow. We solve a spatially uniform version of the Lotka-Volterra model using:

Domain: Mesh 1 of Figure 4, Δt = 0.0005, μ₁ = μ₂ = 1,
Initial densities: M₀ = W₀ = 0.5 localized in the square ([0.45,0.55] × [0.45,0.55]),
Three scenarios: α₁ = α₂ = 0.5 < 1 (moderate competition), α₁ = 0.5 < 1 < α₂ = 1.5 (species 1 dominates) and α₂ = 0.5 < 1 < α₁ = 1.5 (species 2 dominates).

Figure 3 illustrates the outcomes of these simulations in the weak and strong situations. In the first scenario, we observe a coexistence of two species with densities asymptotically approaching $N_{3} (\frac{2}{3}, \frac{2}{3})$ , aligning with theoretical predictions. Conversely, in the second case, species 1 dominates the competition, leading to the exclusion of species 2, and the densities (M, W) converge towards (1,0). Similarly, in the third scenario, species 2 emerges victorious in the competition and species 1 is excluded. These results validate the critical dynamics of species interactions and their dependence on competition coefficients.

6.2

Test 2: Full model with strong competition in a fluid domain

Our extended code based on a robust numerical scheme show numerically through this simulation the exponential convergence of species’ density in the competition regime. We now activate all components of the system: diffusion, chemotaxis, and fluid flow, to investigate the full dynamics under strong competition:

Domain: Mesh 1 of Figure 4,
Initial conditions: Random densities M₀(x, y), W₀(x, y) in [0, 1] and C₀(x, y) = 1 in specific regions,
Fluid : Driven cavity flow; top wall velocity set to (1,0), zero elsewhere,
Δt = 0.0005, μ₁ = μ₂ = 1, α = 0.06, β = 0.08, γ = 0.01, λ = 0.01,
c₁ = 0.001, c₂ = 0.001, c₃ = 0, ${\tilde{D}}_{1} = {\tilde{D}}_{2} = 0.001, {\tilde{D}}_{3} = 10^{- 5}$ ,
c₄ = c₅ = 0.1, v = 0.001, κ = 1, ∇ϕ = (0.1;0.1), α₁ = 4 > 1 > α₂ = 0.01,
Diffusion tensors: Species 1 uses a homogeneous anisotropic tensor $Q = [\begin{array}{l} 2.2 & 0.8 \\ 0.8 & 1.8 \end{array}]$ , whereas species 2 employs a heterogeneous rotational anisotropic tensor $R (x, y) = \frac{1}{x^{2} + y^{2}} [\begin{matrix} y^{2} + 0.01 x^{2} & - (1 - 0.01) x y \\ - (1 - 0.01) x y & x^{2} + 0.01 y^{2} \end{matrix}]$ .

Figures 5 and 6 illustrate the behavior of the species within the fluid. After 20 seconds, Species 2 rapidly dominates due to its low competition coefficient and efficient chemotactic strategy. The solution (M, W) exponentially converges toward (0, 1), showing the robustness of competitive advantage under strong interspecies differences. Table 2 reports convergence errors supporting this exponential decay.

Table 2

Relative errors-Mesh size =h= 0.54 × 10⁻².

	M(t = 1)	W(t = 1)	M(t = 3)	W(t = 3)	M(t = 5)	W(t = 5)
Relative L²-error	3.09 × 10⁻⁶	3.36 × 10⁻⁷	1.21 × 10⁻⁷	1.16 × 10⁻⁷	6.43 × 10⁻⁸	6.80 × 10⁻⁸
Relative L^∞-error	1.76 × 10⁻³	5.80 × 10⁻⁴	3.48 × 10⁻⁴	3.41 × 10⁻⁴	2.54 × 10⁻⁴	2.61 × 10⁻⁴

6.3

Test 3: Coexistence in a domain with obstacles

To examine coexistence under more realistic geometries, we introduce a domain with an internal obstacle. For that, we set:

Domain: Mesh 2 of Figure 4,
Initial conditions: M₀ = W₀ = 0.2 and C₀(x, y) = 1 in some regions,
Fluid: Monodirectional, Δt = 0.0005, μ₁ = μ₂ = 1, α = 0.04, β = 0.08, γ = 0.01, λ = 0.01,
c₁ = c₂ = 0.1, c₃ = 0, ${\tilde{D}}_{1} = {\tilde{D}}_{2} = 0.1, {\tilde{D}}_{3} = 0.001$ , c₄ = c₅ = 0.1
v = 0:1, κ = 1, ∇ϕ = (0.1, 0.1), α₁ = α₂ = 0.01 < 1,
Diffusion tensors: $Q = R = [\begin{matrix} 8 & - 7 \\ - 7 & 20 \end{matrix}]$ .

Despite the complex geometry and internal barrier, Figure 7 shows that both species manage to distribute and survive by securing nutrients from different areas.

6.4

Test 4: Fluid influence on competitive dynamics

This test compares species dynamics with and without fluid motion, using a vertically structured initial configuration. For that, we set:

Domain: Mesh 3 of Figure 4,
Initial conditions: M₀(x, y) = 0.05 in [0.45,0.55] × [0.25,0.35], W₀(x, y) = 0.05 in [0.45, 0.55]×[0.65, 0.75] and C₀(x, y) = 5 in [0.2, 0.3]×[0.45, 0.55]⋃[0.7, 0.8]×[0.45, 0.55],
Fluid : Upwards monodirectional with U₀ = (0,5), Δt = 0.0005, μ₁ = μ₂ = 1,
α = 0.01, β = 0.05, γ = 0.01, λ = 0.01, c₁ = 20, c₂ = 20, c₃ = 0, ${\tilde{D}}_{1} = {\tilde{D}}_{2} = 0.001, {\tilde{D}}_{3} = 0.001$ ,
c₄ = c₅ = 0.1, v = 0.01, κ = 1, ∇ϕ = (0, 100), α₁ = α₂ = 0:01 < 1,
Diffusion tensors: Species 1 uses a homogeneous anisotropic tensor $Q = R = [\begin{matrix} 8 & - 7 \\ - 7 & 20 \end{matrix}]$ .

If the fluid is at rest, Figure 8 shows the dynamics of the species that move slowly and accumulate near a single chemical region due to the anisotropic diffusion. With fluid, species are transported and spread, reaching both chemical-rich zones. Fluid enhances nutrient access and creates directional migration patterns. Despite weak competition, fluid flow substantially changes spatial strategies in Figure 9.

6.5

Test 5: Impact of chemotactic sensitivity on cooperation

We simulate how chemotactic cues and anisotropic diffusion can transform competitive dynamics into collaboration to ensure nutrients or coexistence to avoid predators (modeled via chemo-repulsion). We consider an oblique fluid flow such that U₀ = (0.1, 0.1) with negligible pressure P₀ = 0. This numerical experiment is conducted on the mesh 3 given in Figure 4, with the same parameters as Test 3, is divided into two scenarios:

First scenario: (c₄ = c₅ = –0.15, chemo-repellent)

The initial conditions are M₀(x, y) = W₀(x, y) = 0.2 in [0.45,0.55] × [0.75,0.85] and C₀(x, y) = 5 in [0.45,0.55] × [0.45,0.55]. In Figure 10, species 1 follows a heterogeneous isotropic tensor Q = e^–(x+y)Id, whereas species 2 employs a heterogeneous rotational anisotropic tensor $R (x, y) = \frac{1}{x^{2} + y^{2}} [\begin{matrix} y^{2} + 0.02 x^{2} & - (1 - 0.02) x y \\ - (1 - 0.02) x y & x^{2} + 0.02 y^{2} \end{matrix}] .$

Second scenario: (c₄ = c₅ = 0.15, chemo-attractant)

The initial conditions are M₀ = W₀ = 0.2 in [0.45,0.55] × [0.45,0.55] and C₀(x, y) = 5 in four regions. We consider the anisotropic tensors $Q = R = [\begin{array}{l} 5.5 & 3.2 \\ 3.2 & 3.5 \end{array}]$ .

Starting from the first five seconds of the simulation, many remarkable changes occur in the unidirectional fluid flow. In all scenarios, the species coexist, yet with distinct strategies: in the second one, species in Figure 11 coexist while competing for nutrients, distributing themselves across multiple chemical sources, while in the first one, species in Figure 10 actively collaborate to escape predator zones, guided by repellent gradients. Moreover, we have chosen the first scenario to check the time-evolution of the relative L∞ and L² errors given in Table 3, confirming numerical accuracy. Finally, we note that the velocity field gradually diminishes towards zero, with a decline in chemical concentration levels.

Table 3

Relative errors-Mesh size =h= 0.34 × 10⁻³.

	M(t = 1)	W(t = 1)	M(t = 3)	W(t = 3)	M(t = 5)	W(t = 5)
Relative L²-error	1.82 × 10⁻⁶	1.48 × 10⁻⁷	9.70 × 10⁻⁸	1.39 × 10⁻⁷	6.09 × 10⁻⁸	1.36 × 10⁻⁷
Relative L^∞-error	1.35 × 10⁻³	3.85 × 10⁻⁴	3.11 × 10⁻⁴	3.72 × 10⁻⁴	2.47 × 10⁻⁴	3.69 × 10⁻⁴

Contrary to the anisotropic configurations considered in the previous tests, the isotropic case corresponds to uniform diffusion in all spatial directions, as the diffusion tensors reduce to scalar multiples of the identity matrix. In this configuration, the dispersal mechanism becomes direction-independent, leading to a spatially homogeneous smoothing of the species densities.

Conclusion

This work addressed a nonlinear anisotropic chemotaxis-competition model for two interacting species evolving in a fluid environment, incorporating degenerate density-dependent diffusion and Stokes coupling. From the theoretical perspective, we established the global existence of weak solutions in spatial dimensions d ≤ 3, along with a uniqueness result in the Stokes regime under additional regularity assumptions. On the numerical side, a generalized hybrid finite volume-finite element scheme was developed to discretize the fully coupled system on general meshes while preserving confinement properties of cell densities. Numerical experiments highlighted the impact of heterogeneous anisotropic diffusion on spatial segregation dynamics. Future research directions include the extension of the analysis to Navier-Stokes coupling, the incorporation of nonlinear Holling-type functional responses, and the study of optimal control strategies in an anisotropic tumor microenvironments or ecological competitive systems. The proposed numerical framework can also be adapted to multi-species chemotaxis models and porous media transport mechanisms.

Declarations

Directional migration and competition in fluid media: Global existence, uniqueness, and robust simulation of chemotaxis and tumor growth dynamics

Full Article

Paradigm

My account