# Full discretization of the time dependent Navier-Stokes equations with anisotropic slip boundary condition coupled with the convection-diffusion-reaction equation

Rim Aldbaissy, Nancy Chalhoub, Jules K Djoko, Toni Sayah

► To cite this version:

Rim Aldbaissy, Nancy Chalhoub, Jules K Djoko, Toni Sayah. Full discretization of the time dependent Navier-Stokes equations with anisotropic slip boundary condition coupled with the convection-diffusion-reaction equation. SeMA Journal: Boletín de la Sociedad Española de Matemática Aplicada, 2024, 82, pp.1-30. 10.1007/s40324-024-00355-7. http://dx.doi.org/10.1007/s40324-024-00355-7

HAL Id: hal-05055427

https://hal.science/hal-05055427v1

Submitted on 5 May 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Full discretization of the time dependent Navier-Stokes equations with anisotropic slip boundary condition coupled with the convection-diffusion-reaction equation

Rim Aldbaissy<fnref n="1" /> · Nancy Chalhoub<fnref n="1" /> · Jules K. Djoko<fnref n="2" /> · Toni Sayah<fnref n="3" />

Received: 30 June 2023 / Accepted: 4 March 2024

© The Author(s), under exclusive licence to Sociedad Española de Matemática Aplicada 2024

## Abstract

In this article, we consider the time dependent convection-diffusion-reaction equation coupled with the Navier-Stokes system subject partially to anisotropic slip boundary condition. We formulate the weak problem, and construct the weak solutions by using the approximation approach combined with some compactness results. The discrete solution is constructed by using the finite element method in space and the implicit Euler scheme in time. The convergence of the discrete solution is examined thanks to a priori errors estimation, and the rate of convergence is obtained. We propose an iterative scheme and investigate its convergence using the weak convergence approach. Finally, numerical investigations are performed to validate the theoretical results presented.

**Keywords** Power law slip boundary condition · Time dependent Navier-Stokes equations · Heat equation · Space-time discretization · Error estimates

**Mathematics Subject Classification** 65N30, 76M10, 35J85

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

Toni Sayah
toni.sayah@usj.edu.lb

Rim Aldbaissy
rim.dbaissy1@usj.edu.lb

Nancy Chalhoub
nancy.chalhoub@usj.edu.lb

Jules K. Djoko
Jules.DjokoKamdem@nwu.ac.za

<footnote n="1">1 Faculté des Ingénieurs, Université Saint-Joseph de Beyrouth, Beirut, Lebanon</footnote>

<footnote n="2">2 School of Mathematical Sciences, University of North West, Pochefstroom, South Africa</footnote>

<footnote n="3">3 Laboratoire de Mathématiques et Applications, Unité de recherche Mathématiques et Modélisation, Faculté des Sciences, Université Saint-Joseph de Beyrouth, Riad El Solh, B.P 11-514,, Beyrouth 1107 2050 Beirut, Lebanon</footnote>
# 1 Introduction

Let $\Omega$ be a bounded, smooth domain in $\mathbf{R}^d$, $d = 2, 3$. Let $[0, T]$ be an interval of $\mathbf{R}$. We assume that the boundary $\partial\Omega$ is Lipschitz-continuous. We are interested in the numerical solution of the following system of equations

$$ egin{cases} 
\mathbf{u}_t - 2\nu \operatorname{div} \mathbf{D}\mathbf{u} + [\mathbf{u} \cdot \nabla]\mathbf{u} + \nabla p = \mathbf{f}(C) & \text{in } \Omega \times (0, T), \\
C_t + (\mathbf{u} \cdot \nabla)C - \alpha \Delta C + r_0 C = g & \text{in } \Omega \times (0, T), \\
\operatorname{div} \mathbf{u} = 0 & \text{in } \Omega \times (0, T) \\
C = 0 & \text{on } \partial\Omega \times (0, T), \\
\mathbf{u}(\mathbf{x}, 0) = \mathbf{u}_0 & \text{in } \Omega \times \{0\} \\
C(\mathbf{x}, 0) = C_0 & \text{in } \Omega \times \{0\},
\end{cases} 
\quad (1.1) $$

where the unknowns are; the velocity $\mathbf{u}$, the pressure $p$ and the concentration $C$. The function $f$ (given as $f(\mathbf{x}, t, C(\mathbf{x}, t))$) represents the external force that is a function of the concentration $C$, while $g$ represents an external concentration source. $\nu$ is the viscosity of the fluid and $\alpha$ is the diffusivity. Both parameters are non negative and play very important role in the modeling of such phenomena. $\mathbf{u}_0$ and $C_0$ denote the initial velocity and concentration, respectively, and for now we assume that $\operatorname{div} \mathbf{u}_0 = 0$. We recall that the Cauchy stress tensor is given by $\mathbf{T} = -p\mathbf{I} + 2\nu D\mathbf{u}$, where $\mathbf{I}$ is the identity matrix in $\mathbf{R}^{d \times d}$, and the symmetric part of the velocity gradient is $2D\mathbf{u} = \nabla\mathbf{u} + (\nabla\mathbf{u})^T$. We are interested in Problem (1.1) when considering the position and the direction of the slip boundary condition (see [21, 22]). We then assume that the boundary $\partial\Omega$ of $\Omega$ is divided into two components $S$ and $\Gamma$, such that $\partial\Omega = S \cup \Gamma$, with $S \cap \Gamma = \emptyset$. We assume that the solution verifies a homogeneous Dirichlet condition on $\Gamma$, given by:

$$ \mathbf{u} = \mathbf{0} \quad \text{on } \Gamma. \quad (1.2) $$

Thus $\Gamma$ is the porous or artificial boundary where the fluid is prescribed. On $S$, we assume that the solution verifies an impermeability condition given by

$$ \mathbf{u} \cdot \mathbf{n} = 0 \quad \text{on } S, \quad (1.3) $$

where $\mathbf{n}: S \to \mathbf{R}^d$ is the normal outward unit vector to $S$. $S$ is an impermeable solid surface along which the fluid may slip. The power law slip boundary condition is given as follows (cf. [13])

$$ (\mathbf{T}\mathbf{n})_{\tau} + |K\mathbf{u}_{\tau}|^{s-2} K^2 \mathbf{u}_{\tau} = 0 \quad \text{on } S \times (0, T), \quad (1.4) $$

where $|\mathbf{v}|^2 = \mathbf{v} \cdot \mathbf{v}$ is the Euclidean norm and $(\mathbf{T}\mathbf{n})_{\tau}$ denotes the projection of the normal component of stress into the corresponding tangent plane. $K$ is an anisotropic tensor, uniformly positive definite, symmetric, and bounded. $s$ is a real, positive number representing the flow behavior index. This boundary condition may arise when the contact surface is lubricated with a thin layer of non-Newtonian fluid. It is evident that for $s = 2$ and $K = \mathbf{I}$, the classical Navier's slip condition is recovered. We are interested in Problem (1.1) with the boundary conditions given by, (1.2), (1.3), and (1.4). Our objective is to study the discretization of the problem using the Euler method in time and the finite element method in space. For interpretations and derivation of (1.4), we refer to the works [8, 11, 19–22, 25]. For the Stokes and Navier–Stokes with slip boundary condition we refer also to [23, 30–32, 34]. Regarding the coupling of the time-dependent Navier–Stokes system with the heat (or in general the convection-diffusion–reaction) equation, relevant references include [1, 2, 4, 12, 14] and the references therein. In [1, 4], the authors discuss the continuous and numerical solution
of a coupling with Dirichlet boundary condition. The approximation method presented in [2] is based on spectral discretization, while [4] employs the finite element method. In [3], we specifically address the time dependent Navier-Stokes system subject to (1.4). In this work, we extend the model treated in [4] by coupling it with the convection-diffusion-reaction problem. We show the existence of a solution, formulate the finite element approximation and establish its convergence. Then we propose an iterative scheme and provide corresponding numerical simulations. The outline of the paper is as follows:

* Sect. 2 introduces the classical notations and functional spaces relevant to the mathematical formulation and study of the boundary value problem (1.1), (1.2), (1.3) and (1.4).

* Sect. 3 presents a weak formulation for the problem (1.1), (1.2), (1.3), and (1.4) and discusses some mathematical properties.

* In Sect. 4, we introduce the discrete problem, recall their main properties, analyze a priori error estimates, and derive convergence. We also formulate an iterative scheme for the practical implementation of the nonlinear discrete problem.

* Sect. 5 presents the numerical results obtained from the simulations.

## 2 Preliminaries

In this section, we recall the main notations and the results that will be used later. We denote by $[L^p(\Omega)]^d$ the space of measurable vector functions $\mathbf{v}$ such that $|\mathbf{v}|^p$ is integrable. For $\mathbf{v} \in [L^p(\Omega)]^d$, the norm is defined by

$$ \| \mathbf{v} \|_{[L^p(\Omega)]^d} = \left( \int_{\Omega} |\mathbf{v}(\mathbf{x})|^p d\mathbf{x} \right)^{1/p} . $$

For $p = 2$, we denote the norm $||.||_{L^2(\Omega)}^d$ by $||.||$. We introduce the Sobolev space

$$ W^{m,r}(\Omega)^d = \{ \mathbf{v} \in [L^r(\Omega)]^d; \partial^k \mathbf{v} \in [L^r(\Omega)]^d, \forall |k| \le m \}, $$

where $k = (k_1, \dots, k_d)$ is a vector of non-negative integers, $|k| = k_1 + \dots + k_d$, and

$$ \partial^k \mathbf{v} = \frac{\partial^{|k|} \mathbf{v}}{\partial x_1^{k_1} \dots \partial x_d^{k_d}}. $$

This space is equipped with the semi-norm

$$ |\mathbf{v}|_{m,r,\Omega} = \left( \sum_{|k|=m} \int_{\Omega} |\partial^k \mathbf{v}|^r d\mathbf{x} \right)^{1/r}, $$

and is a Banach space for the norm

$$ \| \mathbf{v} \|_{m,r,\Omega} = \left( \sum_{\ell=0}^{m} |\mathbf{v}|_{\ell,r,\Omega}^r d\mathbf{x} \right)^{1/r} . $$

When $r = 2$, this space becomes the Hilbert space $H^m(\Omega)^d$. In particular, we consider the following spaces:

$$ H_0^1(\Omega)^d = \{ \mathbf{v} \in H^1(\Omega)^d, \mathbf{v}|_{\partial\Omega} = 0 \}, $$
equipped with the norm

$$|\mathbf{v}|_{H_0^1(\Omega)^d} = |\mathbf{v}|_{1,\Omega} = \left( \int_{\Omega} |\nabla \mathbf{v}|^2 d\mathbf{x} \right)^{1/2}.$$

The dual of $H_0^1(\Omega)^d$ is denoted by $H^{-1}(\Omega)^d$. We also introduce

$$L_0^2(\Omega) = \{q \in L^2(\Omega); \int_{\Omega} q(\mathbf{x}) d\mathbf{x} = 0\}.$$

For time-dependent problems, it is convenient to consider functions defined on a time interval $[a, b]$ with values in a separable functional space $W$ equipped with a norm $\| \cdot \|_W$. For all $r \ge 1$ we introduce the space

$$L^r(a, b; W) = \left\{ \mathbf{f} \text{ is measurable on } ]a, b[ \text{ and } \int_a^b \| \mathbf{f}(t) \|_W^r dt < \infty \right\},$$

equipped with the norm

$$\| \mathbf{f} \|_{L^r(a,b;W)} = \left( \int_a^b \| \mathbf{f}(t) \|_W^r dt \right)^{1/r}.$$

If $r = \infty$, then

$$L^\infty(a, b; W) = \left\{ \mathbf{f} \text{ is measurable on } (a, b) \text{ and } \sup_{t \in [a,b]} \| \mathbf{f}(t) \|_W < \infty \right\},$$

equipped with the norm

$$\| \mathbf{f} \|_{L^\infty(a,b;W)} = \sup_{t \in [a,b]} \| \mathbf{f}(t) \|_W.$$

**Remark 2.1** $L^r(a, b; W)$, $r \ge 1$, is Banach space if $W$ is a Banach space. In addition, we define $C^j(0, T; W)$, $j \ge 0$, as the space of functions that are of class $C^j$ in time and with values in $W$. We consider the following spaces:

$$X = \{ \mathbf{v} \in H^1(\Omega)^d; \text{ with } \mathbf{v}_{\Gamma} = \mathbf{0} \text{ and } (\mathbf{v} \cdot \mathbf{n})_S = \mathbf{0} \}, \quad M = L_0^2(\Omega), \quad Y = H_0^1(\Omega),$$

$$V = \{ \mathbf{v} \in X; \text{ div } \mathbf{v} = 0 \text{ in } \Omega \} \quad \text{and} \quad V^{\perp} = \{ \mathbf{v} \in X; \forall \mathbf{w} \in V, (\nabla \mathbf{v}, \nabla \mathbf{w}) = 0 \}.$$

**Lemma 2.2** Let $p$ be an integer such that $p \ge 1$ if $d = 1$ or $2$, or $1 \le p \le \frac{2d}{d-2}$ if $d \ge 3$. There exist two positive constants $S_p$ and $S_p^0$ such that (see [7])

$$\forall \mathbf{v} \in X, \quad \| \mathbf{v} \|_{L^p(\Omega)^d} \le S_p^0 |\mathbf{v}|_{1,\Omega},$$

and

$$\forall \mathbf{v} \in H^1(\Omega)^d, \quad \| \mathbf{v} \|_{L^p(\Omega)^d} \le S_p \| \mathbf{v} \|_{1,\Omega}.$$

**Lemma 2.3** For $d = 2$, Ladyzhenskaya's inequality states that (see [29]):

$$\| \mathbf{u} \|_{L^4(\Omega)^2} \le c(\Omega) \| \nabla \mathbf{u} \|^{1/2} \| \mathbf{u} \|^{1/2}. \qquad (2.1)$$

For $d = 3$, then Gagliardo-Nirenberg's inequality reads

$$\| \mathbf{u} \|_{L^4(\Omega)^3} \le c(\Omega) \| \mathbf{u} \|^{1/4} \| \nabla \mathbf{u} \|^{3/4}. \qquad (2.2)$$
**Lemma 2.4** [Korn's inequality (see [9])] There exists a positive constant c such that

$$ \forall \mathbf{v} \in X, \quad \int_{\Omega} |D(\mathbf{u})|^2 d\mathbf{x} \geq \mathbf{c} \int_{\Omega} |\nabla \mathbf{u}|^2 d\mathbf{x}. $$

Henceforth, we suppose the following hypothesis:

**Assumption 2.5** We assume that the data **f**, **g** and **v** satisfy the following conditions:

(i) **f** can be written as follows

$$ \mathbf{f}(x, t, C(x, t)) = \mathbf{f}_0(x, t) + \mathbf{f}_1(x, C(x, t)), $$

where $\mathbf{f}_0 \in C^0(0, T; L^2(\Omega)^d)$ and $\mathbf{f}_1$ is $c_{\mathbf{f}_1}^*$-Lipschitz with respect to its second argument from $\mathbf{R}$ with value in $\mathbf{R}^d$. In addition, we suppose that

$$ \forall x \in \Omega, \forall \xi_1, \xi_2 \in \mathbf{R}, \quad |\mathbf{f}_1(x, \xi_1) - \mathbf{f}_1(x, \xi_2)| \leq c_{\mathbf{f}_1} |\xi_1 - \xi_2|, \quad (2.3) $$

where $c_{\mathbf{f}_1}$ is a positive constant.

(ii) $g \in C^0(0, T; L^2(\Omega))$,

(iii) $\mathbf{u}_0 \in L^2(\Omega)^d$, div $\mathbf{u}_0 = 0$, $\mathbf{u}_0 \cdot \mathbf{n}_{\partial\Omega} = 0$, $C_0 \in L^2(\Omega)$ and $C_0|_{\partial\Omega} = 0$.

**Lemma 2.6** (Discrete Gronwall Lemma) [33, p. 294]. Let $(y_n)_n$, $(\tilde{f}_n)_n$ and $(\tilde{g}_n)_n$ be three positive sequences that verify:

$$ \forall n \in \mathbb{N}, \quad y_n \leq \tilde{f}_n + \sum_{k=0}^{n-1} \tilde{g}_k y_k. $$

We have:

$$ \forall n \in \mathbb{N}, \quad y_n \leq \tilde{f}_n + \sum_{k=0}^{n-1} \tilde{f}_k \tilde{g}_k \exp \left( \sum_{j=k}^{n-1} \tilde{g}_j \right). $$

In this work, if necessary, we will use the notation $\psi(t)$ for the function $\mathbf{x} \mapsto \psi(\mathbf{x}, \mathbf{t})$.

# 3 Analysis of the continuous problem

This section is devoted to the study of the continuous problem (1.1)-(1.4). We begin by introducing the corresponding weak formulation, and then construct a weak solution following the method presented in [24, 29].

## 3.1 Variational formulation

The weak formulation associated with (1.1), (1.2), (1.3) and (1.4) is standard and reads as follows:
Find $(\mathbf{u}, p, C) \in L^2(0, T; X) \times L^2(0, T; M) \times L^2(0, T; Y)$ such that, for a.e. $t$, $0 \le t \le T$,

$$ 
\left\{
\begin{array}{l}
\text{and for all } (\mathbf{v}, q, r) \in X \times M \times Y \\
\frac{d}{dt}(\mathbf{u}(t), \mathbf{v}) + a(\mathbf{u}(t), \mathbf{v}) + c_{\mathbf{u}}(\mathbf{u}(t), \mathbf{u}(t), \mathbf{v}) + \int_S |K\mathbf{u}_{\tau}(t)|^{s-2} K\mathbf{u}_{\tau}(t) \cdot K\mathbf{v}_{\tau} d\sigma \\
-(p(t), \operatorname{div}\mathbf{v}) = (\mathbf{f}(t, C(t)), \mathbf{v}), \\
\frac{d}{dt}(C(t), r) + c_C(\mathbf{u}(t), C(t), r) + \alpha(\nabla C(t), \nabla r) + r_0(C(t), r) = (g(t), r), \\
(\operatorname{div}\mathbf{u}(t), q) = 0, \\
\mathbf{u}(0) = \mathbf{u}_0 \text{ in } \Omega, \\
C(0) = C_0 \text{ in } \Omega.
\end{array}
\right.
\quad (3.1) $$

Here, we have defined the following terms:

$$ a(\mathbf{u}, \mathbf{v}) = 2\nu \int_{\Omega} D\mathbf{u} : D\mathbf{v} \, d\mathbf{x}, \quad c_{\mathbf{u}}(\mathbf{u}, \mathbf{u}, \mathbf{v}) = ((\mathbf{u} \cdot \nabla)\mathbf{u}, \mathbf{v}), \quad c_C(\mathbf{u}, \mathbf{C}, \mathbf{r}) = ((\mathbf{u} \cdot \nabla)\mathbf{C}, \mathbf{r}) $$

with $\mathbf{A}: \mathbf{B} = \sum_{1 \le i, j \le d} A_{ij} B_{ij}$. It is worth mentioning that by deriving (3.1), we make use of the identity:

$$ \sum_{1 \le i, j \le d} D_{ij} \mathbf{u} \frac{\partial v_i}{\partial x_j} = \sum_{1 \le i, j \le d} D_{ij} \mathbf{u} D_{ij} \mathbf{v}. $$

**Remark 3.1** It is straightforward to deduce that Problem (1.1)–(1.4) is equivalent to (3.1) in the distribution sense.

To study (3.1), it is convenient to recall the following monotonicity and continuity properties (see [16, 27]): there exists a constant $c$ such that for $1 \le s < 2$,

$$ (|\mathbf{y}| + |\mathbf{x}|)^{2-s} (|\mathbf{y}|^{s-2}\mathbf{y} - |\mathbf{x}|^{s-2}\mathbf{x}, \mathbf{y} - \mathbf{x}) \ge c|\mathbf{x} - \mathbf{y}|^2 \quad \forall \mathbf{x}, \mathbf{y} \in \mathbf{R}^d \quad (3.2) $$

and

$$ |\mathbf{x}|^{s-2}\mathbf{x} - |\mathbf{y}|^{s-2}\mathbf{y}| \le c|\mathbf{x} - \mathbf{y}|^{s-1}. \quad (3.3) $$

Now, we introduce the mapping $\mathbf{v} \mapsto \mathcal{A}(\mathbf{v})$ defined as follows:

$$ \forall \mathbf{w} \in V, \quad \langle \mathcal{A}(\mathbf{v}), \mathbf{w} \rangle = a(\mathbf{v}, \mathbf{w}) + \int_S |K\mathbf{v}_{\tau}|^{s-2} K\mathbf{v}_{\tau} \cdot K\mathbf{w}_{\tau} d\sigma. $$

We note that the map $\mathcal{A}$ does not contain the convection part $[\mathbf{u} \cdot \nabla]\mathbf{u}$.

From (3.2) and (3.3), we have the following Lemma (see [13]):

**Lemma 3.2** For $1 \le s < 2$, we have that

(a) $\mathcal{A}$ maps $V$ into its dual $V'$, and is bounded on all bounded subsets of $V$.

(b) For all $\mathbf{v}$, $\mathbf{u}$ elements of $V$

$$ \|\mathcal{A}(\mathbf{v}) - \mathcal{A}(\mathbf{u})\|_{V'} \le 2\nu \|\mathbf{v} - \mathbf{u}\|_1 + c \|K\|_{L^\infty(S)}^s \|\mathbf{u} - \mathbf{v}\|_1^{s-1}. $$

(c) The mapping $\mathbf{v} \to \mathcal{A}\mathbf{v}$ is strictly monotone from $V$ into $V'$:

$$ \text{For all } \mathbf{v}, \mathbf{w} \in V, \quad \langle \mathcal{A}(\mathbf{v}) - \mathcal{A}(\mathbf{w}), \mathbf{v} - \mathbf{w} \rangle \ge 0. $$

(d) $\mathcal{A}$ is hemi-continuous in $V$, i.e. for all $\mathbf{u}, \mathbf{v}$ in $V$, the mapping $t \mapsto \mathcal{A}(\mathbf{u} + t\mathbf{v})$ is continuous on $\mathbf{R}$ in $\mathbf{R}$.
## 3.2 Existence of solution

In this paragraph, we construct a solution $(\mathbf{u}, p, C)$ of (1.1) (or (3.1)) by adopting the approach employed in references [24, 29]. Our goal is to construct the solution by starting with semi-discretization in time and then taking the limit as the time step approaches zero. To achieve this, we introduce the following space:

$$H = \{ \mathbf{v} \in L^2(\Omega)^d : \operatorname{div} \mathbf{v} = 0 \text{ in } \Omega, (\mathbf{v} \cdot \mathbf{n})_{\partial\Omega} = \mathbf{0} \}.$$

We have that $V \subset H$, and the injection is compact.

Let $N > 1$ be an integer. We define the time step $k$ as $k = \frac{T}{N}$ and the subdivision points $t_n = nk$. For each $n \ge 1$, we approximate $g(t_n)$ and $\mathbf{f}(t_n, C(t_n))$ by their average defined almost everywhere in $\Omega$ by

$$g^n = g(t_n) \quad \text{and} \quad \mathbf{f}^n(C^n) = \mathbf{f}_0^n + \mathbf{f}_1^n(C^n),$$

where $\mathbf{f}_0^n = \mathbf{f}_0(t_n)$ and $\mathbf{f}_1^n(C^n) = \mathbf{f}_1(t_n, C^n)$. We set $\mathbf{u}^0 = \mathbf{u}_0$ and $C^0 = C_0$, and we introduce the following semi-discrete problem, which is exact in space and discrete in time:
Find sequence $(\mathbf{u}^n, p^n, C^n) \in X \times M \times Y$ such that

$$\left\{
\begin{array}{l}
\text{for all } (\mathbf{v}, r, q) \in X \times Y \times M, \\
\frac{1}{k}(\mathbf{u}^n - \mathbf{u}^{n-1}, \mathbf{v}) + a(\mathbf{u}^n, \mathbf{v}) + c_{\mathbf{u}}(\mathbf{u}^n, \mathbf{u}^n, \mathbf{v}) + \int_S |K\mathbf{u}_{\tau}^n|^{s-2} K\mathbf{u}_{\tau}^n \cdot K\mathbf{v}_{\tau} d\sigma \\
-(p^n, \operatorname{div} \mathbf{v}) = (\mathbf{f}^n(C^n), \mathbf{v}), \\
\frac{1}{k}(C^n - C^{n-1}, r) + \alpha(\nabla C^n, \nabla r) + c_C(\mathbf{u}^n, C^n, r) + r_0(C^n, r) = (g^n, r), \\
(\operatorname{div} \mathbf{u}^n, q) = 0.
\end{array}
\right.
\quad (3.4)$$

Given $\mathbf{u}^{n-1} \in X$ and $C^{n-1} \in Y$, (3.4) is a steady Navier-Stokes problem.

Problem (3.4) can be written equivalently as the following:

Find sequences $(\mathbf{u}^n, C^n, p^n) \in V \times Y$ such that:

$$\left\{
\begin{array}{l}
\displaystyle \frac{1}{k}(\mathbf{u}^n - \mathbf{u}^{n-1}, \mathbf{v}) + a(\mathbf{u}^n, \mathbf{v}) + c_{\mathbf{u}}(\mathbf{u}^n, \mathbf{u}^n, \mathbf{v}) + \int_S |K\mathbf{u}_{\tau}^n|^{s-2} K\mathbf{u}_{\tau}^n \cdot K\mathbf{v}_{\tau} d\sigma \\
\displaystyle = (\mathbf{f}^n(C^n), \mathbf{v}) \quad \forall \mathbf{v} \in V, \\
\displaystyle \frac{1}{k}(C^n - C^{n-1}, r) + \alpha(\nabla C^n, \nabla r) + c_C(\mathbf{u}^n, C^n, r) + r_0(C^n, r) = (g^n, r) \quad \forall r \in Y.
\end{array}
\right.
\quad (3.5)$$

**Theorem 3.3** At each time step $n$ and for a given $\mathbf{u}^{n-1} \in X$ and $C^{n-1} \in Y$, Problem (3.4) admits at least one solution $(\mathbf{u}^n, p^n, C^n) \in X \times M \times Y$.

*Proof* Let $(\mathbf{u}^{n-1}, C^{n-1})$ be a given element of $X \times Y$. We will show that Problem (3.5) admits at least one solution $(\mathbf{u}^n, C^n) \in V \times Y$. To obtain the pressure $p^n$, we use the following inf-sup condition:

There exists $\beta > 0$ such that

$$\inf_{q \in M} \sup_{0 \neq v \in X} \frac{b(\mathbf{v}, q)}{\|\mathbf{v}\|_X \|q\|} \geq \beta.$$

Since $V$ and $Y$ are separable, there exist an increasing sequence $(V_m)_m$ of finite dimensional subspaces of $V$ and an increasing sequence $(Y_m)_m$ of finite dimensional subspaces of $Y$ such
that $\bigcup_{m \in \mathbb{N}} V_m \times Y_m$ is dense in $V \times Y$. Thus, we can approximate problem (3.5) with the following:

Find $(\mathbf{w}_m, \xi_m) \in V_m \times Y_m$ such that

$$ egin{cases} 
\frac{1}{k}(\mathbf{w}_m - \mathbf{u}^{n-1}, \mathbf{v}) + a(\mathbf{w}_m, \mathbf{v}) + c_{\mathbf{u}}(\mathbf{w}_m, \mathbf{w}_m, \mathbf{v}) + \int_S |K \mathbf{w}_{m\tau}|^{s-2} K \mathbf{w}_{m\tau} \cdot K \mathbf{v}_{\tau} d\sigma \\
= (\mathbf{f}^n(\xi_m), \mathbf{v}) \quad \forall \mathbf{v} \in V_m, \\
\frac{1}{k}(\xi_m - C^{n-1}, r) + \alpha(\nabla\xi_m, \nabla r) + c_C(\mathbf{w}_m, \xi_m, r) + r_0(\xi_m, r) = (g^n, r) \quad \forall r \in Y_m.
\end{cases} 
\quad (3.6) $$

We consider the first equation of the last system for a given $\xi_m \in V_m$. It can be written as: $\forall \mathbf{v} \in V_m$,

$$ egin{aligned} & \frac{1}{k}(\mathbf{w}_m, \mathbf{v}) + a(\mathbf{w}_m, \mathbf{v}) + c_{\mathbf{u}}(\mathbf{w}_m, \mathbf{w}_m, \mathbf{v}) + \int_S |K \mathbf{w}_{m\tau}|^{s-2} K \mathbf{w}_{m\tau} \cdot K \mathbf{v}_{\tau} d\sigma \\ & = (\mathbf{f}^n(\xi_m), \mathbf{v}) + \frac{1}{k}(\mathbf{u}^{n-1}, \mathbf{v}). \end{aligned} 
\quad (3.7) $$

We note that for all $\mathbf{v}_m \in V_m$, we have $c_{\mathbf{u}}(\mathbf{v}_m, \mathbf{v}_m, \mathbf{v}_m) = 0$ due to the homogeneous slip condition on $S$ (and the no-slip condition on $\Gamma$). By referring to the proof of Proposition 2.2 in [13] (based on Brouwer's Fixed-Point Theorem) we can deduce that (3.7) has at least one solution $\mathbf{w}(\xi_m) \in V_m$. Hence, the second equation of System (3.6) can be written as follows: $\forall r \in Y_m$,

$$ \frac{1}{k}(\xi_m, r) + \alpha(\nabla\xi_m, \nabla r) + c_C(\mathbf{w}(\xi_m), \xi_m, r) + r_0(\xi_m, r) = (g^n, r) + \frac{1}{k}(C^{n-1}, r). \quad (3.8) $$

We observe that the only unknown variable in the last equation is $\xi_m$. Since the non-linear form $c_C(\mathbf{w}(\xi_m), \xi_m, r) = ((\mathbf{w}(\xi_m) \cdot \nabla)\xi_m, r)$ satisfies $c_C(\mathbf{w}(\xi_m), \xi_m, \xi_m) = 0$, we can apply Brouwer's Fixed-Point Theorem to conclude that Eq. (3.8) admits at least one solution $\xi_m \in Y_m$. Finally we deduce that System (3.6) has at least one solution $(\mathbf{w}_m, \xi_m) \in V_m \times Y_m$. Furthermore, by setting $r = \xi_m$ in the second equation and $\mathbf{v} = \mathbf{w}_m$ in the first equation of System (3.6) and applying the relation $ab \le \frac{2\epsilon}{a} + \frac{\epsilon}{2}b^2$, we obtain the following bounds:

$$ \frac{1}{2k} \|\xi_m\|^2 + \alpha \|\nabla\xi_m\|^2 + r_0 \|\xi_m\|^2 \le k \|g^n\|^2 + \frac{1}{k} \|C^{n-1}\|^2 $$

and

$$ \frac{1}{2k} \|\mathbf{w}_m\|^2 + c \|\nabla \mathbf{w}_m\|^2 + \int_S |K \mathbf{w}_{m\tau}|^s d\sigma \le k \|\mathbf{f}_0^n\|^2 + kc_1^2 \|\xi_m\|^2. $$

Hence, there exist subsequences, still labelled by $\mathbf{w}_m$ and $\xi_m$, such that

$$ \mathbf{w}_m \to \mathbf{u}^n \quad \text{weakly in } H^1(\Omega)^d, $$

$$ \xi_m \to \xi^n \quad \text{weakly in } H^1(\Omega). $$

Now, we proceed to take the limit in the system (3.6) as $m$ tends towards infinity. The limit of all the terms of the first equation of system (3.6), except for the term $(\mathbf{f}^n(\xi_m), \mathbf{v})$, are treated in the proof of Proposition 2.2 [13]. By applying the properties of $\mathbf{f}$, we can easily deduce that this term $(\mathbf{f}^n(\xi^n), \mathbf{v})$ converges for all $\mathbf{v} \in Y$. The same arguments used for the convergence of the first equation of system (3.6) can be applied to deduce the convergence of the second. Therefore, we conclude that $(\mathbf{u}^n, C^n)$ satisfies the system (3.5) and this completes the proof of the theorem. □
We next derive bounds on the expression $\mathbf{u}^n$. The following proposition provides basic uniform a priori estimates for each solution of (3.4).

**Proposition 3.4** Each solution $(\mathbf{u}^n, p^n, C^n)$ of Problem (3.4) satisfies the following uniform a priori estimates: For each $m \ge 1$,

$$ \begin{aligned} & \| C^m \|^2 + \sum_{n=1}^{m} \| C^n - C^{n-1} \|^2 + \sum_{n=1}^{m} k |C^n|_X^2 + r_0 \sum_{n=1}^{m} k \|C^n\|^2 \\ & \le c_1 \left( \sum_{n=1}^{m} k \| g^n \|^2 + \| C^0 \|^2 \right), \end{aligned} \quad (3.9) $$ 
and 
$$ \begin{aligned} & \| \mathbf{u}^m \|^2 + \sum_{n=1}^{m} \| \mathbf{u}^n - \mathbf{u}^{n-1} \|^2 + \sum_{n=1}^{m} k |\mathbf{u}^n|_X^2 + \sum_{n=1}^{m} k \int_S |\mathbf{u}_\tau^n|^s d\sigma \\ & \le c_2 \left( \sum_{n=1}^{m} k \| \mathbf{f}_0^n \|^2 + \sum_{n=1}^{m} k \| g^n \|^2 + \| \mathbf{u}^0 \|^2 + \| C^0 \|^2 \right), \end{aligned} $$

where $c_1$ and $c_2$ are positive constants independent of the time and $k$.

**Proof** We consider the second equation of system (3.4) multiplied by the time step $k$, take $r = C^n$, use the relation $(a - b, a) = \frac{1}{2}|a|^2 - \frac{1}{2}|b|^2 + \frac{1}{2}|a - b|^2$ and the relation $ab \le \frac{1}{2\epsilon}a^2 + \frac{\epsilon}{2}b^2$ for the right-hand side $(g^n, r)$ with a suitable choice of the parameter $\epsilon$. Finally, summing over $n = 1, \dots, m$ yields the first bound of (3.9).

Next, we consider the first equation of system (3.4) multiplied by the time step $k$, take $\mathbf{v} = \mathbf{u}^n$, use the fact that the tensor $K$ is bounded, apply the above relations for the right-hand side $(\mathbf{f}^n(C^n), \mathbf{v})$ use relation (2.3), and then use the first bound of (3.9) to obtain the second inequality of (3.9) after summing over $n$ as described above. □

We define the piecewise linear functions in time as follows:

$$ \begin{aligned} \forall t \in [t_{n-1}, t_n], \quad \mathbf{u}_k(t) &= \mathbf{u}_{n-1} + \frac{t - t_{n-1}}{k} (\mathbf{u}_n - \mathbf{u}_{n-1}), \\ C_k(t) &= C_{n-1} + \frac{t - t_{n-1}}{k} (C_n - C_{n-1}), \quad 0 < n \le N. \end{aligned} \quad (3.10) $$

We also define the step functions

$$ \forall t \in [t_{n-1}, t_n], \quad g_k(t) = g^n, \quad \mathbf{f}_k(t) = \mathbf{f}^n, \quad 0 < n \le N, $$

and

$$ \forall t \in [t_{n-1}, t_n], \quad \mathbf{w}_k(t) = \mathbf{u}^n, \quad \xi_k(t) = C^n, \quad 0 < n \le N. $$

We have the following convergence theorem:

**Proposition 3.5** There exist functions $\mathbf{u} \in L^2(0, T; V) \cap L^\infty(0, T, H)$ and $C \in L^2(0, T; Y) \cap L^\infty(0, T, L^2(\Omega))$ such that subsequences of $k$, still labelled by $k$, satisfy:

$$ \begin{aligned} \lim_{k \to 0} \mathbf{u}_k &= \lim_{k \to 0} \mathbf{w}_k = \mathbf{u} \quad \text{weakly} * \text{ in } L^\infty(0, T; H), \\ \lim_{k \to 0} C_k &= \lim_{k \to 0} \xi_k = C \quad \text{weakly} * \text{ in } L^\infty(0, T; L^2(\Omega)) \end{aligned} $$
and

$$ \lim_{k \to 0} \mathbf{u}_k = \lim_{k \to 0} \mathbf{w}_k = \mathbf{u} \text{ weakly in } L^2(0, T; V), $$

$$ \lim_{k \to 0} C_k = \lim_{k \to 0} \xi_k = C \text{ weakly in } L^2(0, T; Y). $$

**Proof** Relations (3.9) allow us to deduce that $(\mathbf{u}_k)_k$ and $(C_k)$ are respectively uniformly bounded in $L^\infty(0, T; H) \cap L^2(0, T; V)$ and in $L^\infty(0, T; L^2(\Omega)) \cap L^2(0, T; Y)$, since we have:

$$ \lim_{k \to 0} \mathbf{u}_k = \mathbf{u} \text{ weakly } \ast \text{ in } L^\infty(0, T; H), $$

$$ \lim_{k \to 0} C_k = C \text{ weakly } \ast \text{ in } L^\infty(0, T; L^2(\Omega)), $$

$$ \lim_{k \to 0} \mathbf{u}_k = \mathbf{u} \text{ weakly in } L^2(0, T; V), $$

$$ \lim_{k \to 0} C_k = C \text{ weakly in } L^2(0, T; Y), $$

$$ \lim_{k \to 0} \mathbf{w}_k = \mathbf{w} \text{ weakly } \ast \text{ in } L^\infty(0, T; H), $$

$$ \lim_{k \to 0} \xi_k = \xi \text{ weakly } \ast \text{ in } L^\infty(0, T; L^2(\Omega)), $$

$$ \lim_{k \to 0} \mathbf{w}_k = \mathbf{w} \text{ weakly in } L^2(0, T; V), $$

$$ \lim_{k \to 0} \xi_k = \xi \text{ weakly in } L^2(0, T; Y). $$

As far as the function $\mathbf{w}$ is concerned, observe that

$$ \forall t \in ]t_{n-1}, t_n], \quad \mathbf{w}_k(t) - \mathbf{u}_k(t) = \frac{t_n - t}{k}(\mathbf{u}_n - \mathbf{u}_{n-1}), \\ \xi_k(t) - C_k(t) = \frac{t_n - t}{k}(C_n - C_{n-1}), \quad 0 < n \le N. $$

Therefore, we have

$$ \| \mathbf{w}_k - \mathbf{u}_k \|_{L^2(0,T;L^2(\Omega)^d)}^2 = \frac{k}{3} \sum_{n=1}^{N} \| \mathbf{u}^n - \mathbf{u}^{n-1} \|_{L^2(\Omega)^d}, $$

$$ \| \xi_k - C_k \|_{L^2(0,T;L^2(\Omega))}^2 = \frac{k}{3} \sum_{n=1}^{N} \| \mathbf{u}^n - \mathbf{u}^{n-1} \|_{L^2(\Omega)^d}. $$

Then, the relation (3.9) implies that the terms $\sum_{1 \le n \le N} \|\mathbf{u}^n - \mathbf{u}^{n-1}\|_{L^2(\Omega)^d}$ and $\sum_{1 \le n \le N} \|C^n - C^{n-1}\|_{L^2(\Omega)}$ are uniformly bounded, and we obtain:

$$ \lim_{k \to 0} \| \mathbf{w}_k - \mathbf{u}_k \|_{L^2(0,T;L^2(\Omega)^d)}^2 = 0, \\ \lim_{k \to 0} \| \xi_k - C_k \|_{L^2(0,T;L^2(\Omega))}^2 = 0. \qquad (3.11) $$

The uniqueness of the limit allows us to obtain $\mathbf{u} = \mathbf{w}$ and $\xi = C$, thereby deducing the convergence results. □

To pass to the limit to the system (3.4), we need to prove the strong convergence of $(\mathbf{u}_k)$ and $(\mathbf{w}_k)$. We need also to introduce the following spaces (see [29, page 320]):

$$ \mathcal{V} = \{ \mathbf{v} \in \mathcal{D}^d(\Omega); \operatorname{div} \mathbf{v} = 0 \}. $$

Let $V_s$ be the closure of $\mathcal{V}$ in $X \cap H^{3/2}(\Omega)^d$, and let $Y_s$ be the closure of $\mathcal{D}(\Omega)$ in $Y \cap H^{3/2}(\Omega)$. $V_s$ and $Y_s$ are respectively endowed with the usual Hilbert norm of $H^{3/2}(\Omega)^d$ and $H^{3/2}(\Omega)$. We note that $V_s \subset V$.
**Theorem 3.6** Under the assumptions of Proposition 3.5, there exists $\xi \in L^2(0, T; V')$ and subsequences, still labelled by $k$, such that

$$ \lim_{k \to 0} \frac{d\mathbf{u}_k}{dt} = \frac{d\mathbf{u}}{dt} \quad \text{weakly in } L^2(0, T; V'_s), \quad (3.12) $$

$$ \lim_{k \to 0} \frac{dC_k}{dt} = \frac{dC}{dt} \quad \text{weakly in } L^2(0, T; Y'_s), \quad (3.13) $$

$$ \lim_{k \to 0} \mathbf{u}_k = \mathbf{u} \quad \text{strongly in } L^2(0, T; H), \quad (3.14) $$

$$ \lim_{k \to 0} C_k = C \quad \text{strongly in } L^2(0, T; L^2(\Omega)), \quad (3.15) $$

$$ \lim_{k \to 0} \mathbf{w}_k = \mathbf{u} \quad \text{strongly in } L^2(0, T; H), \quad (3.16) $$

$$ \lim_{k \to 0} \xi_k = C \quad \text{strongly in } L^2(0, T; L^2(\Omega)), \quad (3.17) $$

and

$$ \lim_{k \to 0} \mathcal{A}(\mathbf{w}_k) = \xi \quad \text{weakly in } L^2(0, T; V'). \quad (3.18) $$

**Proof** The system (3.4) can be written as follows: $\forall \mathbf{v} \in V$ and $\forall r \in Y$,

$$ \left( \frac{d\mathbf{u}_k}{dt}, \mathbf{v} \right) = (\mathbf{f}_k(\xi_k), \mathbf{v}) - a(\mathbf{w}_k, \mathbf{v}) - c_{\mathbf{u}}(\mathbf{w}_k, \mathbf{w}_k, \mathbf{v}) - \int_S |K\mathbf{w}_{k, \tau}|^{s-2} K\mathbf{w}_{k, \tau} \cdot K\mathbf{v}_{\tau} d\sigma \quad (3.19) $$

and

$$ \left( \frac{dC_k}{dt}, r \right) = (g_k, r) - \alpha(\nabla \xi_k, \nabla r) - c_C(\mathbf{w}_k, \xi_k, r) - r_0(\xi_k, r). \quad (3.20) $$

Using the fact that $\mathcal{A}$ maps $V$ into its dual $V'$ (see Lemma 3.2) and [29, Lemma 4.6, page 327], the following result is obtained: the terms

$$ \sum_{n=1}^{N} k \left\| \frac{\mathbf{u}^n - \mathbf{u}^{n-1}}{k} \right\|_{V'_s}^2 \quad \text{and} \quad \sum_{n=1}^{N} k \left\| \frac{C^n - C^{n-1}}{k} \right\|_{Y'_s}^2 \quad (3.21) $$

are bounded independently of $k$. In fact, we will give the steps to bound the first term of (3.21) and the second one can be bounded similarly. We recall first the following property (see [29, Lemma 4.1, page 320]):

$$ |c_{\mathbf{u}}(\mathbf{u}^n, \mathbf{u}^n, \mathbf{v})| \le c \| \mathbf{u}^n \| \| \mathbf{u}^n \|_{1,\Omega} \| \mathbf{v} \|_{V'_s}. \quad (3.22) $$

Then, we consider relation (3.19) with $\mathbf{v} \in V_s$, replace $\mathbf{u}_k$ by (3.10), use the first property of Lemma 3.2 and the relation (3.9) and (3.22) to get the following (as $V_s \subset V$): for $0 < n \le N$,

$$ \| \frac{\mathbf{u}^n - \mathbf{u}^{n-1}}{k} \|_{V'_s} \le c \left( (S_2^0 \| \mathbf{f}_0^n \| + c_{\mathbf{f}_1} ||C^n||) + \| \mathcal{A}(\mathbf{u}^n) \|_{V'} + \| \mathbf{u}^n \|_X \right). $$

By raising to the power two, multiplying by $k$ and summing over $n$ from 1 to $N$, we get the following bound:

$$ \sum_{n=1}^{N} k \left\| \frac{\mathbf{u}^n - \mathbf{u}^{n-1}}{k} \right\|_{V'_s}^2 \le c \sum_{n=1}^{N} k \left( \| \mathbf{f}_0^n \|^2 + ||C^n||^2 + \| \mathcal{A}(\mathbf{u}^n) \|_{V'}^2 + \| \mathbf{u}^n \|_X^2 \right). \quad (3.23) $$
The first two terms of the last inequality can be bounded by using the fact that $f_0 \in C^0(0, T; L^2(\Omega)^d)$ and relation (3.9). The last term can be bounded by using (3.9). It remains to bound the third term of the right hand side of (3.23). By using Lemma 3.2 (see [13]) we get by using the Young-inequality $ab \le \frac{a^p}{p} + \frac{b^q}{q}$, $p = 1/(s-1)$, $q = p/(p-1)$,

$$ \begin{aligned} \| A(\mathbf{u}^n) \|_{V'}^2 &\le c ( \| \mathbf{u}^n \|_V^2 + \| K \|_{L^\infty(\Gamma)}^{2s} \| \mathbf{u}^n \|_V^{2(s-1)} ) \\ &\le c_1 ( \| \mathbf{u}^n \|_V^2 + \| K \|_{L^\infty(\Gamma)}^{2s} \left( (s-1) \| \mathbf{u}^n \|_V^2 + \frac{p-1}{p} \right) ) \end{aligned} \quad (3.24) $$

By regrouping all the above bounds, relation (3.23) gives that $\sum_{n=1}^N k \| \frac{\mathbf{u}^n - \mathbf{u}^{n-1}}{k} \|_{V'}^2$ is bounded independently of $k$. Similarly we treat the term $\sum_{n=1}^N k \| \frac{C^n - C^{n-1}}{k} \|_{Y_s'}^2$ and the steps are more simple, the main idea is based on relation (3.9) and the following relation (similar to relation (3.22))

$$ c_C(\mathbf{u}^n, C^n, r) \le \| \mathbf{u}^n \| \| C^n \|_{1,\Omega} \| r \|_{Y_s'} $$

Consequently, it follows that $\frac{d\mathbf{u}_k}{dt}$ and $\frac{dC_k}{dt}$ are respectively bounded in $L^2(0, T; V_s')$ and $L^2(0, T; Y_s')$ and there exists a subsequence of $k$, still labelled by $k$ such that

$$ \lim_{k \to 0} \frac{dC_k}{dt} = \frac{dC_1}{dt} \quad \text{weakly in } L^2(0, T; Y_s') $$

and

$$ \lim_{k \to 0} \frac{d\mathbf{u}_k}{dt} = \frac{d\mathbf{u}_1}{dt} \quad \text{weakly in } L^2(0, T; Y_s'). $$

Relation (3.11) implies that $\mathbf{u}_1 = \mathbf{u}$, which in turn leads to (3.12), and also $C_1 = C$, resulting in (3.13). To prove (3.14) and (3.15), it is sufficient to apply [29, Theorem 2.1, page 271]. Furthermore, the relation (3.16) and (3.17) are simple consequence of (3.11). Finally, Lemma 3.2 and relation (3.24) allows us to deduce that there exists a subsequence of $k$, still labelled by $k$ such that $A(\mathbf{w}_k)$ satisfies 3.18 and converge weakly to $\xi$ in $L^2(0, T; V')$. $\square$

**Remark 3.7** We have the following properties: (see [29])

1. Using Lemma 2.3, Relation (3.22) can be replaced by

$$ |c_{\mathbf{u}}(\mathbf{u}^n, \mathbf{u}^n, \mathbf{v})| = |c_{\mathbf{u}}(\mathbf{u}^n, \mathbf{v}, \mathbf{u}^n)| \le c \| \nabla \mathbf{u}^n \|^2 \| \nabla \mathbf{v} \|_V, $$

and we get by using the relations (3.19), (3.9) and (3.24) that

$$ \lim_{k \to 0} \frac{d\mathbf{u}_k}{dt} = \frac{d\mathbf{u}}{dt} \quad \text{weakly in } L^1(0, T; V'). $$

Similarly the following relation

$$ \begin{aligned} |c_C(\mathbf{u}^n, C^n, r)| &= |c_C(\mathbf{u}^n, r, C^n)| \le c \| \nabla \mathbf{u}^n \| \| \nabla C^n \| \| \nabla r \|_Y \\ &\le \frac{c}{2} ( \| \nabla \mathbf{u}^n \|^2 + \| \nabla C^n \|^2 ) \| \nabla r \|_Y \end{aligned} $$

allows us to get

$$ \lim_{k \to 0} \frac{dC_k}{dt} = \frac{dC}{dt} \quad \text{weakly in } L^1(0, T; Y'). $$
2. in 2D: by using Lemma 2.3 and Relations (3.9) we have

$$ 
\begin{aligned}
|c_{\mathbf{u}}(\mathbf{u}^n, \mathbf{u}^n, \mathbf{v})| &= |c_{\mathbf{u}}(\mathbf{u}^n, \mathbf{v}, \mathbf{u}^n)| \\
&\le c \| \mathbf{u}^n \|_{L^4(\Omega)^d}^2 \| \mathbf{v}^n \|_{1, \Omega} \\
&\le c \| \mathbf{u}^n \| \| \nabla \mathbf{u}^n \| \| \mathbf{v}^n \|_V \\
&\le c_1 \| \nabla \mathbf{u}^n \| \| \mathbf{v}^n \|_V,
\end{aligned}
$$

which allows us, by using relation (3.9), (3.19) and (3.24), to deduce that relation (3.12) can be replaced by the following

$$ \lim_{k \to 0} \frac{d\mathbf{u}_k}{dt} = \frac{d\mathbf{u}}{dt} \quad \text{weakly in } L^2(0, T; V'). $$

We can easily show by the same way that

$$ \lim_{k \to 0} \frac{dC_k}{dt} = \frac{dC}{dt} \quad \text{weakly in } L^2(0, T; Y'). $$

3. in 3D: by using Lemma 2.3 and Relations (3.9) we have

$$ 
\begin{aligned}
|c_{\mathbf{u}}(\mathbf{u}^n, \mathbf{u}^n, \mathbf{v})| &= |c_{\mathbf{u}}(\mathbf{u}^n, \mathbf{v}, \mathbf{u}^n)| \\
&\le c \| \mathbf{u}^n \|_{L^4(\Omega)^d}^2 \| \mathbf{v}^n \|_V \\
&\le c \| \mathbf{u}^n \|^{1/2} \| \nabla \mathbf{u}^n \|^{3/2} \| \mathbf{v}^n \|_V \\
&\le c_1 \| \nabla \mathbf{u}^n \|^{3/2} \| \mathbf{v}^n \|_{1, \Omega}
\end{aligned}
$$

which allows us, by using relation (3.9), (3.19) and (3.24), to deduce the following bound:

$$ \lim_{k \to 0} \frac{d\mathbf{u}_k}{dt} = \frac{d\mathbf{u}}{dt} \quad \text{weakly in } L^{4/3}(0, T; V'). $$

Similar steps give the following

$$ \lim_{k \to 0} \frac{dC_k}{dt} = \frac{dC}{dt} \quad \text{weakly in } L^{4/3}(0, T; Y'). $$

**Theorem 3.8** Under assumption 2.5, The Navier-Stokes system (1.1) dmits at least one solution in the following sense: find $(\mathbf{u}, C)$ satisfying,

$$ \mathbf{u} \in L^2(0, T; V) \cap L^\infty(0, T; H), \quad C \in L^2(0, T; Y) \cap L^\infty(0, T; L^2(\Omega)), $$

$$ 
\left\{
\begin{array}{ll}
\displaystyle \frac{d}{dt}(\mathbf{u}, \mathbf{v}) + a(\mathbf{u}, \mathbf{v}) + c_{\mathbf{u}}(\mathbf{u}, \mathbf{u}, \mathbf{v}) + \int_S |K\mathbf{u}_{\tau}|^{s-2} K\mathbf{u}_{\tau} \cdot K\mathbf{v}_{\tau} d\sigma = (f(C), \mathbf{v}), & \forall \mathbf{v} \in V_s \\
\\
\displaystyle \frac{d}{dt}(C, r) + c_C(\mathbf{u}, C, r) + \alpha(\nabla C, \nabla r) + r_0(C, r) = (g, r), & \forall r \in Y_s \\
\\
\mathbf{u}(0) = \mathbf{u}_0, \\
C(0) = C_0.
\end{array}
\right.
\qquad (3.25) $$

**Proof** To show the existence of a solution to the Navier-Stokes system (1.1), we use the results obtained thus far to pass to the limit with the system (3.4), which can be written as:

$$ \left( \frac{d\mathbf{u}_k}{dt}, \mathbf{v} \right) + \langle A(\mathbf{w}_k), \mathbf{v} \rangle + c_{\mathbf{u}}(\mathbf{w}_k, \mathbf{w}_k, \mathbf{v}) = (f_k(\xi_k), \mathbf{v}), \quad \forall \mathbf{v}_s \in V_s, \qquad (3.26) $$

$$ \left( \frac{dC_k}{dt}, r \right) + \alpha(\nabla \xi_k, \nabla r) + c_C(\mathbf{w}_k, \xi_k, r) + r_0(\xi_k, r) = (g_k, r), \quad \forall r \in Y_s. \qquad (3.27) $$

The convergence of Eq. (3.26) is proved in [29, page 328]) and [24, page 159] (see also [3]), where the function $f$ is independent of the concentration. Here, the fact that $f$ satisfies Assumption 2.5 allows us to easily establish the convergence of the right-hand side of Eq. (3.26). The convergence of Eq. (3.27) is relatively straightforward and relies on similar arguments. Hence we get the existence of at least one solution of problem (3.1). □
**Theorem 3.9** If the Assumption 2.5 holds, every solution of (3.1) verifies the bounds

$$ \| C \|_{L^{\infty}(0,T;L^2(\Omega))} + \| C \|_{L^2(0,T;Y)} \leq c_1 (\| g \|_{L^2(0,T;L^2(\Omega))} + \| C_0 \|) \quad (3.28) $$

and

$$ \begin{aligned} & \| \mathbf{u} \|_{L^{\infty}(0,T;L^2(\Omega)^d)} + \| \mathbf{u} \|_{L^2(0,T;X)} \\ & \leq c_2 ( \| \mathbf{f}_0 \|_{L^2(0,T;L^2(\Omega)^d)} + \| \mathbf{u}_0 \|_{L^2(\Omega)^d} + \| g \|_{L^2(0,T;L^2(\Omega))} + \| C_0 \| ) , \end{aligned} \quad (3.29) $$

where $c_1$ and $c_2$ are positive constants independent of $(\mathbf{u}, p, C)$.

**Proof** Let $(\mathbf{u}, p, C)$ be a solution of (3.1). To prove the bound (3.28), we choose $r = C$ in the second equation of problem (3.1). We apply inequality $ab \leq \frac{\varepsilon}{2}a^2 + \frac{1}{2\varepsilon}b^2$ to obtain

$$ \begin{aligned} & \frac{1}{2} \frac{d}{dt} \| C(t) \|^2 + \alpha |C(t)|_Y^2 + r_0 \|C(t)\|^2 \leq S_2^0 \| g(t) \| |C(t)|_Y \\ & \leq \frac{(S_2^0)^2 \varepsilon}{2} \| g(t) \|^2 + \frac{1}{2\varepsilon} |C(t)|_Y^2. \end{aligned} $$

For $\varepsilon = \frac{1}{\alpha}$, we obtain:

$$ \frac{d}{dt} \| C(t) \|^2 + \alpha |C(t)|_Y^2 \leq \frac{(S_2^0)^2}{\alpha} \| g(t) \|^2 . \quad (3.30) $$

We integrate with respect to $t$ between 0 and $T$ to derive the following bound:

$$ \|C\|_{L^{\infty}(0,T;L^2(\Omega))}^2 + \alpha \|C\|_{L^2(0,T;Y)}^2 \leq \frac{(S_2^0)^2}{\alpha} \|g\|_{L^2(0,T;L^2(\Omega))}^2 + \|C_0\|^2. \quad (3.31) $$

We note that relations (3.31) is reformulation of (3.28), whereas the relation (3.29) can be obtained by following the same steps as above and using Assumption 2.5. $\square$

We study next, the circumstances under which (3.1) has a unique solution. We state that

**Theorem 3.10** Let $\Omega$ be a bounded, connected domain of $\mathbb{R}^d$ with a Lipschitz-continuous boundary. Assume that Assumption 2.5 holds. Assume that the solution $(\mathbf{u}, p, C)$ of (3.1) is such that $\nabla C \in L^8(0, T; L^2(\Omega)^d)$ and $\nabla \mathbf{u} \in L^4(0, T; L^2)$. Then the solution of (3.1) is unique.

**Proof** Let $(\mathbf{u}_1, p_1, C_1)$, $(\mathbf{u}_2, p_2, C_2)$ be two solutions of (3.1) and set

$$ \mathbf{u} = \mathbf{u}_1 - \mathbf{u}_2, \quad p = p_1 - p_2, \quad C = C_1 - C_2. $$

Then

$$ \left\{ \begin{array}{l} \displaystyle \frac{1}{2} \frac{d}{dt} \| \mathbf{u} \|^2 + a(\mathbf{u}, \mathbf{u}) + \underbrace{\int_S (|K \mathbf{u}_{1,\tau}|^{s-2} K \mathbf{u}_{1,\tau} - |K \mathbf{u}_{2,\tau}|^{s-2} K \mathbf{u}_{2,\tau}) \cdot K \mathbf{u}_{\tau} d\sigma}_{\geq 0} \\ \\ = -c_{\mathbf{u}}(\mathbf{u}, \mathbf{u}_2, \mathbf{u}) + (\mathbf{f}_1(C_1) - \mathbf{f}_1(C_2), \mathbf{u}), \\ \\ \displaystyle \frac{1}{2} \frac{d}{dt} \| C \|^2 + \alpha \| \nabla C \|^2 + r_0 \| C \|^2 = -c_C(\mathbf{u}, C_2, C). \end{array} \right. \quad (3.32) $$
From Cauchy-Schawrz's inequality, Hölder inequality, Lemma 2.4, and Assumption 2.5, (3.32) implies that

$$ \left\{ \frac{1}{2} \frac{d}{dt} \|\mathbf{u}\|^2 + cv \|\nabla \mathbf{u}\|^2 \le \|\mathbf{u}\|_{L^4}^2 \|\nabla \mathbf{u}_2\| + c \|C\| \|\mathbf{u}\|, \frac{1}{2} \frac{d}{dt} \|C\|^2 + \alpha \|\nabla C\|^2 + r_0 \|C\|^2 \le \|\mathbf{u}\|_{L^4} \|\nabla C_2\| \|C\|_{L^4} \right. \tag{3.33} $$

For $d=2$, we have (see (2.1))

$$ \|\mathbf{u}\|_{L^4}^2 \|\nabla \mathbf{u}_2\| + \|\mathbf{u}\|_{L^4} \|\nabla C_2\| \|C\|_{L^4} \le c \|\nabla \mathbf{u}\| \|\mathbf{u}\| \|\nabla \mathbf{u}_2\| + c \|\mathbf{u}\|^{1/2} \|\nabla \mathbf{u}\|^{1/2} \|\nabla C_2\| \|C\|^{1/2} \|\nabla C\|^{1/2}. $$

For $d=3$, we have (see (2.2))

$$ \|\mathbf{u}\|_{L^4}^2 \|\nabla \mathbf{u}_2\| + \|\mathbf{u}\|_{L^4} \|\nabla C_2\| \|C\|_{L^4} \le c \|\nabla \mathbf{u}\|^{3/2} \|\mathbf{u}\|^{1/2} \|\nabla \mathbf{u}_2\| + c \|\mathbf{u}\|^{1/4} \|\nabla \mathbf{u}\|^{3/4} \|\nabla C_2\| \|C\|^{1/4} \|\nabla C\|^{3/4}. $$

Thus

$$ \begin{aligned} & \|\mathbf{u}\|_{L^4}^2 \|\nabla \mathbf{u}_2\| + \|\mathbf{u}\|_{L^4} \|\nabla C_2\| \|C\|_{L^4} \\ & \le c \|\nabla \mathbf{u}\| \|\mathbf{u}\| \|\nabla \mathbf{u}_2\| + c \|\mathbf{u}\|^{1/2} \|\nabla \mathbf{u}\|^{1/2} \|\nabla C_2\| \|C\|^{1/2} \|\nabla C\|^{1/2} \\ & \qquad + c \|\nabla \mathbf{u}\|^{3/2} \|\mathbf{u}\|^{1/2} \|\nabla \mathbf{u}_2\| + c \|\mathbf{u}\|^{1/4} \|\nabla \mathbf{u}\|^{3/4} \|\nabla C_2\| \|C\|^{1/4} \|\nabla C\|^{3/4}. \end{aligned} \tag{3.34} $$

Adding the above inequalities in (3.33), using Lemma 2.2 and the Young's Inequality, we obtain

$$ \begin{aligned} & \frac{1}{2} \frac{d}{dt} [\|\mathbf{u}\|^2 + \|C\|^2] + \min(cv, \alpha) [\|\nabla \mathbf{u}\|^2 + \|\nabla C\|^2] + r_0 \|C\|^2 \\ & \le c (\|C\| \|\mathbf{u}\| + \|\nabla \mathbf{u}\| \|\mathbf{u}\| \|\nabla \mathbf{u}_2\| + \|\mathbf{u}\|^{1/2} \|\nabla \mathbf{u}\|^{1/2} \|\nabla C_2\| \|C\|^{1/2} \|\nabla C\|^{1/2} \\ & \qquad + \|\nabla \mathbf{u}\|^{3/2} \|\mathbf{u}\|^{1/2} \|\nabla \mathbf{u}_2\| + \|\mathbf{u}\|^{1/4} \|\nabla \mathbf{u}\|^{3/4} \|\nabla C_2\| \|C\|^{1/4} \|\nabla C\|^{3/4}) \\ & \le c_1 \left[ \frac{1}{2} \|C\|^2 + \frac{1}{2} \|\mathbf{u}\|^2 + \varepsilon_2 \|\nabla \mathbf{u}\|^2 + \frac{1}{\varepsilon_2} \|\mathbf{u}\|^2 \|\nabla \mathbf{u}_2\|^2 \\ & \qquad + \varepsilon_3 \|\nabla \mathbf{u}\|^2 + \frac{1}{\varepsilon_3} \left( \frac{1}{\varepsilon_4} \|\nabla C_2\|^4 \|C\|^2 + \varepsilon_4 \|\nabla C\|^2 \right) \\ & \qquad + \varepsilon_5 \|\nabla \mathbf{u}\|^2 + \frac{1}{\varepsilon_5} \|\mathbf{u}\|^2 \|\nabla \mathbf{u}_2\|^4 + \varepsilon_6 \|\nabla \mathbf{u}\|^2 + \frac{1}{\varepsilon_6} \left( \frac{1}{\varepsilon_7} \|\nabla C_2\|^8 \|C\|^2 + \varepsilon_7 \|\nabla C\|^2 \right) \right]. \end{aligned} $$

Denoting $\hat{c}_1 = \min(cv, \alpha)$ and Choosing $\varepsilon_2 = \varepsilon_3 = \varepsilon_5 = \varepsilon_6 = \frac{\hat{c}_1}{16}$ and $\varepsilon_4 = \varepsilon_7 = \frac{\hat{c}_1^2}{64}$, we obtain

$$ \frac{d}{dt} [\|\mathbf{u}\|^2 + \|C\|^2] \le c_2 (1 + \|\nabla \mathbf{u}_2\|^2 + \|\nabla C_2\|^4 + \|\nabla \mathbf{u}_2\|^4 + \|\nabla C_2\|^8) (\|\mathbf{u}\|^2 + \|C\|^2). $$

By integrating between 0 and $t$, we obtain

$$ \|\mathbf{u}\|^2 + \|C\|^2 \le \exp \left( c_3 \int_0^t (1 + \|\nabla \mathbf{u}_2\|^2 + \|\nabla C_2\|^4 + \|\nabla \mathbf{u}_2\|^4 + \|\nabla C_2\|^8) \right) (\|\mathbf{u}(0)\|^2 + \|C(0)\|^2) = 0. $$

This ends the proof.

□
## 4 Space-time discretization and a priori errors

In this section, we propose a time and space discretization of the problem (3.1) and provide an a priori error analysis. We use the semi-implicit Euler method for the time discretization and the finite element method for the space discretization. For the time discretization, we introduce a partition of the interval $[0, T]$ into $N$ subintervals $[t_{n-1}, t_n]$ of length $k$, where $k$ represents the time step. For the space discretization, we suppose that the domain $\Omega$ is a polygon (respectively polyhedron) for $d = 2$ (respectively $d = 3$). We choose a discretization parameter $h > 0$ in space and for each $h$, we define $\mathcal{T}_h$ as a regular (or non-degenerate) family of triangles (respectively tetrahedra) for $d = 2$ (respectively for $d = 3$), satisfying the following properties:

* $\overline{\Omega}$ is the union of all elements of $\mathcal{T}_h$;

* the intersection of two different elements of $\mathcal{T}_h$, if not empty, it is either a vertex or a whole edge (or a whole face of both for $d = 3$);

* the ratio of the diameter of an element $\kappa$ in $\mathcal{T}_h$ to the diameter of its inscribed sphere is bounded by a constant independent of $h$.

Let $X_h \subset X$ and $M_h \subset M$ be a "stable" pair of finite element spaces for discretizing the velocity $\mathbf{u}$ and the pressure $p$. The pair is considered stable in the sense that it satisfies a uniform discrete inf-sup condition. Specifically, there exists a constant $\beta^* \ge 0$, independent of $h$, such that,

$$ \forall q_h \in M_h, \sup_{\mathbf{v}_h \in X_h} \frac{1}{|\mathbf{v}_h|^2 H^1(\Omega)^d} \int_{\Omega} q_h \mathrm{div} \, \mathbf{v}_h dx \ge \beta^* \| q_h \|_{L^2(\Omega)}, \quad (4.1) $$

where $h$ denotes the maximal diameter of the elements of $\mathcal{T}_h$. Moreover, we consider $Y_h \subset Y$ the space of discretization of the concentration $C$.

Let $\mathbf{P}_\ell$ denote the space of polynomials with total degree less than or equal to $\ell$. We choose the "mini-element" (see D. Arnold, F. Brezzi and M. Fortin in [5]) discretization, where in each element $\kappa$, the pressure $p$ is a polynomial of $\mathbf{P}_1$, and each component of the velocity $\mathbf{u}$ is the sum of a polynomial of $\mathbf{P}_1$ and a "bubble" function $b_\kappa$. The bubble function, for each element $\kappa$, is determined as the product of the barycentric coordinates associated with the vertices of $\kappa$).

Therefore, the finite-element spaces for the velocity and the pressure are defined as follows:

$$ X_h = \{ \mathbf{v}_h \in C^0(\overline{\Omega})^d \cap X; \forall \kappa \in \mathcal{T}_h, \mathbf{v}_{h|_\kappa} \in \mathcal{P}(\kappa) \}, $$

$$ Y_h = \{ q_h \in C^0(\overline{\Omega}) \cap Y; \forall \kappa \in \mathcal{T}_h, \mathbf{v}_{h|_\kappa} \in \mathbf{P}_1(\kappa) \}, $$

and

$$ M_h = \{ q_h \in C^0(\overline{\Omega}) \cap M; \forall \kappa \in \mathcal{T}_h, \mathbf{v}_{h|_\kappa} \in \mathbf{P}_1(\kappa) \}, $$

where

$$ \mathcal{P}(\kappa) = [\mathbb{P}_1 \oplus \mathrm{Span}(b_\kappa)]^d. $$

We introduce the following discrete space:

$$ V_h = \{ \mathbf{v}_h \in X_h; \forall q_h \in M_h, \int_{\Omega} q_h(x) \mathrm{div} \, \mathbf{v}_h(x) dx = 0 \}. $$
There exists an approximation operator $P_h \in \mathcal{L}(H_0^1(\Omega)^d; X_h)$ such that (see V. Girault and P.-A. Raviart in [15]):

$$ \forall v \in H_0^1(\Omega)^d, \quad \forall q_h \in M_h, \quad \int_{\Omega} q_h \operatorname{div}(P_h(v) - v)d\mathbf{x} = 0, $$

and for $\ell = 0$ or 1,

$$ \forall v \in [H^{1+\ell}(\Omega) \cap H_0^1(\Omega)]^d, \quad \| P_h(v) - v \|_{L^2(\Omega)^d} \le C_1 h^{1+\ell} |v|_{H^{1+\ell}(\Omega)^d}, $$

and for all $r \ge 2$, $\ell = 0$ or 1,

$$ \forall v \in [W^{1+\ell,r}(\Omega) \cap H_0^1(\Omega)]^d, \quad |P_h(v) - v|_{W^{1,r}(\Omega)^d} \le C_2 h^\ell |v|_{W^{1+\ell,r}(\Omega)^d}. $$

Furthermore, there exists an approximation operator (when $d = 2$, see Bernardi and Girault [6] or Clément [10]; when $d = 2$ or $d = 3$, see Scott and Zhang [28]), $R_h$ in $\mathcal{L}(W^{1,p}(\Omega); Z_h)$ and in $\mathcal{L}(W^{1,p}(\Omega) \cap H_0^1(\Omega); X_h)$ such that for all $m = 0, 1, l = 0, 1$, and all $p \ge 2$,

$$ \forall S \in W^{l+1,p}(\Omega), \quad |S - R_h(S)|_{W^{m,p}(\Omega)} \le C(p, m, l) h^{l+1-m} |S|_{W^{l+1,p}(\Omega)}. $$

In order to introduce the discrete scheme, we define the following forms:

1. For all $\mathbf{u}_h, \mathbf{v}_h, \mathbf{w}_h \in X_h$,

$$ d_{\mathbf{u}}(\mathbf{u}_h, \mathbf{v}_h, \mathbf{w}_h) = c_{\mathbf{u}}(\mathbf{u}_h, \mathbf{v}_h, \mathbf{w}_h) + \frac{1}{2}(\operatorname{div}(\mathbf{u}_h) \mathbf{v}_h, \mathbf{w}_h). $$

2. For all $\mathbf{v}_h \in X_h, \xi_h, \eta_h \in Y_h$,

$$ d_C(\mathbf{v}_h, \xi_h, \eta_h) = c_C(\mathbf{v}_h, \xi_h, \eta_h) + \frac{1}{2}(\operatorname{div}(\mathbf{v}_h) \xi_h, \eta_h). $$

**Remark 4.1** For all $\mathbf{u}_h, \mathbf{v}_h$ in $X_h$, the form $d_{\mathbf{u}}$ verifies the following stability relation

$$ d_{\mathbf{u}}(\mathbf{u}_h, \mathbf{v}_h, \mathbf{v}_h) = 0. $$

Furthermore, for all $\mathbf{v}_h \in X_h, \xi_h \in Y_h$, the form $d_C$ verifies the following stability relation

$$ d_C(\mathbf{v}_h, \xi_h, \xi_h) = 0. $$

## 4.1 Non linear discrete scheme

In this section, we introduce the following space time discretization of (3.1) and analyse its well-posedness and the convergence by means of estimating the expressions $\mathbf{u}(t_n) - \mathbf{u}_h^n, p(t_n) - p_h^n$ and $C(t_n) - C_h^n$.

For each $n \in \{1, ..., N\}$, knowing $\mathbf{u}_h^{n-1} \in X_h$ and $C_h^{n-1} \in Y_h$, we compute $(\mathbf{u}_h^n, p_h^n, C_h^n) \in X_h \times M_h \times Y_h$ such that for all $(\mathbf{v}_h, q_h, r_h) \in X_h \times M_h \times Y_h$

$$ \left\{ \begin{array}{l} \displaystyle \frac{1}{k}(\mathbf{u}_h^n - \mathbf{u}_h^{n-1}, \mathbf{v}_h) + a(\mathbf{u}_h^n, \mathbf{v}_h) + d_{\mathbf{u}}(\mathbf{u}_h^n, \mathbf{u}_h^n, \mathbf{v}_h) - (p_h^n, \operatorname{div} \mathbf{v}_h) \\ \displaystyle + \int_S |K \mathbf{u}_h^n|^{s-2} K \mathbf{u}_h^n \cdot K \mathbf{v}_h d\sigma = (f^n(C_h^n), \mathbf{v}_h), \quad (4.2) \\ \displaystyle \frac{1}{k}(C_h^n - C_h^{n-1}, r_h) + \alpha(\nabla C_h^n, \nabla r_h) + d_C(\mathbf{u}_h^n, C_h^n, r_h) + r_0(C_h^n, r_h) = (g^n, r_h), \\ (q_h, \operatorname{div} \mathbf{u}_h^n) = 0, \end{array} \right. $$

where $\mathbf{u}_h^0$ and $C_h^0$ are approximations of $\mathbf{u}_0$ and $C_0$, respectively. We construct next, the solution of discrete problem (4.2).
**Remark 4.2** We note that in the numerical non-linear scheme (4.2), we use the forms $d_u$ and $d_C$ instead of $c_u$ and $c_C$ for stability reasons (in this discrete case), as explained in Remark 4.1, which ensures the existence of a corresponding solution. In contrast, the analogous problem (3.4) in paragraph 3.2 uses the forms $c_u$ and $c_C$ since they satisfy the stability properties, given the continuous nature of the problem.

**Theorem 4.3** (Existence of discrete solutions) At each time step $n$ and for a given $\mathbf{u}_h^{n-1} \in X_h$ and $C_h^{n-1} \in Y_h$, problem (4.2) admits at least one solution $(\mathbf{u}_h^n, p_h^n, C_h^n) \in X_h \times M_h \times Y_h$.

*Proof* Given $\mathbf{u}_h^{n-1} \in X_h$ and $C_h^{n-1} \in Y_h$. We observe that (4.2) has the same structure as the problem discussed in paragraph 3.2, (see the problem (3.4)). This the conclusion of Theorem 3.3 is applicable. □

We address next, the stability of the solution of (4.2). We claim that

**Theorem 4.4** (Stability of the discrete problem) Each solution of Problem (4.2) verifies, for $m = 1, ..., N$, the following bounds

$$ ||C_h^m||^2 + \sum_{n=1}^{m} k |C_h^n|_{1,\Omega}^2 \le c_1 \left( \sum_{n=1}^{N} k ||g^n||^2 + ||C_h^0||^2 \right) $$

and

$$ ||\mathbf{u}_h^m||^2 + \sum_{n=1}^{m} k |\mathbf{u}_h^n|_{1,\Omega}^2 + \sum_{n=1}^{m} k \int_S |K \mathbf{u}_h^n \tau|^s d\sigma \\ \le c_2 \left( \sum_{n=1}^{N} k ||\mathbf{f}_0^n||^2 + ||\mathbf{u}_h^0||^2 + \sum_{n=1}^{N} k ||g^n||^2 + ||C_h^0||^2 \right), $$

where $c_1$ and $c_2$ are positive constants independent of $h$ and $k$.

*Proof* The proof is identical to that of Proposition 3.4. □

**Theorem 4.5** The solution of system (4.2) is unique if the time step $k$ is sufficiently small.

*Proof* Let $(\mathbf{u}_{h1}^n, p_{h1}^n, C_{h1}^n)$ and $(\mathbf{u}_{h2}^n, p_{h2}^n, C_{h2}^n)$ be two solutions of Problem (4.2). We denote by $\mathbf{w}_h^n = \mathbf{u}_{h1}^n - \mathbf{u}_{h2}^n, \xi_h^n = p_{h1}^n - p_{h2}^n$ and $\theta_h^n = C_{h1}^n - C_{h2}^n$. By taking $\mathbf{v}_h = \mathbf{w}_h^n, q_h = \xi_h^n$ and $r_h = \theta_h^n$ in System (4.2) we obtain

$$ \begin{aligned} & \frac{1}{k}(\mathbf{w}_h^n, \mathbf{w}_h^n) + a(\mathbf{w}_h^n, \mathbf{w}_h^n) + \int_S (|K\mathbf{u}_{h1\tau}^n|^{s-2} K\mathbf{u}_{h1\tau}^n - |K\mathbf{u}_{h2\tau}^n|^{s-2} K\mathbf{u}_{h2\tau}^n) \\ & \cdot (K\mathbf{u}_{h1\tau}^n - K\mathbf{u}_{h2\tau}^n)d\sigma \\ & = -d_u(\mathbf{w}_h^n, \mathbf{u}_{h2}^n, \mathbf{w}_h^n) + (\mathbf{f}_1^n(C_{h1}^n) - \mathbf{f}_1^n(C_{h2}^n), \mathbf{u}_h) \end{aligned} \quad (4.3) $$

and

$$ \frac{1}{k}(\theta_h^n, \theta_h^n) + \alpha(\nabla\theta_h^n, \nabla\theta_h^n) + r_0(\theta_h^n, \theta_h^n) = -d_C(\mathbf{u}_h, C_{h2}^n, \theta_h^n). \quad (4.4) $$

Remark that the third term of Eq. (4.3) is positive (see Lemma 3.2), and using relation (2.3) and Lemma 2.4, we get

$$ \begin{aligned} & \| \mathbf{w}_h^n \|^2 + 2cvk \| \nabla \mathbf{w}_h^n \|^2 \le k \| \mathbf{w}_h^n \|_{L^4(\Omega)^d}^2 \| \nabla \mathbf{u}_{h2}^n \| + c_{f_1} k \| \theta_h^n \| \| \mathbf{w}_h^n \| \\ & \le k \| \mathbf{w}_h^n \|_{L^4(\Omega)^d}^2 \| \nabla \mathbf{u}_{h2}^n \| + \frac{c_{f_1}}{2} k \| \theta_h^n \|^2 + \frac{c_{f_1}}{2} k \| \mathbf{w}_h^n \|^2 \end{aligned} $$
and

$$ \| \theta_h^n \|^2 + \alpha k \| \nabla \theta_h^n \|^2 \leq k \| \mathbf{w}_h^n \|_{L^4(\Omega)^d} \| \nabla C_{h2}^n \| \| \theta_h^n \|_{L^4(\Omega)} $$

We use the fact that in finite dimension spaces all the norms are equivalent (for a fixed mesh step h), sum the above inequalities from n = 0 to m, remark that $\mathbf{w}_h^0 = \mathbf{0}$ and $\theta_h^0 = 0$, and use relations (3.9) to get after summing the obtained inequalities the following bound:

$$ \| \mathbf{w}_h^m \|^2 + \| \theta_h^m \|^2 + \sum_{n=0}^{m} k (\| \nabla \mathbf{w}_h^n \|^2 + \| \nabla \theta_h^n \|^2) \leq C \sum_{n=0}^{m} k (\| \mathbf{w}_h^n \|^2 + \| \theta_h^n \|^2), $$

where C is a constant depending on h but does not depend on k.

For a time step k sufficiently small ($k \le \frac{1}{C}$), the last bound gives

$$ \| \mathbf{w}_h^m \|^2 + \| \theta_h^m \|^2 + \sum_{n=0}^{m} k (\| \nabla \mathbf{w}_h^n \|^2 + \| \nabla \theta_h^n \|^2) \leq C_1 \sum_{n=0}^{m-1} k (\| \mathbf{w}_h^n \|^2 + \| \theta_h^n \|^2). $$

Applying Gronwall's lemma, we obtain $\mathbf{w}_h^n = \mathbf{0}$ and $\theta_h^n = 0$, and then $\mathbf{u}_{h1}^n = \mathbf{u}_{h2}^n$ and $C_{h1}^n = C_{h2}^n$. The discrete inf-sup condition (4.1) deduces that $p_{h1}^n = p_{h2}^n$. Hence, we get the desired result. □

Next, we will establish an a priori error estimate between the exact and the numerical solutions.

**Theorem 4.6** Let $(\mathbf{u}, p)$ be the solution of Problem (3.1) and $(\mathbf{u}_h^n, p_h^n)$ be the solution of Problem (4.2). If $\mathbf{u} \in L^\infty(0, T; H^2(\Omega)^d)$, $\mathbf{u}_\tau \in L^\infty(0, T; H^2(S)^d)$, $\mathbf{u}' \in L^2(0, T; H^1(\Omega)^d)$, $\mathbf{u}'' \in L^\infty(0, T; L^2(\Omega)^d)$, $p \in L^\infty(0, T; H^1(\Omega))$, $C \in L^\infty(0, T; H^2(\Omega))$, $C' \in L^2(0, T; H^1(\Omega))$ and $C'' \in L^\infty(0, T; L^2(\Omega))$, there exists positive constants $C_1$ and $C_2$ independent of h and k such that for all $m \in \{1, \dots, N\}$, the following error estimates hold:

$$ \begin{aligned} & \| \mathbf{u}_h^m - P_h \mathbf{u}(t_m) \|_{L^2(\Omega)^d}^2 + \sum_{n=1}^{m} k |\mathbf{u}_h^n - P_h \mathbf{u}(t_n)|_{H^1(\Omega)^d}^2 && (4.5) \\ & \leq C_1(h^2 + k^2 + ||C_h^0 - R_h C(t_0)||^2 + ||\mathbf{u}_h^0 - P_h \mathbf{u}(t_0)||^2), \end{aligned} $$

and

$$ \begin{aligned} & \| C_h^m - R_h C(t_m) \|_{L^2(\Omega)}^2 + \sum_{n=1}^{m} k |C_h^n - R_h C(t_n)|_{H^1(\Omega)}^2 && (4.6) \\ & \leq C_2(h^2 + k^2 + ||C_h^0 - R_h C(t_0)||^2 + ||\mathbf{u}_h^0 - P_h \mathbf{u}(t_0)||^2). \end{aligned} $$

**Proof** We choose the test functions $\mathbf{v} = \mathbf{v}_h = \mathbf{u}_h^n - P_h \mathbf{u}(t_n)$ (denoted by $\mathbf{v}_h^n$) in the first equations of (3.1) and (4.2) multiplied by k. By setting $t = t_n$ in the first equation of (3.1) and replacing

$$ k \mathbf{u}'(t_n) = \mathbf{u}(t_{n+1}) - \mathbf{u}(t_n) - k^2 \mathbf{u}''(\xi), \quad \text{for } \xi \in [t_n, t_{n+1}]. $$
Then we subtract the obtained continuous and discrete equations, and insert $\pm P_h u(t_{n-1})$ and $\pm P_h u(t_n)$ to obtain:

$$ 
\begin{aligned}
& \frac{1}{2} \| \mathbf{v}_h^n \|^2 - \frac{1}{2} \| \mathbf{v}_h^n \|^2 + \frac{1}{2} \| \mathbf{v}_h^n - \mathbf{v}_h^{n-1} \|^2 + 2\nu k (D(\mathbf{v}_h^n), D(\mathbf{v}_h^n)) \\
& + k \int_S \left( |K \mathbf{u}_{h\tau}^n|^{s-2} k \mathbf{u}_{h\tau}^n - |K \mathbf{u}_\tau(t_n)|^{s-2} K \mathbf{u}_\tau(t_n) \right) K \mathbf{v}_h^n d\sigma \\
& = k(\mathbf{f}_1^n(C_h^n) - \mathbf{f}_1(C(t_n)), \mathbf{v}_h^n) + k(p_h^n - p(t_n), \text{div}(\mathbf{v}_h^n)) - k^2 \left( \frac{d^2\mathbf{u}}{dt^2}(\xi), \mathbf{v}_h^n \right) \\
& - 2\nu k (D(P_h \mathbf{u}(t_n) - \mathbf{u}(t_n)), D\mathbf{v}_h^n) - \int_{t_{n-1}}^{t_n} (P_h \mathbf{u}'(t) - \mathbf{u}'(t), \mathbf{v}_h^n) dt \\
& - k(\mathbf{u}_h^n \nabla) \cdot \mathbf{u}_h^n - (\mathbf{u}(t_n)\nabla) \cdot \mathbf{u}(t_n), \mathbf{v}_h^n - \frac{1}{2} k (\text{div}(\mathbf{u}_h^n)\mathbf{u}_h^n, \mathbf{v}_h^n).
\end{aligned} 
\quad (4.7) $$

The last term on the left-hand side of (4.7) can be written as follows:

$$ 
\begin{aligned}
& k \int_S \left( |K \mathbf{u}_{h\tau}^n|^{s-2} k \mathbf{u}_{h\tau}^n - |K \mathbf{u}_\tau(t_n)|^{s-2} K \mathbf{u}_\tau(t_n) \right) (K \mathbf{u}_{h\tau}^n - K P_h \mathbf{u}_\tau(t_n)) d\sigma \\
& = k \int_S \left( |K \mathbf{u}_{h\tau}^n|^{s-2} K \mathbf{u}_{h\tau}^n - |K \mathbf{u}_\tau(t_n)|^{s-2} K \mathbf{u}_\tau(t_n) \right) (K \mathbf{u}_{h\tau}^n - K \mathbf{u}_\tau(t_n)) d\sigma \\
& + k \int_S \left( |K \mathbf{u}_{h\tau}^n|^{s-2} K \mathbf{u}_{h\tau}^n - |K \mathbf{u}_\tau(t_n)|^{s-2} K \mathbf{u}_\tau(t_n) \right) (K \mathbf{u}_\tau(t_n) - K P_h \mathbf{u}_\tau(t_n)) d\sigma.
\end{aligned} 
\quad (4.8) $$

We replace (4.8) in (4.7) to obtain:

$$ 
\begin{aligned}
& \frac{1}{2} \| \mathbf{v}_h^n \|^2 - \frac{1}{2} \| \mathbf{v}_h^{n-1} \|^2 + \frac{1}{2} \| \mathbf{v}_h^n - \mathbf{v}_h^{n-1} \|^2 + 2\nu k (D(\mathbf{v}_h^n), D(\mathbf{v}_h^n)) \\
& + k \int_S \left( |K \mathbf{u}_{h\tau}^n|^{s-2} K \mathbf{u}_{h\tau}^n - |K \mathbf{u}_\tau(t_n)|^{s-2} K \mathbf{u}_\tau(t_n) \right) (K \mathbf{u}_{h\tau}^n - K \mathbf{u}_\tau(t_n)) d\sigma \\
& = k(\mathbf{f}_1^n(C_h^n) - \mathbf{f}_1(C(t_n)), \mathbf{v}_h^n) + k(p_h^n - p(t_n), \text{div}(\mathbf{v}_h^n)) - k^2 \left( \frac{d^2\mathbf{u}}{dt^2}(\xi), \mathbf{v}_h^n \right) \\
& - 2\nu k (D(P_h \mathbf{u}(t_n) - \mathbf{u}(t_n)), D\mathbf{v}_h^n) - \int_{t_{n-1}}^{t_n} (P_h \mathbf{u}'(t) - \mathbf{u}'(t), \mathbf{v}_h^n) dt \\
& - k(\mathbf{u}_h^n \nabla) \cdot \mathbf{u}_h^n - (\mathbf{u}(t_n)\nabla) \cdot \mathbf{u}(t_n), \mathbf{v}_h^n - \frac{1}{2} k (\text{div}(\mathbf{u}_h^n)\mathbf{u}_h^n, \mathbf{v}_h^n) \\
& - k \int_S \left( |K \mathbf{u}_{h\tau}^n|^{s-2} K \mathbf{u}_{h\tau}^n - |K \mathbf{u}_\tau(t_n)|^{s-2} K \mathbf{u}_\tau(t_n) \right) (K \mathbf{u}_\tau(t_n) - K P_h \mathbf{u}_\tau(t_n)) d\sigma.
\end{aligned} 
\quad (4.9) $$

From (3.2), we deduce that the third term on the left-hand side satisfies the following inequality:

$$ k \int_S \left( |K \mathbf{u}_{h\tau}^n|^{s-2} K \mathbf{u}_{h\tau}^n - |K \mathbf{u}_\tau(t_n)|^{s-2} K \mathbf{u}_\tau(t_n) \right) (K \mathbf{u}_{h\tau}^n - K \mathbf{u}_\tau(t_n)) d\sigma \geq 0. $$

All the terms of the right-hand side of (4.9) can be treated as in the proof of Theorem 4.5 in [3] except for the first one which can be bounded as follows:

$$ 
\begin{aligned}
& k |(\mathbf{f}_1^n(C_h^n) - \mathbf{f}_1(C(t_n)), \mathbf{v}_h^n)| \\
& \leq c_1^* k ||C_h^n - C(t_n)|| ||\mathbf{v}_h^n|| \\
& \leq \frac{c_1^2}{2\epsilon} k ||C_h^n - C(t_n)||^2 + \frac{\epsilon}{2} k |\mathbf{v}_h^n|^{H^1(\Omega)^d}.
\end{aligned} 
$$
Therefore, by following the same steps in the proof of Theorem 4.5 of [3] we obtain the following inequality: for $m = 1, \dots, N$,

$$ \begin{aligned} \| \mathbf{v}_h^m \|^2 + \sum_{n=1}^{m} k |\mathbf{v}_h^n|_{1,\Omega}^2 &\le C_1(k^2 + h^2 + \|\mathbf{u}_h^0 - P_h \mathbf{u}(t_0)\|^2) \\ &+ C_2 \sum_{n=1}^{m} k \|C_h^n - C(t_n)\|^2 + C_3 \sum_{n=1}^{m} k \|\mathbf{v}_h^n\|^2. \end{aligned} $$

By applying Lemma 2.6 (a discrete Gronwall's Lemma), we have

$$ \| \mathbf{v}_h^m \|^2 + \sum_{n=1}^{m} k |\mathbf{v}_h^n|_{1,\Omega}^2 \le c_1(k^2 + h^2 + \|\mathbf{u}_h^0 - P_h \mathbf{u}(t_0)\|^2) + C_2 \sum_{n=1}^{m} k \|C_h^n - C(t_n)\|^2. $$

Let us now treat by the same way the concentration equations. We consider the second equation of (3.1) and (4.2) multiplied by $k$, and choose $r = r_h = C_h^n - R_h C(t_n)$ (denoted by $r_h^n$). By setting $t = t_n$ in the second equation of (3.1) and replacing

$$ kC'(t_n) = C(t_{n+1}) - C(t_n) - k^2 C''(\xi), \quad \text{for } \xi \in [t_n, t_{n+1}]. $$

We obtain

$$ \left\{ \begin{aligned} & \frac{1}{2} \| r_h^n \|^2 - \frac{1}{2} \| r_h^n \|^2 + \frac{1}{2} \| r_h^n - r_h^{n-1} \|^2 + \alpha k (\nabla r_h^n, \nabla r_h^n) + r_0(r_h^n, r_h^n) \\ & = -k^2 \left( \frac{d^2 C}{dt^2}(\xi), r_h^n \right) - \alpha k (\nabla(R_h C(t_n) - C(t_n)), \nabla r_h^n) - r_0(R_h C(t_n) - C(t_n), r_h^n) \\ & \quad - \int_{t_{n-1}}^{t_n} (R_h C'(t) - C'(t), r_h^n) dt \\ & \quad - k \left( (\mathbf{u}_h^n \cdot \nabla) C_h^n - (\mathbf{u}(t_n) \cdot \nabla) C(t_n), r_h^n \right) - \frac{1}{2} k (\text{div } (\mathbf{u}_h^n) C_h^n, r_h^n). \end{aligned} \right. \quad (4.11) $$

We can treat the first term of the right-hand side of Eq. (4.11) as follows

$$ k^2 \left| \left( \frac{d^2 C}{dt^2}(\xi), r_h^n \right) \right| \le \frac{c_2^2}{2\epsilon_1} k^3 \|C''\|^2_{L^\infty(0,T;L^2(\Omega))} + \frac{\epsilon_2}{2} k \|r_h^n\|^2_{H^1(\Omega)}. $$

The second and third terms can be treated as follows

$$ \begin{aligned} & \alpha k |(\nabla(R_h C(t_n) - C(t_n)), \nabla r_h^n) + r_0(R_h C(t_n) - C(t_n), r_h^n)| \\ & \le \frac{c_3^2}{2\epsilon_3} h^2 k \|C\|^2_{L^\infty(0,T;H^2(\Omega))} + \frac{\epsilon_3}{2} k \|r_h^n\|^2_{H^1(\Omega)}. \end{aligned} $$

It is easy to check that the fourth term satisfies the following inequality:

$$ \left| \int_{t_{n-1}}^{t_n} (R_h C'(t) - C'(t), r_h^n) dt \right| \le \frac{c_4^2}{2\epsilon_4} h^2 \|C'\|^2_{L^2(t_{n-1}, t_n; H^2(\Omega))} + \frac{\epsilon_4}{2} k \|r_h^n\|^2_{H^1(\Omega)}. $$

Let us now treat the last two terms of the right-hand side of Eq. (4.11) denoted by $L_n$. By using Remark 4.1, inserting $\pm R_h C(t_n)$, $\pm C(t_n)$, $\pm R_h C(t_{n-1})$, we apply the Green formula
and obtain

$$
\begin{align*}
L_n &= -k((\mathbf{u}_h^n \nabla) \cdot C_h^n - (\mathbf{u}(t_n)\nabla) \cdot C(t_n), r_h^n) - \frac{k}{2} (\operatorname{div}(\mathbf{u}_h^n) C_h^n, r_h^n) \\
&= -k\left[(\mathbf{u}_h^n \nabla) \cdot (R_h C(t_n) - C(t_n)), r_h^n\right] - \frac{k}{2} (\operatorname{div}(\mathbf{u}_h^n)(R_h C(t_n) - C(t_n)), r_h^n) \right] \\
&\quad -k\left[(\mathbf{u}_h^n - \mathbf{u}(t_n))\nabla) \cdot C(t_n), r_h^n\right] - \frac{k}{2} (\operatorname{div}(\mathbf{u}_h^n - \mathbf{u}(t_n))C(t_n), r_h^n) \right] \\
&= -k\left[(\mathbf{u}_h^n \nabla) \cdot (R_h C(t_n) - C(t_n)), r_h^n\right] - \frac{k}{2} (\operatorname{div}(\mathbf{u}_h^n)(R_h C(t_n) - C(t_n)), r_h^n) \right] \\
&\quad -\frac{k}{2}\left[(\mathbf{u}_h^n - \mathbf{u}(t_n))\nabla) \cdot C(t_n), r_h^n\right] + \frac{k}{2}\left[(\mathbf{u}_h^n - \mathbf{u}(t_n))\nabla) \cdot r_h^n, C(t_n)\right].
\end{align*}
$$

Thus, as $C \in L^\infty(0, T; H^2(\Omega))$, $L_n$ can be bounded as follows

$$
|L_n| \le c_5 h^2 \left( \frac{1}{2\epsilon_5} k ||\mathbf{u}_h^n||_{1,\Omega}^2 + \frac{\epsilon_5}{2} k ||r_h^n||_{1,\Omega}^2 \right) \\
+ c_6 \left( \frac{1}{2\epsilon_6} k ||\mathbf{u}_h^n - \mathbf{u}(t_n)||_{L^2(\Omega)^d}^2 + \frac{\epsilon_6}{2} k ||r_h^n||_{1,\Omega}^2 \right).
$$

Using Eq. (4.11) along with the above inequalities and Theorem 4.4 we obtain

$$
\| r_h^m \|^2 + \sum_{n=1}^{m} k |r_h^n|_{1,\Omega}^2 \le C_4 (k^2 + h^2 + \|C_h^0 - R_h C(t_0)\|^2) + C_5 \sum_{n=1}^{m} k \|\mathbf{u}_h^n - \mathbf{u}(t_n)\|^2. \quad (4.12)
$$

Next, we insert $P_h u(t_n)$ in the last term of (4.12), make use of triangle inequality, (4.10) to obtain:

$$
\| r_h^m \|^2 + \sum_{n=1}^{m} k |r_h^n|_{1,\Omega}^2 \le C_6 (k^2 + h^2 + \|C_h^0 - R_h C(t_0)\|^2 \\
\qquad \qquad \qquad + \|\mathbf{u}_h^0 - P_h \mathbf{u}(t_0)\|^2) + C_7 \sum_{n=1}^{m} k \|C_h^n - C(t_n)\|^2.
$$

We then apply the discrete Gronwall's Lemma (Lemma 2.6) and get the asserted estimate
((4.6)). The estimate (4.5) is a simple deduction of the last inequality and (4.6). $\square$

Next, we observe after application of the triangle's inequality that the following a priori error
estimation:

**Corollary 4.7** Under the assumption of Theorem 4.6, the exact solution $(\mathbf{u}, p)$ of Problem (3.1) and the discrete solution $(\mathbf{u}_h^n, p_h^n)$ of Problem (4.2) satisfies the following a priori error estimates:

$$
\begin{align*}
\| \mathbf{u}_h^m - \mathbf{u}(t_m) \|_{L^2(\Omega)^2}^2 &+ \sum_{n=1}^m k |\mathbf{u}_h^n - \mathbf{u}(t_n)|_{H^1(\Omega)^2}^2 + \| C_h^m - C(t_m) \|_{L^2(\Omega)^2}^2 \\
&\quad + \sum_{n=1}^m k |C_h^n - C(t_n)|_{H^1(\Omega)^2}^2 \\
&\le c(h^2 + k^2),
\end{align*}
$$

where c is a positive constant independent of h and k.
## 4.2 Iterative scheme

To approximate the solution of the non-linear discrete system (4.2), we introduce an iterative problem consisting of two steps:

For every $n \in \{1, ..., N\}$, knowing $\mathbf{u}_h^{n-1} \in X_h$ and $C_h^{n-1} \in Y_h$, we compute $(\mathbf{u}_h^n, p_h^n, C_h^n) \in (X_h, M_h, Y_h)$ by solving the following iterative problem (over $i = 1, 2, ...$):

For each $\mathbf{w}_h^{i-1} \in X_h$, compute $(\mathbf{w}_h^i, \xi_h^i, \eta_h^i) \in X_h \times M_h \times Y_h$ such that for all $(\mathbf{v}_h, q_h, r_h) \in (X_h, M_h, Y_h)$

1. Step 1: we compute $(\mathbf{w}_h^i, \xi_h^i) \in X_h \times M_h$ such that,

$$ egin{cases} 
\frac{1}{k}(\mathbf{w}_h^i - \mathbf{u}_h^{n-1}, \mathbf{v}_h) + a(\mathbf{w}_h^i, \mathbf{v}_h) + d_{\mathbf{u}}(\mathbf{w}_h^{i-1}, \mathbf{w}_h^i, \mathbf{v}_h) - (\xi_h^i, \operatorname{div} \mathbf{v}_h) \\ 
+ \int_S |K \mathbf{w}_h^{i-1}|^{s-2} K \mathbf{w}_h^i \cdot K \mathbf{v}_h \tau \, d\sigma = (\mathbf{f}^n(\eta_h^{i-1}), \mathbf{v}_h), \\ 
(q_h, \operatorname{div} \mathbf{w}_h^i) = 0, 
\end{cases} 
\quad (4.13) $$

where $\mathbf{w}_h^0 = \mathbf{u}_h^{n-1}$ and $\mathbf{f}^n(C_h^{i-1}) = \mathbf{f}_0^n + \mathbf{f}_1^n(C_h^{i-1})$.

2. Step 2: Using $\mathbf{w}_h^i$ obtained from the first step, we compute $\eta_h^i$ by solving the following equation:

$$ \frac{1}{k}(\eta_h^i - C_h^{n-1}, r_h) + \alpha(\nabla\eta_h^i, \nabla r_h) + d_C(\mathbf{w}_h^i, \eta_h^i, r_h) + r_0(\eta_h^i, r_h) = (\mathbf{g}^n, r_h). \quad (4.14) $$

We repeat the iterative process until convergence is reached. Hence, yielding the numerical solution $(\mathbf{u}_h^n, p_h^n, C_h^n) \in X_h \times M_h \times Y_h$ which represents the limit of $(\mathbf{w}_h^i, \xi_h^i, \eta_h^i)$ as the iterative index $i$ tends to $+\infty$.

**Remark 4.8** At each step $n$ and given $(\mathbf{u}_h^{n-1}, C_h^{n-1})$, Problem (4.13) and (4.14) is an iterative one starting from the initial guess $\mathbf{w}_h^0 = \mathbf{u}_h^{n-1}$ and at each iteration $i$, knowing $\mathbf{w}_h^{i-1}$, we compute $(\mathbf{w}_h^i, \xi_h^i, \eta_h^i)$. As Problem (4.13) and (4.14) is a linearized coupling Navier-Stokes system with the concentration equation, the existence and uniqueness of the solution $(\mathbf{w}_h^i, \xi_h^i, \eta_h^i)$ is straightforward. In fact it suffices to show the uniqueness as we are with finite dimensional problem.

**Theorem 4.9** For $i = 1, 2, ..., $, the solution $(\mathbf{w}_h^i, p_h^i, C_h^i)$ of Problem (4.13)–(4.14) satisfies the following bounds:

$$ \|\eta_h^i\|^2 + \|\eta_h^i\|_{1,\Omega}^2 \le C_1 \|C_h^{n-1}\|^2, \quad (4.15) $$

$$ \|\mathbf{w}_h^i\|^2 + \|\mathbf{w}_h^i\|_{1,\Omega}^2 \le C_2 (\|\mathbf{u}_h^{n-1}\|^2 + \|C_h^{n-1}\|^2), \quad (4.16) $$

and

$$ \|p_h^i\| \le C_3 (\|\mathbf{w}_h^i\| + \|\mathbf{u}_h^{n-1}\| + \|\mathbf{w}_h^i\|_1 + \|\mathbf{w}_h^{i-1}\|_1 \|\mathbf{w}_h^i\|_1) \quad (4.17) $$

where $C_1, C_2$ and $C_3$ are positive constants independent of $i$ but depend on the time step $k$.

**Proof** By setting $\mathbf{v}_h = \mathbf{w}_h^i$ in (4.13) and $r_h = \eta_h^i$ in (4.14), and by using Remark 4.1 and the positivity of the term on $\Gamma$, we obtain the relations

$$ \left(\frac{1}{k} + r_0\right) \|\eta_h^i\|^2 + \|\eta_h^i\|_{1,\Omega}^2 \le \frac{c_1}{k} \|C_h^{n-1}\|^2, $$
and

$$ \frac{1}{k} \|\mathbf{w}_h^i\|^2 + \|\mathbf{w}_h^i\|_{1,\Omega}^2 \le c_2 \left( \frac{1}{k} \|\mathbf{u}_h^{n-1}\|^2 + \|C_h^{n-1}\|^2 \right), $$

where $c_1$ and $c_2$ are positive constants independent of $i$. Thus, we get the relations (4.15) and (4.16).

For the pressure, we use the following discrete inf-sup condition (a particular form of (4.1)):

$$ \forall q_h \in M_h, \quad \|q_h\|_{L^2(\Omega)} \le \frac{1}{\beta^*} \sup_{\mathbf{v}_h \in X_h \cap \{\mathbf{v}_h|\Gamma=0\}} \frac{1}{|\mathbf{v}_h|_{H^1(\Omega)^d}} \int_{\Omega} q_h \operatorname{div} \mathbf{v}_h dx, $$

and we get immediately (4.17). □

**Theorem 4.10** For a given $(\mathbf{u}_h^{n-1}, C_h^{n-1}) \in X_h \times Y_h$, the iterative equations (4.13)–(4.14) are consistent with Problem (4.2) in the sense that if $(\mathbf{w}_h^0, \eta_h^0) = (\mathbf{u}_h^n, C_h^n)$ with $(\mathbf{u}_h^n, p_h^n, C_h^n)$ being a solution of Eq. (4.2), then for all $i \ge 1$, $(\mathbf{w}_h^i, \xi_h^i, \eta_h^i) = (\mathbf{u}_h^n, p_h^n, C_h^n)$.

*Proof* For a given $(\mathbf{u}_h^{n-1}, C_h^{n-1}) \in X_h \times Y_h$, let $(\mathbf{u}_h^n, p_h^n, C_h^n)$ be the solution of (4.2). The proof of the theorem is done by induction on the index $i$ for system (4.13) and (4.14). We assume that $\mathbf{w}_h^{i-1} = \mathbf{u}_h^n$ and we will show that $(\mathbf{w}_h^i, \xi_h^i, \eta_h^i) = (\mathbf{u}_h^n, p_h^n, C_h^n)$.

We take the difference between the second equation of (4.2) and (4.14) and use Remark 4.1, this yields the following for $r_h = C_h^n - \eta_h^i$:

$$ \frac{1}{k} \|C_h^n - \eta_h^i\|^2 + \alpha(\nabla(C_h^n - \eta_h^i), \nabla(C_h^n - \eta_h^i)) = 0. $$

We deduce that $\eta_h^i = C_h^n$.

Next, we take the difference between the first equations of (4.2) and (4.13) and use Remark 4.1. This gives us, for $\mathbf{v}_h = \mathbf{u}_h^n - \mathbf{w}_h^i$:

$$ \frac{1}{k} \|\mathbf{u}_h^n - \mathbf{w}_h^i\|^2 + a(\mathbf{u}_h^n - \mathbf{w}_h^i, \mathbf{u}_h^n - \mathbf{w}_h^i)) + \int_S |K \mathbf{u}_h^n h \boldsymbol{\tau}|^{s-2} |K(\mathbf{u}_h^i \boldsymbol{\tau} - \mathbf{w}_h^i \boldsymbol{\tau})|^2 d\sigma = 0. $$

All the terms of the left-hand side of the last equation are positive and we deduce that $\mathbf{w}_h^i = \mathbf{u}_h^n$. The inf-sup condition (4.1) concludes that $\xi_h^i = p_h^n$.

**Theorem 4.11** We assume that the solution $(\mathbf{u}_h^n, p_h^n, C_h^n)$ of Problem (4.2) is unique. Let $(\mathbf{w}_h^i, \xi_h^i, \eta_h^i)$ the iterative solution of Problem (4.13)–(4.14). There exists subsequences of $(\mathbf{w}_h^i)_i$, $(\xi_h^i)_i$ and $(\eta_h^i)_i$ denoted also with the index $i$, such that (as we are in finite dimensional spaces):

$$ \lim_{i \to +\infty} \mathbf{w}_h^i = \mathbf{u}_h^n \quad \text{in } X_h, $$

$$ \lim_{i \to +\infty} \xi_h^i = p_h^n \quad \text{in } M_h, $$

and

$$ \lim_{i \to +\infty} \eta_h^i = C_h^n \quad \text{in } Y_h. $$
**Proof** The bounds (4.15) and (4.16) allows us to deduce that subsequences of $(\mathbf{w}_h^i)_i$, $(\xi_h^i)_i$ and $(\eta_h^i)_i$ denoted also with the index $i$, such that (as we are in finite dimensional spaces):

$$ \lim_{i \to +\infty} \mathbf{w}_h^i = \tilde{\mathbf{u}}_h^n \text{ in } X_h, $$

$$ \lim_{i \to +\infty} \xi_h^i = \tilde{p}_h^n \text{ in } M_h, $$

and

$$ \lim_{i \to +\infty} \eta_h^i = \tilde{C}_h^n \text{ in } Y_h. $$

Let us now show that $\tilde{\mathbf{u}}_h^n = \mathbf{u}_h^n, \tilde{p}_h^n$ and $\tilde{C}_h^n = C_h^n$ are the solutions of Problem (4.2). We pass to the limit with $i$ in System (4.13)-(4.14). The difficulty lies at the level of the third and fifth terms of the left-hand side of Eq. (4.13) and its right-hand side term, and also the third term of Eq. (4.14). Let us check the convergence of every term:

The continuity of $f$ and the convergence of $\eta_i$ with respect to $i$ ensure the convergence of the right-hand side of Eq. (4.13). The convergence of the third and fifth terms of the left-hand side of (4.13) can be treated exactly as in the proof of Theorem 5.1 of [13]. Let us first treat the convergence of the term $d_u(\mathbf{w}_h^{i-1}, \mathbf{w}_h^i, \mathbf{v}_h)$, which can be written by using the definition of $d_u$ as

$$ \begin{aligned} & d_u(\mathbf{w}_h^{i-1}, \mathbf{w}_h^i, \mathbf{v}_h) - d_u(\tilde{\mathbf{u}}_h^n, \tilde{\mathbf{u}}_h^n, \mathbf{v}_h) \\ & = d_u(\mathbf{w}_h^{i-1} - \tilde{\mathbf{u}}_h^n, \mathbf{w}_h^i, \mathbf{v}_h) + d_u(\tilde{\mathbf{u}}_h^n, \mathbf{w}_h^i - \tilde{\mathbf{u}}_h^n, \mathbf{v}_h). \end{aligned} $$

We obtain the following bound

$$ |d_u(\mathbf{w}_h^{i-1}, \mathbf{w}_h^i, \mathbf{v}_h) - d_u(\tilde{\mathbf{u}}_h^n, \tilde{\mathbf{u}}_h^n, \mathbf{v}_h)| \le c \left( \| \mathbf{w}_h^{i-1} - \tilde{\mathbf{u}}_h^n \|_{1,\Omega} \| \mathbf{w}_h^i \|_{1,\Omega} \| \mathbf{v}_h \|_{1,\Omega} \right. \\ \left. + \| \mathbf{w}_h^i - \tilde{\mathbf{u}}_h^n \|_{1,\Omega} \| \tilde{\mathbf{u}}_h^n \|_{1,\Omega} \| \mathbf{v}_h \|_{1,\Omega} \right). $$

This bound, combined with the convergence of $\mathbf{w}_h^i$ to $\tilde{\mathbf{u}}_h^n$ in $X_h$, ensures the convergence of the third term of the left-hand side of (4.13).

Let us now verify the convergence of the fifth term of the left-hand side of (4.13). From the convergence of $\mathbf{w}_h^i$ to $\tilde{\mathbf{u}}_h^n$ in $X_h$ and the continuity of the trace theorem on $S$, we deduce that

$$ \mathbf{w}_h^i|_S \xrightarrow{i \to +\infty} \tilde{\mathbf{u}}_h^n|_S \quad \text{in } L^2(S). $$

So, we deduce that (see [24, chapter 2])

$$ \lim_{i \to +\infty} (|K \mathbf{w}_{h\tau}^{i-1}|^{s-2} K \mathbf{w}_{h\tau}^i, K \mathbf{v}_{h\tau})_S = (|K \tilde{\mathbf{u}}_{h\tau}^n|^{s-2} K \tilde{\mathbf{u}}_{h\tau}^n, K \mathbf{v}_{h\tau})_S. $$

Finally, the third term of (4.14) can be treated as the third term of Eq. (4.13). We conclude that $(\tilde{\mathbf{u}}_h^n, \tilde{p}_h^n, \tilde{C}_h^n)$ satisfies system (4.2), yielding the desired result. $\square$

**Remark 4.12** We note that a simple idea is to replace the iterative problem (4.13) ... (4.14) by the following one:

for every $n \in \{1, ..., N\}$, knowing $(\mathbf{u}_h^{n-1}, C_h^{n-1}) \in X_h \times Y_h$, we compute $(\mathbf{u}_h^n, p_h^n, C_h^n) \in (X_h, M_h, Y_h)$ such that for $(\mathbf{v}_h^n, q_h^n, r_h) \in (X_h, M_h, Y_h)$

1. Step 1: Compute first $(\mathbf{u}_h^n, p_h^n)$ such that:

$$ \begin{aligned} & \frac{1}{k}(\mathbf{u}_h^n - \mathbf{u}_h^{n-1}, \mathbf{v}_h) + a(\mathbf{u}_h^n, \mathbf{v}_h) + d_u(\mathbf{u}_h^{n-1}, \mathbf{u}_h^n, \mathbf{v}_h) - (p_h^n, \text{div } \mathbf{v}_h) \\ & + \int_S |K \mathbf{u}_{h\tau}^{n-1}|^{s-2} K \mathbf{u}_{h\tau}^n \cdot K \mathbf{v}_{\tau} d\sigma = (\mathbf{f}^n(C_h^{n-1}), \mathbf{v}_h), \quad (4.18) \\ & (q_h, \text{div } \mathbf{u}_h^n) = 0, \end{aligned} $$
2. step 2: knowing $\mathbf{u}_h^n$ obtained in the first step, we compute $C_h^n$ such that:

$$ \frac{1}{k}(C_h^n - C_h^{n-1}, r_h) + \alpha(\nabla C_h^n, \nabla v_h) + d_C(\mathbf{u}_h^n, C_h^n, r_h) + r_0(C_h^n, r_h) = (\mathbf{g}^n, r_h). $$

We note that one of the difficulties in the scheme (4.18) is the derivation of the corresponding a priori error bounds because the term $\int_S |K \mathbf{u}_h^{n-1}|^{s-2} K \mathbf{u}_h^n \cdot K v_\tau d\sigma$ is difficult to manipulate. That is why we introduced the method (4.13).

# 5 Numerical simulations and concluding remarks

In this section, we perform several numerical simulations using the FreeFem++ code (see [18]) in two-dimensions. We perform two validation tests. The first one is a simple case where we analyse the numerical results obtained using the scheme (4.2). It is important to note that the solution of the non-linear discrete scheme (4.2) is approximated using the iterative scheme (4.13). Subsequently, we consider a more complex case, the "Lid driven cavity flow" and show the corresponding numerical results.

For the numerical approximation of the solution of (3.1), we use the iterative scheme (4.13) and (4.14) as follows: for a given mesh step $h$ and a specific time step $n$, we compute the numerical solution $(\mathbf{u}_h^n, p_h^n, C_h^n) \in X_h \times M_h \times Y_h$ by solving the iterative scheme (4.13) and (4.14) with $\mathbf{w}_h^0 = \mathbf{u}_h^{n-1}$. To ensure the convergence of the iterative scheme, we employ a stopping criterion based on a predefined tolerance (see [3]):

$$ \frac{\|\mathbf{w}_h^i - \mathbf{w}_h^{i-1}\|_1}{\|\mathbf{w}_h^i\|_1} \le C \min(h, k), \quad (5.1) $$

where $C$ is a positive constant (but in this work, we take $C = 1$).

## 5.1 Coupling in a square domain

We consider the square $\Omega = (0, 1)^2$ and the time interval $[0, T]$ where $T = 1$. To discretize the domain, each edge of $\partial\Omega$ is divided into $M$ equal segments, resulting in a mesh with $2M^2$ triangles. Therefore, the mesh step size is given by $h = \frac{1}{M}$. For the time discretization, we choose the time step size to be equal to the mesh step size, $k = h$. As a result, the a priori error bound given in Theorem 4.6 can be bounded by $h$, since $k = h$ in this case.

We set $K = I$, $\nu = 1$ and

$$ \mathbf{f}_0 = (1 + (x^2 + y^2) \sin(t(x^2 + y^2)), txy \cos(txy))^T, $$

$$ \mathbf{f}_1(C) = (2 + \sin(C), 1 + e^C) $$

and

$$ g = x^2(x - 1)^2 y^2(y - 1)^2. $$

We impose the anisotropic slip boundary condition (1.4), on right and top boundaries of $\Omega$, that is on $S = \{(x, y) | x = 1 \text{ or } y = 1\} = \bar{\Gamma}_2 \cup \bar{\Gamma}_3$. For the remaining boundaries, i.e, the bottom segment and the left segment $\Gamma = \bar{\Gamma}_1 \cup \bar{\Gamma}_4$, we impose the homogeneous Dirichlet boundary condition $\mathbf{u} = \mathbf{0}$. Furthermore, we set $C = 0$ on $\Gamma$. As for the initial conditions, we set $\mathbf{u}_0 = \mathbf{0}$ and $C_0 = 0$.

To analyze the convergence rate of the finite element solution $(\mathbf{u}_h^n, p_h^n, C_h^n)$ given by Theorem 4.6, we need to compare it with an approximate reference solution. Since we don't
<figure><img src="image_3.png" /><figcaption>**Fig. 1** The a priori error estimates with respect to the mesh step h in logarithmic scale for s = 1., 1.2, 1.5, 1.8, 2</figcaption></figure>

have the exact solution for this particular case, we will approximate it by computing the numerical solution of the iterative scheme (4.13) for M = 300, denoted as (u_r, p_r, C_r). The reference solution depends on the parameter s ∈ [1, 2].

Next, we compute the numerical solution of (4.2) approximated by (4.13) and (4.14), where M takes the values 5, 8, 10, 12, 15, 20, 25 and for k = h = 1/M. For all the numerical tests treated in this part, Algorithm (4.13) and (4.14) is stopped for i ≤ 2 by using the stopping criterion (5.1). Figure 1 shows the curves corresponding to the following error:

$$E_{rr}^2 = \frac{\sum_{n=1}^{M} k|\mathbf{u}_h^n - \mathbf{u}_r(t_n)|_{H^1(\Omega)}^2 + \sum_{n=1}^{M} k|p_h^n - p_r(t_n)|_{L^2(\Omega)}^2 + \sum_{n=1}^{M} k|C_h^n - C_r(t_n)|_{H^1(\Omega)}^2}{\sum_{n=1}^{M} k|\mathbf{u}_r(t_n)|_{H^1(\Omega)}^2 + \sum_{n=1}^{M} k|p_r(t_n)|_{L^2(\Omega)}^2 + \sum_{n=1}^{M} k|C_r(t_n)|_{H^1(\Omega)}^2} \quad (5.2)$$

with respect to the mesh step h = 1/M in logarithmic scale for s = 1., 1.2, 1.5, 1.8, 2. and with the stopping criteria (5.1). The computed slopes of the corresponding lines are 1.11 (for s = 1.), 1.1 (for s = 1.2), 1.07 (for s = 1.5 and 1.8), and 1.08 (for s = 2.), which are in agreement with the theoretical finding.

## 5.2 Lid driven cavity

In this paragraph we focus on the numerical simulation of the Lid Driven cavity problem, which is a well-known example studied extensively in various works (see for instance [17, 26]). The fluid is confined in the domain $\Omega = (0, 1)^2$ with the velocity u satisfying $u = 0$ on $\Gamma = \{(x, y) : x = 0 \text{ or } y = 0\}$, the relation (1.4) on the right boundary (part S of $\partial\Omega$), and the following Dirichlet condition on the top boundary:

$$\mathbf{u}_{\tau} = (800(1 + th(-t))x^2(x - 1)^2, 0)^T \quad \text{on } y = 1, 0 < x < 1. \quad (5.3)$$

One observes that for $t \in [0, T]$, $\mathbf{u}_{\tau}$ attains its maximum $(50(1 + th(-T)), 0)^T$ at the center $x = 1/2$. The boundary condition (5.3) has the specific purpose of preventing local singularities at the top-right and top-left corners of the domain. By imposing this boundary condition, both the velocity and velocity gradient are guaranteed to vanish at these corners. For the heat equation, we take $C = 1$ on the right boundary of $\Omega$ (corresponding to $x = 1$) and homogeneous Dirichlet condition $C = 0$ on the other part of $\partial\Omega$. For the numerical results considered in this example, we set $N = 100$, $s = 1.5$, $K = 1$, $f = 0$, $g = 0$ and $v = \alpha = r_0 = 1$. Figure 2 displays the distribution of the velocity and the temperature in the domain $\Omega$. Furthermore the dependence of the temperature distribution on the velocity
<figure><img src="image_4.png" /><figcaption>**Fig. 2** The numerical velocity (left) and concentration (right) for the driven-cavity</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 3** The a priori error estimates with respect to the mesh step h in logarithmic scale for the lid-driven cavity ($s = 1.2, 1.8$)</figcaption></figure>

is clearly highlighted. In order to show the a priori error estimate and as we don't have the exact solution also in this case, we will approximate it by computing the numerical solution of the iterative scheme (4.13) and (4.14) for $M = 300$, and then consider $M = 5, 8, 10, 12, 15, 20, 25$ and $k = h = 1/M$ to compute the error $E_{rr}^2$ given by (5.2). Figure 3 shows the a priori error estimate ($E_{rr}^2$) with respect to the mesh step $h$ in logarithmic scale for the cases $s = 1.2$ and $s = 1.8$. The corresponding slope are equal to 1.12 which is close to the theoretical one.

# 6 Conclusion

In this paper we treated the fully discrete time dependent Navier-Stokes equations under anisotropic slip boundary condition coupled with the convection-diffusion-reaction equation. We formulated a weak formulation and construct the weak solutions by using the “approximation approach”. Moreover, we introduce a discrete scheme and establish the convergence through the a priori error estimates. Next, we propose an iterative scheme and study the convergence of the discrete solution. Finally, we perform several numerical simulations using the FreeFem++ software to validate the theoretical results.

**Data Availability** not applicable.
# Declarations

**Conflict of interest** The authors declare that they have no conflict of interest.

# References

1. Agroum, R., Aouadi, S.M., Bernardi, C., Satouri, J.: Spectral discretization of the Navier-Stokes equations coupled with the heat equation. ESAIM. Math. Model. Numer. Anal. **49**(3), 621–639 (2015)

2. Agroum, R., Bernardi, C., Satouri, J.: Spectral discretization of the time-dependent Navier-Stokes problem coupled with the heat equation. Appl. Math. Comput. **49**, 59–82 (2015)

3. Aldbaissy, R., Chalhoub, N., Djoko, J.K., Sayah, T.: Full discretization of the time dependent Navier-Stokes equations with anisotropic slip boundary condition. Int. J. Numer. Anal. Model. **20**(4), 497–517 (2023)

4. Aldbaissy, R., Hecht, F., Mansour, G., Sayah, T.: A full discretization of the time-dependent Boussinesq (buoyancy) model with nonlinear viscosity. Calcolo **55**(44), 285 (2018)

5. Arnold, D., Brezzi, F., Fortin, M.: A stable finite element for the Stokes equations. Calcolo **21**, 337–344 (1984)

6. Bernardi, B., Girault, V.: A local regularisation operation for triangular and quadrilateral finite elements. SIAM J. Numer. Anal. **35**, 1893–1916 (1998)

7. Brezis, H.: Functional Analysis. Sobolev Spaces and Partial Differential Equations. Springer, New York (2010)

8. Busse, A., Sandham, N.D.: Influence of an anisotropic slip-length boundary condition on turbulent channel flow. Phys. Fluids **24**, 055111 (2012)

9. Ciarlet, P.G.: Three Dimensional Elasticity. North Holland, Amsterdam (1988)

10. Clément, P.: Approximation by finite element functions using local regularization. RAIRO Anal. Numér. **9**, 77–84 (1975)

11. Cooper, A.J., Harris, J.H., Garrett, S.J., Ozkan, M., Thomas, P.J.: The effect of anisotropic and isotropic roughness on the convective stability of the rotating disk boundary layer. Phys. Fluids **27**, 014107 (2015)

12. Dakroub, J., Faddoul, J., Omnes, P., Sayah, T.: A posteriori error estimates for the time-dependent Navier-Stokes system coupled with the convection-diffusion-reaction equation. Adv. Comput. Math. **49**, 67 (2023)

13. Djoko, J.K., Koko, J., Mbehou, M., Sayah, T.: Stokes and Navier-Stokes equations under power law slip boundary condition: numerical analysis. Comput. Math. Appl. **128**, 198–213 (2022)

14. Gerstner, P., Heuveline, V.: Finite element approximation of dielectrophoretic force driven flow problems. ESAIM: M2AN **57**(3), 1691–1729 (2023)

15. Girault, V., Raviart, P.-A.: Finite element methods for the Navier-Stokes equations. In: Theory and Algorithms in Springer Series in Computational Mathematics, vol. 5. Springer, Berlin (1986)

16. Glowinski, R., Morrocco, A.: Sur l'approximation par éléments finis d ordre un et la résolution par penalisation-dualite d une classe de problèmes de Dirichlet non线aires. Rairo Ser. Rouge-Anal. Numer. **9**, 41–76 (1975)

17. Glowinski, R., Guidoboni, G., Pan, T.-W.: Wall-driven incompressible viscous flow in a two-dimensional semi-circular cavity. J. Comput. Phys. **216**, 76–91 (2006)

18. Hecht, H.: New development in FreeFem++. J. Numer. Math. **20**, 251–266 (2012)

19. Jao, H.-C., Chang, K.-M., Chu, L.-M., Li, W.-L.: A lubrication theory for anisotropic slips and flow rheology. Tribol. Trans. **59**(2), 252–266 (2016)

20. Khan, N.A., Sohail, A., Sultan, F.: Effect of anisotropic slip and magnetic field on the flow and heat transfer of Eyring Powell fluid over an infinite rotating disk. Int. J. Fluid Mech. Res. **44**(3), 257–273 (2017)

21. Le Roux, C.: Flows of incompressible viscous liquids with anisotropic wall slip. J. Math. Anal. Appl. **465**, 723–730 (2018)

22. Le Roux, C.: On the Navier-Stokes equations with anisotropic wall slip conditions. Appl. Math. (2022). https://doi.org/10.21136/AM.2021.0079-21

23. Liakos, A.: Discretization of the Navier-Stokes equations with slip boundary condition. Numer. Methods Partial Differ. Equ. Int. J. **17**, 26–42 (2001)

24. Lions, J-L.: Quelques méthodes de résolution des problèmes aux limites non linéaires. Paris, Dunod (1968)
25. Lu, J.G., Jang, H.K., Lee, S.B., Hwang, W.R.: Characterization on the anisotropic slip for flows over unidirectional fibrous porous media for advanced composites manufacturing. Compos. Part A Appl. Sci. Manuf. **100**, 9–19 (2017)

26. Owens, R.G.: Computational Rheology. Imperial College Press, London (2002)

27. Sandri, D.A.: Sur l'approximation des ecoulements numeriques quasi-Newtoniens dont la viscosite obeit a la loi de puissance ou de Carreau. M2AN **27**, 131–155 (1993)

28. Scott, L.R., Zhang, S.: Finite element interpolation of nonsmooth functions satisfying boundary conditions. Math. Comput. **54**, 483–493 (1990)

29. Temam, R.: Finite element interpolation of nonsmooth functions satisfying boundary conditions. Math. Comput. **54**, 483–493 (1990)

30. Urquiza, J.M., Garon, A., Farinas, M.-I.: Weak imposition of the slip boundary condition on curved boundaries for Stokes flow. J. Comput. Phys. **256**, 748–767 (2014)

31. Verfürth, R.: Finite element approximation on incompressible Navier–Stokes equations with slip boundary condition. Numer. Math. **50**, 697–721 (1986)

32. Verfürth, R.: Finite Element Approximation of incompressible Navier–Stokes equations with slip boundary condition II. Numer. Math. **59**, 615–636 (1991)

33. Vidyasagar, M.: Nonlinear Systems Analysis, 2nd edn. Prentice Hall, Englewood Cliffs (1993)

34. Zhou, G., Kashiwabara, T., Oikawa, I.: Penalty method for the stationary Navier–Stokes problems under the slip boundary condition. J. Sci. Comput. **68**, 339–374 (2016)

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.