# Mortar finite element discretization of the time dependent nonlinear Darcy's equations

Karima Amoura, Christine Bernardi, Samira Saadi

► To cite this version:

Karima Amoura, Christine Bernardi, Samira Saadi. Mortar finite element discretization of the time dependent nonlinear Darcy's equations. Calcolo, 2016, 53, pp.597-619. hal-00834721

HAL Id: hal-00834721

https://hal.science/hal-00834721

Submitted on 18 Jun 2013

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.

Public Domain
# Mortar finite element discretization of the time dependent nonlinear Darcy's equations

by Karima Amoura<fnref n="1" />, Christine Bernardi<fnref n="2" />, and Samira Saadi<fnref n="3" />

**Abstract:** We consider the non stationary flow of a viscous incompressible fluid in a rigid homogeneous porous medium provided with mixed boundary conditions. Since the medium is nonhomogeneous, its permeability is only piecewise continuous. We are thus led to use the mortar method to handle these discontinuities. We propose a space and time discretization of the full system. We prove optimal a priori error estimates, which confirms the interest of the discretization.

**Résumé:** Nous considérons l'écoulement instationnaire d'un fluide visqueux incompressible dans un milieu poreux rigide avec conditions aux limites mixtes. Comme le milieu n'est pas homogène, sa perméabilité est seulement continue par morceaux, ce qui nous amène à utiliser la méthode d'éléments finis avec joints pour traiter ces discontinuités. Nous proposons une discrétisation en temps et en espace du système complet. Nous prouvons des estimations d'erreur a priori qui s'avèrent optimales, ce qui confirme l'intérêt de la discrétisation.

<footnote n="1">¹ Laboratoire de Modélisation Mathématique et Simulation Numérique,
Université Badji-Mokhtar, Faculté des Sciences, Département de Mathématiques,
B.P. 12, 23000 Annaba, Algérie.</footnote>

<footnote>e-mail address: amouradz@yahoo.fr</footnote>

<footnote n="2">² Laboratoire Jacques-Louis Lions, C.N.R.S. & Université Pierre et Marie Curie,
B.C. 187, 4 place Jussieu, 75252 Paris Cedex 05, France.</footnote>

<footnote>e-mail address: bernardi@ann.jussieu.fr</footnote>

<footnote n="3">³ Laboratoire d'Analyse Numérique, Optimisation et Statistique,
Université Badji-Mokhtar, Faculté des Sciences, Département de Mathématiques,
B.P. 12, 23000 Annaba, Algérie.</footnote>

<footnote>e-mail address: Samira.Saadi@univ-annaba.org</footnote>
<figure><img src="image_1.png" /></figure>
# 1. Introduction.

Let $\Omega$ be a bounded connected domain in $\mathbb{R}^d$, $d = 2$ or $3$, with a Lipschitz-continuous boundary $\partial\Omega$. We assume that this boundary is divided into two disjoint parts $\Gamma_1$ and $\Gamma_2$ such that $\partial\Gamma_1$ and $\partial\Gamma_2$ are Lipschitz-continuous submanifolds of $\partial\Omega$. Let also $T$ be a positive real number. We are interested in studying the following model, first suggested by K.R. Rajagopal [14], for the flow of a viscous incompressible fluid in a rigid saturated porous medium

$$ 
\left\{
\begin{array}{ll}
\partial_t u + \alpha(\cdot, p) u + \text{grad} p = f & \text{in } \Omega \times ]0, T[, \\
\text{div } u = 0 & \text{in } \Omega \times ]0, T[, \\
p = p_1 & \text{on } \Gamma_1 \times ]0, T[, \\
u \cdot n = g_2 & \text{on } \Gamma_2 \times ]0, T[, \\
u|_{t=0} = u_0 & \text{in } \Omega,
\end{array}
\right.
\quad (1.1) $$

where the unknowns are the velocity $u$ and the pressure $p$ of the fluid. Indeed this system models the flow of a viscous incompressible fluid in a rigid saturated porous medium. It seems rather realistic from a mechanical point of view, at least when the coefficient $\alpha$ which is linked to the permeability of the medium and the viscosity of the fluid satisfies the following conditions:

* When $p$ is fixed, the function: $x \mapsto \alpha(x, p)$ is piecewise constant. Indeed, the previous equations model the flow in a porous medium, but underground porous media are most often nonhomogeneous, made of several components such as clay or calcalenite, and the values of $\alpha$ for these different components are very different (the ratio of the maximal value to the minimal one is often of order $10^6$, see [6, Table 1] for instance);

* Where the pressure presents high variations, for instance induced by the boundary values, the coefficient $\alpha$ depends on these values in an exponential way. We refer to [14] for details on the way of deriving this model.

The data are a density of body forces $f$ (which most often represents the gravity action), the boundary pressure $p_1$ and the boundary flux $g_2$, and also the initial condition $u_0$.

We refer to [3] for a first work on this nonlinear problem in the steady case and for a simple three-dimensional geometry. On the other hand, piecewise constant coefficients $\alpha$ with possible high variations have been handled in [1] and [6], while time-dependent Darcy's equations have been studied in [5] and more recently in [10] for a piecewise constant coefficient $\alpha$. However, up to our knowledge, combining the three difficulties of nonlinearity, nonhomogeneity and time dependence has not been considered up to now. We also recall that, even in the basic case of a constant coefficient $\alpha$ (and also for the stationary problem), system (1.1) admits several variational formulations. We have chosen to work with one of them in view of the discretization and according to the approach proposed in [2]. We first write this space variational formulation which is equivalent to problem (1.1) when the partition of $\partial\Omega$ into $\Gamma_1$ and $\Gamma_2$ is smooth enough. Next we prove the existence of a solution by handling separately the boundary conditions, uncoupling the unknowns and applying the Cauchy-Lipschitz theorem for the velocity, the Brouwer fixed point theorem for the pressure. We also establish a uniqueness result.
The time discretization that we propose relies on the Euler's implicit scheme, however we have decided to treat the nonlinear term in an explicit way to make the implementation both simpler and less expensive. For the space discretization, we combine the mortar element method introduced in [7] with the use of nonconforming finite elements proposed in [2, §4], which involves the elements due to M. Crouzeix and P.-A. Raviart [9]. Such a discretization has already been studied in [6] in the simpler case of a linear steady problem with piecewise constant coefficient $\alpha$. We describe the corresponding discrete problem and prove the existence of a solution. Next, we prove a priori error estimates for this problem which turn out to be fully optimal.

For simplicity, we work with a fixed time step and a triangulation which does not depend on the time. Note however that a very similar mortar finite element discretization has been studied in [6] in the linear stationary case and that, by combining the arguments in this paper with those in [4], we could perform the a posteriori analysis of the discretization for variable time steps and triangulations, see also [5] for similar results. However the adaptation of either the time step or the mesh seems useless since there is neither birth nor transport of singularities of the solution. So we omit this a posteriori analysis.

An outline of the paper is as follows.

• In Section 2, we write the variational formulation of problem (1.1) and investigate its wellposedness in appropriate Sobolev spaces.

• Section 3 is devoted to the description of the time semi-discrete problem and of the full discrete problem. We also prove the existence of a solution.

• In Section 4, we perform the a priori analysis of the discretization and prove optimal error estimates.

**Acknowledgement:** We thank Professor Frédéric Hecht for his clever comments on this paper which allow us to greatly improve it.
## 2. Variational formulation and well-posedness.

From now on, we assume that $\Gamma_1$ has a positive measure in $\partial\Omega$. Before writing the variational formulation of problem (1.1), we make precise the assumptions that we make on the coefficient $\alpha$: For each $\xi$ in $\mathbb{R}$, the mapping: $x \mapsto \alpha(x, \xi)$ is piecewise constant, more precisely equal to a constant $\alpha_i(\xi)$ on each subdomain $\Omega_i$, $1 \le i \le I$, of a partition of $\Omega$ without overlap. We make the further and realistic assumption (the medium is rigid) that this partition is independent of $\xi$. Next, we assume that each mapping: $\xi \mapsto \alpha_i(\xi)$, $1 \le i \le I$, is

(i) bounded and positive on $\mathbb{R}$: There exist positive constants $\alpha_{ib}$ and $\alpha_{i\sharp}$ such that

$$ \forall x \in \Omega_i, \forall \xi \in \mathbb{R}, \quad \alpha_{ib} \le \alpha_i(x, \xi) \le \alpha_{i\sharp}; \qquad (2.1) $$

(ii) Lipschitz-continuous on $\mathbb{R}$: There exists a constant $\alpha_{i*}$ such that

$$ \forall x \in \Omega_i, \forall \xi_1 \in \mathbb{R}, \forall \xi_2 \in \mathbb{R}, \quad |\alpha_i(x, \xi_1) - \alpha_i(x, \xi_2)| \le \alpha_{i*} |\xi_1 - \xi_2|. \qquad (2.2) $$

We denote by $\alpha_b$ and $\alpha_{\sharp}$, the min of the $\alpha_{ib}$ and the max of the $\alpha_{i\sharp}$, $1 \le i \le I$, respectively. The main idea for the introduction of this notation is that all ratios $\alpha_{i\sharp}/\alpha_{ib}$ are of reasonable size and will not be taken into account in the next estimates while the ratio $\alpha_{\sharp}/\alpha_b$ can be very large.

In what follows, the scalar product defined on $L^2(\Omega)$ or $L^2(\Omega)^d$ is denoted by $(\cdot, \cdot)$. We use the whole scale of Sobolev spaces $H^s(\Omega)$, $s \ge 0$, equipped with the norm $\|\cdot\|_{H^s(\Omega)}$ and seminorm $|\cdot|_{H^s(\Omega)}$, and their subspaces $H_0^s(\Omega)$. In view of the boundary conditions in (1.1), we introduce the space

$$ H_{\triangle}^{1}(\Omega) = \{q \in H^{1}(\Omega); q = 0 \text{ on } \Gamma_{1}\}. \qquad (2.3) $$

Since the traces on $\Gamma_2$ of functions in $H_{\triangle}^1(\Omega)$ belong to $H_{00}^{1/2}(\Gamma_2)$ (see [13, Chap. 1, Th. 11.7] for the definition of this last space), we introduce its dual space $H_{00}^{1/2}(\Gamma_2)'$ and denote by $(\cdot, \cdot)_{\Gamma_2}$ the duality pairing between $H_{00}^{1/2}(\Gamma_2)'$ and $H_{00}^{1/2}(\Gamma_2)$.

For any separable Banach space $E$ equipped with the norm $\|\cdot\|_E$, we denote by $C^0(0, T; E)$ the space of continuous functions from $[0, T]$ with values in $E$. For each $s \ge 0$, we also introduce the space $H^s(0, T; E)$ in the following way: When $s$ is an integer, it is the space of measurable functions on $]0, T[$ with values in $E$ such that the mappings: $v \mapsto \|\partial_t^\ell v\|_E$, $0 \le \ell \le s$, are square-integrable on $]0, T[$; otherwise, it is defined by interpolation between $H^{\lfloor s \rfloor+1}(0, T; E)$ and $H^{\lfloor s \rfloor}(0, T; E)$, where $\lfloor s \rfloor$ stands for the integer part of $s$.

Assuming the data $(f, u_0, p_1, g_2)$ sufficient smooth for the next equations to make sense, we now consider the following variational problem:

Find $(u, p)$ in $C^0(0, T; L^2(\Omega)^d) \times L^2(0, T; H^1(\Omega))$ such that

$$ \begin{array}{ll} u(x, 0) = u_0(x) & \text{for a.e. } x \text{ in } \Omega \\ \text{and} & p(x, t) = p_1(x, t) \text{ for a.e. } (x, t) \text{ in } \Gamma_1 \times ]0, T[, \end{array} \qquad (2.4) $$
and, for a.e. $t$ in $]0, T[$,

$$ 
\begin{aligned}
& \forall \mathbf{v} \in L^2(\Omega)^d, \quad (\partial_t \mathbf{u}, \mathbf{v}) + (\alpha(\cdot, p) \mathbf{u}, \mathbf{v}) + (\mathbf{v}, \mathrm{grad} \, p) = (\mathbf{f}, \mathbf{v}), \\
& \forall q \in H_{\triangle}^1(\Omega), \quad (\mathbf{u}, \mathrm{grad} \, q) = \langle g_2, q \rangle_{\Gamma_2}.
\end{aligned} 
\quad (2.5) $$

Indeed standard density arguments yield the equivalence of this problem with system (1.1).

**Proposition 2.1.** Assume that the partition of $\partial\Omega$ into $\Gamma_1$ and $\Gamma_2$ is sufficiently smooth for $\mathcal{D}(\Omega \cup \Gamma_2)$ to be dense in $H_{\triangle}^1(\Omega)$. Then, problems (1.1) and (2.4) – (2.5) are equivalent in the sense that a pair $(\mathbf{u}, p)$ in $\mathcal{C}^0(0, T; L^2(\Omega)^d) \times L^2(0, T; H^1(\Omega))$ is a solution of (1.1) (in the distribution sense) if and only if it is a solution of (2.4) – (2.5).

To prove the well-posedness of problem (2.4) – (2.5), we begin with a priori estimates on the solution. To do this, we consider the norm on $L^2(\Omega)$

$$ \|v\|_{\alpha} = \left( \sum_{i=1}^{I} \alpha_{ib} \int_{\Omega_i} v^2(\mathbf{x}) \, d\mathbf{x} \right)^{\frac{1}{2}}. \quad (2.6) $$

Its extension to $L^2(\Omega)^d$ is obvious. We also define its “dual” norm

$$ \|v\|_{\alpha^{-1}} = \left( \sum_{i=1}^{I} \alpha_{ib}^{-1} \int_{\Omega_i} v^2(\mathbf{x}) \, d\mathbf{x} \right)^{\frac{1}{2}}. \quad (2.7) $$

Analogous norms were introduced in [1] in the simpler case of a piecewise constant coefficient $\alpha$, see also [6, §2]. It can also be noted that, owing to the Poincaré–Friedrichs inequality, the mapping: $q \mapsto \|\mathrm{grad} \, q\|_{\alpha^{-1}}$ is a norm on $H_{\triangle}^1(\Omega)$.

**Lemma 2.2.** The following inf-sup condition holds

$$ \forall q \in H_{\triangle}^{1}(\Omega), \quad \sup_{v \in L^{2}(\Omega)^{d}} \frac{(v, \mathrm{grad} \, q)}{\|v\|_{\alpha}} \geq \|\mathrm{grad} \, q\|_{\alpha^{-1}}. \quad (2.8) $$

**Proof:** When taking $v$ equal to $\alpha_{ib}^{-1} \mathrm{grad} \, q$ on each $\Omega_i$, $1 \le i \le I$, we obtain

$$ (v, \mathrm{grad} \, q) = \|\mathrm{grad} \, q\|_{\alpha^{-1}}^2 \quad \text{and} \quad \|v\|_{\alpha} = \|\mathrm{grad} \, q\|_{\alpha^{-1}}, $$

which yields the desired condition.

The next statement requires the introduction of the kernel

$$ \mathcal{V}(\Omega) = \{v \in L^2(\Omega)^d; \forall q \in H_{\triangle}^1(\Omega), (v, \mathrm{grad} \, q) = 0\}. \quad (2.9) $$

It is readily checked that this kernel admits the following characterization

$$ \mathcal{V}(\Omega) = \{v \in L^2(\Omega)^d; \operatorname{div} v = 0 \text{ in } \Omega \text{ and } v \cdot n = 0 \text{ on } \Gamma_2\}. \quad (2.10) $$
**Lemma 2.3.** For any data

$$f \in L^2(0, T; L^2(\Omega)^d), \quad u_0 \in L^2(\Omega)^d, \quad \text{and} \quad g_2 \in \mathcal{C}^0(0, T; H_{00}^{1/2}(\Gamma_2)') \quad (2.11)$$

the following a priori estimate holds for any solution $(u, p)$ of problem (2.4) – (2.5)

$$\begin{align}
\sup_{0 \le t \le T} \|\mathbf{u}(t)\|_{L^2(\Omega)} + \left( \int_0^T \|\mathbf{u}(t)\|_{\alpha}^2 dt \right)^{\frac{1}{2}} & \tag{2.12} \\
& \le 3 \left( \int_0^T \|\mathbf{f}(t)\|_{\alpha^{-1}}^2 dt \right)^{\frac{1}{2}} + \|\mathbf{u}_0\|_{L^2(\Omega)^d} + c(\alpha) \|g_2\|_{\mathcal{C}^0(0,T; H_{00}^{1/2}(\Gamma_2)')} \nonumber
\end{align}$$

where the constant $c(\alpha)$ only depends on the coefficient $\alpha$.

**Proof:** Owing to the inf-sup condition (2.8) and the continuity of the trace operator from $H_{\triangle}^1(\Omega)$ into $H_{00}^{1/2}(\Gamma_2)$, for all $t$ in $[0, T]$, there exists (see [12, Chap. I, Lemma 4.1]) an isomorphism $B$ from $H_{00}^{1/2}(\Gamma_2)'$ into the orthogonal of $V(\Omega)$ (for the scalar product associated with the norm $\| \cdot \|_{\alpha}$) such that, when setting $\mathbf{w}(t) = Bg_2(t)$,

$$\forall q \in H_{\triangle}^1(\Omega), \quad (\mathbf{w}(t), \mathrm{grad} \, q) = \langle g_2(t), q \rangle_{\Gamma_2} \quad (2.13)$$

and also

$$\|\mathbf{w}(t)\|_{\alpha} \le c(\alpha) \|g_2(t)\|_{H_{00}^{1/2}(\Gamma_2)'} \quad (2.14)$$

Next, we observe that the function $\mathbf{u}^* = \mathbf{u} - \mathbf{w}$ belongs to $V(\Omega)$ and satisfies

$$\forall \mathbf{v} \in V(\Omega), \quad (\partial_t \mathbf{u}^*, \mathbf{v}) + (\alpha(\cdot, p) \mathbf{u}^*, \mathbf{v}) = (\mathbf{f}, \mathbf{v}) - (\partial_t \mathbf{w}, \mathbf{v}) - (\alpha(\cdot, p) \mathbf{w}, \mathbf{v}). \quad (2.15)$$

Taking $\mathbf{v}$ equal to $\mathbf{u}^*$ in this equation gives, with $c = \max_{1 \le i \le I} \frac{\alpha_{i\bar{i}}}{\alpha_{i\bar{b}}}$,

$$\begin{align*}
\frac{1}{2} d_t \|\mathbf{u}^*(t)\|_{L^2(\Omega)^d}^2 + \|\mathbf{u}^*(t)\|_{\alpha}^2 &\le (\|\mathbf{f}\|_{\alpha^{-1}} + \alpha_{\bar{b}}^{-\frac{1}{2}} \|\partial_t \mathbf{w}(t)\|_{L^2(\Omega)^d} + c \|\mathbf{w}(t)\|_{\alpha}) \|\mathbf{u}^*(t)\|_{\alpha} \\
&\le \frac{1}{2} (\|\mathbf{f}\|_{\alpha^{-1}} + \alpha_{\bar{b}}^{-\frac{1}{2}} \|\partial_t \mathbf{w}(t)\|_{L^2(\Omega)^d} + c \|\mathbf{w}(t)\|_{\alpha})^2 + \frac{1}{2} \|\mathbf{u}^*(t)\|_{\alpha}^2,
\end{align*}$$

whence, by integrating with respect to $t$,

$$\begin{align*}
\|\mathbf{u}^*(t)\|_{L^2(\Omega)^d}^2 - \|\mathbf{u}_0\|_{L^2(\Omega)^d}^2 + \int_0^t \|\mathbf{u}^*(s)\|_{\alpha}^2 ds \\
\le 3 \left( \int_0^t \|\mathbf{f}(s)\|_{\alpha^{-1}}^2 ds + \alpha_{\bar{b}}^{-1} \int_0^t \|\partial_t \mathbf{w}(s)\|_{L^2(\Omega)^d}^2 ds + c^2 \int_0^t \|\mathbf{w}(s)\|_{\alpha}^2 ds \right).
\end{align*}$$

When combined with (2.14), this yields the desired estimate.

We have not made precise the dependence of the constant $c(\alpha)$ in (2.12) with respect to the different $\alpha_i$ first for simplicity but also because, in practical situations, $\Gamma_2$ intersects a few number of subdomains $\bar{\Omega}_i$.

We are now in a position to derive the main result of this section, namely the existence of a solution to problem (2.4) – (2.5). However its rather technical proof requires some
preliminary results. Assuming that the data satisfy (2.11), we need the function **w** introduced in the proof of Lemma 2.3, which satisfies (2.13) and (2.14), and also the harmonic lifting $\tilde{p}_1$ of any extension of $p_1$ in $H^{1/2}(\partial\Omega)$, which satisfies for a.e. $t$ in $(0, T)$,

$$ \|\tilde{p}_1(t)\|_{H^1(\Omega)} \le c \|p_1(t)\|_{H^{1/2}(\Gamma_1)}. \quad (2.16) $$

Next, we set: $u^* = \mathbf{u} - \mathbf{w}$, $p^* = p - \tilde{p}_1$, and we observe that the pair $(\mathbf{u}^*, p^*)$ belongs to $\mathcal{C}^0(0, T; \mathbb{V}(\Omega)) \times L^2(0, T; H^1_\Delta(\Omega))$ and satisfies, for a.e. $t$ in $]0, T[$,

$$ \forall \mathbf{v} \in L^2(\Omega)^d, \quad (\partial_t \mathbf{u}^*, \mathbf{v}) + (\alpha^*(\cdot, p^*) \mathbf{u}^*, \mathbf{v}) + (\alpha^*(\cdot, p^*) \mathbf{w}, \mathbf{v}) + (\mathbf{v}, \mathrm{grad} \, p^*) = (\mathbf{f}^*, \mathbf{v}), \quad (2.17) $$

where obviously the new function $\mathbf{f}^*$ is equal to $\mathbf{f} - \partial_t \mathbf{w} - \mathrm{grad} \, \tilde{p}_1$ and the function $\alpha^*$ is defined by

$$ \forall (\mathbf{x}, t) \in \Omega \times ]0, T[, \quad \forall \xi \in \mathbb{R}, \quad \alpha^*(\mathbf{x}, t, \xi) = \alpha(\mathbf{x}, \xi + \tilde{p}_1(\mathbf{x}, t)). \quad (2.18) $$

It is readily checked that the function $\alpha^*$ satisfies the same properties (2.1) and (2.2) as the coefficient $\alpha$, with the same constants $\alpha_{ib}$, $\alpha_{i\ddag}$, and $\alpha_{i*}$. Finally, it follows from the definition (2.9) of $\mathbb{V}(\Omega)$ that $L^2(\Omega)^d$ is the orthogonal sum of $\mathbb{V}(\Omega)$ and $\mathrm{grad} \, H^1_\Delta(\Omega)$, so that problem (2.17) can be written as the sum of two coupled equations:

$$ \forall q \in H^1_\Delta(\Omega), \quad (\alpha^*(\cdot, p^*) \mathbf{u}^*, \mathrm{grad} \, q) + (\alpha^*(\cdot, p^*) \mathbf{w}, \mathrm{grad} \, q) + (\mathrm{grad} \, p^*, \mathrm{grad} \, q) = (\mathbf{f}^*, \mathrm{grad} \, q), \quad (2.19) $$

and

$$ \forall \mathbf{v} \in \mathbb{V}(\Omega), \quad (\partial_t \mathbf{u}^*, \mathbf{v}) + (\alpha^*(\cdot, p^*) \mathbf{u}^*, \mathbf{v}) + (\alpha^*(\cdot, p^*) \mathbf{w}, \mathbf{v}) = (\mathbf{f}^*, \mathbf{v}). \quad (2.20) $$

We begin with a preliminary lemma. We skip its proof since the result is a direct consequence of the Cauchy–Lipschitz theorem and the separability of $\mathbb{V}(\Omega)$ (see [5, Thm 2.4] for a similar result).

**Lemma 2.4.** For any data

$$ \mathbf{f} \in L^2(0, T; L^2(\Omega)^d), \quad \mathbf{u}_0 \in L^2(\Omega)^d, \quad p_1 \in L^2(0, T; H^{1/2}(\Gamma_1)), \\ \text{and} \quad g_2 \in H^1(0, T; H^{1/2}_{00}(\Gamma_2)') \quad (2.21) $$

and for any function $p^*$ in $L^2(0, T; L^2(\Omega))$, problem (2.20) provided with the initial condition $\mathbf{u}^*(\cdot, 0) = \mathbf{u}_0 - \mathbf{w}(\cdot, 0)$ has a unique solution $\mathbf{u}^*$ in $\mathcal{C}^0(0, T; \mathbb{V}(\Omega))$. Moreover this solution satisfies

$$ \sup_{0 \le t \le T} \|\mathbf{u}^*(t)\|_{L^2(\Omega)} + \left(\int_0^T \|\mathbf{u}^*(t)\|_{\alpha}^2 dt\right)^{\frac{1}{2}} \\ \le 3 \left(\int_0^T \|\mathbf{f}(t)\|_{\alpha^{-1}}^2 dt\right)^{\frac{1}{2}} + \|\mathbf{u}_0\|_{L^2(\Omega)^d} \\ + c(\alpha) \|p_1\|_{L^2(0,T; H^{1/2}(\Gamma_1))} + c(\alpha) \|g_2\|_{H^1(0,T; H^{1/2}_{00}(\Gamma_2))}. \quad (2.22) $$
We call $\mathcal{F}$ the mapping: $p^* \mapsto u^*$, where $u^*$ is the solution exhibited in Lemma 2.4. It is readily checked from this lemma that $\mathcal{F}$ is continuous from $L^2(0, T; L^2(\Omega))$ into $\mathcal{C}^0(0, T; L^2(\Omega)^d)$ and also, thanks to (2.22), that its range is bounded.

**Theorem 2.5.** Assume that the partition of $\partial\Omega$ into $\Gamma_1$ and $\Gamma_2$ is sufficiently smooth for $\mathcal{D}(\Omega \cup \Gamma_2)$ to be dense in $H^1_\Delta(\Omega)$. For any data $(f, u_0, p_1, g_2)$ satisfying (2.21), problem (2.4)-(2.5) has a solution $(u, p)$ in $\mathcal{C}^0(0, T; L^2(\Omega)^d) \times L^2(0, T; H^1(\Omega))$. Moreover, its part $u$ satisfies (2.12) and its part $p$ satisfies

$$ \left( \int_0^T \|p\|_{H^1(\Omega)}^2 dt \right)^{\frac{1}{2}} \le c \alpha_{\#} \left( \|f\|_{L^2(0,T;L^2(\Omega)^d)} + \|u_0\|_{L^2(\Omega)^d} \right. \\ \left. \qquad + c(\alpha) \left( \|p_1\|_{L^2(0,T;H^{1/2}(\Gamma_1))} + \|g_2\|_{\mathcal{C}^0(0,T;H_{00}^{1/2}(\Gamma_2)')} \right) \right). \quad (2.23) $$

**Proof:** The main idea is to apply Brouwer's fixed point theorem to problem (2.19). Several steps are needed for that.

1) For a.e. $t$ in $(0, T)$, we first define a mapping $\Phi$ from $H^1_\Delta(\Omega)$ into its dual space by

$$ \forall p \in H^1_\Delta(\Omega), \forall q \in H^1_\Delta(\Omega), \\ \langle \Phi(p), q \rangle = (\alpha^*(\cdot, p) \mathcal{F}(p), \mathbf{grad} \, q) + (\alpha^*(\cdot, p) \mathbf{w}, \mathbf{grad} \, q) + (\mathbf{grad} \, p, \mathbf{grad} \, q) \\ - (\mathbf{f}^*, \mathbf{grad} \, q), $$

where $\mathbf{w}$ is still the lifting introduced in the proof of Lemma 2.3. This mapping is continuous and moreover satisfies

$$ \langle \Phi(p), p \rangle \ge |p|_{H^1(\Omega)}^2 - c_0(t) |p|_{H^1(\Omega)}, $$

where $c_0(t)$ is given by

$$ c_0(t) = 2c_1 \alpha_{\#} + \|\mathbf{f}^*(t)\|_{L^2(\Omega)^d}, $$

and $c_1$ stands for the right-hand side of (2.22). Thus, $\langle \Phi(p), p \rangle$ is nonnegative on the sphere with radius $c_0(t)$.

2) It follows from the density assumption that there exists an increasing sequence $(\mathbb{H}_n)_n$ of finite-dimensional subspaces of $H^1_\Delta(\Omega)$ such that $\cup_{n \in \mathbb{N}} \mathbb{H}_n = H^1_\Delta(\Omega)$. For each $n$, the restriction of the mapping $\Phi$ to $\mathbb{H}_n$ satisfies exactly the same properties as previously, so that applying Brouwer's fixed point theorem (see [12, Chap. IV, Cor. 1.1] for instance) yields for each $n$ the existence of a function $p_n$ in $\mathbb{H}_n$ so that

$$ \forall q_n \in \mathbb{H}_n, \quad \langle \Phi(p_n), q_n \rangle = 0 \quad \text{and} \quad |p_n|_{H^1(\Omega)} \le c_0(t). $$

3) Owing to the Poincaré-Friedrichs inequality, the seminorm $|\cdot|_{H^1(\Omega)}$ is a norm on $H^1_\Delta(\Omega)$. Since the sequence $(p_n)_n$ is bounded in this norm, there exists a subsequence, still denoted by $(p_n)_n$ for simplicity, which converges to a function $p^*$ weakly in $H^1(\Omega)$ and strongly in $L^2(\Omega)$ (we use here the compactness of the imbedding of $H^1(\Omega)$ into $L^2(\Omega)$). It thus follows from the weak lower semi-continuity of the norm that

$$ |p^*(t)|_{H^1(\Omega)} \le c_0(t). $$
Then, using a Poincaré–Friedrichs inequality, combining this with (2.16) and integrating the square of the result with respect to t gives estimate (2.23).

4) The function $p_n$ satisfies for all $q_m$ in $\mathbb{H}_m$, $m \le n$,

$$ (\alpha^*(\cdot, p_n) \mathcal{F}(p_n), \mathbf{grad} \, q_m) + (\alpha^*(\cdot, p_n) \mathbf{w}, \mathbf{grad} \, q_m) + (\mathbf{grad} \, p_n, \mathbf{grad} \, q_m) = (\mathbf{f}^*, \mathbf{grad} \, q_m). $$

Thanks the continuity of the mapping $\mathcal{F}$ on $L^2(0, T; L^2(\Omega))$, the sequence $(\mathcal{F}(p_n))_n$ converges to $\mathcal{F}(p^*)$ strongly in $L^2(0, T; L^2(\Omega))$. We now set: $\mathbf{u}^* = \mathcal{F}(p^*)$. On the other hand, since $(p_n(\mathbf{x}, t))$ converges to $p^*$ a.e. in $\Omega \times ]0, T[$, the sequence $(\alpha^*(\cdot, p_n) \mathbf{grad} \, q_m)_n$ converges to $\alpha^*(\cdot, p^*) \mathbf{grad} \, q_m$ a.e. in $\Omega \times ]0, T[$ and, for a.e. $t$, is bounded by $\alpha_{\#} \mathbf{grad} \, q_m$ which is square-integrable; thus applying the Lebesgue dominated convergence theorem implies that $(\alpha^*(\cdot, p_n) \mathbf{grad} \, q_m)_n$ converges to $\alpha^*(\cdot, p^*) \mathbf{grad} \, q_m$ in $L^2(\Omega)$. By combining these convergence properties, we observe that $p^*$ satisfies for all $q_m$ in $\mathbb{H}_m$,

$$ (\alpha^*(\cdot, p^*) \mathbf{u}^*, \mathbf{grad} \, q_m) + (\alpha^*(\cdot, p^*) \mathbf{w}, \mathbf{grad} \, q_m) + (\mathbf{grad} \, p^*, \mathbf{grad} \, q_m) = (\mathbf{f}^*, \mathbf{grad} \, q_m). $$

Finally, the density of $\cup_{m \in \mathbb{N}} \mathbb{H}_m$ in $H^1_{\Delta}(\Omega)$ implies that $p^*$ is a solution of problem (2.19). To conclude, the pair $(\mathbf{u} = \mathbf{u}^* + \mathbf{w}, p = p^* + \hat{p}_1)$ is a solution of problem (2.4) – (2.5) and satisfies (2.12) and (2.23).

**Remark 2.6.** Some further regularity properties can be derived on the solution $(\mathbf{u}, p)$ if the data are smooth enough, for instance:

(i) By formally taking $\mathbf{v}$ equal to $\partial_t \mathbf{u}^*$ in (2.20) and with a little more regularity on the data, we obtain that the part $\mathbf{u}$ belongs to $H^1(0, T; L^2(\Omega)^d)$.

(ii) By using an extension of [5, Prop. 2.5] to the case of mixed boundary conditions, it can be checked that, if the data are smooth enough and when the coefficient $\alpha$ is constant, the solution $(\mathbf{u}, p)$ belongs to $H^1(0, T; H^s(\Omega)^d) \times L^2(0, T; H^{s+1}(\Omega))$ for a real number $s > 0$ only depending on the geometry of $\Omega$. However, in the case of pressure dependent permeability $\alpha$, this property seems much more difficult to establish.

We conclude with a uniqueness result.

**Proposition 2.7.** If problem (2.4) – (2.5) admits a solution $(\mathbf{u}, p)$ such that its part $\mathbf{u}$ belongs to $L^\rho(\Omega)^d$ with $\rho > 2$ in dimension $d = 2$ and $\rho \ge 3$ in dimension $d = 3$ and satisfies for an appropriate constant $\kappa$

$$ \left( \max_{1 \le i \le I} \alpha_{i*} \right) \sup_{0 \le t \le T} \| \mathbf{u}(t) \|_{L^\rho(\Omega)^d} \le \kappa, \qquad (2.24) $$

then this solution is unique.

**Proof:** Let $(\mathbf{u}_1, p_1)$ and $(\mathbf{u}_2, p_2)$ be two solutions of problem (2.4)–(2.5), with $\mathbf{u}_1$ satisfying (2.24). It is readily checked that the pair $(\mathbf{u}^\dagger = \mathbf{u}_1 - \mathbf{u}_2, p^\dagger = p_1 - p_2)$ belongs to $\mathcal{C}^0(0, T; \mathbb{V}(\Omega)) \times L^2(0, T; H^1_{\Delta}(\Omega))$ and satisfies (note that this requires the definition of $\mathbb{V}(\Omega)$)

$$ \forall q \in H^1_{\Delta}(\Omega), \quad (\alpha(p_2) \mathbf{u}^\dagger, \mathbf{grad} \, q) + (\mathbf{grad} \, p^\dagger, \mathbf{grad} \, q) = ((\alpha(p_2) - \alpha(p_1))\mathbf{u}_1, \mathbf{grad} \, q), \quad (2.25) $$

$$ \forall \mathbf{v} \in \mathbb{V}(\Omega), \quad (\partial_t \mathbf{u}^\dagger, \mathbf{v}) + (\alpha(p_2) \mathbf{u}^\dagger, \mathbf{v}) = ((\alpha(p_2) - \alpha(p_1))\mathbf{u}_1, \mathbf{v}). $$
1) Taking $q$ equal to $p^\dagger$ in the first equation (2.25) yields

$$|p^\dagger|_{H^1(\Omega)} \le \alpha_\# \|\mathbf{u}^\dagger\|_{L^2(\Omega)^d} + \left(\max_{1 \le i \le I} \alpha_{i*}\right) \|p^\dagger\|_{L^{\rho^*}(\Omega)} \|\mathbf{u}_1\|_{L^\rho(\Omega)^d},$$

with $\frac{1}{\rho} + \frac{1}{\rho^*} = \frac{1}{2}$. Using the imbedding of $H^1(\Omega)$ into $L^{\rho^*}(\Omega)$ and the Poincaré–Friedrichs inequality thus yields

$$|p^\dagger|_{H^1(\Omega)} \le \alpha_\# \|\mathbf{u}^\dagger\|_{L^2(\Omega)^d} + c \left( \max_{1 \le i \le I} \alpha_{i*} \right) |p^\dagger|_{H^1(\Omega)} \|\mathbf{u}_1\|_{L^\rho(\Omega)^d}.$$

Thus, choosing the constant $\kappa$ in (2.24) smaller than $1/2c$ yields

$$|p^\dagger|_{H^1(\Omega)} \le 2\alpha_\# \|\mathbf{u}^\dagger\|_{L^2(\Omega)^d}. \qquad (2.26)$$

2) On the other hand, by taking $\mathbf{v}$ equal to $\mathbf{u}^\dagger$ in the second equation (2.25), we obtain

$$\frac{1}{2} d_t \|\mathbf{u}^\dagger\|_{L^2(\Omega)^d}^2 \le \left( \max_{1 \le i \le I} \alpha_{i*} \right) \|p^\dagger\|_{L^{\rho^*}(\Omega)} \|\mathbf{u}_1\|_{L^\rho(\Omega)^d} \|\mathbf{u}^\dagger\|_{L^2(\Omega)^d},$$

whence, by the same arguments as in the first part of the proof,

$$d_t \|\mathbf{u}^\dagger\|_{L^2(\Omega)^d} \le c \left( \max_{1 \le i \le I} \alpha_{i*} \right) |p^\dagger|_{H^1(\Omega)} \|\mathbf{u}_1\|_{L^\rho(\Omega)^d}.$$

Combining this with (2.26) yields

$$d_t \|\mathbf{u}^\dagger\|_{L^2(\Omega)^d} \le \alpha_\# \|\mathbf{u}^\dagger\|_{L^2(\Omega)^d}.$$

Since $\mathbf{u}^\dagger(\cdot, 0)$ is zero, integrating this inequality between 0 and $t$ and applying Grönwall's lemma (see [11, Chap. V, Lemma 1.8] for instance) imply that $\|\mathbf{u}^\dagger\|_{L^2(\Omega)^d}$ is zero. Thus, $\mathbf{u}_1$ and $\mathbf{u}_2$ coincide.

3) Finally, combining the previous result with (2.26) gives that $p_1$ and $p_2$ are equal.

Even if condition (2.24) is not too restrictive, we try to avoid it in what follows.
### 3. The semi-discrete and discrete problems.

We first describe the time semi-discrete problem constructed from the Euler's scheme. Next, we consider the fully discrete problem obtained by combining this scheme with the mortar finite element method used for the space discretization.

#### The time semi-discrete problem

For already explained reasons, we have decided to work with a fixed time step, that we denote by $\delta t$, and with an explicit treatment of the nonlinear term $\alpha(\cdot, p)$. For this, we need an initialization step.

**INITIALIZATION STEP:** For each $i$, $1 \le i \le I$, we fix a constant $\bar{\alpha}_i$ satisfying $\alpha_{ib} \le \bar{\alpha}_i \le \alpha_{i\ddot{}}$ and we denote by $\bar{\alpha}$ the function equal to $\bar{\alpha}_i$ on each $\Omega_i$. Next, assuming that the data $f$, $p_1$ and $g_2$ belong to $\mathcal{C}^0(0, T; L^2(\Omega))$, $\mathcal{C}^0(0, T; H^{\frac{1}{2}}(\Gamma_1))$ and $\mathcal{C}^0(0, T; H_{00}^{\frac{1}{2}}(\Gamma_2))'$, respectively, we consider the linear steady Darcy problem

$$ \text{Find } (\tilde{u}^0, p^0) \text{ in } L^2(\Omega)^d \times H^1(\Omega) \text{ such that} \\ p^0 = p_1(\cdot, 0) \quad \text{on } \Gamma_1, \qquad (3.1) $$

and

$$ \begin{gathered} \forall v \in L^2(\Omega)^d, (\bar{\alpha}(\cdot) \tilde{u}^0, v) + (v, \text{grad } p^0) = (f(\cdot, 0), v), \\ \forall q \in H_{\Delta}^1(\Omega), \quad (\tilde{u}^0, \text{grad } q) = \langle g_2(\cdot, 0), q \rangle_{\Gamma_2}. \end{gathered} \qquad (3.2) $$

Standard arguments (see the inf-sup condition (2.8)) yield that problem (3.1) – (3.2) has a unique solution. Moreover, when applying an Uzawa-like algorithm, we observe that this problem is equivalent to the Laplace equation with piecewise constant diffusion coefficient and mixed boundary conditions

Find $p^0$ in $H^1(\Omega)$ satisfying (3.1) and such that

$$ \forall q \in H_{\Delta}^1(\Omega), \quad \left( \frac{1}{\bar{\alpha}(\cdot)} \text{grad } p^0, \text{grad } q \right) = \left( \frac{1}{\bar{\alpha}(\cdot)} f(\cdot, 0), \text{grad } q \right) - \langle g_2(\cdot, 0), q \rangle_{\Gamma_2}. \qquad (3.3) $$

Solving this problem is not expensive. Moreover, for what follows, we only need $p^0$. Indeed, we initialize the sequence $(u^n, p^n)$ by $(u_0, p^0)$.

**ITERATION STEP:** At each step $n$, $1 \le n \le N$ (with $T = N \delta t$), assuming that $(u^{n-1}, p^{n-1})$ is known, we consider the following problem

$$ \text{Find } (u^n, p^n) \text{ in } L^2(\Omega)^d \times H^1(\Omega) \text{ such that} \\ p^n = p_1(\cdot, n\delta t) \quad \text{on } \Gamma_1, \qquad (3.4) $$

and

$$ \begin{gathered} \forall v \in L^2(\Omega)^d, \quad \left( \frac{u^n - u^{n-1}}{\delta t}, v \right) + (\alpha(\cdot, p^{n-1}) u^n, v) + (v, \text{grad } p^n) \\ = (f(\cdot, n\delta t), v), \end{gathered} \qquad (3.5) $$

$$ \forall q \in H_{\Delta}^1(\Omega), \quad (u^n, \text{grad } q) = \langle g_2(\cdot, n\delta t), q \rangle_{\Gamma_2}. $$
It can be noted that, since $(\mathbf{u}^{n-1}, p^{n-1})$ is known, this problem can be written as a standard steady Darcy equation. Thus, the next statement is easily derived by using the inf-sup condition (2.8).

**Proposition 3.1.** For any data

$$ \mathbf{f} \in \mathcal{C}^0(0, T; L^2(\Omega)^d), \quad \mathbf{u}_0 \in L^2(\Omega)^d, \quad p_1 \in \mathcal{C}^0(0, T; H^{1/2}(\Gamma_1)), \\ \text{and} \quad g_2 \in \mathcal{C}^0(0, T; H_{00}^{1/2}(\Gamma_2)') \tag{3.6} $$

problem (3.4) – (3.5) for $1 \le n \le N$ has a unique solution $(\mathbf{u}^n, p^n)$ in $L^2(\Omega)^d \times H^1(\Omega)$.

## The fully discrete problem

From now on, we assume that each $\Omega_i$ is a polygon ($d = 2$) or a polyhedron with a Lipschitz-continuous boundary ($d = 3$). For each $i$, $1 \le i \le I$, let $(\mathcal{T}_h^i)_{h_i}$ be a regular family of triangulations of $\Omega_i$ (by triangles or tetrahedra), in the usual sense that:

* For each $h_i$, $\bar{\Omega}_i$ is the union of all elements of $\mathcal{T}_h^i$;

* The intersection of two different elements of $\mathcal{T}_h^i$, if not empty, is a vertex or a whole edge or a whole face of both of them;

* The ratio of the diameter $h_K$ of any element $K$ of $\mathcal{T}_h^i$ to the diameter of its inscribed circle or sphere is smaller than a constant $\sigma$ independent of $h_i$.

As usual, $h_i$ stands for the maximum of the diameters $h_K$, $K \in \mathcal{T}_h^i$. We also introduce the triangulation $\mathcal{T}_h$ equal to the union of the $\mathcal{T}_h^i$ and we denote by $\mathbf{h}$ the $I$-tuple $(h_1, h_2, \dots, h_I)$. In what follows, $c, c', \dots$ stand for generic constants which may vary from line to line but are always independent of $\mathbf{h}$. We make the further (and non restrictive) assumption that both $\bar{\Gamma}_1$ and $\bar{\Gamma}_2$ are the unions of whole edges ($d = 2$) or whole faces ($d = 3$) of elements of $\mathcal{T}_h$.

We first introduce the local discrete spaces. For each $k \ge 0$, $\mathcal{P}_k(K)$ stands for the space of restrictions to $K$ of polynomials with $d$ variables and total degree $\le k$. Thus, the discrete space of velocities on each $\Omega_i$ is defined as

$$ \mathbb{X}_h^i = \{ \mathbf{v}_h \in L^2(\Omega_i)^d; \forall K \in \mathcal{T}_h^i, \mathbf{v}_h|_K \in \mathcal{P}_0(K)^d \}. \tag{3.7} $$

On the other hand, the discrete space of pressures $\mathbb{M}_h^i$ is the space of functions $q_h$ in $L^2(\Omega_i)$ such that their restrictions to each $K$ in $\mathcal{T}_h^i$ belongs to $\mathcal{P}_1(K)$ and which are continuous at the midpoint of each edge ($d = 2$) or at the barycenter of each face ($d = 3$) of all elements of $\mathcal{T}_h^i$ (see [9] for more details on the corresponding finite element and [2] for its first use for Darcy's equations).

The global discrete space of velocities is then constructed in an obvious way:

$$ \mathbb{X}_h = \{ \mathbf{v}_h \in L^2(\Omega)^d; \mathbf{v}_h|_{\Omega_i} \in \mathbb{X}_h^i, 1 \le i \le I. \} \tag{3.8} $$

However, defining the global discrete space of pressures requires some further notation linked to the mortar method.
Let $\mathcal{S}$ stand for the skeleton of the decomposition

$$ \mathcal{S} = \bigcup_{i=1}^{I} \partial \Omega_i \setminus \partial \Omega. \qquad (3.9) $$

As usual for the mortar element method, see [7, Section 2], we introduce a further decomposition of $\mathcal{S}$ into disjoint open parts

$$ \bar{\mathcal{S}} = \bigcup_{m=1}^{M} \bar{\gamma}_{m}^{+} \quad \text{and} \quad \gamma_{m}^{+} \cap \gamma_{m'}^{+} = \emptyset, \quad 1 \le m < m' \le M, $$

where each $\gamma_m^+$ is

* a part of an edge of an $\Omega_i$, denoted by $\Omega_m^+$, in dimension $d=2$,

* a part of a face of an $\Omega_i$, denoted by $\Omega_m^+$, in dimension $d=3$.

The choice of the $\gamma_m^+$ is not unique but is made independently of the discretization. We also define $\mathcal{E}_{mh}^+$ as the set of edges ($d=2$) or faces ($d=3$) $\gamma_m^+ \cap \partial K$, for all $K$ in $\mathcal{T}_h^i$ for the $i$ such that $\Omega_i = \Omega_m^+$. In order to enforce the matching conditions through $\mathcal{S}$, we need the space

$$ \mathbb{W}_{mh}^{+} = \{ \varphi_h \in L^2(\gamma_m^{+}); \forall e \in \mathcal{E}_{mh}^{+}, \varphi_{h|e} \in \mathcal{P}_0(e) \}, \qquad (3.10) $$

see [6, §3] for this choice. We associate with each piecewise regular function $q$ its mortar function $\Phi_m(q)$: On each $\gamma_m^+$ and for each $\Omega_i \neq \Omega_m^+$, the restriction of $\Phi_m(q)$ to $\gamma_m^+ \cap \partial \Omega_i$ is equal to the trace of $q|_{\Omega_i}$. The discrete space of pressures is the space $\mathbb{M}_h$ of functions $q_h$ in $L^2(\Omega)$ such that

(i) their restriction to each $\Omega_i$, $1 \le i \le I$, belongs to $\mathbb{M}_h^i$,

(ii) the following matching condition holds on each $\gamma_m^+$, $1 \le m \le M$,

$$ \forall \varphi_h \in \mathbb{W}_{mh}^{+}, \quad \int_{\gamma_m^{+}} (q_{h|\Omega_m^{+}} - \Phi_m(q_h))(\boldsymbol{\tau})\varphi_h(\boldsymbol{\tau}) \, d\boldsymbol{\tau} = 0, \qquad (3.11) $$

where $\boldsymbol{\tau}$ denotes the tangential coordinate(s) on $\gamma_m^+$.

The space $\mathbb{M}_h^{\Delta}$ is the space of functions in $\mathbb{M}_h$ which vanish on $\Gamma_1$.

We also introduce the discrete products, for all functions $v$ and $w$ in $L^2(\Omega)$,

$$ (v, w)_h^i = \sum_{K \in \mathcal{T}_h^i} \int_K v(\boldsymbol{x})w(\boldsymbol{x}) \, d\boldsymbol{x}, \quad (v, w)_h = \sum_{i=1}^{I} (v, w)_h^i, $$

with obvious extension to vector-valued functions. We are thus in a position to define the discrete problems.

**INITIALIZATION STEP:** Let $\mathcal{E}_h^{ik}$, $k=1,2$, be the set of edges ($d=2$) or faces ($d=3$) of elements of $\mathcal{T}_h^i$ which are contained in $\Gamma_k$. We thus introduce the approximation $p_{1h}^0$ of $p_1(\cdot, 0)$ defined by

$$ \forall e \in \mathcal{E}_h^{i1}, \quad p_{1h}^0|_e = \frac{1}{\text{meas}(e)} \int_e p_1(\boldsymbol{\tau}, 0) \, d\boldsymbol{\tau}, \quad 1 \le i \le I. \qquad (3.12) $$
We also denote by $\mathcal{A}_1$ the set of midpoints of each edge ($d=2$) or barycenters of each face ($d=3$) of elements of $\mathcal{T}_h$ that belong to $\Gamma_1$. In view of problem (3.1) – (3.3), assuming from now on that the function $g_2$ belongs to $\mathcal{C}^0(0, T; L^2(\Gamma_2))$, we consider the discrete problem:

Find $p_h^0$ in $M_h$ such that

$$ \forall \mathbf{a} \in \mathcal{A}_1, \quad p_h^0(\mathbf{a}) = p_{1h}^0(\mathbf{a}), \qquad (3.13) $$

and

$$ \forall q_h \in \mathbb{M}_h^{\Delta}, \quad \left( \frac{1}{\bar{\alpha}(\cdot)} \mathbf{grad} \, p_h^0, \mathbf{grad} \, q_h \right)_h \\ = \left( \frac{1}{\bar{\alpha}(\cdot)} \mathbf{f}(\cdot, 0), \mathbf{grad} \, q_h \right)_h - \sum_{i=1}^{I} \sum_{e \in \mathcal{E}_h^{i2}} \int_e g_2(\tau, 0) q_h(\tau) \, d\tau. \qquad (3.14) $$

A further (non restrictive) assumption is needed for proving the well-posedness of this problem.

**Assumption 3.2.** For each $i$, $1 \le i \le I$, and each $m$, $1 \le m \le M$, such that $\partial\Omega_i \cap \gamma_m^+$ has a positive measure in $S$, $\partial\Omega_i$ contains an element of $\mathcal{E}_{mh}^+$.

**Proposition 3.3.** If Assumption 3.2 holds, for any data

$$ \mathbf{f} \in \mathcal{C}^0(0, T; L^2(\Omega)^d), \quad \mathbf{u}_0 \in L^2(\Omega)^d, \quad p_1 \in \mathcal{C}^0(0, T; H^{1/2}(\Gamma_1)), \\ \text{and} \quad g_2 \in \mathcal{C}^0(0, T; L^2(\Gamma_2)), \qquad (3.15) $$

problem (3.13) – (3.14) has a unique solution $p_h^0$ in $M_h$.

**Proof:** Since problem (3.13) – (3.14) results into a square linear system, it suffices to check that its solution for zero data is zero. When $p_{1h}^0$, $\mathbf{f}(\cdot, 0)$ and $g_2(\cdot, 0)$ are equal to zero, taking $q_h$ equal to $p_h^0$ (which now belongs to $\mathbb{M}_h^\Delta$) yields that $\mathbf{grad} \, p_h^0$ is zero on each $K$ in $\mathcal{T}_h$, so that $p_h^0$ is constant on each $K$. Thus, the continuity conditions which appear in the definition of each $\mathbb{M}_h^i$ imply that $p_h^0$ is constant on each $\Omega_i$. It follows from the definition (3.10) of the $W_{mh}^+$ that the matching conditions (3.11) are local on each $e$ in $\mathcal{E}_{mh}^+$ so that, due to Assumption 3.2, $p_h^0$ is constant on $\Omega$. Thus the boundary conditions (3.13) imply that it is zero, whence the desired result.

Finally, we introduce the orthogonal projection operator $\Pi_h$ from $L^2(\Omega)^d$ onto $X_h$ and we initialize the sequence $(\mathbf{u}_h^n, p_h^n)$ by $(\Pi_h \mathbf{u}_0, p_h^0)$.

**ITERATION STEP:** At each step $n$, $1 \le n \le N$, we introduce the approximation $p_{1h}^n$ of $p_1(\cdot, n\delta t)$ defined by

$$ \forall e \in \mathcal{E}_h^{i1}, \quad p_{1h}^n|_e = \frac{1}{\text{meas}(e)} \int_e p_1(\tau, n\delta t) \, d\tau, \quad 1 \le i \le I. \qquad (3.16) $$

Next, assuming that $(\mathbf{u}_h^{n-1}, p_h^{n-1})$ is known, we consider the following problem
Find $(\mathbf{u}_h^n, p_h^n)$ in $X_h \times \mathbb{M}_h$ such that

$$ \forall \mathbf{a} \in \mathcal{A}_1, \quad p_h^n(\mathbf{a}) = p_{1h}^n(\mathbf{a}), \qquad (3.17) $$

and

$$ \forall \mathbf{v}_h \in X_h, \quad \left( \frac{\mathbf{u}_h^n - \mathbf{u}_h^{n-1}}{\delta t}, \mathbf{v}_h \right) + \left( \alpha(\cdot, p_h^{n-1}) \mathbf{u}_h^n, \mathbf{v}_h \right) + \left( \mathbf{v}_h, \mathbf{grad} \, p_h^n \right)_h \\ = (\mathbf{f}(\cdot, n\delta t), \mathbf{v}_h), \qquad (3.18) $$

$$ \forall q_h \in \mathbb{M}_h^{\Delta}, \quad (\mathbf{u}_h^n, \mathbf{grad} \, q_h)_h = \sum_{i=1}^{I} \sum_{e \in \mathcal{E}_h^{i2}} \int_e g_2(\tau, n\delta t) q_h(\tau) \, d\tau. $$

**Remark 3.4.** In the implementation of this problem, the term $\alpha(\cdot, p_h^{n-1})$ is most often replaced by its Lagrange interpolate $\alpha_h(\cdot, p_h^{n-1})$: On each $K$ in $\mathcal{T}_h$, $\alpha_h(\cdot, p_h^{n-1})$ belongs to $\mathcal{P}_1(K)$ and is equal to $\alpha(\mathbf{a}, p_h^{n-1}(\mathbf{a}))$ at the $d+1$ vertices $\mathbf{a}$ of $K$. Then, all integrals in the previous problem are easily computed by the trapeze formula. We do not take this modification into account for simplicity.

There also, this problem is well-posed, as stated in the next proposition.

**Proposition 3.5.** If Assumption 3.2 holds, for any data $\mathbf{f}$, $\mathbf{u}_0$, $p_1$, and $g_2$ satisfying (3.15), problem (3.17) – (3.18) has a unique solution $(\mathbf{u}_h^n, p_h^n)$ in $X_h \times \mathbb{M}_h$.

**Proof:** There also, problem (3.17) – (3.18) results into a square linear system, so that it suffices to check that the only solution of the next problem: Find $(\mathbf{u}_h^n, p_h^n)$ in $X_h \times \mathbb{M}_h^\Delta$ such that

$$ \forall \mathbf{v}_h \in X_h, \quad \left( \frac{\mathbf{u}_h^n}{\delta t}, \mathbf{v}_h \right) + \left( \alpha(\cdot, p_h^{n-1}) \mathbf{u}_h^n, \mathbf{v}_h \right) + \left( \mathbf{v}_h, \mathbf{grad} \, p_h^n \right)_h = 0, \qquad (3.19) $$

$$ \forall q_h \in \mathbb{M}_h^{\Delta}, \quad (\mathbf{u}_h^n, \mathbf{grad} \, q_h)_h = 0, $$

is zero. First, taking $\mathbf{v}_h$ equal to $\mathbf{u}_h^n$ and using (2.1) yields that $\mathbf{u}_h^n$ is zero. Second, defining $\mathbf{v}_h$ as equal to the gradient of $p_h^n$ on each $K$ in $\mathcal{T}_h$ and using the same arguments as in the proof of Proposition 3.3 yields that $p_h^n$ is zero. This concludes the proof.

It must be noted that Assumption 3.2 is not at all restrictive: It always holds when the $h_i$ are small enough. Thus, we have exhibited discrete problems which are well-posed when the data satisfy very weak regularity conditions.
## 4. Error estimates.

As now standard for elliptic or parabolic problems, the a priori analysis of the discretization relies on the approach proposed by F. Brezzi, J. Rappaz and P.-A. Raviart [8]. We first write another formulation of the continuous problem.

Let $\bar{\alpha}$ the function introduced in Section 3, i.e. the function which is piecewise constant, equal to $\bar{\alpha}_i$ on each $\Omega_i$. We denote by $\mathcal{T}$ the operator which associates with any data $(f, u_0, p_1, g_2)$ satisfying (2.21) the solution $(u, p)$ of the problem: Find $(u, p)$ in $\mathcal{C}^0(0, T; L^2(\Omega)^d) \times L^2(0, T; H^1(\Omega))$ satisfying (2.4) and such that, for a.e. $t$ in $]0, T[$,

$$ 
\begin{aligned}
\forall v \in L^2(\Omega)^d, \quad (\partial_t u, v) + (\bar{\alpha}(\cdot) u, v) + (v, \mathbf{grad} \, p) &= (f, v), \\
\forall q \in H^1_{\Delta}(\Omega), \quad (u, \mathbf{grad} \, q) &= \langle g_2, q \rangle_{\Gamma_2}.
\end{aligned} 
\quad (4.1) $$

It is readily checked that this problem has a unique solution (see [5, Thm 2.1] for the case of boundary conditions on the pressure only). Moreover, the operator $\mathcal{F}$ is continuous from the space of data introduced in (2.21) with values in $\mathcal{C}^0(0, T; L^2(\Omega)^d) \times L^2(0, T; H^1(\Omega))$. With the notation $U = (u, p)$, problem (2.4) – (2.5) can equivalently be written

$$ \mathcal{F}(U) = U - \mathcal{T}\mathcal{G}(U) = 0, \quad \text{with} \quad \mathcal{G}(U) = (f + (\bar{\alpha}(\cdot) - \alpha(\cdot, p))u, u_0, p_1, g_2). \quad (4.2) $$

Similarly, with any data $(f, u_0, p_1, g_2)$ satisfying (3.15), we associate the sequence $(u_h^n, p_h^n)_{1 \le n \le N}$ defined as follows:

(i) $u_h^0$ is equal to $\Pi_h u_0$;

(ii) for $1 \le n \le N$, $(u_h^n, p_h^n)$ is the solution of the problem: Find $(u_h^n, p_h^n)$ in $X_h \times M_h$ satisfying (3.17) and such that

$$ 
\begin{aligned}
\forall v_h \in X_h, \quad & \left( \frac{u_h^n - u_h^{n-1}}{\delta t}, v_h \right) + (\bar{\alpha}(\cdot) u_h^n, v_h) + (v_h, \mathbf{grad} \, p_h^n)_h = (f(\cdot, n\delta t), v_h), \\
\forall q_h \in M_h^{\Delta}, \quad & (u_h^n, \mathbf{grad} \, q_h)_h = \sum_{i=1}^{I} \sum_{e \in \mathcal{E}_h^{i2}} \int_{e} g_2(\tau, n\delta t) q_h(\tau) \, d\tau,
\end{aligned} 
\quad (4.3) $$

where $p_h^n$ is defined by (3.16).

With each sequence $(v^n)_{0 \le n \le N}$, we associate the function $i_{\delta t}v$ which is affine on each interval $[(n-1)\delta t, n\delta t]$, $1 \le n \le N$, and equal to $v^n$ in $n\delta t$, $0 \le n \le N$ ($i_{\delta t}$ stands for the time Lagrange interpolation operator and can also be applied to continuous functions on $[0, T]$). With each sequence $(q^n)_{0 \le n \le N}$, we associate the functions $\pi_{\delta t}^+q$ and $\pi_{\delta t}^-q$, which are constant, equal to $q^n$ and $q^{n-1}$, respectively, on each interval $[(n-1)\delta t, n\delta t]$, $1 \le n \le N$. Thus, denoting by $\delta$ the discretization parameter $(\delta t, h)$, we introduce the operator $\mathcal{T}_\delta$ which associates with the data $(f, u_0, p_1, g_2)$ the pair $U_\delta = (i_{\delta t}u_h, \pi_{\delta t}^+p_h)$. Thus, the discrete problems (3.13) – (3.14) and (3.17) – (3.18) can equivalently be written

$$ 
\begin{aligned}
\mathcal{F}_{\delta}(U_{\delta}) = U_{\delta} - \mathcal{T}_{\delta}\mathcal{G}_{\delta}(U_{\delta}) = 0, \\
\text{with} \quad \mathcal{G}_{\delta}(U) = (f + (\bar{\alpha}(\cdot) - \alpha(\cdot, \pi_{\delta t}^{-}p))i_{\delta t}u, u_0, p_1, g_2).
\end{aligned} 
\quad (4.4) $$
## Some properties of the operators $\mathcal{T}$ and $\mathcal{T}_\delta$

To go further, we make an assumption that is slightly stronger than Assumption 3.2 but still non-restrictive.

**Assumption 4.1.** For each $m$, $1 \le m \le M$, let $\mathcal{E}_{mh}^-$ be the set of the edges ($d=2$) or faces ($d=3$) of all triangulations $\mathcal{T}_h^i$ such that $\Omega_i$ does not coincide with $\Omega_m^+$ and that the intersection $e \cap \gamma_m^+$ has a positive measure. There exists a mapping $\Psi_m$ from $\mathcal{E}_{mh}^-$ into $\mathcal{E}_{mh}^+$ such that the matrix $A_m$ with coefficients

$$ \mathrm{meas}(\Psi_m(e) \cap e'), \quad e \in \mathcal{E}_{mh}^-, e' \in \mathcal{E}_{mh}^-, $$

is invertible.

We now introduce the quantity, defined on all functions $q_h$ in $\mathbb{M}_h$,

$$ \|q_h\|_{\alpha^{-1}*h} = \left( \sum_{i=1} \alpha_{ib}^{-1} \sum_{K \in \mathcal{T}_h^i} |q_h|_K^2 \right)^{\frac{1}{2}}. \quad (4.5) $$

Indeed, the following lemma states that it is a norm on $\mathbb{M}_h$. We refer to [6, Lemma 10] for its proof.

**Lemma 4.2.** If Assumption 4.1 holds, there exists a constant $c$ independent of $h$ and $\alpha$ such that

$$ \forall q_h \in \mathbb{M}_h, \quad \|q_h\|_{L^2(\Omega)} \le c \alpha_{\#}^{\frac{1}{2}} \|q_h\|_{\alpha^{-1}*h}. \quad (4.6) $$

Exactly the same argument as for Lemma 2.2 leads to the discrete inf-sup condition.

**Lemma 4.3.** The following inf-sup condition holds

$$ \forall q_h \in \mathbb{M}_h, \quad \sup_{v_h \in \mathbb{X}_h} \frac{(v_h, \mathrm{grad} \, q_h)_h}{\|v_h\|_{\alpha}} \ge \|q_h\|_{\alpha^{-1}*h}. \quad (4.7) $$

In view of the next estimates, we introduce the space

$$ \mathcal{Z} = C^0(0, T; L^2(\Omega)^d) \times L^2(0, T; H^1(\Omega)), \quad (4.8) $$

provided with the norm

$$ \|V\|_{\mathcal{Z}} = \sup_{0 \le t \le T} \|v(t)\|_{L^2(\Omega)} + \left( \int_0^T \|v(t)\|_{\alpha}^2 dt + \int_0^T \|q\|_{H^1(\Omega)}^2 dt \right)^{\frac{1}{2}}. \quad (4.9) $$

Note moreover that this norm is fully equivalent on all $V = (i_{\delta t}v, \pi_{\delta t}^+q)$ to

$$ \|V\|_{\mathcal{Z}, \delta t} = \max_{0 \le n \le N} \|v^n\|_{L^2(\Omega)} + \left( \delta t \sum_{n=0}^{N} \|v^n\|_{\alpha}^2 + \delta t \sum_{n=1}^{N} \|q^n\|_{H^1(\Omega)}^2 \right)^{\frac{1}{2}}, \quad (4.10) $$
and on all $V = (i_{\delta t}v_h, \pi_{\delta t}^+q_h)$ with $v_h$ in $X_h$ and $q_h$ in $M_h$, to

$$ \|V\|_{Z,\delta} = \max_{0 \le n \le N} \|v_h^n\|_{L^2(\Omega)} + \left( \delta t \sum_{n=0}^{N} \|v_h^n\|_{\alpha}^2 + \delta t \sum_{n=1}^{N} \sum_{K \in \mathcal{T}_h} \|q_h^n\|_{H^1(K)}^2 \right)^{\frac{1}{2}}. \quad (4.11) $$

**Lemma 4.4.** The following stability property holds for any $\mathbf{f}$ in $C^0(0, T; L^2(\Omega)^d)$ and $\mathbf{u}_0$ in $L^2(\Omega)^d$

$$ \|\mathcal{T}_{\delta}(\mathbf{f}, \mathbf{u}_0, 0, 0)\|_{Z,\delta} \le c \left( \left( \delta t \sum_{n=1}^{N} \|\mathbf{f}(\cdot, n\delta t)\|_{\alpha^{-1}}^2 \right)^{\frac{1}{2}} + \|\mathbf{u}_0\|_{L^2(\Omega)^d} \right). \quad (4.12) $$

**Proof:** We set: $(\mathbf{u}_h, p_h)_{1 \le n \le N} = \mathcal{T}_{\delta}(\mathbf{f}, \mathbf{u}_0, 0, 0)$. By taking $v_h$ equal to $\mathbf{u}_h^n$ in (4.3) and noting that $p_h^n$ belongs to $M_h^{\Delta}$, we obtain

$$ \|\mathbf{u}_h^n\|_{L^2(\Omega)}^2 + \delta t \|\mathbf{u}_h^n\|_{\alpha}^2 \le (\mathbf{u}_h^{n-1}, \mathbf{u}_h^n) + \delta t (\mathbf{f}(\cdot, n\delta t), \mathbf{u}_h^n), $$

whence, by using twice the inequality $ab \le \frac{1}{2}a^2 + \frac{1}{2}b^2$,

$$ \|\mathbf{u}_h^n\|_{L^2(\Omega)^d}^2 + \delta t \|\mathbf{u}_h^n\|_{\alpha}^2 \le \|\mathbf{u}_h^{n-1}\|_{L^2(\Omega)^d}^2 + \delta t \|\mathbf{f}(\cdot, n\delta t)\|_{\alpha^{-1}}^2. $$

Thus, summing on $n$ and using the definition of $\Pi_h$ yields the bound for the first two terms of $\|\mathcal{T}_{\delta}(\mathbf{f}, \mathbf{u}_0, 0, 0)\|_{Z,\delta}$. To bound the third one, we take $v_h$ equal to $\mathbf{grad} \, p_h^n$ in (4.3) and observe that both $(\mathbf{u}_h^n, \mathbf{grad} \, p_h^n)$ and $(\mathbf{u}_h^{n-1}, \mathbf{grad} \, p_h^n)$ vanish, which leads to

$$ \left( \sum_{K \in \mathcal{T}_h} \| \mathbf{grad} \, p_h^n \|_{L^2(K)^d}^2 \right)^{\frac{1}{2}} \le c \left( \| \mathbf{u}_h^n \|_{\alpha} + \| \mathbf{f}(\cdot, n\delta t) \|_{\alpha^{-1}} \right). $$

We conclude by taking the square of this inequality, multiplying by $\delta t$ and summing on the $n$, finally using the generalized Poincaré-Friedrichs inequality (4.6).

We skip the proof of the next lemma since it relies on very standard arguments for the Euler's scheme applied to a linear problem (see [5] for Darcy's equations) and to the same arguments as in [6, Section 4] for the mortar finite element discretization. Let $h$ stand for the maximum of the $h_i$, $1 \le i \le I$.

**Lemma 4.5.** Assume that there exists a positive constant $\lambda$ independent of $h$ such that

$$ \sup_{1 \le m \le M} \sup_{\text{meas}(\partial \Omega_{\ell} \cap \gamma_m^{+}) > 0} \frac{\alpha_{\ell} h_{\ell}}{\alpha_m^{+} \underline{h}_m^{+}} < \lambda, \quad \text{with} \quad \underline{h}_m^{+} = \min_{K \in \mathcal{T}_h^{m+}, \text{meas}(\partial K \cap \gamma_m^{+}) > 0} h_K. \quad (4.13) $$

For any data $\mathbf{f}$, $\mathbf{u}_0$, $p_1$ and $g_2$ satisfying (3.15) such that $\mathcal{T}(\mathbf{f}, \mathbf{u}_0, p_1, g_2)$ belongs to the space

$$ \mathcal{Z}_+ = (C^1(0, T; L^2(\Omega)^d) \cap C^0(0, T; H^1(\Omega)^d)) \cap \mathcal{L}^2(0, T; H^2(\Omega)), \quad (4.14) $$

the following error estimate holds

$$ \|(\mathcal{T} - \mathcal{T}_{\delta})(\mathbf{f}, \mathbf{u}_0, p_1, g_2)\|_{Z,\delta} \le C(\mathbf{f}, \mathbf{u}_0, p_1, g_2) (\delta t + h), \quad (4.15) $$
where the constant $C(f, u_0, p_1, g_2)$ only depends on the norm of $\mathcal{T}(f, u_0, p_1, g_2)$ in $\mathcal{Z}_+$.

The next result is derived from Lemmas 4.4 and 4.5 in a standard way.

**Corollary 4.6.** If condition (4.13) is satisfied, the following convergence property holds for any $f$ in $\mathcal{C}^0(0, T; L^2(\Omega)^d)$ and $u_0$ in $L^2(\Omega)^d$,

$$ \lim_{\delta \to (0,0)} \|(\mathcal{T} - \mathcal{T}_\delta)(f, u_0, 0, 0)\|_{\mathcal{Z}, \delta} = 0. \quad (4.16) $$

## Preliminary lemmas

As standard for applying the theory of [8], we are led to make some assumptions on the solution $(u, p)$ that we want to approximate. We denote by $D$ the differential operator with respect to $U$.

**Assumption 4.7.** The solution $U = (u, p)$ of problem (2.4) – (2.5)

(i) belongs to the space $\mathcal{Z}_{++}$ defined by

$$ \mathcal{Z}_{++} = (\mathcal{C}^1(0, T; L^2(\Omega)^d) \cap \mathcal{C}^0(0, T; H^1(\Omega)^d)) \times (H^1(0, T; H^1(\Omega)) \cap L^2(0, T; H^2(\Omega))), \quad (4.17) $$

(ii) is such that $\mathcal{DF}(U)$ is an isomorphism of the space $\mathcal{C}^0(0, T; L^2(\Omega)^d) \times L^2(0, T; H^1(\Omega))$.

It can be noted that part (ii) is much weaker than the uniqueness result established in Proposition 2.7 since it only yields the local uniqueness of the solution. To go further, we introduce:

1) The space $\mathcal{X}_\delta$ of continuous functions which are affine on each interval $](n-1)\delta t, n\delta t[$, $1 \le n \le N$, and such that their value at each point $t_n$ belongs to $X_h$;

2) The space $\mathcal{M}_\delta$ of functions which are constant, equal to a function of $M_h$, on each interval $](n-1)\delta t, n\delta t[$, $1 \le n \le N$.

Due to part (i) of Assumption 4.7, it is readily checked that there exists a function $u_\delta^*$ in $\mathcal{X}_\delta$ which satisfies

$$ \sup_{0 \le t \le T} \|u(\cdot, t) - u_\delta^*(\cdot, t)\|_{L^2(\Omega)^d} \le c(u) (\delta t + h); \quad (4.18) $$

and a function $p_\delta^*$ in $\mathcal{M}_\delta$ such that

$$ \|p - p_\delta^*\|_{L^2(0,T;H^1_*(\Omega))} \le c(p)(\delta t + h). \quad (4.19) $$

Here, $H^1_*(\Omega)$ means the space of piecewise functions in $H^1(\Omega)$, its norm is the Hilbertian sum of the norms of $H^1(K)$ on all $K$ in $\mathcal{T}_h$.

We set $U_\delta^* = (u_\delta^*, p_\delta^*)$.

Let finally $\mathcal{E}_\delta$ stand for the space of endomorphisms of $\mathcal{X}_\delta \times \mathcal{M}_\delta$.

**Lemma 4.8.** If Assumption 4.7 and condition (4.13) hold, there exists a positive real number $\delta_0$ such that, if $\delta t + h \le \delta_0$, the operator $\mathcal{DF}_\delta(U_\delta^*)$ is an isomorphism of the space $\mathcal{X}_\delta \times \mathcal{M}_\delta$. Moreover the norm of its inverse is bounded independently of $\delta$.
**Proof:** We use the expansion

$$ D\mathcal{F}_{\delta}(U) = D\mathcal{F}(U) + (\mathcal{T} - \mathcal{T}_{\delta})D\mathcal{G}(U) + \mathcal{T}_{\delta}(D\mathcal{G}(U) - D\mathcal{G}(U_{\delta}^{*})) + \mathcal{T}_{\delta}(D\mathcal{G}(U_{\delta}^{*}) - D\mathcal{G}_{\delta}(U_{\delta}^{*})). $$

Owing to part (ii) of Assumption 4.7, it suffices to check that the last three terms tend to zero when $\delta$ tend to $(0, 0)$.

1) We have, for $W = (w, r)$ in $\mathcal{X}_{\delta} \times \mathcal{M}_{\delta}$,

$$ D\mathcal{G}(U).W = ((\bar{\alpha}(\cdot) - \alpha(\cdot, p))w - \alpha'(p)r \mathbf{u}, 0, 0, 0). $$

When $W$ runs through the unit sphere of $\mathcal{X}_{\delta} \times \mathcal{M}_{\delta}$ and since this space is finite-dimensional, $D\mathcal{G}(U).W$ runs through a compact set of $\mathcal{C}^0(0, T; L^2(\Omega)^d)$. Thus, it follows from Corollary 4.6 that

$$ \lim_{\delta \to (0,0)} \| (\mathcal{T} - \mathcal{T}_{\delta}) D\mathcal{G}(U) \|_{\varepsilon_{\delta}} = 0. \qquad (4.20) $$

2) By combining Lemma 4.4 and the continuity of the operator $D\mathcal{G}$, we deduce from (4.18) and (4.19) that

$$ \lim_{\delta \to (0,0)} \| \mathcal{T}_{\delta}(D\mathcal{G}(U) - D\mathcal{G}(U_{\delta}^{*})) \|_{\varepsilon_{\delta}} = 0. \qquad (4.21) $$

3) Since $i_{\delta t}$ and $\pi_{\delta t}^{-}$ are equal to the identity on $\mathcal{X}_{\delta}$ and $\mathcal{M}_{\delta}$, respectively, we have, for $W = (w, r)$ in $\mathcal{X}_{\delta} \times \mathcal{M}_{\delta}$,

$$ D\mathcal{G}_{\delta}(U_{\delta}^{*}).W = ((\bar{\alpha}(\cdot) - \alpha(\cdot, p_{\delta}^{*}))w - \alpha'(p_{\delta}^{*})r \mathbf{u}, 0, 0, 0) = D\mathcal{G}(U_{\delta}^{*}).W, $$

whence

$$ \mathcal{T}_{\delta}(D\mathcal{G}(U_{\delta}^{*}) - D\mathcal{G}_{\delta}(U_{\delta}^{*})) = 0. \qquad (4.22) $$

We obtain the desired result from (4.20), (4.21) and (4.22).

We skip the proof of the next lemma since it is a direct consequence of Lemma 4.4, the definition of $D\mathcal{G}_{\delta}$ and the Lipschitz continuity of $\alpha$.

**Lemma 4.9.** There exists a neighbourhood of $U_{\delta}^{*}$ in $\mathcal{X}_{\delta} \times \mathcal{M}_{\delta}$ such that the following Lipschitz property holds for any $V_{\delta} = (v_{\delta}, q_{\delta})$ in this neighbourhood

$$ \| \mathcal{T}_{\delta}(D\mathcal{G}_{\delta}(U_{\delta}^{*}) - D\mathcal{G}_{\delta}(V_{\delta})) \|_{\varepsilon_{\delta}} \\ \leq c(\mathbf{u}, p) \left( \sup_{0 \leq t \leq T} \| \mathbf{u}_{\delta}^{*}(\cdot, t) - v_{\delta}(\cdot, t) \|_{L^2(\Omega)^d} + \| p_{\delta}^{*} - q_{\delta} \|_{L^2(0,T; H_{*}^{1}(\Omega))} \right). \qquad (4.23) $$

The last lemma provides an estimate of the quantity $\|\mathcal{F}_{\delta}(U_{\delta}^{*})\|_{\mathcal{Z}}$.

**Lemma 4.10.** If Assumption 4.7 and condition (4.13) hold, the following estimate for the quantity $\varepsilon_{\delta} = \|\mathcal{F}_{\delta}(U_{\delta}^{*})\|_{\mathcal{Z}}$ is satisfied:

$$ \varepsilon_{\delta} \leq c(\mathbf{u}, p)(\delta t + h). \qquad (4.24) $$
**Proof:** Recalling that $\mathcal{F}(U) = 0$, we have

$$ \varepsilon_{\delta} \le \|(T - \mathcal{T}_{\delta})\mathcal{G}(U)\|_{\mathcal{Z}} + \|\mathcal{T}_{\delta}(\mathcal{G}(U) - \mathcal{G}(U_{\delta}^{*}))\|_{\mathcal{Z}} + \|\mathcal{T}_{\delta}(\mathcal{G}(U_{\delta}^{*}) - \mathcal{G}_{\delta}(U_{\delta}^{*}))\|_{\mathcal{Z}}. $$

The estimate for the first term is an immediate consequence of Lemma 4.5, the estimate for the second term is easily derived from Lemma 4.4, the continuity of $\mathcal{G}$, (4.18) and (4.19). Finally, for the same reasons as in the proof of Lemma 4.8, the last term is zero. This concludes the proof.

## The final result

Owing to Lemmas 4.8 to 4.10, all assumptions needed for applying [8, Thm 1] (see also [12, Chap. IV, Thm 3.1]) are satisfied. This leads to the final a priori error estimate.

**Theorem 4.11.** If Assumption 4.7 and condition (4.13) hold, there exist a positive real number $\delta_0^*$ and a neighbourhood of $U_\delta^*$ in $\mathcal{X}_\delta \times \mathcal{M}_\delta$ such that, if $\delta t + h \le \delta_0^*$, problems (3.13) – (3.14) has a unique solution such that the associated pair $(u_\delta, p_\delta)$ belongs to this neighbourhood. Moreover, the following error estimate holds

$$ \sup_{0 \le t \le T} \|u(\cdot, t) - u_\delta(\cdot, t)\|_{L^2(\Omega)^d} + \|p - p_\delta\|_{L^2(0,T; H_\delta^1(\Omega))} \le c(u, p)(\delta t + h). \quad (4.25) $$

Estimate (4.25) is fully optimal. Moreover, condition (4.13) is often verified in practical situations and the convergence of the method can be proved with a weaker assumption on the regularity of $(u, p)$ than required in Assumption 4.7. So the mortar finite element method that we propose seems very efficient for the discretization of this rather complex problem.
References

[1] Y. Achdou, C. Bernardi — Un schéma de volumes ou éléments finis adaptatif pour les équations de Darcy à perméabilité variable, C.R. Acad. Sci. Paris 333 série I (2001), 693–698.

[2] Y. Achdou, C. Bernardi, F. Coquel — A priori and a posteriori analysis of finite volume discretizations of Darcy’s equations, Numer. Math. 96 (2003), 17–42.

[3] M. Azaïez, F. Ben Belgacem, C. Bernardi, N. Chorfi — Spectral discretization of Darcy’s equations with pressure dependent porosity, Appl. Math. Comput. 217 (2010), 1838–1856.

[4] A. Bergam, C. Bernardi, Z. Mghazli — A posteriori analysis of the finite element discretization of some parabolic equations, Math. Comput. 74 (2005), 1117–1138.

[5] C. Bernardi, V. Girault, K.R. Rajagopal — Discretization of an unsteady flow through a porous solid modeled by Darcy’s equations, Math. Models Methods Appl. Sci. 18 (2008), 2087–2123.

[6] C. Bernardi, F. Hecht, Z. Mghazli — Mortar finite element discretization for the flow in a non homogeneous porous medium, Comput. Methods Appl. Mech. Engrg. 196 (2007), 1554–1573.

[7] C. Bernardi, Y. Maday, A.T. Patera — A new nonconforming approach to domain decomposition: the mortar element method, Collège de France Seminar XI, H. Brezis & J.-L. Lions eds., Pitman (1994), 13–51.

[8] F. Brezzi, J. Rappaz, P.-A. Raviart — Finite dimensional approximation of nonlinear problems, Part I: Branches of nonsingular solutions, Numer. Math. 36 (1980), 1–25.

[9] M. Crouzeix, P.-A. Raviart — Conforming and nonconforming finite element methods for solving the stationary Stokes equations, R.A.I.R.O. Anal. Numér. 7 R3 (1973), 33–76.

[10] M. Daadaa — Discrétisation spectrale et par éléments spectraux des équations de Darcy, PhD Thesis, Université Pierre et Marie Curie, Paris (2009).

[11] V. Girault, P.-A. Raviart — Finite Element Approximation of the Navier–Stokes Equations, Lecture Notes in Mathematics 749, Springer-Verlag (1979).

[12] V. Girault, P.-A. Raviart — Finite Element Methods for Navier–Stokes Equations, Theory and Algorithms, Springer–Verlag (1986).

[13] J.-L. Lions, E. Magenes — Problèmes aux limites non homogènes et applications, Vol. I, Dunod, Paris (1968).

[14] K.R. Rajagopal — On a hierarchy of approximate models for flows of incompressible fluids through porous solid, Math. Models Methods Appl. Sci. 17 (2007), 215–252.