Numerical analysis for a chemotaxis-Navier-Stokes system

TL;DR

This paper develops and analyzes a finite element numerical scheme for simulating a chemotaxis-Navier-Stokes system, capturing cellular movement in fluids with error estimates and convergence proofs.

Contribution

It introduces a novel splitting-based finite element scheme for the chemotaxis-Navier-Stokes system with proven convergence and error estimates.

Findings

The scheme is well-posed and mass-conservative.

Numerical simulations confirm the scheme's accuracy and stability.

Error estimates align with theoretical predictions.

Abstract

In this paper we develop a numerical scheme for approximating a $d$-dimensional chemotaxis-Navier-Stokes system, $d=2,3$, modeling cellular swimming in incompressible fluids. This model describes the chemotaxis-fluid interaction in cases where the chemical signal is consumed with a rate proportional to the amount of organisms. We construct numerical approximations based on the Finite Element method and analyze some error estimates and convergence towards weak solutions. In order to construct the numerical scheme, we use a splitting technique to deal with the chemo-attraction term in the cell-density equation, leading to introduce a new variable given by the gradient of the chemical concentration. Having the equivalent model, we consider a fully discrete Finite Element approximation which is well-posed and it is mass-conservative. We obtain uniform estimates and analyze the convergence of the scheme. Finally, we present some numerical simulations to verify the good behavior of our scheme, as well as to check numerically the error estimates proved in our theoretical analysis.

Equations

$$\left\{\begin{array}[]{lc}\eta_{t}+\mathbf{u}\cdot\nabla\eta=D_{n}\Delta\eta-\chi\nabla\cdot(\eta\nabla c),&\\ c_{t}+\mathbf{u}\cdot\nabla c=D_{c}\Delta c-\gamma\eta c,&\\ \rho\left(\mathbf{u}_{t}+(\mathbf{u}\cdot\nabla)\mathbf{u}\right)=D_{\mathbf{u}}\Delta\mathbf{u}-\nabla\pi+\eta\nabla\phi,&\\ \nabla\cdot\mathbf{u}=0,&\end{array}\right.$$

$$\left\{\begin{array}[]{lc}\left[\eta(\textit{{x}},0),c(\textit{{x}},0),\mathbf{u}(\textit{{x}},0)\right]=\left[\eta_{0}(\textit{{x}}),c_{0}(\textit{{x}}),\mathbf{u}_{0}(\textit{{x}})\right],\ \textit{{x}}\in\Omega,\\[8.5359pt] \frac{\partial\eta(\textit{{x}},t)}{\partial\bm{\nu}}=\frac{\partial c(\textit{{x}},t)}{\partial\bm{\nu}}=0,\quad\mathbf{u}(\textit{{x}},t)=0,\quad\textit{{x}}\in\partial\Omega,\quad t\in(0,T).\end{array}\right.$$

$$\mathbf{V}$$

$${\mathbf{H}}_{\sigma }^1(\Omega )$$

$$L_0^2(\Omega )$$

$${\hat{H}}^1(\Omega )$$

$$\Vert u{\Vert }_{H^1}^2=\Vert \nabla u{\Vert }_0^2+{\left({\int }_{\Omega }u\right)}^2,\forall u\in H^1(\Omega ),$$

$$\Vert \mathbf{\sigma }{\Vert }_{H^1}^2=\Vert \mathbf{\sigma }{\Vert }_0^2+\Vert \text{mbox}rot\mathbf{\sigma }{\Vert }_0^2+\Vert \nabla \cdot \mathbf{\sigma }{\Vert }_0^2,\forall \mathbf{\sigma }\in {\mathbf{H}}_{\sigma }^1(\Omega ).$$

$$\Vert \mathbf{v}{\Vert }_{H^1}\le C_P\Vert \nabla \mathbf{v}{\Vert }_{L^2},\forall \mathbf{v}\in {\mathbf{H}}_0^1(\Omega ),$$

$$\Vert u{\Vert }_{L^4}\le C\Vert u{\Vert }_{L^2}^{1/2}\Vert u{\Vert }_{H^1}^{1/2}\forall u\in H^1(\Omega )\text{mbox}(in2Ddomains),$$

$$\Vert u{\Vert }_{L^3}\le C\Vert u{\Vert }_{L^2}^{1/2}\Vert u{\Vert }_{L^6}^{1/2}\forall u\in H^1(\Omega )\text{mbox}(in3Ddomains).$$

$$\left\{\begin{array}[]{lc}\eta_{0}\in C^{0}(\overline{\Omega}),\ \ n_{0}>0\ \ \text{in }\>\overline{\Omega},&\\ c_{0}\in W^{1,q}(\Omega)\>\text{ for some }q>2,\ \ c_{0}>0\ \ \text{in }\>\overline{\Omega},&\\ \mathbf{u}_{0}\in D(A^{\alpha})\ \ \text{ for some }\alpha\in(\frac{1}{2},1),\nabla\phi\in C^{1}(\overline{\Omega}),&\end{array}\right.$$

$$\left\{\begin{array}[]{lc}\eta\in C([0,T);L^{2}(\Omega))\cap L^{\infty}((0,T);C(\bar{\Omega}))\cap C^{2,1}(\bar{\Omega}\times(0,T)),\\ c\in C([0,T);L^{2}(\Omega))\cap L^{\infty}((0,T);W^{1,q}({\Omega}))\cap C^{2,1}(\bar{\Omega}\times(0,T)),\\ \mathbf{u}\in C([0,T);L^{2}(\Omega))\cap L^{\infty}((0,T);D(A^{\alpha}))\cap C^{2,1}(\bar{\Omega}\times(0,T)),\\ \pi\in L^{1}((0,T);W^{1,2}).\end{array}\right.$$

$$-{\int }_0^T{\int }_{\Omega }\eta {\phi }_t-{\int }_{\Omega }{\eta }_0\phi (\cdot ,0)-{\int }_0^T{\int }_{\Omega }\eta \mathbf{u}\cdot \nabla \phi =-{\int }_0^T{\int }_{\Omega }\nabla \eta \cdot \nabla \phi +\chi {\int }_0^T{\int }_{\Omega }\eta \nabla c\cdot \phi ,$$

$$-{\int }_0^T{\int }_{\Omega }c{\phi }_t-{\int }_{\Omega }{c}_0\phi (\cdot ,0)-{\int }_0^T{\int }_{\Omega }c\mathbf{u}\cdot \nabla \phi +{\int }_0^T{\int }_{\Omega }\nabla c\cdot \nabla \phi =-\gamma {\int }_0^T{\int }_{\Omega }\eta c\phi ,$$

$$-{\int }_0^T{\int }_{\Omega }\mathbf{u}{\psi }_t-{\int }_{\Omega }{\mathbf{u}}_0\psi (\cdot ,0)-{\int }_0^T{\int }_{\Omega }\mathbf{u}\otimes \mathbf{u}\cdot \nabla \psi +{\int }_0^T{\int }_{\Omega }\nabla \mathbf{u}\cdot \nabla \psi =-{\int }_0^T{\int }_{\Omega }\eta \nabla \varphi \cdot \nabla \psi ,$$

$$\left\{\begin{array}[]{lc}\eta_{0}\in L\log L(\Omega)\ \mbox{is nonnegative with}\ \eta_{0}\neq 0,\\ c_{0}\in L^{\infty}(\Omega)\ \mbox{is nonnegative with}\ \sqrt{c_{0}}\in W^{1,2}(\Omega),\\ {\mathbf{u}}_{0}\in L^{2}(\Omega),\ \nabla\cdot\mathbf{u}_{0}=0,\end{array}\right.$$

$${\int }_{\Omega }\eta (\cdot ,t)d\text{x}={\int }_{\Omega }{\eta }_{0}d\text{x}.$$

$$\left\{\begin{array}[]{lc}\vspace{0.1 cm}(n_{t},\bar{n})+D_{n}(\nabla n,\nabla\bar{n})+(\mathbf{u}\cdot\nabla n,\bar{n})=\chi((n+\alpha_{0})\bm{\sigma},\nabla\bar{n}),\quad\forall\bar{n}\in\widehat{H}^{1}(\Omega),&\\ (\bm{\sigma},\bar{\bm{\sigma}})+(c,\nabla\cdot\bar{\bm{\sigma}})=0,\quad\forall\bar{\bm{\sigma}}\in\mathbf{H}_{\sigma}^{1}(\Omega),&\\ \displaystyle(c_{t},\bar{c})-D_{c}(\nabla\cdot\bm{\sigma},\bar{c})+(\mathbf{u}\cdot\bm{\sigma},\bar{c})=-\gamma((n+\alpha_{0})c,\bar{c}),\quad\forall\bar{c}\in L^{2}(\Omega),&\\ \rho(\mathbf{u}_{t},\bar{\mathbf{u}})+D_{\mathbf{u}}(\nabla\mathbf{u},\nabla\bar{\mathbf{u}})+\rho((\mathbf{u}\cdot\nabla)\mathbf{u},\bar{\mathbf{u}})-(\pi,\nabla\cdot\bar{\mathbf{u}})=((n+\alpha_{0})\nabla\phi,\bar{\mathbf{u}}),\ \ \forall\bar{\mathbf{u}}\in\mathbf{H}_{0}^{1}(\Omega),&\\ (\bar{\pi},\nabla\cdot{\mathbf{u}})=0,\quad\forall\bar{\pi}\in{L}_{0}^{2}(\Omega).&\end{array}\right.$$

$$({\mathbf{\sigma }}_t,\bar{\mathbf{\sigma }})+(c_t,\nabla \cdot \bar{\mathbf{\sigma }})=0,\forall \bar{\mathbf{\sigma }}\in {\mathbf{H}}_{\sigma }^1(\Omega ).$$

$$\left\{\begin{array}[]{lc}\vspace{0.1 cm}(n_{t},\bar{n})+D_{n}(\nabla n,\nabla\bar{n})+(\mathbf{u}\cdot\nabla n,\bar{n})=\chi((n+\alpha_{0})\bm{\sigma},\nabla\bar{n}),&\\ (\bm{\sigma}_{t},\bar{\bm{\sigma}})+D_{c}(\nabla\cdot\bm{\sigma},\nabla\cdot\bar{\bm{\sigma}})+D_{c}(\mbox{rot }{\bm{\sigma}},\mbox{rot }\bar{\bm{\sigma}})-(\mathbf{u}\cdot\bm{\sigma},\nabla\cdot\bar{\bm{\sigma}})=\gamma((n+\alpha_{0})c,\nabla\cdot\bar{\bm{\sigma}}),&\\ (c_{t},\bar{c})+D_{c}(\nabla c,\nabla\bar{c})+(\mathbf{u}\cdot\nabla c,\bar{c})=-\gamma((n+\alpha_{0})c,\bar{c}),&\\ \rho(\mathbf{u}_{t},\bar{\mathbf{u}})+D_{\mathbf{u}}(\nabla\mathbf{u},\nabla\bar{\mathbf{u}})+\rho((\mathbf{u}\cdot\nabla)\mathbf{u},\bar{\mathbf{u}})-(\pi,\nabla\cdot\bar{\mathbf{u}})=((n+\alpha_{0})\nabla\phi,\bar{\mathbf{u}}),&\\ (\bar{\pi},\nabla\cdot{\mathbf{u}})=0,&\end{array}\right.$$

$$\underset{\mathbf{v}\in {\mathcal{X}}_{\text{u}}\setminus \{0\}}{\sup }\frac{-(w,\nabla \cdot \mathbf{v})}{\Vert \mathbf{v}{\Vert }_{{\mathcal{X}}_{\text{u}}}}\ge \beta \Vert w{\Vert }_{{\mathcal{X}}_{\pi }},\forall w\in {\mathcal{X}}_{\pi }.$$

$$\left\{\begin{array}[]{lc}\vspace{0.1 cm}D_{\mathbf{u}}(\nabla(\mathbb{P}_{\mathbf{u}}{\mathbf{u}}-{\mathbf{u}}),\nabla\bar{\mathbf{u}})-(\mathbb{P}_{\pi}\pi-\pi,\nabla\cdot\bar{\mathbf{u}})=0,\quad\forall\bar{\mathbf{u}}\in\mathcal{X}_{\mathbf{u}},&\\ (\nabla\cdot(\mathbb{P}_{\mathbf{u}}{\mathbf{u}}-{\mathbf{u}}),\bar{\pi})=0,\quad\forall\bar{\pi}\in\mathcal{X}_{\pi},&\end{array}\right.$$

$$\Vert [\mathbf{u}-{\mathbb{P}}_{\mathbf{u}}\mathbf{u},\pi -{\mathbb{P}}_{\pi }\pi ]{\Vert }_{H^1\times L^2}+\frac{1}{h}\Vert \mathbf{u}-{\mathbb{P}}_{\mathbf{u}}\mathbf{u}{\Vert }_{L^2}\le K{h}^r\Vert [\mathbf{u},\pi ]{\Vert }_{H^{r+1}\times H^r},$$

$$\Vert [{\mathbb{P}}_{\mathbf{u}}\mathbf{u},{\mathbb{P}}_{\pi }\pi ]{\Vert }_{W^{1,6}\times L^6}\le C\Vert [\mathbf{u},\pi ]{\Vert }_{H^2\times H^1}.$$

$${\mathbb{P}}_n:{\hat{H}}^1(\Omega )\to {\mathcal{X}}_n,{\mathbb{P}}_c:H^1(\Omega )\to {\mathcal{X}}_c,{\mathbb{P}}_{\mathbf{\sigma }}:{\mathbf{H}}_{\sigma }^1(\Omega )\to {\mathcal{X}}_{\mathbf{\sigma }},$$

$$\left\{\begin{array}[]{lc}\vspace{0.1 cm}(\nabla(\mathbb{P}_{n}n-n),\nabla\bar{n})=0,\quad\forall\bar{n}\in\mathcal{X}_{n},&\\ \vspace{0.1 cm}(\nabla(\mathbb{P}_{c}c-c),\nabla\bar{c})+(\mathbb{P}_{c}c-c,\bar{c})=0,\quad\forall\bar{c}\in\mathcal{X}_{c},&\\ \vspace{0.1 cm}(\nabla\cdot(\mathbb{P}_{\bm{\sigma}}{\bm{\sigma}}-{\bm{\sigma}}),\nabla\cdot\bar{\bm{\sigma}})+(\mbox{rot}(\mathbb{P}_{\bm{\sigma}}{\bm{\sigma}}-{\bm{\sigma}}),\mbox{rot }\bar{\bm{\sigma}})+(\mathbb{P}_{\bm{\sigma}}{\bm{\sigma}}-{\bm{\sigma}},\bar{\bm{\sigma}})=0,\ \ \forall\bar{\bm{\sigma}}\in\mathcal{X}_{\bm{\sigma}},&\end{array}\right.$$

$$\left\{\begin{array}[]{lc}\|n-\mathbb{P}_{n}n\|_{L^{2}}+h\|n-\mathbb{P}_{n}n\|_{H^{1}}\leq Kh^{r_{1}+1}\|n\|_{H^{r_{1}+1}},&\forall n\in H^{r_{1}+1}(\Omega),\\[5.69046pt] \|c-\mathbb{P}_{c}c\|_{L^{2}}+h\|c-\mathbb{P}_{c}c\|_{H^{1}}\leq Kh^{r_{2}+1}\|c\|_{H^{r_{2}+1}},&\forall c\in H^{r_{2}+1}(\Omega),\\[5.69046pt] \|{\bm{\sigma}}-\mathbb{P}_{\bm{\sigma}}{\bm{\sigma}}\|_{L^{2}}+h\|{\bm{\sigma}}-\mathbb{P}_{\bm{\sigma}}{\bm{\sigma}}\|_{H^{1}}\leq Ch^{r_{3}+1}\|{\bm{\sigma}}\|_{H^{r_{3}+1}},&\forall{\bm{\sigma}}\in\mathbf{H}^{r_{3}+1}(\Omega).\end{array}\right.$$

$$\Vert [{\mathbb{P}}_{n}n,{\mathbb{P}}_{c}c,{\mathbb{P}}_{\mathbf{\sigma }}\mathbf{\sigma }]{\Vert }_{H^1}\le \Vert [n,c,\mathbf{\sigma }]{\Vert }_{H^1},$$

$$\Vert [{\mathbb{P}}_{n}n,{\mathbb{P}}_{c}c,{\mathbb{P}}_{\mathbf{\sigma }}\mathbf{\sigma }]{\Vert }_{W^{1,6}}\le C\Vert [n,c,\mathbf{\sigma }]{\Vert }_{H^2},$$

$$\Vert [n_h,c_h,{\mathbf{\sigma }}_h]{\Vert }_{W^{1,6}}\le C{h}^{-p}\Vert [n_h,c_h,{\mathbf{\sigma }}_h]{\Vert }_{H^1}\text{mbox}forall[n_h,c_h,{\mathbf{\sigma }}_h]\in {\mathbb{P}}_n\times {\mathbb{P}}_c\times {\mathbb{P}}_{\mathbf{\sigma }},$$

$$B({\mathbf{v}}_1,{\mathbf{v}}_2,{\mathbf{v}}_3)$$

$$A(\mathbf{v},w_1,w_2)$$

Full text

Keywords. Chemotaxis-Navier-Stokes system, finite elements, convergence rates, error estimates.

AMS subject classifications. 35Q35; 35Q92; 92C17; 65M12; 65M15; 65M60

1 Introduction

Chemotaxis is the oriented movement of cells towards the concentration gradient of certain chemicals in their environment. In particular, when the movement of cells is directed towards the increasing concentration of a signal, the phenomenon is known as chemotaxis by atracction. This kind of phenomena, which play an outstanding role in a large range of biological applications, are modeled, in their simplest form, by the Keller-Segel system. However, some experimental studies, as reported in [17, 29], reveal that the chemotactic motion in liquid environments affects substantially the migration of cells. Some examples of these facts are the phenomenon of broadcast spawning, the pattern generation and spontaneous emergence of turbulence in populations of aerobic bacteria suspended in sessile drops of water [28]. This kind of interaction can be modeled through a chemotaxis-Navier-Stokes system. This model is given by the following system of PDEs:

[TABLE]

where $\eta=\eta(\textit{{x}},t),$ $c=c(\textit{{x}},t),$ $\pi(\textit{{x}},t)$ and $\mathbf{u}(\textit{{x}},t)$ denote respectively the cell density, the concentration of an attractive chemical signal, the hydrostatic pressure, and the velocity field of the fluid at position $\textit{{x}}\in\Omega\subset\mathbb{R}^{d},$ $d=2,3,$ and time $t\in(0,T]$. This model describes the interaction between a type of cells (e.g., bacteria), and a chemical signal which is consumed with a rate proportional to the amount of organisms. The cells and chemical substance are transported by a viscous incompressible fluid under the influence of a force due to the aggregation of cells. The equation for the velocity field $\mathbf{u}$ is described by the incompressible Navier-Stokes system with forcing term given by $-\eta\nabla\phi,$ which represents the effects due to density variations caused by cell aggregation. The parameters $\chi,D_{n},D_{c},\rho$ and $D_{\mathbf{u}}$ are positive constants that represent the chemotactic coefficient, the cell diffusion constant, the chemical diffusion constant, the fluid density and the viscosity of fluid, respectively.

System (1.1) is completed with the following initial and boundary data:

[TABLE]

The mathematical understanding of the existence and uniqueness of a solution for (1.1)-(1.2) is quite challenging, due the coupling between the Navier-Stokes equations and the chemotaxis system. However, there are several results of existence, uniqueness, regularity and qualitative properties of the solutions for system (1.1)-(1.2) and related models (see for instance [7, 18, 20, 30, 31, 32, 34] and references therein). In [30], it was proved the existence of global classical solutions in two-dimensional bounded convex domains. The results of [30] were extended to nonconvex domains in [18]. Results of convergence of classical solutions to the corresponding stationary model were analyzed in [32, 34]. In [31], it was proved the existence of global weak solutions in bounded three-dimensional convex domains, considering the Stokes system instead of the Navier-Stokes equations, that is, neglecting the nonlinear term $(\mathbf{u}\cdot\nabla)\mathbf{u}$ in the fluid equation. The existence of global weak solutions for the full three-dimentional Chemotaxis Navier-Stokes model in bounded convex domains was obtained in [31]. In [31] the author also proved that any eventual energy solution becomes smooth after some waiting time. In [20] the author proved the existence of weak solutions for general bounded domains of $\mathbb{R}^{3}$ for the chemotaxis-Navier-Stokes system with logistic source. In the same paper, the author also proved that after some waiting time, the weak solutions become smooth and converge to some steady state. Results of global existence of mild solutions in bounded domains of $\mathbb{R}^{d}$ , $d=2,3,$ for small initial data in $L^{p}$-spaces were obtained in [7]. However, to the best of our knowledge, there is no convergence numerical analysis for the solutions of (1.1)-(1.2). We only know some numerical simulations to investigate the patterns formation and predict numerically the nonlinear dynamic of the chemotaxis-fluid system (see for instance [4, 6, 19, 21, 29]).

This paper is focused on the analysis of a numerical method to approximate the solutions of (1.1)-(1.2). As far as we know, there is no numerical analysis related to the convergence and error estimates of approximations for the weak solutions of (1.1)-(1.2). Even more, there are only a few works about numerical analysis for Keller-Segel system (i.e., with chemo-attraction and linear production) and assuming that there is no interaction with the fluid [8, 10, 22, 25, 26]. In [10] the author analyzed the existence of discrete solutions and the convergence of a finite volume scheme. In [25, 26], was obtained some error estimates for a conservative Finite Element (FE) approximation. In [8] were proved some error estimates for a fully discrete discontinuous FE method. A mixed FE approximation was analyzed in [22]. On the other hand, in the same framework of Keller-Segel system some previous energy stable numerical schemes have also been analyzed (cf. [2]). In addition, unconditionally energy stable time-discrete numerical schemes and fully discrete FE schemes for a chemo-repulsion model with quadratic production has been analyzed in [14, 15]. Some unconditionally energy stable fully discrete schemes for a parabolic repulsive-productive chemotaxis model (with linear production term) was recently analyzed in [13]. In our case, the difficulties to deal with the numerical analysis of (1.1)-(1.2) come from the strong coupling nonlinear term $\nabla\cdot(\eta\nabla c)$ and the transport terms $\mathbf{u}\cdot\nabla\eta,$ $\mathbf{u}\cdot\nabla c$ and $\mathbf{u}\cdot\nabla\mathbf{u}.$ Therefore, in order to deal the chemotaxis term $\nabla\cdot(\eta\nabla c)$ we introduce a new variable given by the gradient of the chemical concentration, that is, $\bm{\sigma}=\nabla c,$ which allow us to control the strong regularity request by the system; on the other hand, the transport terms have the difficulty that the corresponding discrete forms do not preserve the same properties of the continuous case. We concentrate our analysis by introducing a splitting mixed finite element method for approximating the weak solutions of (1.1)-(1.2) and analyze some error estimates. This scheme is well-posed and it is mass-conservative. We obtain uniform estimates and analyze the convergence towards weak solutions. We also present some numerical simulations in order to validity the numerical results.

The plan of this paper is as follows. In Section 2, we first establish some basic notations and recall some existence and uniqueness results of (1.1)-(1.2) in the continuous case. We also define the weak formulation of (1.1)-(1.2), which will be used to construct the numerical approximation. In Section 3, we construct a mass-conservative numerical approximation for the weak solutions of the chemotaxis-fluid system (1.1) with initial and boundary data (1.2), by using FE approximations in space and finite differences in time, and give some preliminary results concerning to the FE spaces. In Section 4, we analyze the well-posedness and the mass-conservation property (see (2.9) below) for the numerical scheme; and following an inductive procedure, in Section 5, we obtain some uniform estimates for any solution of scheme (3.2), subsequently required in the convergence analysis. The analysis of Section 5 is focused on the two-dimensional case. In Section 6, we make some comments about the extensions of these results to the three-dimensional case. Finally, in Section 7, we provide some numerical simulations in agreement with the numerical results.

2 The continuous problem

In this section, we establish a variational formulation of (1.1)-(1.2) which will be used to construct the numerical scheme. In order to deal with the chemotaxis term in the cell density equation, we introduce a splitting variational formulation for the chemoatractant concentration, which yields to consider an auxiliar variable representing the gradient of the chemical signal. In this section, after to establish the basic notation which will be used, we recall some existence results for two and three dimensional bounded domains, and define a variational formulation of the continuous problem.

2.1 Notations

We consider the standard Sobolev and Lebesgue spaces $W^{k,p}(\Omega)$ and $L^{p}(\Omega),$ with respective norms $\|\cdot\|_{W^{k,p}}$ and $\|\cdot\|_{L^{p}}.$ In particular, we denote $W^{k,2}(\Omega)=H^{k}(\Omega).$ Also, $W_{0}^{1,p}(\Omega)$ denotes the elements of $W^{1,p}(\Omega)$ with trace zero on the boundary of $\Omega.$ The $L^{2}(\Omega)$-inner product will be represented by $(\cdot,\cdot).$ Corresponding Sobolev spaces of vector valued functions will be denoted by ${\bf W}^{k,p}(\Omega),$ ${\bf L}^{p}(\Omega),$ and so on. Also, we will use the following function spaces

[TABLE]

where $\bm{\nu}$ denotes the unit outward normal vector to the boundary. We will use the following equivalent norms in $H^{1}(\Omega)$ and ${\bf H}_{\sigma}^{1}(\Omega)$, respectively (see [24] and [1, Corollary 3.5]):

[TABLE]

[TABLE]

We recall the Poincaré inequality

[TABLE]

for some constant $C_{P}>0$ independent of $\mathbf{v}$. Also, we will use the classical interpolations inequalities

[TABLE]

[TABLE]

Finally, the letter $C$ denotes different positive constants (independent of discrete parameters) which may change from line to line (or even within the same line).

2.2 Variational formulation

In this subsection we give a variational formulation for (1.1)-(1.2). We start by recalling a result in [30] which provides the existence and uniqueness of global classical solution for (1.1)-(1.2) in bounded and convex domains of $\mathbb{R}^{2}$ with boundary $\partial\Omega$ smooth enough. Let us assume that the initial data satisfy the following regularity conditions:

[TABLE]

where $A^{\alpha}$ denotes the fractional Stokes operator with domain $D(A^{\alpha}).$ Then, the following theorem holds.

Theorem 2.1**.**

Let $\Omega\subset\mathbb{R}^{2}$ be a bounded domain with smooth boundary, and suppose that $\eta_{0},c_{0},\mathbf{u}_{0},\phi$ satisfy (2.6). Then (1.1)-(1.2) possesses a classical solution which is global in time. This solution is unique, up to addition of constants to the pressure $\pi$, and for all $T\in(0,\infty)$ the solution has the following regularity properties

[TABLE]

Proof.

The proof is given in [30, Theorem 1.1]. Indeed, the proof of [30] requires that $\Omega$ be a bounded convex domain. This condition comes from an integration term on the boundary which appears dealing with the differentiation of $\int_{\Omega}|\nabla\sqrt{c}|^{2}.$ That boundary term takes the form $\int_{\partial\Omega}\frac{\partial}{\partial\nu}|\nabla c|^{2}dS.$ Thus, if $\Omega$ is convex, the last boundary integral turns out to be non-positive and thus, this term can be neglected. However, the convexity condition can be removed by using Lemma 4.2 in [23], as used in [18, 20]. ∎

For the three-dimensional setting of (1.1)-(1.2) and large initial data, Winkler in [31] proved the existence of weak solutions (for bounded convex domains). The notion of weak solution considered in [31] is the following one:

Definition 2.2**.**

(Weak solution) A weak solution of (1.1)-(1.2) is a triple $[n,c,{u}]$ of functions such that $\eta\in{L^{1}([0,T],W^{1,1}}(\Omega)),$ $c\in{L^{1}([0,T],W^{1,1}}(\Omega)),$ $\mathbf{u}\in L^{1}([0,T],W^{1,1}),$ verifying $\eta\geq 0$ and $c\geq 0$ a.e. in $\Omega\times(0,T),$ $nc\in L^{1}(0,T;L^{1}(\Omega)),$ $\mathbf{u}\otimes\mathbf{u}\in L^{1}(0,T;L^{1}(\Omega)),$ $\eta\nabla c,\eta c$ and $\eta\mathbf{u}\in L^{1}(0,T;L^{1}(\Omega)),$ $\nabla\cdot\mathbf{u}=0$ a.e in $\Omega\times(0,T),$ and such that

[TABLE]

for all $\varphi\in C_{0}^{\infty}([0,T]\times\bar{\Omega})$ with $\varphi(T)=0,$ and ${\psi}$ in $C_{0}^{\infty}([0,T]\times\bar{\Omega})^{3}$ with ${\psi}(T)=0.$

By assuming $\phi\in W^{2,\infty}({\Omega})$ and the initial data $[n_{0},c_{0},\mathbf{u}_{0}]$ satisfying

[TABLE]

in [31], the author proves the following theorem of global existence of weak solutions.

Theorem 2.3**.**

Assume that the initial data $[n_{0},c_{0},{\mathbf{u}}_{0}]$ satisfies (2.8) and $\phi\in W^{2,\infty}({\Omega}).$ Then there exist $\eta\in{L^{\infty}(0,T;L^{1}(\Omega))\cap{L^{\frac{5}{4}}([0,T],W^{1,\frac{5}{4}}}(\Omega))},$ $c\in{L^{\infty}(0,T;L^{\infty}(\Omega))\cap L^{4}([0,T];W^{1,4}(\Omega))},$ $\mathbf{u}\in L^{2}([0,T],\bf{V}),$ such that $[n,c,\bf{u}]$ is a weak solution of (1.1)-(1.2) in the sense of Definition 2.2.

Notice that if $[\eta,c,\mathbf{u},\pi]$ is a classical solution of system (1.1)-(1.2), then we get the following basic property:

[TABLE]

Equality (2.9) is obtained integrating the equation (1.1)1 and using the boundary conditions (1.2)2; and it means that the total mass of bacteria is conserved in time and equal to the initial mass.

In order to design the scheme, taking into account the property (2.9), we first consider the auxiliary variable $n:=\eta-\frac{1}{|\Omega|}\int_{\Omega}\eta=\eta-\frac{1}{|\Omega|}\int_{\Omega}\eta_{0}:=\eta-\alpha_{0}$. Also, in order to introduce the weak formulation for the $c$-equation (1.1)2, we consider the auxiliary unknown $\bm{\sigma}=\nabla c.$ Thus, we obtain the following mixed variational form for the variables $\eta,c,\bm{\sigma},\mathbf{u}$ and $\pi$:

[TABLE]

Taking the time derivative in the equation (2.10)2 we have

[TABLE]

Then, choosing $\bar{c}=\nabla\cdot\ \bar{\bm{\sigma}}$ in (2.10)3, substituting it into (2.11), and adding the term $D_{c}(\mbox{rot }{\bm{\sigma}},\mbox{rot }\bar{\bm{\sigma}})$ using that $\mbox{rot }{\bm{\sigma}}=\mbox{rot}(\nabla c)=0$, we get the following variational formulation:

[TABLE]

for all $[\bar{n},\bar{\bm{\sigma}},\bar{c},\bar{\mathbf{u}},\bar{\pi}]\in\widehat{H}^{1}(\Omega)\times\mathbf{H}_{\sigma}^{1}(\Omega)\times H^{1}(\Omega)\times\mathbf{H}_{0}^{1}(\Omega)\times{L}_{0}^{2}(\Omega)$. The integral equations (2.12) define a variational form of (1.1)-(1.2). Also, it is straightforward to check that if we have a smooth solution of (2.12), then $\bm{\sigma}=\nabla c.$

3 Numerical scheme

In this section, we construct a mass-conservative numerical approximation for the weak solutions of the chemotaxis-fluid system (1.1) with initial and boundary data (1.2). The idea is to use finite element approximations in space and finite differences in time (considered for simplicity on a uniform partition of $[0,T]$ with time step $\Delta t=T/N:(t_{n}=n\Delta t)_{n=0}^{n=N}$), combined with splitting ideas to decouple the computation of the fluid part from the chemotaxis one. Moreover, in order to deal with the velocity trilinear form and the nonlinear convective terms, we will use the skew-symmetric forms $A$ and $B$ given in (3.9)-(3.10). Concerning to the space discretization, we consider conformed finite element spaces: $\mathcal{X}_{n}\times\mathcal{X}_{c}\times\mathcal{X}_{\bm{\sigma}}\times\mathcal{X}_{\textbf{u}}\times\mathcal{X}_{\pi}\subset\widehat{H}^{1}(\Omega)\times{H}^{1}(\Omega)\times\mathbf{H}^{1}_{\sigma}(\Omega)\times\textbf{H}^{1}_{0}(\Omega)\times{L}^{2}_{0}(\Omega)$ corresponding to a family of shape-regular and quasi-uniform triangulations of $\overline{\Omega}$, $\{\mathcal{T}_{h}\}_{h>0}$, made up of triangles $K$, so that $\overline{\Omega}=\cup_{K\in\mathcal{T}_{h}}K$, where $h=\max_{K\in\mathcal{T}_{h}}h_{K}$, with $h_{K}$ being the diameter of $K$.

3.1 Hyphotesis and preliminary results

We assume that $\mathcal{X}_{\textbf{u}}$ and $\mathcal{X}_{\pi}$ satisfy the following discrete inf-sup condition: There exists a constant $\beta>0$, independent of $h$, such that

[TABLE]

Remark 3.1**.**

(Some possibilities for the choice of the discrete spaces)* For the spaces $[\mathcal{X}_{\mathbf{u}},\mathcal{X}_{\pi}]$, we can choose the Taylor-Hood approximation $\mathbb{P}_{r}\times\mathbb{P}_{r-1}$ (for $r\geq 2$) ([11, 27]); or the approximation $\mathbb{P}_{1}-bubble\times\mathbb{P}_{1}$ ([11]) (for $r=1$). On the other hand, the spaces $[\mathcal{X}_{n},\mathcal{X}_{c},\mathcal{X}_{\bm{\sigma}}]$ are approximated by $\mathbb{P}_{r_{1}}\times\mathbb{P}_{r_{2}}\times\mathbb{P}_{r_{3}}$- continuous FE, with $r_{i}\geq 1$ ($i=1,2,3$).*

Then, we can define the well-known Stokes operator $(\mathbb{P}_{\mathbf{u}},\mathbb{P}_{\pi}):\mathbf{H}^{1}_{0}(\Omega)\times L^{2}_{0}(\Omega)\rightarrow\mathcal{X}_{\mathbf{u}}\times\mathcal{X}_{\pi}$ such that $[\mathbb{P}_{\mathbf{u}}\mathbf{u},\mathbb{P}_{\pi}\pi]\in\mathcal{X}_{\mathbf{u}}\times\mathcal{X}_{\pi}$ satisfies

[TABLE]

and the following approximation and stability properties hold ([12]):

[TABLE]

[TABLE]

Moreover, we consider the following interpolation operators:

[TABLE]

such that for all $n\in\widehat{H}^{1}(\Omega)$, $c\in H^{1}(\Omega)$ and ${\bm{\sigma}}\in\mathbf{H}^{1}_{\sigma}(\Omega)$, $\mathbb{P}_{n}n\in\mathcal{X}_{n}$, $\mathbb{P}_{c}c\in\mathcal{X}_{c}$ and $\mathbb{P}_{\bm{\sigma}}{\bm{\sigma}}\in\mathcal{X}_{\bm{\sigma}}$ satisfy

[TABLE]

respectively. Observe that from Lax-Milgram Theorem, we have that the interpolation operators $\mathbb{P}_{n}$, $\mathbb{P}_{c}$ and $\mathbb{P}_{\bm{\sigma}}$ are well defined. Moreover, it is well known that the following interpolation errors hold:

[TABLE]

Also, the stability properties

[TABLE]

[TABLE]

hold. Inequality (3.7) can be deduced from (3.5), and (3.8) can be obtained from (3.6) using the inverse inequality

[TABLE]

with $p=2/3$ (in the 2D case) and $p=1$ (in the 3D case), and comparing $\mathbb{P}_{n,c,{\bm{\sigma}}}$ with an average interpolation of Clement or Scott-Zhang type (which are stable in $W^{1,6}$-norm). Finally, we consider the following skew-symmetric trilinear forms which will be used in the formulation of the numerical scheme:

[TABLE]

It is easy to verify that

[TABLE]

3.2 Definition of the scheme

Taking into account (2.12), we consider the following first order in time, linear and semi-coupled scheme:

Initialization: Let $[n^{0}_{h},c^{0}_{h},\bm{\sigma}^{0}_{h},\mathbf{u}^{0}_{h}]=[\mathbb{P}_{n}n_{0},\mathbb{P}_{c}c_{0},\mathbb{P}_{\bm{\sigma}}{\bm{\sigma}}_{0},\mathbb{P}_{\mathbf{u}}\mathbf{u}_{0}]\in\mathcal{X}_{n}\times\mathcal{X}_{c}\times\mathcal{X}_{\bm{\sigma}}\times\mathcal{X}_{\mathbf{u}}$.

Time step $m$: Given the vector $[n^{m-1}_{h},c^{m-1}_{h},\bm{\sigma}^{m-1}_{h},\mathbf{u}^{m-1}_{h}]\in\mathcal{X}_{n}\times\mathcal{X}_{c}\times\mathcal{X}_{\bm{\sigma}}\times\mathcal{X}_{\mathbf{u}}$, compute $[n^{m}_{h},c^{m}_{h},\bm{\sigma}^{m}_{h},\mathbf{u}^{m}_{h},\pi_{h}^{m}]\in\mathcal{X}_{n}\times\mathcal{X}_{c}\times\mathcal{X}_{\bm{\sigma}}\times\mathcal{X}_{\mathbf{u}}\times\mathcal{X}_{\pi}$ such that for each $[\bar{n},\bar{c},\bar{\bm{\sigma}},\bar{\mathbf{u}},\bar{\pi}]\in\mathcal{X}_{n}\times\mathcal{X}_{c}\times\mathcal{X}_{\bm{\sigma}}\times\mathcal{X}_{\mathbf{u}}\times\mathcal{X}_{\pi}$ it holds:

[TABLE]

where, in general, we denote $\delta_{t}a^{m}_{h}=\frac{a^{m}_{h}-a^{m-1}_{h}}{\Delta t}$.

Remark 3.2**.**

The skew-symmetric forms $A$ and $B$ verifying the properties (3.11)-(3.14) will be important in order to get uniform estimates for the discrete solutions and to develop convergence analysis.

4 Well-posedness and mass-conservation

In this section, we analyze the well-posedness of the scheme (3.2) and the mass-conservation property. With this aim, we define $\eta^{m}_{h}:=n^{m}_{h}+\alpha_{0}$. Then $\eta^{m}_{h}$ verifies

[TABLE]

Lemma 4.1**.**

(Mass conservation)* The discrete cell density $\eta^{m}_{h}$ satisfies the mass-conservation property (2.9).*

Proof.

From the construction of scheme (3.2), for each $m\geq 0$ it holds that $\int_{\Omega}n^{m}_{h}=0$ (since $n^{m}_{h}\in\mathcal{X}_{n}\subset\widehat{H}^{1}(\Omega)$); from which we deduce that $\eta^{m}_{h}$ satisfies

[TABLE]

∎

Now, in the next result, we prove the well-posedness of the scheme (3.2).

Theorem 4.2**.**

(Unconditional well-posedness)* The numerical scheme (3.2) is well-posed, that is, there exists a unique $[n^{m}_{h},c^{m}_{h},\bm{\sigma}^{m}_{h},\mathbf{u}^{m}_{h},\pi_{h}^{m}]\in\mathcal{X}_{n}\times\mathcal{X}_{c}\times\mathcal{X}_{\bm{\sigma}}\times\mathcal{X}_{\mathbf{u}}\times\mathcal{X}_{\pi}$ solution of the scheme (3.2).*

Proof.

Taking into account that the scheme (3.2) is an algebraic linear system, it suffices to prove the uniqueness. For that, suppose that there exist $[n^{m}_{h,1},c^{m}_{h,1},\bm{\sigma}^{m}_{h,1},\mathbf{u}^{m}_{h,1},\pi_{h,1}^{m}],[n^{m}_{h,2},c^{m}_{h,2},\bm{\sigma}^{m}_{h,2},\mathbf{u}^{m}_{h,2},\pi_{h,2}^{m}]\in\mathcal{X}_{n}\times\mathcal{X}_{c}\times\mathcal{X}_{\bm{\sigma}}\times\mathcal{X}_{\mathbf{u}}\times\mathcal{X}_{\pi}$ two possible solutions of the scheme (3.2). Then defining $n^{m}_{h}=n^{m}_{h,1}-n^{m}_{h,2}$, $c^{m}_{h}=c^{m}_{h,1}-c^{m}_{h,2}$, $\bm{\sigma}^{m}_{h}=\bm{\sigma}^{m}_{h,1}-\bm{\sigma}^{m}_{h,2}$ $\mathbf{u}^{m}_{h}=\mathbf{u}^{m}_{h,1}-\mathbf{u}^{m}_{h,2}$ and $\pi_{h}^{m}=\pi_{h,1}^{m}-\pi_{h,2}^{m}$, we have that $[n^{m}_{h},c^{m}_{h},\bm{\sigma}^{m}_{h},\mathbf{u}^{m}_{h},\pi_{h}^{m}]\in\mathcal{X}_{n}\times\mathcal{X}_{c}\times\mathcal{X}_{\bm{\sigma}}\times\mathcal{X}_{\mathbf{u}}\times\mathcal{X}_{\pi}$ satisfies

[TABLE]

for all $[\bar{n},\bar{c},\bar{\bm{\sigma}},\bar{\mathbf{u}},\bar{\pi}]\in\mathcal{X}_{n}\times\mathcal{X}_{c}\times\mathcal{X}_{\bm{\sigma}}\times\mathcal{X}_{\mathbf{u}}\times\mathcal{X}_{\pi}$. Taking $[\bar{n},\bar{c},\bar{\bm{\sigma}},\bar{\mathbf{u}},\bar{\pi}]=[n^{m}_{h},c^{m}_{h},\bm{\sigma}^{m}_{h},\mathbf{u}^{m}_{h},\Delta t\pi_{h}^{m}]$ in (4.1), using the properties (3.13)-(3.14), and adding, we obtain

[TABLE]

which implies that $[n^{m}_{h},c^{m}_{h},\bm{\sigma}^{m}_{h},\mathbf{u}^{m}_{h}]=[0,0,\bm{0},\bm{0}]$. Finally, using that $\mathbf{u}^{m}_{h}=\bm{0}$ in (4.1)4, and using the discrete inf-sup condition (3.1) we deduce that $\pi^{m}_{h}=0$. ∎

5 Uniform estimates and convergence

In this section, we will obtain some uniform estimates for any solution of scheme (3.2) that will be used in the convergence analysis. Here, we focus on the two-dimensional case. Later, in subsection 6, we study the three-dimensional case. With this aim, we make the following inductive hypothesis: there exists a positive constant $K>0$, independent of $m$, such that

[TABLE]

After the convergence analysis we verify the vality of (5.1). It is worthwhile to remark that the use of induction hypothesis to deal with the convergence analysis of approximating schemes in nonlinear partial differential equations have been considered by several authors (see for instance, [5, 9, 33, 35]). In [9] the authors analyze a numerical method for a model describing compressible miscible displacement in porous media, and consider an inductive hypothesis on the difference of pressures $\pi=p-p_{h}$ comming from the Darcy model, in norm $W^{1,\infty}.$ In the same spirit, by using inductive hypotheses, a superconvergence estimate of a Combined Mixed Finite Element and Discontinuous Galerkin Method for a Compressible Miscible Displacement Problem was obtained in [33] (see also [5]). In the context of the Keller-Segel system without fluid influence, in [35], an inductive hypothesis on $\|{\bm{\sigma}}^{m}_{h}\|_{W^{1,\infty}}$ is considered. Observe that the norm in (5.1) is in $H^{1}(\Omega)$ instead of $W^{1,\infty}(\Omega).$

Additionally, we will use the following discrete Gronwall lemma:

Lemma 5.1**.**

([16, p. 369]) Assume that $\Delta t>0$ and $B,b^{k},d^{k},g^{k},h^{k}\geq 0$ satisfy:

[TABLE]

Then, it holds

[TABLE]

5.1 Uniform estimates in finite time

We can prove the following uniform estimates in finite time for the discrete cell and chemical variables, $n^{m}_{h}$ and $c^{m}_{h}$, in weak norms.

Lemma 5.2**.**

(Uniform estimates for $n_{h}^{m}$)* Assume the inductive hypothesis (5.1). Then $n^{m}_{h}$ is bounded in $l^{\infty}(L^{2})\cap l^{2}(H^{1})$.*

Proof.

Testing (3.2)a by $\bar{n}=2\Delta tn^{m}_{h}$, using the equality $(a-b,2a)=|a|^{2}-|b|^{2}+|a-b|^{2}$, taking into account the property (3.14) and using the fact that $\int_{\Omega}n^{m}_{h}=0$ (since $n^{m}_{h}\in\mathcal{X}_{n}$), we have

[TABLE]

Using the Hölder and Young inequalities and the 2D interpolation inequality (2.4), we get

[TABLE]

Thus, taking into account the inductive hypothesis (5.1), from (5.2)-(5.3) we arrive at

[TABLE]

Then, summing in (5.4) from $m=1$ to $m=r$, we arrive at

[TABLE]

Therefore, applying Lemma 5.1 to (5.5) we deduce

[TABLE]

where the constant $C>0$ depends on $(\chi,D_{n},\alpha_{0},T,n_{0})$, but is independent of $(\Delta t,h)$ and $r$; thus we conclude the proof. ∎

Lemma 5.3**.**

(Uniform estimates for $c_{h}^{m}$)* Assume the inductive hypothesis (5.1). Then, $c^{m}_{h}$ is bounded in $l^{\infty}(L^{2})\cap l^{2}(H^{1})$.*

Proof.

Testing (3.2)c by $\bar{c}=2\Delta tc^{m}_{h}$, using the equality $(a-b,2a)=|a|^{2}-|b|^{2}+|a-b|^{2}$ and taking into account the property (3.14), we obtain

[TABLE]

Using the Hölder and Young inequalities, we have

[TABLE]

Therefore, from (5.6) and (5.7), we arrive at

[TABLE]

Then, summing in (5.8) from $m=1$ to $m=r$, applying Lemma 5.1 and taking into account that $n^{m}_{h}$ is bounded in $l^{2}(H^{1})$ (see Lemma 5.2), we arrive at

[TABLE]

the constant $C>0$ depends on $(\chi,D_{c},D_{n},\gamma,\alpha_{0},T,n_{0},c_{0})$, but is independent of $(\Delta t,h)$ and $r$. ∎

5.2 Error estimates in weak norms

The aim of this subsection is to obtain optimal error estimates for any solution $[n^{m}_{h},c^{m}_{h},{\bm{\sigma}}^{m}_{h},\mathbf{u}^{m}_{h},\pi^{m}_{h}]$ of the scheme (3.2), with respect to a sufficiently regular solution $[n,c,{\bm{\sigma}},\mathbf{u},\pi]$ of (2.12) in the two dimensional case. We start by introducing the following notations for the errors at $t=t_{m}$: $e_{n}^{m}=n^{m}-n^{m}_{h}$, $e_{c}^{m}=c^{m}-c^{m}_{h}$, $e_{\bm{\sigma}}^{m}={\bm{\sigma}}^{m}-{\bm{\sigma}}^{m}_{h}$, $e_{\mathbf{u}}^{m}=\mathbf{u}^{m}-\mathbf{u}^{m}_{h}$ and $e_{\pi}^{m}=\pi^{m}-\pi^{m}_{h}$, where $w^{m}$ denote, in general, the value of $w$ at time $t_{m}$. Then, subtracting the scheme (3.2) to (2.12) at $t=t_{m}$, we obtain that $[e_{n}^{m},e^{m}_{c},e_{{\bm{\sigma}}}^{m},e^{m}_{\mathbf{u}},e^{m}_{\pi}]$ satisfies

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

for all $[\bar{n},\bar{c},\bar{\bm{\sigma}},\bar{\mathbf{u}},\bar{\pi}]\in\mathcal{X}_{n}\times\mathcal{X}_{c}\times\mathcal{X}_{\bm{\sigma}}\times\mathcal{X}_{\mathbf{u}}\times\mathcal{X}_{\pi}$, where $\omega_{n}^{m},\omega_{c}^{m},\omega^{m}_{\bm{\sigma}},\omega^{m}_{\mathbf{u}}$ are the consistency errors associated to the scheme (3.2), that is, $\omega_{n}^{m}=\delta_{t}n^{m}-(n^{m})_{t}$ and so on.

1. Error estimate for the cell density $n$

Considering the interpolation operators $\mathbb{P}_{n},\mathbb{P}_{\bm{\sigma}},\mathbb{P}_{\mathbf{u}}$ defined in (3.5) and (3.2), we decompose the total errors $e_{n}^{m},e_{\bm{\sigma}}^{m},e_{\mathbf{u}}^{m}$ as follows:

[TABLE]

[TABLE]

[TABLE]

Then, taking into account (3.5)1, from (5.9) and (5.14)-(5.16) we have

[TABLE]

Taking $\bar{n}=\xi_{n}^{m}$ in (5.17), using (3.14) and taking into account that $\int_{\Omega}\xi_{n}^{m}=0$ (since $\xi_{n}^{m}\in\mathcal{X}_{n}$), we get

[TABLE]

Then, using the Hölder and Young inequalities, (3.3)-(3.4), (3.6) and (3.8), we control the terms on the right hand side of (5.18) as follows:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where, in the first and second inequalities of (5.24), the 2D interpolation inequality (2.4) and the inductive hypothesis (5.1) were used. Therefore, from (5.18)-(5.24), we arrive at

[TABLE]

2. Error estimate for the chemical concentration $c$

Considering the interpolation operator $\mathbb{P}_{c}$ defined in (3.5)2, we decompose the total error $e_{c}^{m}$ as follows:

[TABLE]

Then, taking into account (3.5)2, from (5.11), (5.14), (5.16) and (5.26), we have

[TABLE]

Taking $\bar{c}=\xi_{c}^{m}$ in (5.27) and using (3.14), we get

[TABLE]

Then, using the Hölder and Young inequalities, (3.3)-(3.4), (3.6) and (3.8), we control the terms on the right hand side of (5.28) as follows:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Therefore, denoting $\gamma_{1}:=\frac{1}{D_{c}}+\Delta t+1$ and $\gamma_{2}:=\Delta t+1$, from (5.28)-(5.34), we arrive at

[TABLE]

3. Error estimate for the velocity $\mathbf{u}$

Taking into account the interpolation operator $\mathbb{P}_{\pi}$ defined in (3.2), we decompose the total error $e_{\pi}^{m}$ as follows:

[TABLE]

Then, taking into account (3.2), from (5.12)-(5.13), (5.14), (5.16) and (5.36), we have

[TABLE]

[TABLE]

Taking $\bar{\mathbf{u}}=\xi_{\mathbf{u}}^{m}$ in (5.37), $\bar{\pi}=\frac{1}{\rho}\xi^{m}_{\pi}$ in (5.38), using (3.13) and adding the resulting expressions, we obtain

[TABLE]

Then, using the Hölder and Young inequalities, the Poincaré inequality (2.3), (3.3)-(3.4) and (3.6), we control the terms on the right hand side of (5.39) as follows

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Therefore, from (5.39)-(5.44), we arrive at

[TABLE]

4. Error estimate for the flux $\bm{\sigma}$

Then, taking into account (3.5)3, from (5.10), (5.14)-(5.16) and (5.26), we have

[TABLE]

Taking $\bar{\bm{\sigma}}=\xi_{\bm{\sigma}}^{m}$ in (5.46), we arrive at

[TABLE]

Then, using the Hölder and Young inequalities, the equivalent norm in $\mathbf{H}^{1}_{\sigma}(\Omega)$ given in (2.2), as well as (3.3)-(3.4) and (3.6)-(3.8), we control the terms on the right hand side of (5.47) as follows

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where, in the first inequality of (5.51), the 2D interpolation inequality (2.4) was used, and in the second inequality, the inductive hypothesis (5.1) and Lemma 5.3 were used. Therefore, from (5.47)-(5.51), we arrive at

[TABLE]

5. Estimate for the terms $\|n^{m}-n^{m-1}\|_{L^{2}}$, $\|c^{m}-c^{m-1}\|_{L^{2}}$, $\|\mathbf{u}^{m}-\mathbf{u}^{m-1}\|_{L^{2}}$ and $\|{\bm{\sigma}^{m}}-\bm{\sigma}^{m-1}]\|_{L^{2}}$

Observe that the following estimate holds

[TABLE]

In fact,

[TABLE]

where the last inequality was obtained as in (5.19), with the space norm in $L^{2}$ instead of $(H^{1})^{\prime},$ the dual of $H^{1}$. Therefore, we deduce

[TABLE]

Analogously, we obtain the estimate for $c,\mathbf{u}$ and ${\bm{\sigma}}$.

Then, we can prove the following results:

Theorem 5.4**.**

Under Hypothesis of Lemma 5.2, the following error estimate holds

[TABLE]

where the constant $C(T)>0$ is independent of $m,\Delta t$ and $h$.

Proof.

The proof follows adding the four inequalities (5.25), (5.35), (5.45) and (5.52), multiplying the resulting expresion by $\Delta t$, adding from $m=1$ to $m=s$, using (5.53), (2.7) and Lemmas 5.2 - 5.3, and applying the discrete Gronwall Lemma 5.1 (recalling that $[\xi^{0}_{n},\xi^{0}_{c},\xi^{0}_{\mathbf{u}},\xi^{0}_{\bm{\sigma}}]=[0,0,{\bf 0},{\bf 0}]$). It is clear that the error estimates were derived under the inductive hypothesis (5.1). Now we have to check it. We derive (5.1) by using (5.54) recursively. First, observe that $\|\mathbb{P}_{\bm{\sigma}}\bm{\sigma}^{m-1}\|_{H^{1}}\leq\|{\bm{\sigma}}\|_{L^{\infty}(H^{1})}:=C_{0}$ for all $m\geq 1,$, from which we deduce that $\|{\bm{\sigma}}^{0}_{h}\|_{H^{1}}=\|\mathbb{P}_{\bm{\sigma}}{\bm{\sigma}}_{0}\|_{H^{1}}\leq C_{0}\leq C_{0}+1:=K$. Now, notice that

[TABLE]

Then, it is enough to show that $\|\xi^{m-1}_{\bm{\sigma}}\|_{H^{1}}\leq 1,$ for each $m\geq 2.$ For that, we consider two cases: First, if $\frac{(\Delta t)^{1/2}}{h}\leq C,$ then, by using the inverse inequality $\|\xi^{m}_{\bm{\sigma}}\|_{H^{1}}\leq h^{-1}\|\xi^{m}_{\bm{\sigma}}\|_{L^{2}},$ denoting $h_{max}:=\max\{h^{r_{1}+1},h^{r_{2}+1},h^{r_{3}+1},h^{r+1}\},$ recalling that $h_{max}=h^{k},$ (for some $k\geq 2$), and from the calculation of the norm $l^{\infty}(L^{2})$ in (5.54), we get

[TABLE]

On the other hand, if $\frac{(\Delta t)^{1/2}}{h}$ is not bounded, then $\frac{h_{max}}{\Delta t}\leq\frac{h^{2}}{\Delta t}\leq C;$ thus, from the calculation of the norm $l^{2}(H^{1})$ in (5.54), we get

[TABLE]

In any case, assuming $\Delta t$ and $h$ small enough (without any additional restriction relating the discrete parameters $(\Delta t,h)$), we conclude that $\|{\bm{\sigma}}^{1}_{h}\|_{H^{1}}\leq K$. Arguing recursively we conclude that $\|{\bm{\sigma}}^{m-1}_{h}\|_{H^{1}}\leq K$, for all $m\geq 1$. ∎

Theorem 5.5**.**

Under Hypothesis of Lemma 5.2, the following error estimate holds

[TABLE]

where the constant $C(T)>0$ is independent of $m,\Delta t$ and $h$.

Proof.

Applying Lemma 11 of [12] to (5.37)-(5.38) and usind the interpolation inequality (2.5), we have

[TABLE]

where the constant $C$ depends on the data $(\rho,D_{\mathbf{u}})$. Now, testing (5.37) by $\bar{\mathbf{u}}=\delta_{t}\xi_{\mathbf{u}}^{m}$, bounding the terms on the right hand side as in (5.59) (taking into account that the pressure term vanishes), and adding the resulting expresion with $\lambda$(5.59) (for $0<\lambda<1$), we arrive at

[TABLE]

We take $\lambda$ small enough in order to absorb the term $C\lambda\|\delta_{t}\xi_{\mathbf{u}}^{m}\|_{L^{2}}^{2}$ at the right hand side, and $\varepsilon$ small enough with respect to $\lambda$. Moreover, observe that

[TABLE]

In fact, recalling that $h_{max}:=\max\{h^{r_{1}+1},h^{r_{2}+1},h^{r_{3}+1},h^{r+1}\}$, from estimate (5.54) we have in particular that \|\xi^{m}_{\mathbf{u}}\|_{H^{1}}^{2}\leq C(T)\Big{(}\Delta t+\frac{1}{\Delta t}h_{max}^{2}\Big{)}, which implies that

[TABLE]

Therefore, we conclude (5.61) under the hypothesis

[TABLE]

On the other hand, from (5.54) we also have \|\xi^{m}_{\mathbf{u}}\|_{L^{2}}^{2}\leq C(T)\Big{(}\Delta t+h_{max}^{2}\Big{)} for each $m$. Therefore, by using the inverse inequality $\|\xi^{m}_{\mathbf{u}}\|_{H^{1}}\leq h^{-1}\|\xi^{m}_{\mathbf{u}}\|_{L^{2}}$ we obtain

[TABLE]

Therefore, we conclude (5.61) under the hypothesis

[TABLE]

Thus, we conclude (5.61) without imposing any restriction on the discrete parameters $(\Delta t,h)$ because for any choice of $(\Delta t,h)$ either (5.62) or (5.63) holds. Therefore, multiplyng (5.60) by $\Delta t$, adding from $m=1$ to $m=s$, bounding the terms $\|\omega_{\mathbf{u}}^{m}\|_{L^{2}}^{2}$ and $\|\delta_{t}\theta_{\mathbf{u}}^{m}\|_{L^{2}}^{2}$ as in (5.40) and (5.41) respectively, using (2.7), (3.3), (3.4), (3.6)1, (5.53), (5.54) and taking into account that $\xi^{0}_{\mathbf{u}}={\bf 0}$, we conclude (5.58). ∎

Remark 5.6**.**

*From (5.54)-(5.58), in particular we deduce that $\|[\mathbf{u}^{m}_{h},{\bm{\sigma}}^{m}_{h}]\|_{l^{\infty}(L^{2})\cap l^{2}(H^{1})}\leq C(T)$ and

$\|[\mathbf{u}^{m}_{h},\pi^{m}_{h}]\|_{l^{\infty}(H^{1})\times l^{2}(L^{6})}\leq C(T)$, for all $m=1,...,N$.*

Remark 5.7**.**

Theorems 5.4 and 5.5 in particular imply the convergence of the discrete solutions of the scheme (3.2) towards weak solutions of model (1.1), when the parameters $\Delta t$ and $h$ go to [math].

6 Some comments on the three-dimensional case

The results obtained in Subsections 3-5.2 can be proved also in the three-dimensional case. In fact, observe that the unconditional well-posedness of the scheme (3.2) and the mass-conservation property (2.9) were proved independent on the dimension (see Lemma 4.1 and Theorem 4.2, respectively). Now, in order to analyze the convergence in the 3D case, we need to modify the inductive hypothesis (5.1) by the following one: there exists a positive constant $K>0$, independent of $m$, such that

[TABLE]

Then, we can prove error estimates for any solution of the scheme (3.2), with respect to a sufficiently regular solution of (2.12).

Theorem 6.1**.**

Assume that there exists a sufficiently regular solution of (2.12). If the inductive hypothesis (6.1) is satisfied, then the following error estimates hold

[TABLE]

[TABLE]

where the constant $C(T)>0$ is independent of $m,\Delta t$ and $h$.

Proof.

The proof follows as in Theorems 5.4 and 5.5, but in this case, in order to bound the terms in the estimates (5.24) and (5.51), we need to use the 3D interpolation inequality (2.5) and the inductive hypothesis (6.1), as follows:

[TABLE]

and

[TABLE]

Finally, the inductive hypothesis (6.1) can be verified in the same spirit of the two-dimensional case (see the proof of Theorem 5.4). ∎

7 Numerical simulations

In this section, we present two numerical experiments: the first one, is used to verify that our scheme gives a good approximation to chemotaxis phenomena in a liquid environment; and the second has been considered in order to check numerically the error estimates proved in our theoretical analysis. All the numerical results are computed by using the software Freefem++. We have considered the spaces for $\eta$, $c$, ${\bm{\sigma}}$, $\mathbf{u}$ and $\pi$, generated by $\mathbb{P}_{1},\mathbb{P}_{1},\mathbb{P}_{1},\mathbb{P}_{1}-bubble,\mathbb{P}_{1}$-continuous FE, respectively.

Test 1: In this experiment we consider the rectangular domain $\Omega=[0,2]\times[0,1],$ and the initial conditions

[TABLE]

where $s_{1}=0.2$, $s_{2}=0.5$ and $s_{3}=1.2$. The numerical solution is computed with mesh parameter $h=1/40$ and time step $\Delta t=1e-5$. Additionally, we consider the parameters values $\chi=8,$ $D_{c}=5,$ $\gamma=8,$ $D_{\mathbf{u}}=10$, $D_{n}=\rho=1$ and $\phi(x,y)=-1000y$. We show the simulations results for the times $t=0$, $t=12e-5$ and $t=30e-5$. The evolution results for the cell density and chemical signal are presented in Figure 7, while the evolution of the velocity field is shown in Figure 1. Initially, the cells are in two clusters in the upper part of our domain, and we observe that they begin to orient their movement in the direction of greater concentration of the chemical signal (in this case, the center of the domain). We can see how the clusters of organisms generate a kind of bridge between them and after, we see how organisms tend to agglomerate in the center of the rectangle. Moreover, we can see that the attraction signal $c_{h}$ is consumed by the organisms, and some changes in the velocity field are evidenced, influenced by the movement of the cells.

[TABLE]

\figcaption

Cell density vs Chemical concentration.

Test 2 (Convergence rates): In this experiment we take $\Omega=[0,1]\times[0,1]$ and we consider the exact solutions $\eta=e^{-t}(cos(2\pi x)+cos(2\pi y)+3)$, $c=e^{-t}(sin(2\pi y)+cos(2\pi x)-2\pi y+9)$, ${\bm{\sigma}}=\nabla c=2\pi e^{-t}(-sin(2\pi x),cos(2\pi y)-1)$, $\mathbf{u}=e^{-t}(sin(2\pi y)(-cos(2\pi x+\pi)-1),sin(2\pi x)(cos(2\pi y+\pi)+1))$ and $\pi=e^{-t}(cos(2\pi x)+sin(2\pi y))$; and all parameters in (3.2) equal to 1. Note that $\mathbf{u}={\bf 0}$ and $\frac{\partial c}{\partial\bm{\nu}}=\frac{\partial\eta}{\partial\bm{\nu}}=0$ on $\partial\Omega$, $\nabla\cdot\mathbf{u}=0$ in $\Omega$ and $\int_{\Omega}\pi=0$. Moreover, we use a uniform partition with $k+1$ nodes in each direction.

Numerical results are presented in Tables 7-7 for $\Delta t=2e-4$ with respect to time $T=0.01$. We observe the second-order convergence for the total errors in $e^{m}_{\eta},e^{m}_{c},e^{m}_{\mathbf{u}}$ in $l^{\infty}(L^{2})$-norm, and the first-order convergence for $e^{m}_{\eta},e^{m}_{c},e^{m}_{\mathbf{u}}$ in $l^{2}(H^{1})$-norm and $e^{m}_{\mathbf{u}}$ in $l^{\infty}(H^{1})$-norm, which is agreement with our theoretical analysis.

[TABLE]

\tabcaption

Convergence rates in space for $\eta$.

[TABLE]

\tabcaption

Convergence rates in space for $c$.

[TABLE]

\tabcaption

Convergence rates in space for $\mathbf{u}_{1}$ in weak norms.

[TABLE]

\tabcaption

Convergence rates in space for $\mathbf{u}_{2}$ in weak norms.

[TABLE]

\tabcaption

Convergence rates in space for $\mathbf{u}_{1}$ and $\mathbf{u}_{2}$ in strong norms.