# Coupling linear virtual element and non-linear finite volume methods for poroelasticity

Guillaume Enchéry, Léo Agélas

► To cite this version:

Guillaume Enchéry, Léo Agélas. Coupling linear virtual element and non-linear finite volume methods for poroelasticity. Comptes Rendus. Mécanique, 2023, 351 (S1), pp.395-410. 10.5802/crmeca.225. hal-04618593

HAL Id: hal-04618593

https://ifp.hal.science/hal-04618593

Submitted on 20 Jun 2024

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
# Coupling linear virtual element and non-linear finite volume methods for poroelasticity

Guillaume Enchéry, Léo Agélas

► To cite this version:

Guillaume Enchéry, Léo Agélas. Coupling linear virtual element and non-linear finite volume methods for poroelasticity. Comptes Rendus. Mécanique, 2023, 351 (S1), pp.395-410. 10.5802/crmeca.225. hal-04618593

HAL Id: hal-04618593

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

Submitted on 20 Jun 2024

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
# Comptes Rendus
## Mécanique

Guillaume Enchéry and Léo Agélas

**Coupling linear virtual element and non-linear finite volume methods for
poroelasticity**

Volume 351, Special Issue S1 (2023), p. 395-410

Online since: 17 November 2023

**Part of Special Issue:** The scientific legacy of Roland Glowinski

**Guest editors:** Grégoire Allaire (CMAP, Ecole Polytechnique, Institut Polytechnique de Paris, Palaiseau, France),
Jean-Michel Coron (Laboratoire Jacques-Louis Lions, Sorbonne Université)
and Vivette Girault (Laboratoire Jacques-Louis Lions, Sorbonne Université)

https://doi.org/10.5802/crmeca.225

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

©

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

This article is licensed under the
CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.

http://creativecommons.org/licenses/by/4.0/
# Comptes Rendus
## Mécanique

Guillaume Enchéry and Léo Agélas

**Coupling linear virtual element and non-linear finite volume methods for
poroelasticity**

Volume 351, Special Issue S1 (2023), p. 395-410

Online since: 17 November 2023

**Part of Special Issue:** The scientific legacy of Roland Glowinski

**Guest editors:** Grégoire Allaire (CMAP, Ecole Polytechnique, Institut Polytechnique de Paris, Palaiseau, France),
Jean-Michel Coron (Laboratoire Jacques-Louis Lions, Sorbonne Université)
and Vivette Girault (Laboratoire Jacques-Louis Lions, Sorbonne Université)

https://doi.org/10.5802/crmeca.225

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

This article is licensed under the
CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.
http://creativecommons.org/licenses/by/4.0/

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

The Comptes Rendus. Mécanique are a member of the
*Mersenne Center for open scientific publishing*
www.centre-mersenne.org — e-ISSN : 1873-7234
<figure><img src="image_4.png" /></figure>

# The scientific legacy of Roland Glowinski / L'héritage scientifique de Roland Glowinski

## Coupling linear virtual element and non-linear finite volume methods for poroelasticity

## Couplage de méthodes basées sur des éléments virtuels linéaires et des volumes-finis nonlinéaires pour la poroélasticité

Guillaume Enchéry*,a and Léo Agélas a

**Abstract.** Solving poroelastic problems on a single grid is of paramount importance in the applications, especially in geosciences. However, these applications sometimes require to work with meshes made of distorted cells. With such grids, dedicated numerical schemes should be chosen to obtain consistent approximations for the stresses and for the fluxes. In J. Coulet et al. (2020), a coupled scheme based on the joint use of linear virtual elements and a two-point finite volume approximation on the same grid was proposed for Biot's poroelastic problem. This work has also provided an analytic and numerical convergence study of the discrete coupled system. As a continuation, we here propose an extension of this study to more general polyhedral meshes and to heterogeneous anisotropic mobility tensors for the flow equations. At this occasion, a general finite-volume framework, which includes non-linear or sushi finite volume methods for instance, is introduced to extend this analysis.

**Résumé.** La résolution des problèmes de poroélasticité sur une grille unique est d'une importance capitale dans les applications, en particulier dans le domaine des géosciences. Cependant, ces applications nécessitent parfois de travailler avec des grilles composées de cellules déformées. Avec de telles grilles, des schémas numériques spécifiques doivent être choisis pour obtenir des approximations convergentes pour les contraintes et les flux. Dans J. Coulet et al. (2020), un schéma couplé basé sur l'utilisation conjointe d'éléments virtuels linéaires et d'une approximation de volume fini à deux points sur la même grille a été proposé pour le problème poroélastique de Biot. Ce travail a également fourni une étude analytique et numérique de la convergence du système couplé discret. Comme suite, nous proposons ici une extension de cette étude à des maillages polyédriques plus généraux et à des tenseurs de mobilité anisotropes hétérogènes pour les équations d'écoulement. À cette occasion, un cadre général de volumes finis, qui inclut par exemple les méthodes de volumes finis non linéaires ou de sushi, est introduit pour étendre cette analyse.

**Keywords.** finite volume, virtual element, poroelasticity.

**Mots-clés.** volume fini, élément virtuel, poroélasticité.

**Electronic supplementary material.** Supplementary material for this article is supplied as a separate archive available from the journal's website under article's URL or from the author.

<footnote>a IFPEN, 1-4 av. du Bois Préau 92852 Rueil-Malmaison France</footnote>

<footnote>E-mails: guillaume.enchery@ifpen.fr (G. Enchéry), leo.agelas@ifpen.fr (L. Agélas)</footnote>

<footnote>* Corresponding author.</footnote>
<figure><img src="image_4.png" /></figure>

The scientific legacy of Roland Glowinski / L'héritage scientifique
de Roland Glowinski

Coupling linear virtual element and
non-linear finite volume methods for
poroelasticity

Couplage de méthodes basées sur des éléments virtuels
linéaires et des volumes-finis nonlinéaires pour la
poroélasticité

Guillaume Enchéry *, a and Léo Agélas a

**Abstract.** Solving poroelastic problems on a single grid is of paramount importance in the applications, especially in geosciences. However, these applications sometimes require to work with meshes made of distorted cells. With such grids, dedicated numerical schemes should be chosen to obtain consistent approximations for the stresses and for the fluxes. In J. Coulet et al. (2020), a coupled scheme based on the joint use of linear virtual elements and a two-point finite volume approximation on the same grid was proposed for Biot's poroelastic problem. This work has also provided an analytic and numerical convergence study of the discrete coupled system. As a continuation, we here propose an extension of this study to more general polyhedral meshes and to heterogeneous anisotropic mobility tensors for the flow equations. At this occasion, a general finite-volume framework, which includes non-linear or sushi finite volume methods for instance, is introduced to extend this analysis.

**Résumé.** La résolution des problèmes de poroélasticité sur une grille unique est d'une importance capitale dans les applications, en particulier dans le domaine des géosciences. Cependant, ces applications nécessitent parfois de travailler avec des grilles composées de cellules déformées. Avec de telles grilles, des schémas numériques spécifiques doivent être choisis pour obtenir des approximations convergentes pour les contraintes et les flux. Dans J. Coulet et al. (2020), un schéma couplé basé sur l'utilisation conjointe d'éléments virtuels linéaires et d'une approximation de volume fini à deux points sur la même grille a été proposé pour le problème poroélastique de Biot. Ce travail a également fourni une étude analytique et numérique de la convergence du système couplé discret. Comme suite, nous proposons ici une extension de cette étude à des maillages polyédriques plus généraux et à des tenseurs de mobilité anisotropes hétérogènes pour les équations d'écoulement. À cette occasion, un cadre général de volumes finis, qui inclut par exemple les méthodes de volumes finis non linéaires ou de sushi, est introduit pour étendre cette analyse.

**Keywords.** finite volume, virtual element, poroelasticity.

**Mots-clés.** volume fini, élément virtuel, poroélasticité.

**Electronic supplementary material.** Supplementary material for this article is supplied as a separate archive available from the journal's website under article's URL or from the author.

<footnote>a IFPEN, 1-4 av. du Bois Préau 92852 Rueil-Malmaison France</footnote>

<footnote>E-mails: guillaume.enchery@ifpen.fr (G. Enchéry), leo.agelas@ifpen.fr (L. Agélas)</footnote>

<footnote>* Corresponding author.</footnote>
Manuscript received 31 August 2023, accepted 6 September 2023.

# 1. Introduction

The numerical resolution of poroelastic problems is often carried out on two different grids, in particular in geosciences when simulating the flow through deformable porous media. For this kind of application, a first grid is dedicated to mechanics and built in order to guarantee its compatibility with the finite-element method which is often used to discretize that part of the system. A second mesh is also designed in parallel to discretize the flow equations with finite-volume methods. This second family of discretization schemes is very popular to simulate porous-media flows since they provide conservative fluxes in a natural way. They are based on piecewise-constant approximations of the solutions and the flow grid is therefore built according to the shape of the geological layers which can have strongly varying flow properties and complicated geometries due, for example, to the presence of faults or erosion events. A space discretisation of these layers is therefore not an easy task since, in addition, some conditions should be met in order to guarantee the consistency of the finite-volume fluxes.

Recent space discretization methods, such as hybrid finite-volumes [1], virtual elements [2], now enable to solve the elastic problem on the flow grid. We do not claim to make an exhaustive review of all the methods that have been proposed over the past years to discretize this system. We here only mention the ones that will be used and analyzed in this work and invite the reader to refer to the introduction of [3] to have a broader overview of the family of techniques which can be used for this problem.

In this latter work, in particular, linear virtual elements and the two-point finite volume approximation (TPFA) were used and analysed. However, a drawback of the TPFA scheme is that, for practical applications, the TPFA scheme does not provide consistent fluxes on general grids and/or with anisotropic mobility tensors. In the past years, several works have been done to correct this shortcoming (see [4, 5] and the references therein for instance) leading to linear or even non-linear approximations of the fluxes.

In this work, we show that the convergence and existence results established in [3] can be extended to more general finite-volume schemes and in particular to non-linear ones.

This paper is thus organized as follows. The poroelastic problem is introduced in Section 2. Then the coupled discretization scheme based on linear virtual elements and on a generic family of finite volume methods is presented in Section 3. In particular, a reminder on the linear virtual element method is made in Section 3.1 and the generic finite-volume scheme is presented in Section 3.2 along with two examples of flux approximation. This coupled scheme is analysed in Section 4 where the existence of discrete solutions (see Proposition 7) is shown as well as an a priori error estimate (see Theorem 8). To illustrate this estimate, numerical results performed on grids of the FVCA5 benchmark [6] are presented in Section 5.

# 2. Mathematical formulation of the poroelasticity problem

Let $T > 0$ and $\Omega \in \mathbb{R}^d$ be an open subset with a Lipschitz boundary. The coupled system describing the interactions between mechanical deformations and flows in porous media have been established in [7]. The expressions of rock mechanical equilibrium and fluid volume conservation here assume quasistatic strains and a slightly compressible single-phase flow. Under these assumptions, let $\mathbf{u} : \Omega \times (0, T] \to \mathbb{R}^d$ be the rock displacement and $p : \Omega \times (0, T] \to \mathbb{R}$ the pore pressure. In $\Omega \times (0, T]$, the pair $(\mathbf{u}, p)$ is solution to:

$$ -\operatorname{div}(\bar{\sigma}^e(\mathbf{u}) - \alpha p \bar{I}_d) = \mathbf{f}, \quad (1) $$

$$ \partial_t (c_0 p + \alpha \operatorname{div}(\mathbf{u})) - \operatorname{div}(\bar{\Lambda}(\nabla p - \rho \mathbf{g})) = s. \quad (2) $$
Manuscript received 31 August 2023, accepted 6 September 2023.

# 1. Introduction

The numerical resolution of poroelastic problems is often carried out on two different grids, in particular in geosciences when simulating the flow through deformable porous media. For this kind of application, a first grid is dedicated to mechanics and built in order to guarantee its compatibility with the finite-element method which is often used to discretize that part of the system. A second mesh is also designed in parallel to discretize the flow equations with finite-volume methods. This second family of discretization schemes is very popular to simulate porous-media flows since they provide conservative fluxes in a natural way. They are based on piecewise-constant approximations of the solutions and the flow grid is therefore built according to the shape of the geological layers which can have strongly varying flow properties and complicated geometries due, for example, to the presence of faults or erosion events. A space discretisation of these layers is therefore not an easy task since, in addition, some conditions should be met in order to guarantee the consistency of the finite-volume fluxes.

Recent space discretization methods, such as hybrid finite-volumes [1], virtual elements [2], now enable to solve the elastic problem on the flow grid. We do not claim to make an exhaustive review of all the methods that have been proposed over the past years to discretize this system. We here only mention the ones that will be used and analyzed in this work and invite the reader to refer to the introduction of [3] to have a broader overview of the family of techniques which can be used for this problem.

In this latter work, in particular, linear virtual elements and the two-point finite volume approximation (TPFA) were used and analysed. However, a drawback of the TPFA scheme is that, for practical applications, the TPFA scheme does not provide consistent fluxes on general grids and/or with anisotropic mobility tensors. In the past years, several works have been done to correct this shortcoming (see [4, 5] and the references therein for instance) leading to linear or even non-linear approximations of the fluxes.

In this work, we show that the convergence and existence results established in [3] can be extended to more general finite-volume schemes and in particular to non-linear ones.

This paper is thus organized as follows. The poroelastic problem is introduced in Section 2. Then the coupled discretization scheme based on linear virtual elements and on a generic family of finite volume methods is presented in Section 3. In particular, a reminder on the linear virtual element method is made in Section 3.1 and the generic finite-volume scheme is presented in Section 3.2 along with two examples of flux approximation. This coupled scheme is analysed in Section 4 where the existence of discrete solutions (see Proposition 7) is shown as well as an a priori error estimate (see Theorem 8). To illustrate this estimate, numerical results performed on grids of the FVCA5 benchmark [6] are presented in Section 5.

# 2. Mathematical formulation of the poroelasticity problem

Let $T > 0$ and $\Omega \in \mathbb{R}^d$ be an open subset with a Lipschitz boundary. The coupled system describing the interactions between mechanical deformations and flows in porous media have been established in [7]. The expressions of rock mechanical equilibrium and fluid volume conservation here assume quasistatic strains and a slightly compressible single-phase flow. Under these assumptions, let $u: \Omega \times (0, T] \to \mathbb{R}^d$ be the rock displacement and $p: \Omega \times (0, T] \to \mathbb{R}$ the pore pressure. In $\Omega \times (0, T]$, the pair $(u, p)$ is solution to:

$$ -\operatorname{div}(\bar{\sigma}^e(u) - \alpha p \bar{I}_d) = f, \qquad (1) $$

$$ \partial_t (c_0 p + \alpha \operatorname{div}(u)) - \operatorname{div}(\bar{\Lambda}(\nabla p - \rho g)) = s. \qquad (2) $$
In (1), $\bar{I}_d$ denotes the identity tensor, $\bar{\sigma}^e(\mathbf{u}) - \alpha p \bar{I}_d$ corresponds to the total stress, $\bar{\sigma}^e(\mathbf{u})$ to Terzaghi's effective stress and $\alpha$ to Biot-Willis' coefficient. The sign convention for the stresses follows the one of classical continuum mechanics where compressive stresses are considered to be negative. The effective stress is assumed to follow a linear elastic behaviour with respect to the strain $\bar{\varepsilon}(\mathbf{u}) := \frac{1}{2}(\nabla \mathbf{u} + \nabla \mathbf{u}^T)$ and is thus given by $\bar{\sigma}^e(\mathbf{u}) := \bar{C}\bar{\varepsilon}(\mathbf{u})$ where $\bar{C}$ is the fourth-order stiffness tensor. The stress tensor can be defined thanks to the Lamé constants $(\mu, \lambda) \in [\mu_1, \mu_2] \times (0, \infty)$, $0 < \mu_1 < \mu_2$. In that case

$$ \bar{\sigma}^e(\mathbf{u}) = 2\mu\bar{\varepsilon}(\mathbf{u}) + \lambda \operatorname{div}(\mathbf{u})\bar{I}_d. \quad (3) $$

$f$ denotes a volumetric force. In equation (2), $c_0 > 0$ is the product of porosity by the fluid compressibility, $\bar{\Lambda}$ the mobility tensor defined as the ratio between the rock permeability and the fluid dynamic viscosity, $\rho g$ is the volumetric weight and $s$ a volumetric source or sink term. For $\bar{\Lambda}$, it is assumed that there exist $(\alpha_0, \beta_0) \in (0, \infty)$ s.t. the spectrum of $\bar{\Lambda}$ is included in $[\alpha_0, \beta_0]$. Using the index $\alpha \in \{d, p\}$ to refer to the displacement or pressure variable, in both cases, the boundary $\partial\Omega$ is divided into two complementary areas where either Dirichlet or Neumann conditions are applied. Thus $\partial\Omega = \Gamma_{D_a} \cup \Gamma_{N_a}$, with $\Gamma_{D_a} \cap \Gamma_{N_a} = \emptyset$ and $m_{\Gamma_{D_a}} > 0$. On these boundaries, the conditions are given by:

$$ \mathbf{u} = 0 \text{ on } \Gamma_{D_d}, \quad (4) $$

$$ (\bar{\sigma}^e(\mathbf{u}) - \alpha p \bar{I}_d) \cdot \mathbf{n} = \mathbf{t} \text{ on } \Gamma_{N_d}, \quad (5) $$

$$ p = 0 \text{ on } \Gamma_{D_p}, \quad (6) $$

$$ \bar{\Lambda}(\nabla p - \rho g) \cdot \mathbf{n} = \phi \text{ on } \Gamma_{N_p}, \quad (7) $$

where $\mathbf{n}$ stands for the unit normal oriented towards the outside of $\Omega$. At initial time, we have

$$ p(\mathbf{x}, 0) = p_0(\mathbf{x}) \text{ in } \Omega \quad (8) $$

and $\mathbf{u}_0 = \mathbf{u}(\mathbf{x}, 0)$ is solution to (1)-(4)-(5) with $p = p_0$. Under assumptions of regularity on the data, the existence and uniqueness of the continuous solution to (1)–(8) have been established in [8].

### 3. Coupled discretization scheme

We now detail the discretization of the problem (1)–(8). From then on, for simplicity's sake, gravity effects are neglected in the following. For scalar functions defined on any domain $\omega$, we use the norms:

$$ \|v\|_{L^2(\omega)} = \left(\int_{\Omega} v^2\right)^{\frac{1}{2}}, \quad \|v\|_{H^1(\omega)} = \left(\|v\|_{L^2(\omega)}^2 + \int_{\omega} \nabla v \cdot \nabla v\right)^{\frac{1}{2}} $$

and

$$ \|\mathbf{v}\|_{L^2(\omega)} = \left(\int_{\omega} \mathbf{v} \cdot \mathbf{v}\right)^{\frac{1}{2}}, \quad \|\mathbf{v}\|_{H^1(\omega)} = \left(\|\mathbf{v}\|_{L^2(\omega)}^2 + \int_{\omega} \nabla \mathbf{v} : \nabla \mathbf{v}\right)^{\frac{1}{2}} $$

for the vector ones. The aim of this section is to define a coupled scheme to approximate a weak solution to (1)–(8). Weak solutions $(\tilde{\mathbf{u}}, \tilde{\mathbf{p}})$ of the previous problem are sought in the spaces

$$ V_0 = \{ \mathbf{v} \in H^1(\Omega)^d : \mathbf{v}_{\Gamma_{D_d}} = \mathbf{0} \}, $$

$$ Q_0 = \{ q \in H^1(\Omega) : q_{\Gamma_{D_p}} = 0 \}. $$

We assume that the time interval is discretized with a constant time step $\Delta t$ such that

$$ [0, T] = \bigcup_{n=0}^{N} [t^n, t^{n+1}] = \bigcup_{n=0}^{N} [n\Delta t, (n+1)\Delta t]. $$
In (1), $\overline{I}_d$ denotes the identity tensor, $\overline{\sigma}^e(\mathbf{u}) - \alpha p \overline{I}_d$ corresponds to the total stress, $\overline{\sigma}^e(\mathbf{u})$ to Terzaghi's effective stress and $\alpha$ to Biot-Willis' coefficient. The sign convention for the stresses follows the one of classical continuum mechanics where compressive stresses are considered to be negative. The effective stress is assumed to follow a linear elastic behaviour with respect to the strain $\overline{\varepsilon}(\mathbf{u}) := \frac{1}{2}(\nabla \mathbf{u} + \nabla \mathbf{u}^T)$ and is thus given by $\overline{\sigma}^e(\mathbf{u}) := \overline{C\varepsilon}(\mathbf{u})$ where $\overline{C}$ is the fourth-order stiffness tensor. The stress tensor can be defined thanks to the Lamé constants $(\mu, \lambda) \in [\mu_1, \mu_2] \times (0, \infty)$, $0 < \mu_1 < \mu_2$. In that case

$$ \overline{\sigma}^e(\mathbf{u}) = 2\mu\overline{\varepsilon}(\mathbf{u}) + \lambda \operatorname{div}(\mathbf{u})\overline{I}_d. \quad (3) $$

$f$ denotes a volumetric force. In equation (2), $c_0 > 0$ is the product of porosity by the fluid compressibility, $\overline{\Lambda}$ the mobility tensor defined as the ratio between the rock permeability and the fluid dynamic viscosity, $\rho g$ is the volumetric weight and $s$ a volumetric source or sink term. For $\overline{\Lambda}$, it is assumed that there exist $(\alpha_0, \beta_0) \in (0, \infty)$ s.t. the spectrum of $\overline{\Lambda}$ is included in $[\alpha_0, \beta_0]$. Using the index $\alpha \in \{d, p\}$ to refer to the displacement or pressure variable, in both cases, the boundary $\partial\Omega$ is divided into two complementary areas where either Dirichlet or Neumann conditions are applied. Thus $\partial\Omega = \Gamma_{D_a} \cup \Gamma_{N_a}$, with $\Gamma_{D_a} \cap \Gamma_{N_a} = \emptyset$ and $m_{\Gamma_{D_a}} > 0$. On these boundaries, the conditions are given by:

$$ \mathbf{u} = 0 \text{ on } \Gamma_{D_d}, \quad (4) $$

$$ (\overline{\sigma}^e(\mathbf{u}) - \alpha p \overline{I}_d) \cdot \mathbf{n} = \mathbf{t} \text{ on } \Gamma_{N_d}, \quad (5) $$

$$ p = 0 \text{ on } \Gamma_{D_p}, \quad (6) $$

$$ \overline{\Lambda}(\nabla p - \rho g) \cdot \mathbf{n} = \phi \text{ on } \Gamma_{N_p}, \quad (7) $$

where $\mathbf{n}$ stands for the unit normal oriented towards the outside of $\Omega$. At initial time, we have

$$ p(\mathbf{x}, 0) = p_0(\mathbf{x}) \text{ in } \Omega \quad (8) $$

and $\mathbf{u}_0 = \mathbf{u}(\mathbf{x}, 0)$ is solution to (1)-(4)-(5) with $p = p_0$. Under assumptions of regularity on the data, the existence and uniqueness of the continuous solution to (1)–(8) have been established in [8].

### 3. Coupled discretization scheme

We now detail the discretization of the problem (1)–(8). From then on, for simplicity's sake, gravity effects are neglected in the following. For scalar functions defined on any domain $\omega$, we use the norms:

$$ \|v\|_{L^2(\omega)} = \left(\int_{\Omega} v^2\right)^{\frac{1}{2}}, \quad \|v\|_{H^1(\omega)} = \left(\|v\|_{L^2(\omega)}^2 + \int_{\omega} \nabla v \cdot \nabla v\right)^{\frac{1}{2}} $$

and

$$ \|\mathbf{v}\|_{L^2(\omega)} = \left(\int_{\omega} \mathbf{v} \cdot \mathbf{v}\right)^{\frac{1}{2}}, \quad \|\mathbf{v}\|_{H^1(\omega)} = \left(\|\mathbf{v}\|_{L^2(\omega)}^2 + \int_{\omega} \nabla \mathbf{v} : \nabla \mathbf{v}\right)^{\frac{1}{2}} $$

for the vector ones. The aim of this section is to define a coupled scheme to approximate a weak solution to (1)–(8). Weak solutions $(\tilde{\mathbf{u}}, \tilde{\mathbf{p}})$ of the previous problem are sought in the spaces

$$ V_0 = \{ \mathbf{v} \in H^1(\Omega)^d : \mathbf{v}_{\Gamma_{D_d}} = \mathbf{0} \}, $$

$$ Q_0 = \{ q \in H^1(\Omega) : q_{\Gamma_{D_p}} = 0 \}. $$

We assume that the time interval is discretized with a constant time step $\Delta t$ such that

$$ [0, T] = \bigcup_{n=0}^{N} [t^n, t^{n+1}] = \bigcup_{n=0}^{N} [n\Delta t, (n+1)\Delta t]. $$
If $f \in H^1((0, T) \times L^2(\Omega))^d$ and $t \in H^1((0, T) \times L^2(\Gamma_{N_d}))^d$, an integration by part of (1) over $\Omega \times (t^n, t^{n+1})$, for all $n \in \{0, \dots, N\}$ and the use of an Euler time-implicit scheme lead to a variational form of the elasticity problem which reads, for any test function $v \in V_0$:

$$a(\mathbf{u}^{n+1}, \mathbf{v}) - \int_{\Omega} \alpha p^{n+1} \mathrm{div}(\mathbf{v}) = \int_{\Omega} \mathbf{f}^{n+1} \cdot \mathbf{v} + \int_{\Gamma_{N_d}} \mathbf{t}^{n+1} \cdot \mathbf{v} \quad (9)$$

with

$$a(\mathbf{u}, \mathbf{v}) = \int_{\Omega} \bar{\sigma}^{e}(\mathbf{u}) : \bar{\varepsilon}(\mathbf{v}). \quad (10)$$

Since $m_{\Gamma_{D_d}} > 0$, a corollary of Korn's inequality used with (3) (see [9, Chapter 11]) gives the coercivity of $a$, that is, for all $v \in V_0$,

$$a(v, v) \geq \alpha_a \|v\|_{H^1(\Omega)}^2. \quad (11)$$

We also define the norm $\beta_a$ of $a$ by

$$\beta_a = \sup_{v \in H^1(\Omega)^d, \|v\|_{H^1(\Omega)}^d = 1} a(v, v). \quad (12)$$

In the same way as for (1), assuming that $s \in L^2(\Omega \times (0, T))$ and $\phi \in L^2(\Gamma_{N_p} \times (0, T))$, a multiplication of (2) by any test function $q \in Q_0$ followed by an integration over $\Omega$ and the use of an Euler time-implicit scheme give, for all $n \in \{0, \dots, N\}$:

$$\int_{\Omega} (D^{n+1} (c_0 p + \alpha \mathrm{div}(\mathbf{u})) q) - \Delta t \int_{\Omega} \mathrm{div}(\bar{\Lambda} \nabla p^{n+1}) q = \Delta t \int_{\Omega} s^{n+1} q \quad (13)$$

where $D^{n+1}q = q^{n+1} - q^n$. The next two subsections will detail the virtual-element and finite-volume approximations but an admissible space discretization for both methods should be first defined. Following assumptions H0 in [2] and [5, Definition 1], we consider the following definition.

**Definition 1 (Admissible virtual-element/finite-volume discretization).** An admissible discretization $\mathcal{D}$ is a triplet $\mathcal{D} = (\mathcal{T}, \mathcal{E}, \mathcal{P})$, where

(i) $\mathcal{T}$ is a finite family of non-empty connected open disjoint subsets of $\Omega$ (the cells) s.t. $\bar{\Omega} = \cup_{K \in \mathcal{T}} \bar{K}$. For all $K \in \mathcal{T}$, we denote by $m_K > 0$ its $d$-dimensional measure (the volume) and set $\partial K \stackrel{\mathrm{def}}{=} \bar{K} \setminus K$;

(ii) $\mathcal{E}$ is a subset of $\bar{\Omega}$ (the faces) s.t., for all $\sigma \in \mathcal{E}$, $\sigma$ is a non-empty closed subset of a hyperplane of $\mathbb{R}^d$ with $(d-1)$-dimensional measure $m_\sigma > 0$ (the area), and s.t. the intersection of two different faces has zero $(d-1)$-dimensional measure. We assume that there exist, for all $K \in \mathcal{T}$, a subset $\mathcal{E}_K$ of $\mathcal{E}$ such that $\partial K = \cup_{\sigma \in \mathcal{E}_K} \sigma$ and $\rho_1 < +\infty$ independent of $n$ s.t.

$$|\mathcal{E}_K| \leq \rho_1. \quad (14)$$

For a given $\sigma \in \mathcal{E}$, either $\mathcal{T}_\sigma \stackrel{\mathrm{def}}{=} \{K \in \mathcal{T} | \sigma \in \mathcal{E}_K\}$ has exactly one element and then $\sigma \subset \partial\Omega$ (boundary face) or $\mathcal{T}_\sigma$ has exactly two elements (inner face); the sets of inner and boundary faces are denoted by $\mathcal{E}_{\mathrm{int}}$ and $\mathcal{E}_{\mathrm{ext}}$ respectively; for all $\alpha \in \{d, p\}$ and $\sigma \in \mathcal{E}_{\mathrm{ext}}$, $\sigma$ lies either completely in $\Gamma_{D_a}$ or in $\Gamma_{N_a}$; for all $K \in \mathcal{T}_n$ and $\sigma \in \mathcal{E}_K$, $\mathbf{n}_{K,\sigma}$ denotes the unit vector that is normal to $\sigma$ and outward to $K$;

(iii) $\mathcal{P} = \{\mathbf{x}_K\}_{K \in \mathcal{T}}$ s.t. $\mathbf{x}_K$ is the barycenter of $K$. $K$ is star-shaped with respect to a ball centered on $\mathbf{x}_K$ and with a radius $r h_K$ where $r > 0$ does not depend on $K$ or on $\mathcal{D}$. $r$ is smaller than the ratio between the shortest edge of $K$ and its diameter. For all $K \in \mathcal{T}$ and for all $\sigma \in \mathcal{E}_K$ we denote by $d_{K,\sigma}$ the Euclidean distance between $\mathbf{x}_K$ and the hyperplane supporting $\sigma$.
If $f \in H^1((0, T) \times L^2(\Omega))^d$ and $t \in H^1((0, T) \times L^2(\Gamma_{N_d}))^d$, an integration by part of (1) over $\Omega \times (t^n, t^{n+1})$, for all $n \in \{0, \dots, N\}$ and the use of an Euler time-implicit scheme lead to a variational form of the elasticity problem which reads, for any test function $v \in V_0$:

$$a(\mathbf{u}^{n+1}, \mathbf{v}) - \int_{\Omega} \alpha p^{n+1} \mathrm{div}(\mathbf{v}) = \int_{\Omega} \mathbf{f}^{n+1} \cdot \mathbf{v} + \int_{\Gamma_{N_d}} \mathbf{t}^{n+1} \cdot \mathbf{v} \quad (9)$$

with

$$a(\mathbf{u}, \mathbf{v}) = \int_{\Omega} \bar{\sigma}^{e}(\mathbf{u}) : \bar{\varepsilon}(\mathbf{v}). \quad (10)$$

Since $m_{\Gamma_{D_d}} > 0$, a corollary of Korn's inequality used with (3) (see [9, Chapter 11]) gives the coercivity of $a$, that is, for all $v \in V_0$,

$$a(v, v) \geq \alpha_a \|v\|_{H^1(\Omega)}^2. \quad (11)$$

We also define the norm $\beta_a$ of $a$ by

$$\beta_a = \sup_{v \in H^1(\Omega)^d, \|v\|_{H^1(\Omega)}^d = 1} a(v, v). \quad (12)$$

In the same way as for (1), assuming that $s \in L^2(\Omega \times (0, T))$ and $\phi \in L^2(\Gamma_{N_p} \times (0, T))$, a multiplication of (2) by any test function $q \in Q_0$ followed by an integration over $\Omega$ and the use of an Euler time-implicit scheme give, for all $n \in \{0, \dots, N\}$:

$$\int_{\Omega} (D^{n+1} (c_0 p + \alpha \mathrm{div}(\mathbf{u})) q) - \Delta t \int_{\Omega} \mathrm{div}(\bar{\Lambda} \nabla p^{n+1}) q = \Delta t \int_{\Omega} s^{n+1} q \quad (13)$$

where $D^{n+1}q = q^{n+1} - q^n$. The next two subsections will detail the virtual-element and finite-volume approximations but an admissible space discretization for both methods should be first defined. Following assumptions H0 in [2] and [5, Definition 1], we consider the following definition.

**Definition 1 (Admissible virtual-element/finite-volume discretization).** An admissible discretization $\mathcal{D}$ is a triplet $\mathcal{D} = (\mathcal{T}, \mathcal{E}, \mathcal{P})$, where

(i) $\mathcal{T}$ is a finite family of non-empty connected open disjoint subsets of $\Omega$ (the cells) s.t. $\bar{\Omega} = \cup_{K \in \mathcal{T}} \bar{K}$. For all $K \in \mathcal{T}$, we denote by $m_K > 0$ its $d$-dimensional measure (the volume) and set $\partial K \stackrel{\mathrm{def}}{=} \bar{K} \setminus K$;

(ii) $\mathcal{E}$ is a subset of $\bar{\Omega}$ (the faces) s.t., for all $\sigma \in \mathcal{E}$, $\sigma$ is a non-empty closed subset of a hyperplane of $\mathbb{R}^d$ with $(d-1)$-dimensional measure $m_\sigma > 0$ (the area), and s.t. the intersection of two different faces has zero $(d-1)$-dimensional measure. We assume that there exist, for all $K \in \mathcal{T}$, a subset $\mathcal{E}_K$ of $\mathcal{E}$ such that $\partial K = \cup_{\sigma \in \mathcal{E}_K} \sigma$ and $\rho_1 < +\infty$ independent of $n$ s.t.

$$|\mathcal{E}_K| \le \rho_1. \quad (14)$$

For a given $\sigma \in \mathcal{E}$, either $\mathcal{T}_\sigma \stackrel{\mathrm{def}}{=} \{K \in \mathcal{T} \mid \sigma \in \mathcal{E}_K\}$ has exactly one element and then $\sigma \subset \partial\Omega$ (boundary face) or $\mathcal{T}_\sigma$ has exactly two elements (inner face); the sets of inner and boundary faces are denoted by $\mathcal{E}_{\mathrm{int}}$ and $\mathcal{E}_{\mathrm{ext}}$ respectively; for all $\alpha \in \{d, p\}$ and $\sigma \in \mathcal{E}_{\mathrm{ext}}$, $\sigma$ lies either completely in $\Gamma_{D_a}$ or in $\Gamma_{N_a}$; for all $K \in \mathcal{T}_n$ and $\sigma \in \mathcal{E}_K$, $\mathbf{n}_{K,\sigma}$ denotes the unit vector that is normal to $\sigma$ and outward to $K$;

(iii) $\mathcal{P} = \{\mathbf{x}_K\}_{K \in \mathcal{T}}$ s.t. $\mathbf{x}_K$ is the barycenter of $K$. $K$ is star-shaped with respect to a ball centered on $\mathbf{x}_K$ and with a radius $r h_K$ where $r > 0$ does not depend on $K$ or on $\mathcal{D}$. $r$ is smaller than the ratio between the shortest edge of $K$ and its diameter. For all $K \in \mathcal{T}$ and for all $\sigma \in \mathcal{E}_K$ we denote by $d_{K,\sigma}$ the Euclidean distance between $\mathbf{x}_K$ and the hyperplane supporting $\sigma$.
When considering a sequence of admissible discretizations $(\mathcal{D}_n)_{n \in \mathbb{N}}$, we denote by $h_{\mathcal{D}_n}$ the size of the discretization defined by $h_{\mathcal{D}_n} \stackrel{\text{def}}{=} \sup_{K \in \mathcal{T}_n} \text{diam}(K)$. In that case, we suppose that there exist $0 < \rho_2, \rho_3, \rho_4, \rho_5 < +\infty$ independent of $n$ s.t.

$$ \begin{gathered} \min_{K \in \mathcal{T}_n, \sigma \in \mathcal{E}_{n,K}} \frac{d_{K,\sigma}}{\text{diam}(K)} \ge \rho_2, \quad \min_{\sigma \in \mathcal{E}_{n,\text{int}}, \mathcal{T}_{n,\sigma} = \{K,L\}} \frac{\min(d_{K,\sigma}, d_{L,\sigma})}{\max(d_{K,\sigma}, d_{L,\sigma})} \ge \rho_3, \quad \min_{K \in \mathcal{T}_n} \frac{\text{diam}(K)}{h_{\mathcal{D}_n}} \ge \rho_4, \\ \rho_5 \le \frac{d_{K,\sigma}}{d_{L,\sigma}} \le \frac{1}{\rho_5}, \quad \forall \sigma \in \mathcal{E}_{n,\text{int}}, \text{ such that } \mathcal{T}_{n,\sigma} = \{K,L\}. \end{gathered} $$

### 3.1. Virtual elements for the elasticity problem

In this work, the displacements are approximated using linear virtual elements. In that case, the approximation space $V_{\mathcal{D}}$ is given by

$$ V_{\mathcal{D}} = \bigcup_{K \in \mathcal{T}} V_K $$

where

$$ V_K = \{ \boldsymbol{\varphi} \in (H^1(K))^d \mid \boldsymbol{\varphi} \cdot \boldsymbol{e}_i \in V_{\text{scal},K}, \forall i = 1, \dots, d \}, $$

$e_i$ denoting one of the vectors of the canonical basis. The main idea of linear virtual elements is to define the spaces $V_{\text{scal},K}$ without giving an explicit definition of its basis functions even if, for simplicial cells, it reduces to $\mathbb{P}_1$ polynomials. Functions of $V_{\text{scal},K}$ are entirely defined by their values on the cell vertices of $K$. However, the exact expression of the elements of $V_{\text{scal},K}$ is not known in the general case in $\check{K}$. But, for each face or cell element $\gamma$, a projection $\pi_\gamma$ on $\mathbb{P}_1(\gamma)$ can however be evaluated setting

$$ \pi_{\gamma}(\varphi(\mathbf{x})) = \langle \nabla \varphi \rangle_{\gamma} \cdot (\mathbf{x} - \bar{\mathbf{x}}) + \bar{\varphi} $$

where $\langle \nabla \varphi \rangle_{\gamma} = \frac{1}{m_{\gamma}} \int_{\gamma} \nabla \varphi$, $\bar{u} = \frac{1}{\text{card}(\mathcal{V}(\gamma))} \sum_{v \in \mathcal{V}(\gamma)} u(\mathbf{x}_v)$, $u \in \{\varphi, \mathbf{x}\}$ and $\mathcal{V}(\gamma)$ is the set of all vertices of $\gamma$. In two dimensions, $V_{\text{scal},K}$ is defined by:

$$ V_{\text{scal},K} = \{ \varphi \in H^1(K) \mid \forall \sigma \in \partial K, \varphi|_{\sigma} \in \mathbb{P}_1(\sigma); \varphi|_{\partial K} \in \mathcal{C}^0(\partial K); \Delta \varphi = 0 \text{ in } K \}. $$

In three dimensions, we first define the approximate spaces for all faces $\sigma \in \mathcal{E}$ by

$$ V_{\text{scal},\sigma} = \left\{ \varphi \in H^1(\sigma) \mid \forall \gamma \in \partial\sigma, \varphi|_{\gamma} \in \mathbb{P}_1(\gamma); \varphi|_{\partial\sigma} \in \mathcal{C}^0(\partial\sigma); \Delta\varphi \in \mathbb{P}_1(\sigma); \int_{\sigma} \pi_{\sigma}(\varphi)q = \int_{\sigma} \varphi q, \forall q \in \mathbb{P}_1(\sigma) \right\} $$

and then, for all $K \in \mathcal{T}$, we take:

$$ V_{\text{scal},K} = \{ \varphi \in H^1(K) \mid \forall \sigma \in \partial K, \varphi|_{\sigma} \in V_{\text{scal},\sigma}; \Delta\varphi = 0 \text{ in } K \}. $$

Let $\mathcal{V}$ denote the set of the cell vertices and $\mathcal{V}_{D_d}$ its subset containing the vertices that are included in $\Gamma_{D_d}$. We then define:

$$ V_{\mathcal{D},0} = \{ \boldsymbol{\varphi} \in (H^1(\Omega))^d \mid \boldsymbol{\varphi} \in V_{\mathcal{D}}; \boldsymbol{\varphi}(\omega) = \mathbf{0} \text{ if } \omega \in \mathcal{V}_{D_d} \}. $$

Assuming we have an approximation for the pressure in $Q_{\mathcal{D}} = \mathbb{P}^0(\mathcal{T})$, which is the space of piecewise constant functions on each cell of the mesh, at a time-step $n \in \{0...N+1\}$, a discrete approximation of (9) in $V_{\mathcal{D}}$ can be computed by solving, for any test function $v_{\mathcal{D}}$ in this latter space, the following problem:

$$ \sum_{K \in \mathcal{T}} \left( a_{\mathcal{D}}^{K}(\mathbf{u}_{\mathcal{D}}^{n}, \mathbf{v}_{\mathcal{D}}) - \int_{K} \alpha p_{\mathcal{D}}^{n} \operatorname{div}(\mathbf{v}_{\mathcal{D}}) - \int_{K} \mathbf{f}^{n} \cdot \mathbf{v}_{\mathcal{D}} - \int_{\partial K \cap \Gamma_{N_d}} \mathbf{t}_{N}^{n} \cdot \mathbf{v}_{\mathcal{D}} \right) = 0 \quad (15) $$
When considering a sequence of admissible discretizations $(\mathcal{D}_n)_{n \in \mathbb{N}}$, we denote by $h_{\mathcal{D}_n}$ the size of the discretization defined by $h_{\mathcal{D}_n} \stackrel{\text{def}}{=} \sup_{K \in \mathcal{T}_n} \text{diam}(K)$. In that case, we suppose that there exist $0 < \rho_2, \rho_3, \rho_4, \rho_5 < +\infty$ independent of $n$ s.t.

$$ \begin{gathered} \min_{K \in \mathcal{T}_n, \sigma \in \mathcal{E}_{n,K}} \frac{d_{K,\sigma}}{\text{diam}(K)} \ge \rho_2, \quad \min_{\sigma \in \mathcal{E}_{n,\text{int}}, \mathcal{T}_{n,\sigma} = \{K,L\}} \frac{\min(d_{K,\sigma}, d_{L,\sigma})}{\max(d_{K,\sigma}, d_{L,\sigma})} \ge \rho_3, \quad \min_{K \in \mathcal{T}_n} \frac{\text{diam}(K)}{h_{\mathcal{D}_n}} \ge \rho_4, \\ \rho_5 \le \frac{d_{K,\sigma}}{d_{L,\sigma}} \le \frac{1}{\rho_5}, \quad \forall \sigma \in \mathcal{E}_{n,\text{int}}, \text{ such that } \mathcal{T}_{n,\sigma} = \{K,L\}. \end{gathered} $$

### 3.1. Virtual elements for the elasticity problem

In this work, the displacements are approximated using linear virtual elements. In that case, the approximation space $V_{\mathcal{D}}$ is given by

$$ V_{\mathcal{D}} = \bigcup_{K \in \mathcal{T}} V_K $$

where

$$ V_K = \{ \boldsymbol{\varphi} \in (H^1(K))^d \mid \boldsymbol{\varphi} \cdot \boldsymbol{e}_i \in V_{\text{scal},K}, \forall i = 1, \dots, d \}, $$

$e_i$ denoting one of the vectors of the canonical basis. The main idea of linear virtual elements is to define the spaces $V_{\text{scal},K}$ without giving an explicit definition of its basis functions even if, for simplicial cells, it reduces to $\mathbb{P}_1$ polynomials. Functions of $V_{\text{scal},K}$ are entirely defined by their values on the cell vertices of $K$. However, the exact expression of the elements of $V_{\text{scal},K}$ is not known in the general case in $\check{K}$. But, for each face or cell element $\gamma$, a projection $\pi_\gamma$ on $\mathbb{P}_1(\gamma)$ can however be evaluated setting

$$ \pi_{\gamma}(\varphi(\mathbf{x})) = \langle \nabla \varphi \rangle_{\gamma} \cdot (\mathbf{x} - \bar{\mathbf{x}}) + \bar{\varphi} $$

where $\langle \nabla \varphi \rangle_\gamma = \frac{1}{m_\gamma} \int_\gamma \nabla \varphi, \bar{u} = \frac{1}{\text{card}(\mathcal{V}(\gamma))} \sum_{v \in \mathcal{V}(\gamma)} u(\mathbf{x}_v), u \in \{\varphi, \mathbf{x}\}$ and $\mathcal{V}(\gamma)$ is the set of all vertices of $\gamma$. In two dimensions, $V_{\text{scal},K}$ is defined by:

$$ V_{\text{scal},K} = \{ \varphi \in H^1(K) \mid \forall \sigma \in \partial K, \varphi|_{\sigma} \in \mathbb{P}_1(\sigma); \varphi|_{\partial K} \in \mathcal{C}^0(\partial K); \Delta\varphi = 0 \text{ in } K \}. $$

In three dimensions, we first define the approximate spaces for all faces $\sigma \in \mathcal{E}$ by

$$ V_{\text{scal},\sigma} = \{ \varphi \in H^1(\sigma) \mid \forall \gamma \in \partial\sigma, \varphi|_{\gamma} \in \mathbb{P}_1(\gamma); \varphi|_{\partial\sigma} \in \mathcal{C}^0(\partial\sigma); \Delta\varphi \in \mathbb{P}_1(\sigma); \int_{\sigma} \pi_{\sigma}(\varphi) q = \int_{\sigma} \varphi q, \forall q \in \mathbb{P}_1(\sigma) \} $$

and then, for all $K \in \mathcal{T}$, we take:

$$ V_{\text{scal},K} = \{ \varphi \in H^1(K) \mid \forall \sigma \in \partial K, \varphi|_{\sigma} \in V_{\text{scal},\sigma}; \Delta\varphi = 0 \text{ in } K \}. $$

Let $\mathcal{V}$ denote the set of the cell vertices and $\mathcal{V}_{D_d}$ its subset containing the vertices that are included in $\Gamma_{D_d}$. We then define:

$$ V_{\mathcal{D},0} = \{ \boldsymbol{\varphi} \in (H^1(\Omega))^d \mid \boldsymbol{\varphi} \in V_{\mathcal{D}}; \boldsymbol{\varphi}(\omega) = \mathbf{0} \text{ if } \omega \in \mathcal{V}_{D_d} \}. $$

Assuming we have an approximation for the pressure in $Q_{\mathcal{D}} = \mathbb{P}^0(\mathcal{T})$, which is the space of piecewise constant functions on each cell of the mesh, at a time-step $n \in \{0...N+1\}$, a discrete approximation of (9) in $V_{\mathcal{D}}$ can be computed by solving, for any test function $v_{\mathcal{D}}$ in this latter space, the following problem:

$$ \sum_{K \in \mathcal{T}} \left( a_{\mathcal{D}}^{K}(\mathbf{u}_{\mathcal{D}}^{n}, \mathbf{v}_{\mathcal{D}}) - \int_{K} \alpha p_{\mathcal{D}}^{n} \operatorname{div}(\mathbf{v}_{\mathcal{D}}) - \int_{K} \mathbf{f}^{n} \cdot \mathbf{v}_{\mathcal{D}} - \int_{\partial K \cap \Gamma_{N_d}} \mathbf{t}_{N}^{n} \cdot \mathbf{v}_{\mathcal{D}} \right) = 0 \quad (15) $$
where

$$
\begin{align}
a_{\mathcal{D}}^{K}(\mathbf{w}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}}) &= a^{K}(\pi_{K}(\mathbf{w}_{\mathcal{D}}), \pi_{K}(\mathbf{v}_{\mathcal{D}})) + s_{K}(\mathbf{w}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}}) \nonumber \\
&= \int_{K} \bar{\sigma}^{e}(\pi_{K}(\mathbf{w}_{\mathcal{D}})) : \bar{\varepsilon}(\pi_{K}(\mathbf{v}_{\mathcal{D}})) + s_{K}(\mathbf{w}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}}), \tag{16}
\end{align}
$$

$s_K$ is a stabilisation term taken, following [10], as:

$$
s_K(\mathbf{w}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}}) = h_{\mathcal{D}}^{d-2} \|\bar{\mathcal{C}}\|_{\infty} \sum_{v \in \mathcal{V}(K)} ((\mathbf{w}_{\mathcal{D}} - \pi_K(\mathbf{w}_{\mathcal{D}})) \cdot (\mathbf{v}_{\mathcal{D}} - \pi_K(\mathbf{v}_{\mathcal{D}}))) (\mathbf{x}_v) \quad (17)
$$

with $\|\bar{\mathcal{C}}\|_{\infty} = \max_{i,j,k,l} \bar{\mathcal{C}}_{i,j,k,l}$. The bilinear form $a_{\mathcal{D}}$ being defined on each cell $K$ by (16)-(17), there exist two positive constants $\alpha_*, \alpha^*$ s.t.

$$
\alpha_* a^K(\mathbf{v}_\mathcal{D}, \mathbf{v}_\mathcal{D}) \le a_\mathcal{D}^K(\mathbf{v}_\mathcal{D}, \mathbf{v}_\mathcal{D}) \le \alpha^* a^K(\mathbf{v}_\mathcal{D}, \mathbf{v}_\mathcal{D}). \quad (18)
$$

Property (18) is, for instance, established in [11] with the assumption made on $r$ in Definition 1.

3.2. A generic finite-volume scheme for the flow problem

The discrete pressures are computed in the space $Q_{\mathcal{D}} = \mathbb{F}^0(\mathcal{T})$. In the following, for all $v \in Q_{\mathcal{D}}$ and $K \in \mathcal{T}$, we use the notation $v_K$ so that

$v_K = v(x), \forall x \in K.$

The space $Q_{\mathcal{D}}$ is equipped with the scalar product defined, for all $(v, w) \in (Q_{\mathcal{D}})^2$, by

$$
[v, w]_{\mathcal{D}} = \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K} \frac{m_{\sigma}}{d_{K, \sigma}} (\gamma_{\sigma} v - v_K) (\gamma_{\sigma} w - w_K) \quad (19)
$$

with

$$
\left\{
\begin{array}{ll}
\gamma_{\sigma} v = \frac{d_{L,\sigma} v_K + d_{K,\sigma} v_L}{d_{L,\sigma} + d_{K,\sigma}} & \text{if } \sigma \in \mathcal{E}_{\text{int}} \text{ with } \mathcal{T}_{\sigma} = \{K, L\}, \\
\gamma_{\sigma} v = 0 & \text{if } \sigma \in \mathcal{E}_{\text{ext}}.
\end{array}
\right.
\qquad (20)
$$

The associated norm is $\|v\|_{\mathcal{D}} = ([v, v]_{\mathcal{D}})^{1/2}$. A finite-volume approximation of (9) consists in choosing a test function $q$ s.t., for all $K \in \mathcal{T}$, $q = 1$ on $K$ and $q = 0$ elsewhere. Thus, using Stokes formula and denoting by $F_{K,\sigma}(u, v)$, a numerical conservative flux function with two entries that is meant to approximate the diffusive flux, $\int_{\sigma} \bar{\Lambda} \nabla p \cdot n_{K,\sigma}$, flowing out of $K$ through $\sigma \in \mathcal{E}_K$, we have

$$
m_K D^{n+1} \left( c_0 p_K + \alpha \int_K \operatorname{div}(\boldsymbol{u}_{\mathcal{D}}) \right) - \Delta t \sum_{\sigma \in \mathcal{E}_K} F_{K,\sigma} (p^{n+1}, p^{n+1}) = \Delta t \int_K s^{n+1}. \quad (21)
$$

**Example 2 (The non-linear two-point flux approximation).** As an example of flux $F_{K,\sigma}$, we here detail a non-linear two-point flux approximation (NLTPFA). This approximation is built to ensure the positivity of the solutions in case of elliptic problems. It satisfies Hypotheses 4 (see below) with exception of the coercivity property (36) which does not hold unconditionally. To obtain the positivity property:

* for each interior face $\sigma \in \mathcal{E}_{\text{int}}$ with $\mathcal{T}_{\sigma} = \{K, L\}$, the approximated flux $F_{K,\sigma}$ is written as a convex combination of two consistent linear fluxes $\tilde{F}_{K,\sigma}(p_{\mathcal{D}})$ and $\tilde{F}_{L,\sigma}(p_{\mathcal{D}})$ whose coefficients depend on $p_{\mathcal{D}}$ too. Thus, we set

$$
\begin{gather}
F_{K,\sigma}(p_{\mathcal{D}}, p_{\mathcal{D}}) = \mu_{K,\sigma}(p_{\mathcal{D}}) \tilde{F}_{K,\sigma}(p_{\mathcal{D}}) - \mu_{L,\sigma}(p_{\mathcal{D}}) \tilde{F}_{L,\sigma}(p_{\mathcal{D}}), \tag{22} \\
\text{with } \mu_{K,\sigma}(p_{\mathcal{D}}) \ge 0, \mu_{L,\sigma}(p_{\mathcal{D}}) \ge 0 \quad \text{and} \quad \mu_{K,\sigma}(p_{\mathcal{D}}) + \mu_{L,\sigma}(p_{\mathcal{D}}) = 1, \nonumber
\end{gather}
$$

* for each outer face $\sigma \in \mathcal{E}_{\text{ext}}$ with $\mathcal{T}_{\sigma} = \{K\}$, the approximated flux $F_{K,\sigma}$ can be chosen to coincide with $\tilde{F}_{K,\sigma}(p_{\mathcal{D}})$:

$$
F_{K,\sigma}(p_{\mathcal{D}}, p_{\mathcal{D}}) = \tilde{F}_{K,\sigma}(p_{\mathcal{D}}). \quad (23)
$$
where

$$ 
\begin{aligned}
a_{\mathcal{D}}^{K}(\mathbf{w}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}}) &= a^{K}(\pi_{K}(\mathbf{w}_{\mathcal{D}}), \pi_{K}(\mathbf{v}_{\mathcal{D}})) + s_{K}(\mathbf{w}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}}) \\
&= \int_{K} \bar{\sigma}^{e}(\pi_{K}(\mathbf{w}_{\mathcal{D}})) : \bar{\varepsilon}(\pi_{K}(\mathbf{v}_{\mathcal{D}})) + s_{K}(\mathbf{w}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}}),
\end{aligned} 
\quad (16) $$

$s_K$ is a stabilisation term taken, following [10], as:

$$ s_K(\mathbf{w}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}}) = h_{\mathcal{D}}^{d-2} \|\bar{\mathcal{C}}\|_{\infty} \sum_{v \in \mathcal{V}(K)} ((\mathbf{w}_{\mathcal{D}} - \pi_K(\mathbf{w}_{\mathcal{D}})) \cdot (\mathbf{v}_{\mathcal{D}} - \pi_K(\mathbf{v}_{\mathcal{D}}))) (\mathbf{x}_v) \quad (17) $$

with $\|\bar{\mathcal{C}}\|_{\infty} = \max_{i,j,k,l} \bar{C}_{i,j,k,l}$. The bilinear form $a_{\mathcal{D}}$ being defined on each cell $K$ by (16)-(17), there exist two positive constants $\alpha_\star$, $\alpha^*$ s.t.

$$ \alpha_{\star} a^{K}(\mathbf{v}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}}) \leq a_{\mathcal{D}}^{K}(\mathbf{v}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}}) \leq \alpha^{*} a^{K}(\mathbf{v}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}}). \quad (18) $$

Property (18) is, for instance, established in [11] with the assumption made on $r$ in Definition 1.

## 3.2. A generic finite-volume scheme for the flow problem

The discrete pressures are computed in the space $Q_{\mathcal{D}} = \mathbb{F}^0(\mathcal{T})$. In the following, for all $v \in Q_{\mathcal{D}}$ and $K \in \mathcal{T}$, we use the notation $v_K$ so that

$$ v_K = v(x), \forall x \in K. $$

The space $Q_{\mathcal{D}}$ is equipped with the scalar product defined, for all $(v, w) \in (Q_{\mathcal{D}})^2$, by

$$ [v, w]_{\mathcal{D}} = \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K} \frac{m_{\sigma}}{d_{K, \sigma}} (\gamma_{\sigma} v - v_K) (\gamma_{\sigma} w - w_K) \quad (19) $$

with

$$ \left\{ \begin{array}{ll} \gamma_{\sigma} v = \frac{d_{L, \sigma} v_K + d_{K, \sigma} v_L}{d_{L, \sigma} + d_{K, \sigma}} & \text{if } \sigma \in \mathcal{E}_{\text{int}} \text{ with } \mathcal{T}_{\sigma} = \{K, L\}, \\ \gamma_{\sigma} v = 0 & \text{if } \sigma \in \mathcal{E}_{\text{ext}}. \end{array} \right. \quad (20) $$

The associated norm is $\|v\|_{\mathcal{D}} = ([v, v]_{\mathcal{D}})^{1/2}$. A finite-volume approximation of (9) consists in choosing a test function $q$ s.t., for all $K \in \mathcal{T}$, $q = 1$ on $K$ and $q = 0$ elsewhere. Thus, using Stokes formula and denoting by $F_{K,\sigma}(u, v)$, a numerical conservative flux function with two entries that is meant to approximate the diffusive flux, $\int_{\sigma} \bar{\Lambda} \nabla p \cdot n_{K,\sigma}$, flowing out of $K$ through $\sigma \in \mathcal{E}_K$, we have

$$ m_K D^{n+1} \left( c_0 p_K + \alpha \int_K \mathrm{div}(\boldsymbol{u}_{\mathcal{D}}) \right) - \Delta t \sum_{\sigma \in \mathcal{E}_K} F_{K,\sigma} (p^{n+1}, p^{n+1}) = \Delta t \int_K s^{n+1}. \quad (21) $$

**Example 2 (The non-linear two-point flux approximation).** As an example of flux $F_{K,\sigma}$, we here detail a non-linear two-point flux approximation (NLTPFA). This approximation is built to ensure the positivity of the solutions in case of elliptic problems. It satisfies Hypotheses 4 (see below) with exception of the coercivity property (36) which does not hold unconditionally. To obtain the positivity property:

* for each interior face $\sigma \in \mathcal{E}_{\text{int}}$ with $\mathcal{T}_{\sigma} = \{K, L\}$, the approximated flux $F_{K,\sigma}$ is written as a convex combination of two consistent linear fluxes $\tilde{F}_{K,\sigma}(p_{\mathcal{D}})$ and $\tilde{F}_{L,\sigma}(p_{\mathcal{D}})$ whose coefficients depend on $p_{\mathcal{D}}$ too. Thus, we set

$$ F_{K,\sigma}(p_{\mathcal{D}}, p_{\mathcal{D}}) = \mu_{K,\sigma}(p_{\mathcal{D}}) \tilde{F}_{K,\sigma}(p_{\mathcal{D}}) - \mu_{L,\sigma}(p_{\mathcal{D}}) \tilde{F}_{L,\sigma}(p_{\mathcal{D}}), \quad (22) $$
with $\mu_{K,\sigma}(p_{\mathcal{D}}) \ge 0, \mu_{L,\sigma}(p_{\mathcal{D}}) \ge 0$ and $\mu_{K,\sigma}(p_{\mathcal{D}}) + \mu_{L,\sigma}(p_{\mathcal{D}}) = 1,$

* for each outer face $\sigma \in \mathcal{E}_{\text{ext}}$ with $\mathcal{T}_{\sigma} = \{K\}$, the approximated flux $F_{K,\sigma}$ can be chosen to coincide with $\tilde{F}_{K,\sigma}(p_{\mathcal{D}})$:

$$ F_{K,\sigma}(p_{\mathcal{D}}, p_{\mathcal{D}}) = \tilde{F}_{K,\sigma}(p_{\mathcal{D}}). \quad (23) $$
The linear fluxes are defined by

$$ \tilde{F}_{K,\sigma}(u) = m_{\sigma} \sum_{\sigma' \in \mathcal{S}_{K,\sigma}} \alpha_{K,\sigma\sigma'} (I_{\sigma'} u - u_K). \quad (24) $$

In (24), $\mathcal{S}_{K,\sigma}$ denotes the stencil of the pair $(K, \sigma)$ which corresponds to a subset of faces included in $\partial K$ where the coefficients $\alpha_{K,\sigma\sigma'}$ are nonzero. These coefficients are computed in order to verify the conormal decomposition:

$$ \bar{\Lambda}_K \mathbf{n}_{K,\sigma} = \sum_{\sigma' \in \mathcal{S}_{K,\sigma}} \alpha_{K,\sigma\sigma'} (\mathbf{x}_{\sigma'} - \mathbf{x}_K), \quad (25) $$

where $\mathbf{x}_{\sigma'}$ is a point depending on the face $\sigma$ and on the trace reconstruction operator that is considered (see below). The coefficients of the conormal decomposition are non-negative if the geometry of the cell and the mobility tensor enable it. References [5, 12] have proposed to compute this decomposition by means of an optimization problem.

The trace reconstruction operator $I \in \mathcal{L}(Q_{\mathcal{D}}; \mathbb{P}^0(\mathcal{E}))$ is defined so that,

* for each face $\sigma \in \mathcal{E}_{\text{int}}$,

$$ I_{\sigma} u = \sum_{M \in \mathcal{I}_{\sigma}} \omega_{M,\sigma} u_{M}, \quad \sum_{M \in \mathcal{I}_{\sigma}} \omega_{M,\sigma} = 1, \quad \omega_{M,\sigma} \ge 0 \text{ and } \forall \varphi \in \mathcal{D}, \quad |I_{\sigma} \varphi_{\mathcal{D}} - \varphi(\mathbf{x}_{\sigma})| \le C_{\mathcal{D}} h_{\mathcal{D}}^2 \quad (26) $$

where $\mathcal{I}_{\sigma}$ is a cell subset of $\mathcal{T}$, $\mathcal{D}$ is a dense subspace of $H^1(\Omega)$ s.t. $\mathcal{D} \subset C_0(\bar{\Omega})$, $\varphi_{\mathcal{D}} \in Q_{\mathcal{D}}$ with $(\varphi_{\mathcal{D}})_K = \varphi(\mathbf{x}_K)$ for all $K \in \mathcal{T}$, and $C_{\mathcal{D}}$ is a real number only depending on the mesh regularity parameters $\varrho_2, \varrho_3, \varrho_4$ (more details on this operator are given in [5]),

* for each face $\sigma \in \mathcal{E}_{\text{ext}}$,

- if $\sigma \subset \Gamma_{D_p}$ then $I_{\sigma} u = 0$

- if $\sigma \subset \Gamma_{N_p}$ then $I_{\sigma} u = u_{\sigma}$, where $u_{\sigma}$ is a face unknown added to ensure the Neumann boundary condition.

To derive a two-point approximation, the different terms in (22) are reordered so as to obtain

$$ F_{K,\sigma}(u) = t_{L,\sigma}(u)u_L - t_{K,\sigma}(u)u_K - \underbrace{\frac{\mu_{L,\sigma}(u)\lambda_{L,\sigma}(u) - \mu_{K,\sigma}(u)\lambda_{K,\sigma}(u)}{\text{def}=R_{K,\sigma}(u)}}_{\text{def}=R_{K,\sigma}(u)} \quad (27) $$

with the transmissibilities

$$ \begin{aligned} t_{K,\sigma}(u) &= m_{\sigma} \left( \mu_{K,\sigma}(u) \sum_{\sigma' \in \mathcal{S}_{K,\sigma}} \sum_{M \in \{\mathcal{I}_{\sigma'} \setminus \{K\}\}} \alpha_{K,\sigma\sigma'} \omega_{M,\sigma'} + \mu_{L,\sigma}(u) \sum_{\sigma' \in \mathcal{S}_{L,\sigma}} \sum_{M \in \{\mathcal{I}_{\sigma'} \cap \{K\}\}} \alpha_{L,\sigma\sigma'} \omega_{M,\sigma'} \right), \quad (28) \\ t_{L,\sigma}(u) &= m_{\sigma} \left( \mu_{L,\sigma}(u) \sum_{\sigma' \in \mathcal{S}_{L,\sigma}} \sum_{M \in \{\mathcal{I}_{\sigma'} \setminus \{L\}\}} \alpha_{L,\sigma\sigma'} \omega_{M,\sigma'} + \mu_{K,\sigma}(u) \sum_{\sigma' \in \mathcal{S}_{K,\sigma}} \sum_{M \in \{\mathcal{I}_{\sigma'} \cap \{L\}\}} \alpha_{K,\sigma\sigma'} \omega_{M,\sigma'} \right), \end{aligned} $$

and

$$ \begin{aligned} \lambda_{K,\sigma}(v) &= m_{\sigma} \sum_{\sigma' \in \mathcal{S}_{K,\sigma}} \sum_{M \in \{\mathcal{I}_{\sigma'} \setminus \{K,L\}\}} \alpha_{K,\sigma\sigma'} \omega_{M,\sigma'} v_M, \\ \lambda_{L,\sigma}(v) &= m_{\sigma} \sum_{\sigma' \in \mathcal{S}_{L,\sigma}} \sum_{M \in \{\mathcal{I}_{\sigma'} \setminus \{K,L\}\}} \alpha_{L,\sigma\sigma'} \omega_{M,\sigma'} v_M. \end{aligned} \quad (29) $$

The idea of the NLTPFA scheme is to choose the weights such that $R_{K,\sigma}(u) = 0$. From a numerical point of view, it is sufficient that $|R_{K,\sigma}(u)| \le \epsilon$. Under the assumption that $\lambda_{K,\sigma}, \lambda_{L,\sigma} \ge 0$, this can be ensured, by taking

$$ \mu_{K,\sigma}(u) = \frac{|\lambda_{L,\sigma}(u)| + \epsilon}{|\lambda_{K,\sigma}(u)| + |\lambda_{L,\sigma}(u)| + 2\epsilon}, \quad \mu_{L,\sigma}(u) = \frac{|\lambda_{K,\sigma}(u)| + \epsilon}{|\lambda_{K,\sigma}(u)| + |\lambda_{L,\sigma}(u)| + 2\epsilon}. \quad (30) $$
The linear fluxes are defined by

$$ \tilde{F}_{K,\sigma}(u) = m_{\sigma} \sum_{\sigma' \in \mathcal{S}_{K,\sigma}} \alpha_{K,\sigma\sigma'} (I_{\sigma'} u - u_K). \quad (24) $$

In (24), $\mathcal{S}_{K,\sigma}$ denotes the stencil of the pair $(K, \sigma)$ which corresponds to a subset of faces included in $\partial K$ where the coefficients $\alpha_{K,\sigma\sigma'}$ are nonzero. These coefficients are computed in order to verify the conormal decomposition:

$$ \bar{\Lambda}_K \mathbf{n}_{K,\sigma} = \sum_{\sigma' \in \mathcal{S}_{K,\sigma}} \alpha_{K,\sigma\sigma'} (\mathbf{x}_{\sigma'} - \mathbf{x}_K), \quad (25) $$

where $\mathbf{x}_{\sigma'}$ is a point depending on the face $\sigma$ and on the trace reconstruction operator that is considered (see below). The coefficients of the conormal decomposition are non-negative if the geometry of the cell and the mobility tensor enable it. References [5, 12] have proposed to compute this decomposition by means of an optimization problem.

The trace reconstruction operator $I \in \mathcal{L}(Q_{\mathcal{D}}; \mathbb{P}^0(\mathcal{E}))$ is defined so that,

* for each face $\sigma \in \mathcal{E}_{\text{int}}$,

$$ I_{\sigma} u = \sum_{M \in \mathcal{I}_{\sigma}} \omega_{M,\sigma} u_{M}, \quad \sum_{M \in \mathcal{I}_{\sigma}} \omega_{M,\sigma} = 1, \quad \omega_{M,\sigma} \ge 0 \text{ and } \forall \varphi \in \mathcal{D}, \quad |I_{\sigma} \varphi_{\mathcal{D}} - \varphi(\mathbf{x}_{\sigma})| \le C_{\mathcal{D}} h_{\mathcal{D}}^2 \quad (26) $$

where $\mathcal{I}_{\sigma}$ is a cell subset of $\mathcal{T}$, $\mathcal{D}$ is a dense subspace of $H^1(\Omega)$ s.t. $\mathcal{D} \subset C_0(\bar{\Omega})$, $\varphi_{\mathcal{D}} \in Q_{\mathcal{D}}$ with $(\varphi_{\mathcal{D}})_K = \varphi(\mathbf{x}_K)$ for all $K \in \mathcal{T}$, and $C_{\mathcal{D}}$ is a real number only depending on the mesh regularity parameters $\varrho_2, \varrho_3, \varrho_4$ (more details on this operator are given in [5]),

* for each face $\sigma \in \mathcal{E}_{\text{ext}}$,

- if $\sigma \subset \Gamma_{D_p}$ then $I_{\sigma} u = 0$

- if $\sigma \subset \Gamma_{N_p}$ then $I_{\sigma} u = u_{\sigma}$, where $u_{\sigma}$ is a face unknown added to ensure the Neumann boundary condition.

To derive a two-point approximation, the different terms in (22) are reordered so as to obtain

$$ F_{K,\sigma}(u) = t_{L,\sigma}(u)u_L - t_{K,\sigma}(u)u_K - \underbrace{\frac{\mu_{L,\sigma}(u)\lambda_{L,\sigma}(u) - \mu_{K,\sigma}(u)\lambda_{K,\sigma}(u)}{\text{def}=R_{K,\sigma}(u)}}_{\text{def}=R_{K,\sigma}(u)} \quad (27) $$

with the transmissibilities

$$ \begin{aligned} t_{K,\sigma}(u) &= m_{\sigma} \left( \mu_{K,\sigma}(u) \sum_{\sigma' \in \mathcal{S}_{K,\sigma}} \sum_{M \in \{\mathcal{I}_{\sigma'} \setminus \{K\}\}} \alpha_{K,\sigma\sigma'} \omega_{M,\sigma'} + \mu_{L,\sigma}(u) \sum_{\sigma' \in \mathcal{S}_{L,\sigma}} \sum_{M \in \{\mathcal{I}_{\sigma'} \cap \{K\}\}} \alpha_{L,\sigma\sigma'} \omega_{M,\sigma'} \right), \\ t_{L,\sigma}(u) &= m_{\sigma} \left( \mu_{L,\sigma}(u) \sum_{\sigma' \in \mathcal{S}_{L,\sigma}} \sum_{M \in \{\mathcal{I}_{\sigma'} \setminus \{L\}\}} \alpha_{L,\sigma\sigma'} \omega_{M,\sigma'} + \mu_{K,\sigma}(u) \sum_{\sigma' \in \mathcal{S}_{K,\sigma}} \sum_{M \in \{\mathcal{I}_{\sigma'} \cap \{L\}\}} \alpha_{K,\sigma\sigma'} \omega_{M,\sigma'} \right), \end{aligned} \quad (28) $$

and

$$ \begin{aligned} \lambda_{K,\sigma}(v) &= m_{\sigma} \sum_{\sigma' \in \mathcal{S}_{K,\sigma}} \sum_{M \in \{\mathcal{I}_{\sigma'} \setminus \{K,L\}\}} \alpha_{K,\sigma\sigma'} \omega_{M,\sigma'} v_M, \\ \lambda_{L,\sigma}(v) &= m_{\sigma} \sum_{\sigma' \in \mathcal{S}_{L,\sigma}} \sum_{M \in \{\mathcal{I}_{\sigma'} \setminus \{K,L\}\}} \alpha_{L,\sigma\sigma'} \omega_{M,\sigma'} v_M. \end{aligned} \quad (29) $$

The idea of the NLTPFA scheme is to choose the weights such that $R_{K,\sigma}(u) = 0$. From a numerical point of view, it is sufficient that $|R_{K,\sigma}(u)| \le \epsilon$. Under the assumption that $\lambda_{K,\sigma}, \lambda_{L,\sigma} \ge 0$, this can be ensured, by taking

$$ \mu_{K,\sigma}(u) = \frac{|\lambda_{L,\sigma}(u)| + \epsilon}{|\lambda_{K,\sigma}(u)| + |\lambda_{L,\sigma}(u)| + 2\epsilon}, \quad \mu_{L,\sigma}(u) = \frac{|\lambda_{K,\sigma}(u)| + \epsilon}{|\lambda_{K,\sigma}(u)| + |\lambda_{L,\sigma}(u)| + 2\epsilon}. \quad (30) $$
With the weights (30), the residual term is given as

$$R_{K,\sigma}(u) = \epsilon \frac{\lambda_{L,\sigma}(u) - \lambda_{K,\sigma}(u)}{|\lambda_{K,\sigma}(u)| + |\lambda_{L,\sigma}(u)| + 2\epsilon},$$ 
$$|R_{K,\sigma}(u)| \le \epsilon. \tag{31}$$

for which it holds that

**Example 3 (The SUSHI flux approximation).** As a second example of flux $F_{K,\sigma}$, we here detail the Scheme Using Stabilisation and Hybrid Interfaces (SUSHI) introduced in [13]. This approximation is built to satisfy Hypotheses 4, in particular the coercivity property (36) unconditionally holds. In that case, for each face $\sigma \in \mathcal{E}$, the approximated flux $F_{K,\sigma}$, we take

$$F_{K,\sigma}(u, u) = \sum_{\sigma' \in \mathcal{E}_K} \alpha_{K,\sigma\sigma'} (I_{\sigma'} u - u_K), \tag{32}$$

where $I \in \mathcal{L}(Q_{\mathcal{D}}; \mathbb{P}^0(\mathcal{E}))$ is a trace reconstruction operator. The coefficients $\alpha_{K,\sigma\sigma'}$ are given by

$$\alpha_{K,\sigma\sigma'} = \sum_{\sigma'' \in \mathcal{E}_K} \mathbf{y}^{\sigma''\sigma} \cdot \bar{\Lambda}_{K,\sigma\sigma''} \mathbf{y}^{\sigma''\sigma'}, \quad \text{and} \quad \bar{\Lambda}_{K,\sigma\sigma''} = \int_{D_{K,\sigma''}} \bar{\Lambda}(x) dx \tag{33}$$

where $D_{K,\sigma''}$ is the cone with vertex $\mathbf{x}_K$ and basis $\sigma''$ and

$$\mathbf{y}^{\sigma\sigma'} = \begin{cases} \frac{m_{\sigma}}{m_K} \mathbf{n}_{K,\sigma} + \frac{\sqrt{d}}{d_{K,\sigma}} \left( 1 - \frac{m_{\sigma}}{m_K} \mathbf{n}_{K,\sigma} \cdot (\mathbf{x}_{\sigma} - \mathbf{x}_K) \right) \mathbf{n}_{K,\sigma} & \text{if } \sigma = \sigma', \\ \frac{m_{\sigma'}}{m_K} \mathbf{n}_{K,\sigma'} - \frac{\sqrt{d}}{d_{K,\sigma}} \mathbf{n}_{K,\sigma'} \cdot (\mathbf{x}_{\sigma} - \mathbf{x}_K) \mathbf{n}_{K,\sigma} & \text{otherwise,} \end{cases} \tag{34}$$

where, for any $\sigma \in \mathcal{E}$, $\mathbf{x}_{\sigma}$ is the barycenter of $\sigma$. The trace reconstruction operator $I_{\sigma}$ is s.t.,

* for each face $\sigma \in \mathcal{E}_{\text{int}}$, with $\mathcal{T}_{\sigma} = \{K, L\}$, $I_{\sigma} u$ ensures the conservation of the fluxes, i.e.

$$F_{K,\sigma}(u, u) + F_{L,\sigma}(u, u) = 0, \tag{35}$$

* for each face $\sigma \in \mathcal{E}_{\text{ext}}$,

- if $\sigma \subset \Gamma_{D_p}$ then $I_{\sigma} u = 0$,

- if $\sigma \subset \Gamma_{N_p}$ then $I_{\sigma} u = u_{\sigma}$, where $u_{\sigma}$ is a face unknown added to ensure the Neumann boundary condition.

## 4. Numerical analysis of the coupled scheme

The aim of this section is to prove the existence of discrete solutions to the coupled discrete system (15)-(21) and to establish an a priori error estimate. These results are shown with the following assumptions.

**Hypotheses 4.** Let $\mathcal{D}$ an admissible discretization in the sense of Definition 1,

* for any $p \in Q_{\mathcal{D}}$, $K \in \mathcal{T}$, $\sigma \in \mathcal{E}$, $p \mapsto F_{K,\sigma}(p, \cdot)$ is a linear form;

* for any $\sigma \in \mathcal{E}_{\text{int}}$, $\mathcal{T}_{\sigma} = \{K, L\}$, $F_{K,\sigma}(p, q) + F_{L,\sigma}(p, q) = 0$ for all $p, q \in Q_{\mathcal{D}} \times Q_{\mathcal{D}}$;

* for all $\varphi \in C_c^{\infty}(\Omega)$, $K \in \mathcal{T}$ and $\sigma \in \mathcal{E}$, $\sigma \subset \partial K$,

$$\max_{q \in Q_{\mathcal{D}}} \left| F_{K,\sigma}(q, \varphi_{\mathcal{D}}) - \frac{m_{\sigma}}{m_K} \int_K \bar{\Lambda} \nabla \varphi \cdot \mathbf{n}_{K,\sigma} \right| \le C m_{\sigma} h_{\mathcal{D}}$$ 

where $\varphi_{\mathcal{D}} \in Q_{\mathcal{D}}$ is defined, for all $K \in \mathcal{T}$, s.t. $\varphi_{\mathcal{D}}(\mathbf{x}) := \frac{1}{m_K} \int_K \varphi(\mathbf{y}) d\mathbf{y}$; 

* there exists a real number $C_{\mathcal{D}} > 0$ depending on $\varrho_2, \varrho_3, \varrho_4$ s.t. for all $p, q \in Q_{\mathcal{D}} \times Q_{\mathcal{D}}$, 
$$ - \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K} F_{K,\sigma}(p, q) q_K \ge C_{\mathcal{D}} \|q\|_{\mathcal{D}}^2. \tag{36} $$
With the weights (30), the residual term is given as

$$R_{K,\sigma}(u) = \epsilon \frac{\lambda_{L,\sigma}(u) - \lambda_{K,\sigma}(u)}{|\lambda_{K,\sigma}(u)| + |\lambda_{L,\sigma}(u)| + 2\epsilon},$$ 
$$|R_{K,\sigma}(u)| \le \epsilon. \tag{31}$$

for which it holds that

**Example 3 (The SUSHI flux approximation).** As a second example of flux $F_{K,\sigma}$, we here detail the Scheme Using Stabilisation and Hybrid Interfaces (SUSHI) introduced in [13]. This approximation is built to satisfy Hypotheses 4, in particular the coercivity property (36) unconditionally holds. In that case, for each face $\sigma \in \mathcal{E}$, the approximated flux $F_{K,\sigma}$, we take

$$F_{K,\sigma}(u, u) = \sum_{\sigma' \in \mathcal{E}_K} \alpha_{K,\sigma\sigma'} (I_{\sigma'} u - u_K), \tag{32}$$

where $I \in \mathcal{L}(Q_{\mathcal{D}}; \mathbb{P}^0(\mathcal{E}))$ is a trace reconstruction operator. The coefficients $\alpha_{K,\sigma\sigma'}$ are given by

$$\alpha_{K,\sigma\sigma'} = \sum_{\sigma'' \in \mathcal{E}_K} \mathbf{y}^{\sigma''\sigma} \cdot \bar{\Lambda}_{K,\sigma\sigma''} \mathbf{y}^{\sigma''\sigma'}, \quad \text{and} \quad \bar{\Lambda}_{K,\sigma\sigma''} = \int_{D_{K,\sigma''}} \bar{\Lambda}(x) dx \tag{33}$$

where $D_{K,\sigma''}$ is the cone with vertex $\mathbf{x}_K$ and basis $\sigma''$ and

$$\mathbf{y}^{\sigma\sigma'} = \begin{cases} \frac{m_{\sigma}}{m_{K}} \mathbf{n}_{K,\sigma} + \frac{\sqrt{d}}{d_{K,\sigma}} \left(1 - \frac{m_{\sigma}}{m_{K}} \mathbf{n}_{K,\sigma} \cdot (\mathbf{x}_{\sigma} - \mathbf{x}_{K})\right) \mathbf{n}_{K,\sigma} & \text{if } \sigma = \sigma', \\ \frac{m_{\sigma'}}{m_{K}} \mathbf{n}_{K,\sigma'} - \frac{\sqrt{d}}{d_{K,\sigma}} \frac{m_{\sigma'}}{m_{K}} \mathbf{n}_{K,\sigma'} \cdot (\mathbf{x}_{\sigma} - \mathbf{x}_{K}) \mathbf{n}_{K,\sigma} & \text{otherwise,} \end{cases} \tag{34}$$

where, for any $\sigma \in \mathcal{E}$, $\mathbf{x}_{\sigma}$ is the barycenter of $\sigma$. The trace reconstruction operator $I_{\sigma}$ is s.t.,

* for each face $\sigma \in \mathcal{E}_{\text{int}}$, with $\mathcal{T}_{\sigma} = \{K, L\}$, $I_{\sigma} u$ ensures the conservation of the fluxes, i.e.

$$F_{K,\sigma}(u, u) + F_{L,\sigma}(u, u) = 0, \tag{35}$$

* for each face $\sigma \in \mathcal{E}_{\text{ext}}$,

- if $\sigma \subset \Gamma_{D_p}$ then $I_{\sigma} u = 0$,

- if $\sigma \subset \Gamma_{N_p}$ then $I_{\sigma} u = u_{\sigma}$, where $u_{\sigma}$ is a face unknown added to ensure the Neumann boundary condition.

## 4. Numerical analysis of the coupled scheme

The aim of this section is to prove the existence of discrete solutions to the coupled discrete system (15)-(21) and to establish an a priori error estimate. These results are shown with the following assumptions.

**Hypotheses 4.** Let $\mathcal{D}$ an admissible discretization in the sense of Definition 1,

* for any $p \in Q_{\mathcal{D}}$, $K \in \mathcal{T}$, $\sigma \in \mathcal{E}$, $p \mapsto F_{K,\sigma}(p, \cdot)$ is a linear form;

* for any $\sigma \in \mathcal{E}_{\text{int}}$, $\mathcal{T}_{\sigma} = \{K, L\}$, $F_{K,\sigma}(p, q) + F_{L,\sigma}(p, q) = 0$ for all $p, q \in Q_{\mathcal{D}} \times Q_{\mathcal{D}}$;

* for all $\varphi \in C_c^{\infty}(\Omega)$, $K \in \mathcal{T}$ and $\sigma \in \mathcal{E}$, $\sigma \subset \partial K$,

$$\max_{q \in Q_{\mathcal{D}}} \left| F_{K,\sigma}(q, \varphi_{\mathcal{D}}) - \frac{m_{\sigma}}{m_{K}} \int_{K} \bar{\Lambda} \nabla \varphi \cdot \mathbf{n}_{K,\sigma} \right| \le C m_{\sigma} h_{\mathcal{D}}$$ 

where $\varphi_{\mathcal{D}} \in Q_{\mathcal{D}}$ is defined, for all $K \in \mathcal{T}$, s.t. $\varphi_{\mathcal{D}}(\mathbf{x}) := \frac{1}{m_{K}} \int_{K} \varphi(\mathbf{y}) d\mathbf{y}$; 

* there exists a real number $C_{\mathcal{D}} > 0$ depending on $\varrho_2, \varrho_3, \varrho_4$ s.t. for all $p, q \in Q_{\mathcal{D}} \times Q_{\mathcal{D}}$, 
$$- \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K} F_{K,\sigma}(p, q) q_K \ge C_{\mathcal{D}} \|q\|_{\mathcal{D}}^2. \tag{36}$$
**Remark 5.** The coercivity property (36) cannot be easily proved for non-linear fluxes. As mentioned before, it holds for the SUSHI scheme but also for the linear two-point flux approximation that results from (25)–(30) when used on admissible meshes in the sense of [14, Definition 11.1] with a harmonic average of both cell values for the trace reconstruction. In Section 5, the value of $C_\mathcal{D}$ has been evaluated for different sequences of grids.

In order to prove the existence of solutions in Proposition 7, we first state and prove the following a priori estimate.

**Proposition 6 (A priori estimate).** Let $\mathcal{D}$ be an admissible discretization of $\Omega$ in the sense of Definition 1. Under Hypotheses 4 and for the sequence of solutions $(\mathbf{u}_{\mathcal{D}}^n, p_{\mathcal{D}}^n)_{n=1,...,N+1}$ to (15)-(21), there exists a constant $\tilde{B}$ s.t.

$$ \|\mathbf{u}_{\mathcal{D}}^{N+1}\|_{H^1(\Omega)}^2 + \|p_{\mathcal{D}}^{N+1}\|_{L^2(\Omega)}^2 + \sum_{n=1}^{N+1} \Delta t \|p_{\mathcal{D}}^n\|_{\mathcal{D}}^2 \le \tilde{B} \exp(T). \quad (37) $$

**Proof.** Taking $v_{\mathcal{D}} = D^{n+1}(\mathbf{u}_{\mathcal{D}})$ in (15), multiplying (21) by $p_K^{n+1}$, summing over $K \in \mathcal{F}$ and $n \in [0, N]$, we obtain

$$ \sum_{n=0}^{N} \sum_{K \in \mathcal{F}} \left( T_{a,K}^{n+1} + T_{c,K}^{n+1} + T_{F,K}^{n+1} \right) \le \sum_{n=0}^{N} \sum_{K \in \mathcal{F}} \left( T_{f,K}^{n+1} + T_{t,K}^{n+1} + T_{s,K}^{n+1} \right) \quad (38) $$

with

$$ \begin{aligned} T_{a,K}^{n+1} &= a_{\mathcal{D}}^{K}(\mathbf{u}_{\mathcal{D}}^{n+1}, D^{n+1}(\mathbf{u}_{\mathcal{D}})), \\ T_{c,K}^{n+1} &= m_K c_0 D^{n+1}(p_K) p_K^{n+1}, \\ T_{F,K}^{n+1} &= -\Delta t \sum_{\sigma \in \partial K} F_{K,\sigma}(p^{n+1}, p^{n+1}) p_K^{n+1}, \end{aligned} $$

$$ \begin{aligned} T_{f,K}^{n+1} &= \int_K \boldsymbol{f}^{n+1} \cdot D^{n+1}(\boldsymbol{u}_{\mathcal{D}}), \\ T_{t,K}^{n+1} &= \int_{\partial K \cap \Gamma_{N_d}} \boldsymbol{t}^{n+1} \cdot D^{n+1}(\boldsymbol{u}_{\mathcal{D}}), \\ T_{s,K}^{n+1} &= \int_{t^n}^{t^{n+1}} \int_K s p_K^{n+1}. \end{aligned} $$

A reordering of the sums in time (see the fourth tool of [3, section A.1]) yields:

$$ \begin{aligned} \sum_{n=0}^{N} \sum_{K \in \mathcal{F}} T_{f,K}^{n+1} &= \sum_{K \in \mathcal{F}} \left( \int_K \boldsymbol{f}^{N+1} \cdot \boldsymbol{u}_{\mathcal{D}}^{N+1} - \int_K \boldsymbol{f}^1 \cdot \boldsymbol{u}_{\mathcal{D}}^0 - \sum_{n=1}^{N} \int_K D^{n+1}(\boldsymbol{f}) \cdot \boldsymbol{u}_{\mathcal{D}}^n \right), \\ \sum_{n=0}^{N} \sum_{K \in \mathcal{F}} T_{t,K}^{n+1} &= \sum_{K \in \mathcal{F}} \left( \int_{\partial K \cap \Gamma_{N_d}} \boldsymbol{t}^{N+1} \cdot \boldsymbol{u}_{\mathcal{D}}^{N+1} - \int_{\partial K \cap \Gamma_{N_d}} \boldsymbol{t}^1 \cdot \boldsymbol{u}_{\mathcal{D}}^0 - \sum_{n=1}^{N} \int_{\partial K \cap \Gamma_{N_d}} D^{n+1}(\boldsymbol{t}) \cdot \boldsymbol{u}_{\mathcal{D}}^n \right). \end{aligned} $$

Several applications of Cauchy-Schwartz and Young inequalities with parameters $(\varepsilon_i)_{i=1...6} \in (0, 1)$ thus lead to

$$ \left| \sum_{n=0}^{N} \sum_{K \in \mathcal{F}} T_{f,K}^{n+1} \right| \le \frac{1}{2\varepsilon_1} \|\boldsymbol{f}^{N+1}\|_{L^2(\Omega)}^2 + \frac{\varepsilon_1}{2} \|\boldsymbol{u}_{\mathcal{D}}^{N+1}\|_{H^1(\Omega)}^2 + \frac{1}{2\varepsilon_2} \|\boldsymbol{f}^1\|_{L^2(\Omega)}^2 + \frac{\varepsilon_2}{2} \|\boldsymbol{u}_{\mathcal{D}}^0\|_{H^1(\Omega)}^2 \\ + \sum_{n=1}^{N} \left( \frac{1}{2\varepsilon_3} \int_{t^n}^{t^{n+1}} \|\boldsymbol{f}'\|_{L^2(\Omega)}^2 + \frac{\Delta t \varepsilon_3}{2} \|\boldsymbol{u}_{\mathcal{D}}^n\|_{H^1(\Omega)}^2 \right) \quad (39) $$

$$ \left| \sum_{n=0}^{N} \sum_{K \in \mathcal{F}} T_{t,K}^{n+1} \right| \le \frac{1}{2\varepsilon_4} \|\boldsymbol{t}^{N+1}\|_{L^2(\Gamma_{N_d})}^2 + \frac{\varepsilon_4}{2} \|\boldsymbol{u}_{\mathcal{D}}^{N+1}\|_{L^2(\Gamma_{N_d})}^2 + \frac{1}{2\varepsilon_5} \|\boldsymbol{t}^1\|_{L^2(\Gamma_{N_d})}^2 + \frac{\varepsilon_5}{2} \|\boldsymbol{u}_{\mathcal{D}}^0\|_{L^2(\Gamma_{N_d})}^2 $$
**Remark 5.** The coercivity property (36) cannot be easily proved for non-linear fluxes. As mentioned before, it holds for the SUSHI scheme but also for the linear two-point flux approximation that results from (25)–(30) when used on admissible meshes in the sense of [14, Definition 11.1] with a harmonic average of both cell values for the trace reconstruction. In Section 5, the value of $C_\mathcal{D}$ has been evaluated for different sequences of grids.

In order to prove the existence of solutions in Proposition 7, we first state and prove the following a priori estimate.

**Proposition 6 (A priori estimate).** Let $\mathcal{D}$ be an admissible discretization of $\Omega$ in the sense of Definition 1. Under Hypotheses 4 and for the sequence of solutions $(\mathbf{u}_{\mathcal{D}}^n, p_{\mathcal{D}}^n)_{n=1,...,N+1}$ to (15)-(21), there exists a constant $\tilde{B}$ s.t.

$$ \|\mathbf{u}_{\mathcal{D}}^{N+1}\|_{H^1(\Omega)}^2 + \|p_{\mathcal{D}}^{N+1}\|_{L^2(\Omega)}^2 + \sum_{n=1}^{N+1} \Delta t \|p_{\mathcal{D}}^n\|_{\mathcal{D}}^2 \le \tilde{B} \exp(T). \quad (37) $$

**Proof.** Taking $v_{\mathcal{D}} = D^{n+1}(\mathbf{u}_{\mathcal{D}})$ in (15), multiplying (21) by $p_K^{n+1}$, summing over $K \in \mathcal{F}$ and $n \in [0, N]$, we obtain

$$ \sum_{n=0}^{N} \sum_{K \in \mathcal{F}} \left( T_{a,K}^{n+1} + T_{c,K}^{n+1} + T_{F,K}^{n+1} \right) \le \sum_{n=0}^{N} \sum_{K \in \mathcal{F}} \left( T_{f,K}^{n+1} + T_{t,K}^{n+1} + T_{s,K}^{n+1} \right) \quad (38) $$

with

$$ \begin{aligned} T_{a,K}^{n+1} &= a_{\mathcal{D}}^{K}(\mathbf{u}_{\mathcal{D}}^{n+1}, D^{n+1}(\mathbf{u}_{\mathcal{D}})), \\ T_{c,K}^{n+1} &= m_K c_0 D^{n+1}(p_K) p_K^{n+1}, \\ T_{F,K}^{n+1} &= -\Delta t \sum_{\sigma \in \partial K} F_{K,\sigma}(p^{n+1}, p^{n+1}) p_K^{n+1}, \end{aligned} $$

$$ \begin{aligned} T_{f,K}^{n+1} &= \int_K \boldsymbol{f}^{n+1} \cdot D^{n+1}(\boldsymbol{u}_{\mathcal{D}}), \\ T_{t,K}^{n+1} &= \int_{\partial K \cap \Gamma_{N_d}} \boldsymbol{t}^{n+1} \cdot D^{n+1}(\boldsymbol{u}_{\mathcal{D}}), \\ T_{s,K}^{n+1} &= \int_{t^n}^{t^{n+1}} \int_K s p_K^{n+1}. \end{aligned} $$

A reordering of the sums in time (see the fourth tool of [3, section A.1]) yields:

$$ \begin{aligned} \sum_{n=0}^{N} \sum_{K \in \mathcal{F}} T_{f,K}^{n+1} &= \sum_{K \in \mathcal{F}} \left( \int_K \boldsymbol{f}^{N+1} \cdot \mathbf{u}_{\mathcal{D}}^{N+1} - \int_K \boldsymbol{f}^1 \cdot \mathbf{u}_{\mathcal{D}}^0 - \sum_{n=1}^{N} \int_K D^{n+1}(\boldsymbol{f}) \cdot \mathbf{u}_{\mathcal{D}}^n \right), \\ \sum_{n=0}^{N} \sum_{K \in \mathcal{F}} T_{t,K}^{n+1} &= \sum_{K \in \mathcal{F}} \left( \int_{\partial K \cap \Gamma_{N_d}} \boldsymbol{t}^{N+1} \cdot \mathbf{u}_{\mathcal{D}}^{N+1} - \int_{\partial K \cap \Gamma_{N_d}} \boldsymbol{t}^1 \cdot \mathbf{u}_{\mathcal{D}}^0 - \sum_{n=1}^{N} \int_{\partial K \cap \Gamma_{N_d}} D^{n+1}(\boldsymbol{t}) \cdot \mathbf{u}_{\mathcal{D}}^n \right). \end{aligned} $$

Several applications of Cauchy-Schwartz and Young inequalities with parameters $(\varepsilon_i)_{i=1...6} \in (0, 1)$ thus lead to

$$ \left| \sum_{n=0}^{N} \sum_{K \in \mathcal{F}} T_{f,K}^{n+1} \right| \le \frac{1}{2\varepsilon_1} \|\boldsymbol{f}^{N+1}\|_{L^2(\Omega)}^2 + \frac{\varepsilon_1}{2} \|\boldsymbol{u}_{\mathcal{D}}^{N+1}\|_{H^1(\Omega)}^2 + \frac{1}{2\varepsilon_2} \|\boldsymbol{f}^1\|_{L^2(\Omega)}^2 + \frac{\varepsilon_2}{2} \|\boldsymbol{u}_{\mathcal{D}}^0\|_{H^1(\Omega)}^2 \\ + \sum_{n=1}^{N} \left( \frac{1}{2\varepsilon_3} \int_{t^n}^{t^{n+1}} \|\boldsymbol{f}'\|_{L^2(\Omega)}^2 + \frac{\Delta t \varepsilon_3}{2} \|\boldsymbol{u}_{\mathcal{D}}^n\|_{H^1(\Omega)}^2 \right) \quad (39) $$

$$ \left| \sum_{n=0}^{N} \sum_{K \in \mathcal{F}} T_{t,K}^{n+1} \right| \le \frac{1}{2\varepsilon_4} \|\boldsymbol{t}^{N+1}\|_{L^2(\Gamma_{N_d})}^2 + \frac{\varepsilon_4}{2} \|\boldsymbol{u}_{\mathcal{D}}^{N+1}\|_{L^2(\Gamma_{N_d})}^2 + \frac{1}{2\varepsilon_5} \|\boldsymbol{t}^1\|_{L^2(\Gamma_{N_d})}^2 + \frac{\varepsilon_5}{2} \|\boldsymbol{u}_{\mathcal{D}}^0\|_{L^2(\Gamma_{N_d})}^2 $$
$$ + \sum_{n=1}^{N} \left( \frac{1}{2\epsilon_6} \int_{t^n}^{t^{n+1}} \|\mathbf{t}'\|_{L^2(\Gamma_{N_d})}^2 + \frac{\Delta t \epsilon_6}{2} \|\mathbf{u}_{\mathcal{D}}^n\|_{L^2(\Gamma_{N_d})}^2 \right) \quad (40) $$

$$ \left| \sum_{n=0}^{N} \sum_{K \in \mathcal{F}} T_{s,K}^{n+1} \right| \le \frac{1}{2\epsilon_7} \|s\|_{L^2(\Omega \times (0,T))}^2 + \frac{\epsilon_7}{2} \sum_{n=0}^{N} \Delta t \|p_{\mathcal{D}}^{n+1}\|_{L^2(\Omega)}^2. \quad (41) $$

(42)

Let us notice that, since $\Omega$ has a Lipschitz boundary, according to [9, Theorem 1.6.6], for any time step $m$, there exists a constant $C$ s.t.

$$ \| \mathbf{u}_{\mathcal{D}}^{m} \|_{L^2(\Gamma_{N_d})}^{2} \le \| \mathbf{u}_{\mathcal{D}}^{m} \|_{L^2(\partial \Omega)}^{2} \le C \| \mathbf{u}_{\mathcal{D}}^{m} \|_{L^2(\Omega)} \| \mathbf{u}_{\mathcal{D}}^{m} \|_{H^1(\Omega)} \le C \| \mathbf{u}_{\mathcal{D}}^{m} \|_{H^1(\Omega)}^{2}. \quad (43) $$

Since

$$ 2 \int_K a_{\mathcal{D}}^K (\mathbf{u}_{\mathcal{D}}^{n+1}, D^{n+1}(\mathbf{u}_{\mathcal{D}})) = \int_K a_{\mathcal{D}}^K (D^{n+1}(\mathbf{u}_{\mathcal{D}}), D^{n+1}(\mathbf{u}_{\mathcal{D}})) + \int_K a_{\mathcal{D}}^K (\mathbf{u}_{\mathcal{D}}^{n+1}, \mathbf{u}_{\mathcal{D}}^{n+1}) - \int_K a_{\mathcal{D}}^K (\mathbf{u}_{\mathcal{D}}^n, \mathbf{u}_{\mathcal{D}}^n), $$

we have

$$ \frac{1}{2} (a_{\mathcal{D}}(\mathbf{u}_{\mathcal{D}}^{N+1}, \mathbf{u}_{\mathcal{D}}^{N+1}) - a_{\mathcal{D}}(\mathbf{u}_{\mathcal{D}}^{0}, \mathbf{u}_{\mathcal{D}}^{0})) \le \sum_{n=0}^{N} \sum_{K \in \mathcal{F}} T_{a,K}^{n+1}. \quad (44) $$

In a similar way:

$$ \frac{c_0}{2} \left( \|p_{\mathcal{D}}^{N+1}\|_{L^2(\Omega)}^2 - \|p_{\mathcal{D}}^{0}\|_{L^2(\Omega)}^2 \right) \le \sum_{n=0}^{N} \sum_{K \in \mathcal{F}} T_{c,K}^{n+1}. \quad (45) $$

From (36), we deduce that:

$$ \sum_{n=0}^{N} C_{\mathcal{D}} \Delta t \| p_{\mathcal{D}}^{n+1} \|_{\mathcal{D}}^2 \le \sum_{n=0}^{N} \sum_{K \in \mathcal{F}} T_{F,K}^{n+1}. \quad (46) $$

Gathering inequalities (38)–(46) and reordering the terms lead to

$$ \begin{aligned} & a_{\mathcal{D}}(\mathbf{u}_{\mathcal{D}}^{N+1}, \mathbf{u}_{\mathcal{D}}^{N+1}) + c_0 \|p_{\mathcal{D}}^{N+1}\|_{L^2(\Omega)}^2 + 2 \sum_{n=0}^{N} C_{\mathcal{D}} \Delta t \|p_{\mathcal{D}}^{n+1}\|_{\mathcal{D}}^2 \\ & \le a_{\mathcal{D}}(\mathbf{u}_{\mathcal{D}}^{0}, \mathbf{u}_{\mathcal{D}}^{0}) + c_0 \|p_{\mathcal{D}}^{0}\|_{L^2(\Omega)}^2 + \frac{1}{\epsilon_2} \|\mathbf{f}^{1}\|_{L^2(\Omega)}^2 + \frac{1}{\epsilon_3} \|\mathbf{f}'\|_{L^2(\Omega \times (0,T))}^2 + \frac{1}{\epsilon_1} \|\mathbf{f}^{N+1}\|_{L^2(\Omega)}^2 \\ & + \frac{1}{\epsilon_5} \|\mathbf{t}^{1}\|_{L^2(\Gamma_{N_d})}^2 + \frac{1}{\epsilon_6} \|\mathbf{t}'\|_{L^2(\Gamma_{N_d} \times (0,T))}^2 + \frac{1}{\epsilon_4} \|\mathbf{t}^{N+1}\|_{L^2(\Gamma_{N_d})}^2 \\ & + (\epsilon_2 + \epsilon_5 C) \|\mathbf{u}_{\mathcal{D}}^{0}\|_{H^1(\Omega)}^2 + (\epsilon_1 + \epsilon_4 C) \|\mathbf{u}_{\mathcal{D}}^{N+1}\|_{H^1(\Omega)}^2 + (\epsilon_3 + \epsilon_6 C) \sum_{n=1}^{N} \Delta t \|\mathbf{u}_{\mathcal{D}}^{n}\|_{H^1(\Omega)}^2 \\ & + \frac{1}{\epsilon_7} \|s\|_{L^2(\Omega \times (0,T))}^2 + \epsilon_7 \sum_{n=0}^{N} \Delta t \|p_{\mathcal{D}}^{n+1}\|_{L^2(\Omega)}^2. \end{aligned} \quad (47) $$

Let us notice that, thanks to (11), (12), (18), for all $\mathbf{v}_{\mathcal{D}} \in V_{\mathcal{D}}$

$$ \alpha_a \alpha_* \|\mathbf{v}_{\mathcal{D}}\|_{H^1(\Omega)}^2 \le a_{\mathcal{D}}(\mathbf{v}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}}) \le \beta_a \alpha_* \|\mathbf{v}_{\mathcal{D}}\|_{H^1(\Omega)}^2. \quad (48) $$

On the other hand, using Lemma 9 given in appendix, we have

$$ \|p_{\mathcal{D}}^{n+1}\|_{L^2(\Omega)} \le \tilde{C}_{\mathcal{D}} \|p_{\mathcal{D}}^{n+1}\|_{\mathcal{D}}. \quad (49) $$
$$ + \sum_{n=1}^{N} \left( \frac{1}{2\varepsilon_6} \int_{t^n}^{t^{n+1}} \|\mathbf{t}'\|_{L^2(\Gamma_{N_d})}^2 + \frac{\Delta t \varepsilon_6}{2} \|\mathbf{u}_{\mathcal{D}}^n\|_{L^2(\Gamma_{N_d})}^2 \right) \quad (40) $$

$$ \left| \sum_{n=0}^{N} \sum_{K \in \mathcal{F}} T_{s,K}^{n+1} \right| \le \frac{1}{2\varepsilon_7} \|s\|_{L^2(\Omega \times (0,T))}^2 + \frac{\varepsilon_7}{2} \sum_{n=0}^{N} \Delta t \|p_{\mathcal{D}}^{n+1}\|_{L^2(\Omega)}^2. \quad (41) $$

(42)

Let us notice that, since $\Omega$ has a Lipschitz boundary, according to [9, Theorem 1.6.6], for any time step $m$, there exists a constant $C$ s.t.

$$ \| \mathbf{u}_{\mathcal{D}}^m \|_{L^2(\Gamma_{N_d})}^2 \le \| \mathbf{u}_{\mathcal{D}}^m \|_{L^2(\partial \Omega)}^2 \le C \| \mathbf{u}_{\mathcal{D}}^m \|_{L^2(\Omega)} \| \mathbf{u}_{\mathcal{D}}^m \|_{H^1(\Omega)} \le C \| \mathbf{u}_{\mathcal{D}}^m \|_{H^1(\Omega)}^2. \quad (43) $$

Since

$$ 2 \int_K a_{\mathcal{D}}^K (\mathbf{u}_{\mathcal{D}}^{n+1}, D^{n+1}(\mathbf{u}_{\mathcal{D}})) = \int_K a_{\mathcal{D}}^K (D^{n+1}(\mathbf{u}_{\mathcal{D}}), D^{n+1}(\mathbf{u}_{\mathcal{D}})) + \int_K a_{\mathcal{D}}^K (\mathbf{u}_{\mathcal{D}}^{n+1}, \mathbf{u}_{\mathcal{D}}^{n+1}) - \int_K a_{\mathcal{D}}^K (\mathbf{u}_{\mathcal{D}}^n, \mathbf{u}_{\mathcal{D}}^n), $$

we have

$$ \frac{1}{2} (a_{\mathcal{D}}(\mathbf{u}_{\mathcal{D}}^{N+1}, \mathbf{u}_{\mathcal{D}}^{N+1}) - a_{\mathcal{D}}(\mathbf{u}_{\mathcal{D}}^{0}, \mathbf{u}_{\mathcal{D}}^{0})) \le \sum_{n=0}^{N} \sum_{K \in \mathcal{F}} T_{a,K}^{n+1}. \quad (44) $$

In a similar way:

$$ \frac{c_0}{2} \left( \| p_{\mathcal{D}}^{N+1} \|_{L^2(\Omega)}^2 - \| p_{\mathcal{D}}^{0} \|_{L^2(\Omega)}^2 \right) \le \sum_{n=0}^{N} \sum_{K \in \mathcal{F}} T_{c,K}^{n+1}. \quad (45) $$

From (36), we deduce that:

$$ \sum_{n=0}^{N} C_{\mathcal{D}} \Delta t \| p_{\mathcal{D}}^{n+1} \|_{\mathcal{D}}^2 \le \sum_{n=0}^{N} \sum_{K \in \mathcal{F}} T_{F,K}^{n+1}. \quad (46) $$

Gathering inequalities (38)–(46) and reordering the terms lead to

$$ \begin{aligned} & a_{\mathcal{D}}(\mathbf{u}_{\mathcal{D}}^{N+1}, \mathbf{u}_{\mathcal{D}}^{N+1}) + c_0 \|p_{\mathcal{D}}^{N+1}\|_{L^2(\Omega)}^2 + 2 \sum_{n=0}^{N} C_{\mathcal{D}} \Delta t \|p_{\mathcal{D}}^{n+1}\|_{\mathcal{D}}^2 \\ & \le a_{\mathcal{D}}(\mathbf{u}_{\mathcal{D}}^{0}, \mathbf{u}_{\mathcal{D}}^{0}) + c_0 \|p_{\mathcal{D}}^{0}\|_{L^2(\Omega)}^2 + \frac{1}{\varepsilon_2} \|\mathbf{f}^{1}\|_{L^2(\Omega)}^2 + \frac{1}{\varepsilon_3} \|\mathbf{f}'\|_{L^2(\Omega \times (0,T))}^2 + \frac{1}{\varepsilon_1} \|\mathbf{f}^{N+1}\|_{L^2(\Omega)}^2 \\ & + \frac{1}{\varepsilon_5} \|\mathbf{t}^{1}\|_{L^2(\Gamma_{N_d})}^2 + \frac{1}{\varepsilon_6} \|\mathbf{t}'\|_{L^2(\Gamma_{N_d} \times (0,T))}^2 + \frac{1}{\varepsilon_4} \|\mathbf{t}^{N+1}\|_{L^2(\Gamma_{N_d})}^2 \\ & + (\varepsilon_2 + \varepsilon_5 C) \|\mathbf{u}_{\mathcal{D}}^{0}\|_{H^1(\Omega)}^2 + (\varepsilon_1 + \varepsilon_4 C) \|\mathbf{u}_{\mathcal{D}}^{N+1}\|_{H^1(\Omega)}^2 + (\varepsilon_3 + \varepsilon_6 C) \sum_{n=1}^{N} \Delta t \|\mathbf{u}_{\mathcal{D}}^{n}\|_{H^1(\Omega)}^2 \\ & + \frac{1}{\varepsilon_7} \|s\|_{L^2(\Omega \times (0,T))}^2 + \varepsilon_7 \sum_{n=0}^{N} \Delta t \|p_{\mathcal{D}}^{n+1}\|_{L^2(\Omega)}^2. \end{aligned} \quad (47) $$

Let us notice that, thanks to (11), (12), (18), for all $\mathbf{v}_{\mathcal{D}} \in V_{\mathcal{D}}$

$$ \alpha_a \alpha_* \|\mathbf{v}_{\mathcal{D}}\|_{H^1(\Omega)}^2 \le a_{\mathcal{D}}(\mathbf{v}_{\mathcal{D}}, \mathbf{v}_{\mathcal{D}}) \le \beta_a \alpha_* \|\mathbf{v}_{\mathcal{D}}\|_{H^1(\Omega)}^2. \quad (48) $$

On the other hand, using Lemma 9 given in appendix, we have

$$ \| p_{\mathcal{D}}^{n+1} \|_{L^2(\Omega)} \le \tilde{C}_{\mathcal{D}} \| p_{\mathcal{D}}^{n+1} \|_{\mathcal{D}}. \quad (49) $$
Thus, setting

$$
\begin{align*}
B ={}& (\beta_a \alpha^* + \varepsilon_2 + \varepsilon_5 C) \|\mathbf{u}_{\mathcal{D}}^0\|_{H^1(\Omega)}^2 + c_0 \|\mathbf{p}_{\mathcal{D}}^0\|_{L^2(\Omega)}^2 \\
& + \frac{1}{\varepsilon_2} \|\mathbf{f}^1\|_{L^2(\Omega)}^2 + \frac{1}{\varepsilon_3} \|\mathbf{f}'\|_{L^2(\Omega \times (0,T))}^2 + \frac{1}{\varepsilon_1} \|\mathbf{f}^{N+1}\|_{L^2(\Omega)}^2 \\
& + \frac{1}{\varepsilon_5} \|\mathbf{t}^1\|_{L^2(\Gamma_{N_d})}^2 + \frac{1}{\varepsilon_6} \|\mathbf{t}'\|_{L^2(\Gamma_{N_d} \times (0,T))}^2 + \frac{1}{\varepsilon_4} \|\mathbf{t}^{N+1}\|_{L^2(\Gamma_{N_d})}^2 + \frac{1}{\varepsilon_7} \|\mathbf{s}\|_{L^2(\Omega \times (0,T))}^2,
\end{align*}
$$

we deduce from (47)-(49)

$$
\begin{equation}
\begin{aligned}
& (\alpha_a \alpha_* - (\varepsilon_1 + \varepsilon_4 C)) \|\mathbf{u}_{\mathcal{D}}^{N+1}\|_{H^1(\Omega)}^2 + c_0 \|\mathbf{p}_{\mathcal{D}}^{N+1}\|_{L^2(\Omega)}^2 \\
& \qquad + (2C_{\mathcal{D}} - \varepsilon_7 \tilde{C}_{\mathcal{D}}) \sum_{n=0}^{N} \Delta t \|\mathbf{p}_{\mathcal{D}}^{n+1}\|_{\mathcal{D}}^2 \le (\varepsilon_3 + \varepsilon_6 C) \sum_{n=0}^{N} \Delta t \|\mathbf{u}_{\mathcal{D}}^n\|_{H^1(\Omega)}^2 + B.
\end{aligned}
\tag{50}
\end{equation}
$$

Note that B depends on $\|\mathbf{u}_{\mathcal{D}}^0\|_{H^1(\Omega)}^2$. In fact, starting from (15) with $v_{\mathcal{D}}^0 = \mathbf{u}_{\mathcal{D}}^0$ and using Cauchy-Schwartz and Young inequalities in the same way as before, we can also establish that there exists a constant $C_{\mathcal{D}}^0$ which only depends on $\alpha, \alpha_a$ and $\alpha_*$ such that:

$$
\| \mathbf{u}_{\mathcal{D}}^0 \|_{H^1(\Omega)}^2 \le C_{\mathcal{D}}^0 \left( \| \mathbf{p}_{\mathcal{D}}^0 \|_{L^2(\Omega)}^2 + \| \mathbf{f}^0 \|_{L^2(\Omega)}^2 + \| \mathbf{t}^0 \|_{L^2(\Gamma_{N_d})}^2 \right).
$$

In (50), by choosing $\varepsilon_1$ and $\varepsilon_4$ s.t. $(\varepsilon_1 + \varepsilon_4 C) < \frac{\alpha_a \alpha_*}{2}$, $\varepsilon_7 = \frac{C_{\mathcal{D}}}{C_{\mathcal{D}}}$, $\varepsilon_3$ and $\varepsilon_6$ s.t. $\varepsilon_3 + \varepsilon_6 C = \alpha_r = \min(\frac{\alpha_a \alpha_*}{2}, c_0, C_{\mathcal{D}})$, setting $\tilde{B} = \frac{B}{\alpha_r}$, we deduce from (50) that

$$
\| \mathbf{u}_{\mathcal{D}}^{N+1} \|_{H^1(\Omega)}^2 + \| \mathbf{p}_{\mathcal{D}}^{N+1} \|_{L^2(\Omega)}^2 + \sum_{n=1}^{N+1} \Delta t \| \mathbf{p}_{\mathcal{D}}^n \|_{\mathcal{D}}^2 \le \sum_{n=0}^{N} \Delta t \left( \| \mathbf{u}_{\mathcal{D}}^n \|_{H^1(\Omega)}^2 + \| \mathbf{p}_{\mathcal{D}}^n \|_{L^2(\Omega)}^2 \right) + \tilde{B}.
$$

This now allows us to apply the Lemma 4 given in [3] which is a discrete version of Grönwall's lemma by setting $a_n = \|\mathbf{u}_{\mathcal{D}}^n\|_{H^1(\Omega)}^2 + \|\mathbf{p}_{\mathcal{D}}^n\|_{L^2(\Omega)}^2$ and $b_n = \sum_{m=1}^n \Delta t \|\mathbf{p}_{\mathcal{D}}^m\|_{\mathcal{D}}^2$, $\lambda = \Delta t$, and $C = \tilde{B}$, which yields (37). $\square$

**Proposition 7 (Existence of discrete solutions).** Under assumptions of Proposition 6, scheme (15)-(21) admits at each time $t^{n_0}$, $n_0 \in \{0...N+1\}$, a discrete solution $(u_{\mathcal{D}}^{n_0}, p_{\mathcal{D}}^{n_0})$.

**Proof.** The proof of (37) remains valid by replacing $N + 1$ with $n_0$. We thus have

$$
\| \mathbf{u}_{\mathcal{D}}^{n_0} \|_{H^1(\Omega)}^2 + \| \mathbf{p}_{\mathcal{D}}^{n_0} \|_{L^2(\Omega)}^2 \leq C.
$$

The fact that these solutions remain bounded allows us to use Brouwer's topological degree and to conclude to their existence. $\square$

We conclude this section with the theorem below which states an a priori error estimate with a convergence order equal to one in space and time for the proposed scheme.

**Theorem 8 (Convergence and a priori error estimate).** Under hypotheses of Proposition 6, assuming the exact solution $(\mathbf{u}, \mathbf{p})$ to be sufficiently smooth and for $\Delta t \le \frac{c_0}{2}$, we have

$$
\left( \| \mathbf{u} - \mathbf{u}_{\mathcal{D}} \|_{L^{\infty}(0,T;H^1(\Omega))}^2 + \| p - p_{\mathcal{D}} \|_{L^{\infty}(0,T;L^2(\Omega))}^2 + \| p - p_{\mathcal{D}} \|_{L^{2}(0,T;Q_{\mathcal{D}})}^2 \right)^{\frac{1}{2}} \le C (\Delta t + h_{\mathcal{F}}) \quad (51)
$$

where $C > 0$ is a positive constant depending on the physical parameters of the problem, $(\rho_i)_{i=1...4}$ and the time derivatives of the true solution $\partial_{tt}\mathbf{u}$ and $\partial_{tt}\mathbf{p}$.

**Proof.** Replacing [3, Eq. 31] where the Two Point Flux approximation was used for the fluxes $F_{K,\sigma}$, by inequality (36), allows us to apply the same arguments as in [3, Proposition 3 and Theorem 5] and to conclude to (51). $\square$
Thus, setting

$$
\begin{align*}
B ={}& (\beta_a \alpha^* + \varepsilon_2 + \varepsilon_5 C) \|\mathbf{u}_{\mathcal{D}}^0\|_{H^1(\Omega)}^2 + c_0 \|\mathbf{p}_{\mathcal{D}}^0\|_{L^2(\Omega)}^2 \\
& + \frac{1}{\varepsilon_2} \|\mathbf{f}^1\|_{L^2(\Omega)}^2 + \frac{1}{\varepsilon_3} \|\mathbf{f}'\|_{L^2(\Omega \times (0,T))}^2 + \frac{1}{\varepsilon_1} \|\mathbf{f}^{N+1}\|_{L^2(\Omega)}^2 \\
& + \frac{1}{\varepsilon_5} \|\mathbf{t}^1\|_{L^2(\Gamma_{N_d})}^2 + \frac{1}{\varepsilon_6} \|\mathbf{t}'\|_{L^2(\Gamma_{N_d} \times (0,T))}^2 + \frac{1}{\varepsilon_4} \|\mathbf{t}^{N+1}\|_{L^2(\Gamma_{N_d})}^2 + \frac{1}{\varepsilon_7} \|\mathbf{s}\|_{L^2(\Omega \times (0,T))}^2,
\end{align*}
$$

we deduce from (47)-(49)

$$
(\alpha_a \alpha_* - (\varepsilon_1 + \varepsilon_4 C)) \|\mathbf{u}_{\mathcal{D}}^{N+1}\|_{H^1(\Omega)}^2 + c_0 \|\mathbf{p}_{\mathcal{D}}^{N+1}\|_{L^2(\Omega)}^2 \\
+ (2C_{\mathcal{D}} - \varepsilon_7 \tilde{C}_{\mathcal{D}}) \sum_{n=0}^{N} \Delta t \|\mathbf{p}_{\mathcal{D}}^{n+1}\|_{\mathcal{D}}^2 \le (\varepsilon_3 + \varepsilon_6 C) \sum_{n=0}^{N} \Delta t \|\mathbf{u}_{\mathcal{D}}^n\|_{H^1(\Omega)}^2 + B. \quad (50)
$$

Note that B depends on $\|\mathbf{u}_{\mathcal{D}}^0\|_{H^1(\Omega)}^2$. In fact, starting from (15) with $\mathbf{v}_{\mathcal{D}}^0 = \mathbf{u}_{\mathcal{D}}^0$ and using Cauchy-Schwartz and Young inequalities in the same way as before, we can also establish that there exists a constant $C_{\mathcal{D}}^0$ which only depends on $\alpha, \alpha_a$ and $\alpha_*$ such that:

$$
\| \mathbf{u}_{\mathcal{D}}^{0} \|_{H^{1}(\Omega)}^{2} \leq C_{\mathcal{D}}^{0} \left( \| \mathbf{p}_{\mathcal{D}}^{0} \|_{L^{2}(\Omega)}^{2} + \| \mathbf{f}^{0} \|_{L^{2}(\Omega)}^{2} + \| \mathbf{t}^{0} \|_{L^{2}(\Gamma_{N_d})}^{2} \right).
$$

In (50), by choosing $\varepsilon_1$ and $\varepsilon_4$ s.t. $(\varepsilon_1 + \varepsilon_4 C) < \frac{\alpha_a \alpha_*}{2}$, $\varepsilon_7 = \frac{C_{\mathcal{D}}}{C_{\mathcal{D}}}$, $\varepsilon_3$ and $\varepsilon_6$ s.t. $\varepsilon_3 + \varepsilon_6 C = \alpha_r = \min(\frac{\alpha_a \alpha_*}{2}, c_0, C_{\mathcal{D}})$, setting $\tilde{B} = \frac{B}{\alpha_r}$, we deduce from (50) that

$$
\| \mathbf{u}_{\mathcal{D}}^{N+1} \|_{H^1(\Omega)}^2 + \| \mathbf{p}_{\mathcal{D}}^{N+1} \|_{L^2(\Omega)}^2 + \sum_{n=1}^{N+1} \Delta t \| \mathbf{p}_{\mathcal{D}}^n \|_{\mathcal{D}}^2 \le \sum_{n=0}^{N} \Delta t \left( \| \mathbf{u}_{\mathcal{D}}^n \|_{H^1(\Omega)}^2 + \| \mathbf{p}_{\mathcal{D}}^n \|_{L^2(\Omega)}^2 \right) + \tilde{B}.
$$

This now allows us to apply the Lemma 4 given in [3] which is a discrete version of Grönwall's lemma by setting $a_n = \|\mathbf{u}_{\mathcal{D}}^n\|_{H^1(\Omega)}^2 + \|\mathbf{p}_{\mathcal{D}}^n\|_{L^2(\Omega)}^2$ and $b_n = \sum_{m=1}^n \Delta t \|\mathbf{p}_{\mathcal{D}}^m\|_{\mathcal{D}}^2$, $\lambda = \Delta t$, and $C = \tilde{B}$, which yields (37). $\square$

**Proposition 7 (Existence of discrete solutions).** Under assumptions of Proposition 6, scheme (15)-(21) admits at each time $t^{n_0}$, $n_0 \in \{0...N+1\}$, a discrete solution $(\mathbf{u}_{\mathcal{D}}^{n_0}, \mathbf{p}_{\mathcal{D}}^{n_0})$.

**Proof.** The proof of (37) remains valid by replacing $N + 1$ with $n_0$. We thus have

$$
\| \mathbf{u}_{\mathcal{D}}^{n_0} \|_{H^1(\Omega)}^2 + \| \mathbf{p}_{\mathcal{D}}^{n_0} \|_{L^2(\Omega)}^2 \leq C.
$$

The fact that these solutions remain bounded allows us to use Brouwer's topological degree and to conclude to their existence. $\square$

We conclude this section with the theorem below which states an a priori error estimate with a convergence order equal to one in space and time for the proposed scheme.

**Theorem 8 (Convergence and a priori error estimate).** Under hypotheses of Proposition 6, assuming the exact solution $(\mathbf{u}, \mathbf{p})$ to be sufficiently smooth and for $\Delta t \le \frac{c_0}{2}$, we have

$$
\left( \| \mathbf{u} - \mathbf{u}_{\mathcal{D}} \|_{L^{\infty}(0,T;H^1(\Omega))}^2 + \| p - p_{\mathcal{D}} \|_{L^{\infty}(0,T;L^2(\Omega))}^2 + \| p - p_{\mathcal{D}} \|_{L^{2}(0,T;Q_{\mathcal{D}})}^2 \right)^{\frac{1}{2}} \le C (\Delta t + h_{\mathcal{F}}) \quad (51)
$$

where $C > 0$ is a positive constant depending on the physical parameters of the problem, $(\rho_i)_{i=1...4}$ and the time derivatives of the true solution $\partial_{tt}\mathbf{u}$ and $\partial_{tt}\mathbf{p}$.

**Proof.** Replacing [3, Eq. 31] where the Two Point Flux approximation was used for the fluxes $F_{K,\sigma}$, by inequality (36), allows us to apply the same arguments as in [3, Proposition 3 and Theorem 5] and to conclude to (51). $\square$
<figure><img src="image_5.png" /><figcaption>Figure 1. Grids of the FVCA5 benchmark used for the numerical convergence study</figcaption></figure>

5. Numerical results

We here illustrate the numerical convergence of the scheme using the academic test case used in [3, Section 5.1]. The domain Ω is the unit square and T = 1. This medium is assumed to be homogeneous with respect to elastic and flow properties. Here the Young modulus E = 2.5, the poisson ratio v = 0.25 and the mobility matrix Λ is the identity matrix, α = 1 and c₀ = 0.5. The displacements and pressures are given by the following analytical functions:

$$
\mathbf{u}(\mathbf{x}, t) = 10^{-2} e^{-t} \begin{pmatrix} x^2 y \\ -xy^2 \end{pmatrix}, \qquad (52)
$$

$$
p(\mathbf{x}, t) = e^{-t} \sin\left(\frac{x}{\sqrt{2}}\right) \sin\left(\frac{y}{\sqrt{2}}\right). \qquad (53)
$$

In this example Γ<sub>D<sub>d</sub></sub> = Γ<sub>D<sub>p</sub></sub> = ∂Ω. *f* and *s* are chosen in order to make (52)-(53) solutions of (1)-(2) for *t* > 0 and of (1)-(4) for *t* = 0. Here, unlike [3] which used a sequence of Voronoï meshes which yield consistent TPFA flux approximations, we consider the grids *mesh1*, *mesh3* et *mesh4_1* of the FVCA5 benchmark [6] which are represented on Figures 1a–1c at their coarsest resolution. These last two sets of grids, in particular, are of practical interest since such grids may be used when one either refines the mesh around the source term or wants to have piecewise-constant material properties following the geometry of the geological heterogeneities for porous media flows. The time steps Δ*t* = 0.2, 0.1, 0.02, 0.01 have been used for each grid of the data set. For this experiment, the non-linear two-point flux approximation (27) has been used for the interior fluxes and the linear flux (24) for the outer ones. The operator γ defined in (20) was used as discrete trace operator *I*. The approximation error

$$
E_A = \left( \| \mathbf{u} - \mathbf{u}_{\mathcal{D}} \|_{L^{\infty}(0,T; H^1(\Omega))}^2 + \| p - p_{\mathcal{D}} \|_{L^{\infty}(0,T; L^2(\Omega))}^2 + \| p - p_{\mathcal{D}} \|_{L^2(0,T; Q_{\mathcal{D}})}^2 \right)^{\frac{1}{2}},
$$

defined in the left hand side of inequality (51), has been used to verify the convergence of the
proposed scheme in space. Figures 2 show the error evolution with respect to *h* for the three
<figure><img src="image_5.png" /><figcaption>Figure 1. Grids of the FVCA5 benchmark used for the numerical convergence study</figcaption></figure>

## 5. Numerical results

We here illustrate the numerical convergence of the scheme using the academic test case used in [3, Section 5.1]. The domain $\Omega$ is the unit square and $T = 1$. This medium is assumed to be homogeneous with respect to elastic and flow properties. Here the Young modulus $E = 2.5$, the poisson ratio $\nu = 0.25$ and the mobility matrix $\bar{\Lambda}$ is the identity matrix, $\alpha = 1$ and $c_0 = 0.5$. The displacements and pressures are given by the following analytical functions:

$$ \mathbf{u}(\mathbf{x}, t) = 10^{-2} e^{-t} \begin{pmatrix} x^2 y \\ -xy^2 \end{pmatrix}, \qquad (52) $$

$$ p(\mathbf{x}, t) = e^{-t} \sin\left(\frac{x}{\sqrt{2}}\right) \sin\left(\frac{y}{\sqrt{2}}\right). \qquad (53) $$

In this example $\Gamma_{D_d} = \Gamma_{D_p} = \partial\Omega$. $f$ and $s$ are chosen in order to make (52)-(53) solutions of (1)-(2) for $t > 0$ and of (1)-(4) for $t = 0$. Here, unlike [3] which used a sequence of Voronoï meshes which yield consistent TPFA flux approximations, we consider the grids `mesh1`, `mesh3` et `mesh4_1` of the FVCA5 benchmark [6] which are represented on Figures 1a–1c at their coarsest resolution. These last two sets of grids, in particular, are of practical interest since such grids may be used when one either refines the mesh around the source term or wants to have piecewise-constant material properties following the geometry of the geological heterogeneities for porous media flows. The time steps $\Delta t = 0.2, 0.1, 0.02, 0.01$ have been used for each grid of the data set. For this experiment, the non-linear two-point flux approximation (27) has been used for the interior fluxes and the linear flux (24) for the outer ones. The operator $\gamma$ defined in (20) was used as discrete trace operator $I$. The approximation error

$$ E_A = \left( \|\mathbf{u} - \mathbf{u}_{\mathcal{D}}\|_{L^{\infty}(0,T;H^1(\Omega))}^2 + \|p - p_{\mathcal{D}}\|_{L^{\infty}(0,T;L^2(\Omega))}^2 + \|p - p_{\mathcal{D}}\|_{L^2(0,T;Q_{\mathcal{D}})}^2 \right)^{\frac{1}{2}}, $$

defined in the left hand side of inequality (51), has been used to verify the convergence of the proposed scheme in space. Figures 2 show the error evolution with respect to $h$ for the three
cases. As expected a slope of one is obtained for all of them and for sufficiently small values of the time steps. Figure 3 displays the time convergence for the third case. We observe that reducing the time step does not improve the accuracy when h is too large but a convergence order of one can be retrieved at smaller space resolutions.

Although inequality (36) was assumed to be satisfied for the numerical analysis of the scheme, we have numerically checked that it holds. With that purpose, we have evaluated the following ratio

$$c_{\mathcal{D}}(p) = \inf_{t \in (0, T)} \left( \frac{-\sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{S}_K} F_{K, \sigma}(p, p) p_K}{\|p\|_{\mathcal{D}}^2} \right) \quad (54)$$

for each tested grid and with $\Delta t = 0.01$. Results are shown in Figures 4.

<figure><img src="image_6.png" /><figcaption>Figure 2. Space convergences</figcaption></figure>
cases. As expected a slope of one is obtained for all of them and for sufficiently small values of the time steps. Figure 3 displays the time convergence for the third case. We observe that reducing the time step does not improve the accuracy when h is too large but a convergence order of one can be retrieved at smaller space resolutions.

Although inequality (36) was assumed to be satisfied for the numerical analysis of the scheme, we have numerically checked that it holds. With that purpose, we have evaluated the following ratio

$$c_{\mathcal{D}}(p) = \inf_{t \in (0,T)} \left( \frac{-\sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{S}_K} F_{K,\sigma}(p,p) p_K}{\|p\|_{\mathcal{D}}^2} \right) \quad (54)$$

for each tested grid and with $\Delta t = 0.01$. Results are shown in Figures 4.

<figure><img src="image_6.png" /><figcaption>Figure 2. Space convergences</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 3. Time convergences for the mesh3 cases</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 4. Values of (54) observed for the three cases</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 3. Time convergences for the mesh3 cases</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 4. Values of (54) observed for the three cases</figcaption></figure>
## Appendix A.

**Lemma 9 (A discrete Poincaré inequality).** Let $\mathcal{D}$ be an admissible discretization of $\Omega$ in the sense of Definition 1 and $u \in Q_{\mathcal{D}}$. Then there exists a real number $\tilde{C}_{\mathcal{D}} > 0$ depending only on $\Omega$ and $\rho_5$ such that

$$ \|u_{\mathcal{D}}\|_{L^2(\Omega)} \le \tilde{C}_{\mathcal{D}} \|u_{\mathcal{D}}\|_{\mathcal{D}}. \qquad (55) $$

**Proof.** The proof of (55) is in fact a consequence of results established in [13]. Indeed, let us first notice that (see [13, equations (4.5)-(4.6)])

$$ \|u_{\mathcal{D}}\|_{1,2,\mathcal{D}}^2 = \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K} m_{\sigma} d_{K,\sigma} \left( \frac{D_{\sigma} u}{d_{\sigma}} \right)^2 \le \|u_{\mathcal{D}}\|_{\mathcal{D}}^2 \qquad (56) $$

with $d_{\sigma} = d_{K,\sigma} + d_{L,\sigma}$, $D_{\sigma}u = u_K - u_L$. From [13, Lemma 5.4], we deduce that there exists $C_A$ which only depends on $\Omega$ and $\rho_5$ such that

$$ \|u_{\mathcal{D}}\|_{L^4(\Omega)} \le C_A \|u_{\mathcal{D}}\|_{1,2,\mathcal{D}} $$

for $d=1$ or $d=2$ and

$$ \|u_{\mathcal{D}}\|_{L^6(\Omega)} \le C_A \|u_{\mathcal{D}}\|_{1,2,\mathcal{D}} $$

for $d=3$. On the other hand, Hölder inequality, for an integer value $q = 2q' \ge 2$, leads to

$$ \|u_{\mathcal{D}}\|_{L^2(\Omega)} \le \|u_{\mathcal{D}}\|_{L^q(\Omega)} |\Omega|^{\frac{q'-1}{2q'}} $$

where $|\Omega|$ denotes the measure of $\Omega$, which concludes the proof. □

## Declaration of interests

The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.

## References

[1] D. A. Di Pietro, R. Eymard, S. Lemaire, R. Masson, "Hybrid finite volume discretization of linear elasticity models on general meshes", in *Finite Volumes for Complex Applications VI – Problems & Perspectives* (J. Fořt, J. Fürst, J. Halama, R. Herbin, F. Hubert, eds.), Springer Proceedings in Mathematics, vol. 4, Springer, 2011, p. 331-339.

[2] L. Beirão Da Veiga, F. Brezzi, L. D. Marini, "Virtual elements for linear elasticity problems", *SIAM J. Numer. Anal.* **51** (2013), no. 2, p. 794-812.

[3] J. Coulet, I. Faille, V. Girault, N. Guy, F. Nataf, "A fully coupled scheme using virtual element method and finite volume for poroelasticity", *Comput. Geosci.* **24** (2020), p. 381-403.

[4] J. Droniou, "Finite volume schemes for diffusion equations: introduction to and review of modern methods", *Math. Models Methods Appl. Sci.* **24** (2014), no. 08, p. 1575-1619.

[5] M. Schneider, L. Agélas, G. Enchéry, B. Flemisch, "Convergence of nonlinear finite volume schemes for heterogeneous anisotropic diffusion on general meshes", *J. Comput. Phys.* **351** (2017), no. 15, p. 80-107.

[6] R. Herbin, F. Hubert, "Benchmark on discretization schemes for anisotropic diffusion problems on general grids", https://www.i2m.univ-amu.fr/fvca5/, 2008.

[7] M. A. Biot, "General theory of three-dimensional consolidation", *J. Appl. Phys.* **12** (1941), p. 155-164.

[8] R. E. Showalter, "Diffusion in porous-elastic media", *J. Math. Anal. Appl.* **251** (2000), no. 1, p. 310-340.

[9] S. C. Brenner, L. R. Scott, *The mathematical theory of finite element methods*, Texts in Applied Mathematics, Springer, 2008.

[10] L. Beirão Da Veiga, C. Lovadina, D. Mora, "A Virtual Element Method for elastic and inelastic problems on polytope meshes", *Comput. Methods Appl. Mech. Eng.* **295** (2015), p. 327-346.

[11] L. Beirão Da Veiga, F. Brezzi, A. Cangiani, G. Manzini, L. D. Marini, A. Russo, "Basic principles of virtual element methods", *Math. Models Methods Appl. Sci.* **23** (2013), no. 1, p. 199-214.
## Appendix A.

**Lemma 9 (A discrete Poincaré inequality).** Let $\mathcal{D}$ be an admissible discretization of $\Omega$ in the sense of Definition 1 and $u \in Q_{\mathcal{D}}$. Then there exists a real number $\tilde{C}_{\mathcal{D}} > 0$ depending only on $\Omega$ and $\rho_5$ such that

$$ \|u_{\mathcal{D}}\|_{L^2(\Omega)} \le \tilde{C}_{\mathcal{D}} \|u_{\mathcal{D}}\|_{\mathcal{D}}. \qquad (55) $$

**Proof.** The proof of (55) is in fact a consequence of results established in [13]. Indeed, let us first notice that (see [13, equations (4.5)-(4.6)])

$$ \|u_{\mathcal{D}}\|_{1,2,\mathcal{D}}^2 = \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K} m_{\sigma} d_{K,\sigma} \left( \frac{D_{\sigma} u}{d_{\sigma}} \right)^2 \le \|u_{\mathcal{D}}\|_{\mathcal{D}}^2 \qquad (56) $$

with $d_{\sigma} = d_{K,\sigma} + d_{L,\sigma}$, $D_{\sigma}u = u_K - u_L$. From [13, Lemma 5.4], we deduce that there exists $C_A$ which only depends on $\Omega$ and $\rho_5$ such that

$$ \|u_{\mathcal{D}}\|_{L^4(\Omega)} \le C_A \|u_{\mathcal{D}}\|_{1,2,\mathcal{D}} $$

for $d=1$ or $d=2$ and

$$ \|u_{\mathcal{D}}\|_{L^6(\Omega)} \le C_A \|u_{\mathcal{D}}\|_{1,2,\mathcal{D}} $$

for $d=3$. On the other hand, Hölder inequality, for an integer value $q = 2q' \ge 2$, leads to

$$ \|u_{\mathcal{D}}\|_{L^2(\Omega)} \le \|u_{\mathcal{D}}\|_{L^q(\Omega)} |\Omega|^{\frac{q'-1}{2q'}} $$

where $|\Omega|$ denotes the measure of $\Omega$, which concludes the proof. □

## Declaration of interests

The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.

## References

[1] D. A. Di Pietro, R. Eymard, S. Lemaire, R. Masson, "Hybrid finite volume discretization of linear elasticity models on general meshes", in *Finite Volumes for Complex Applications VI – Problems & Perspectives* (J. Fořt, J. Fürst, J. Halama, R. Herbin, F. Hubert, eds.), Springer Proceedings in Mathematics, vol. 4, Springer, 2011, p. 331-339.

[2] L. Beirão Da Veiga, F. Brezzi, L. D. Marini, "Virtual elements for linear elasticity problems", *SIAM J. Numer. Anal.* **51** (2013), no. 2, p. 794-812.

[3] J. Coulet, I. Faille, V. Girault, N. Guy, F. Nataf, "A fully coupled scheme using virtual element method and finite volume for poroelasticity", *Comput. Geosci.* **24** (2020), p. 381-403.

[4] J. Droniou, "Finite volume schemes for diffusion equations: introduction to and review of modern methods", *Math. Models Methods Appl. Sci.* **24** (2014), no. 08, p. 1575-1619.

[5] M. Schneider, L. Agélas, G. Enchéry, B. Flemisch, "Convergence of nonlinear finite volume schemes for heterogeneous anisotropic diffusion on general meshes", *J. Comput. Phys.* **351** (2017), no. 15, p. 80-107.

[6] R. Herbin, F. Hubert, "Benchmark on discretization schemes for anisotropic diffusion problems on general grids", https://www.i2m.univ-amu.fr/fvca5/, 2008.

[7] M. A. Biot, "General theory of three-dimensional consolidation", *J. Appl. Phys.* **12** (1941), p. 155-164.

[8] R. E. Showalter, "Diffusion in poro-elastic media", *J. Math. Anal. Appl.* **251** (2000), no. 1, p. 310-340.

[9] S. C. Brenner, L. R. Scott, *The mathematical theory of finite element methods*, Texts in Applied Mathematics, Springer, 2008.

[10] L. Beirão Da Veiga, C. Lovadina, D. Mora, "A Virtual Element Method for elastic and inelastic problems on polytope meshes", *Comput. Methods Appl. Mech. Eng.* **295** (2015), p. 327-346.

[11] L. Beirão Da Veiga, F. Brezzi, A. Cangiani, G. Manzini, L. D. Marini, A. Russo, "Basic principles of virtual element methods", *Math. Models Methods Appl. Sci.* **23** (2013), no. 1, p. 199-214.
[12] M. Schneider, B. Flemisch, R. Helmig, K. Terekhov, H. Tchelepi, "Monotone nonlinear finite-volume method for challenging grids", *Comput. Geosci.* **22** (2018), p. 565-586.

[13] R. Eymard, T. Gallouët, R. Herbin, "Discretisation of heterogeneous and anisotropic diffusion problems on general nonconforming meshes SUSHI: a scheme using stabilisation and hybrid interfaces", *IMA J. Numer. Anal.* **30** (2010), no. 4, p. 1009-1043.

[14] R. Eymard, T. Gallouët, R. Herbin, "Finite volume methods", in *Techniques of scientific computing (Part 3)*, Handbook of Numerical Analysis, vol. 7, Elsevier, 2000, p. 713-718.
[12] M. Schneider, B. Flemisch, R. Helmig, K. Terekhov, H. Tchelepi, "Monotone nonlinear finite-volume method for challenging grids", *Comput. Geosci.* **22** (2018), p. 565-586.

[13] R. Eymard, T. Gallouët, R. Herbin, "Discretisation of heterogeneous and anisotropic diffusion problems on general nonconforming meshes SUSHI: a scheme using stabilisation and hybrid interfaces", *IMA J. Numer. Anal.* **30** (2010), no. 4, p. 1009-1043.

[14] R. Eymard, T. Gallouët, R. Herbin, "Finite volume methods", in *Techniques of scientific computing (Part 3)*, Handbook of Numerical Analysis, vol. 7, Elsevier, 2000, p. 713-718.