OPTIMAL ERROR BOUNDS FOR THE TWO POINT
FLUX APPROXIMATION FINITE VOLUME
SCHEME

Robert Eymard, Thierry Gallouët, Raphaèle Herbin

► To cite this version:

Robert Eymard, Thierry Gallouët, Raphaele Herbin. OPTIMAL ERROR BOUNDS FOR THE TWO POINT FLUX APPROXIMATION FINITE VOLUME SCHEME. 2024. hal-04586601

HAL Id: hal-04586601

https://hal.science/hal-04586601

Preprint submitted on 27 May 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
OPTIMAL ERROR BOUNDS FOR THE TWO POINT
FLUX APPROXIMATION FINITE VOLUME
SCHEME

Robert Eymard, Thierry Gallouët, Raphaèle Herbin

► To cite this version:

Robert Eymard, Thierry Gallouët, Raphaele Herbin. OPTIMAL ERROR BOUNDS FOR THE TWO POINT FLUX APPROXIMATION FINITE VOLUME SCHEME. 2024. hal-04586601

HAL Id: hal-04586601

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

Preprint submitted on 27 May 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
# OPTIMAL ERROR BOUNDS FOR THE TWO POINT FLUX APPROXIMATION FINITE VOLUME SCHEME

R. EYMARD, T. GALLOUËT, AND R. HERBIN

ABSTRACT. We consider a finite volume scheme with two-point flux approximation (TPFA) to approximate a Laplace problem when the solution exhibits no more regularity than belonging to $H_0^1(\Omega)$. We establish in this case some error bounds for both the solution and the approximation of the gradient component orthogonal to the mesh faces. This estimate is optimal, in the sense that the approximation error has the same order as that of the sum of the interpolation error and a conformity error. A numerical example illustrates the error estimate in the context of a solution with minimal regularity. This result is extended to evolution problems discretized via the implicit Euler scheme in an appendix.

## CONTENTS

<table><tr><td>1. Introduction</td><td>1</td></tr><tr><td>2. The finite volume scheme for the Laplace problem</td><td>4</td></tr><tr><td>3. Convergence analysis</td><td>7</td></tr><tr><td>3.1. An optimal error estimate</td><td>7</td></tr><tr><td>3.2. Convergence of the interpolation and conformity errors</td><td>12</td></tr><tr><td>3.3. Error estimate in the $H^2(\Omega)$ case</td><td>14</td></tr><tr><td>4. Numerical example of a minimal regularity problem</td><td>16</td></tr><tr><td>5. Conclusion and perspectives</td><td>19</td></tr><tr><td>Appendix A. The transient case</td><td>20</td></tr><tr><td>A.1. Description of the implicit Euler scheme</td><td>21</td></tr><tr><td>A.2. An optimal error estimate</td><td>22</td></tr><tr><td>References</td><td>27</td></tr></table>

## 1. INTRODUCTION

Finite volume methods are standardly used for the approximation of elliptic and parabolic problems in several physical or engineering frameworks such as fluid mechanics, reservoir simulation, heat and mass transfer, biology, biomedical research... Here we are more specifically interested in the two-point flow approximation (TPFA) finite volume scheme for the approximation of the Laplace operator; this is a very popular scheme that has been used since the 60's in the various above mentioned applications, see e.g. [3],[6] to cite only a few of these.

<footnote n="2010">2010 Mathematics Subject Classification. 65N30,35K15,47A07.</footnote>

<footnote>Key words and phrases. linear elliptic problem with minimal regularity, optimal error estimate, finite volume method, linear parabolic problem.</footnote>
# OPTIMAL ERROR BOUNDS FOR THE TWO POINT FLUX APPROXIMATION FINITE VOLUME SCHEME

R. EYMARD, T. GALLOUËT, AND R. HERBIN

ABSTRACT. We consider a finite volume scheme with two-point flux approximation (TPFA) to approximate a Laplace problem when the solution exhibits no more regularity than belonging to $H_0^1(\Omega)$. We establish in this case some error bounds for both the solution and the approximation of the gradient component orthogonal to the mesh faces. This estimate is optimal, in the sense that the approximation error has the same order as that of the sum of the interpolation error and a conformity error. A numerical example illustrates the error estimate in the context of a solution with minimal regularity. This result is extended to evolution problems discretized via the implicit Euler scheme in an appendix.

## CONTENTS

<table><tr><td>1. Introduction</td><td>1</td></tr><tr><td>2. The finite volume scheme for the Laplace problem</td><td>4</td></tr><tr><td>3. Convergence analysis</td><td>7</td></tr><tr><td>3.1. An optimal error estimate</td><td>7</td></tr><tr><td>3.2. Convergence of the interpolation and conformity errors</td><td>12</td></tr><tr><td>3.3. Error estimate in the $H^2(\Omega)$ case</td><td>14</td></tr><tr><td>4. Numerical example of a minimal regularity problem</td><td>16</td></tr><tr><td>5. Conclusion and perspectives</td><td>19</td></tr><tr><td>Appendix A. The transient case</td><td>20</td></tr><tr><td>A.1. Description of the implicit Euler scheme</td><td>21</td></tr><tr><td>A.2. An optimal error estimate</td><td>22</td></tr><tr><td>References</td><td>27</td></tr></table>

## 1. INTRODUCTION

Finite volume methods are standardly used for the approximation of elliptic and parabolic problems in several physical or engineering frameworks such as fluid mechanics, reservoir simulation, heat and mass transfer, biology, biomedical research... Here we are more specifically interested in the two-point flow approximation (TPFA) finite volume scheme for the approximation of the Laplace operator; this is a very popular scheme that has been used since the 60's in the various above mentioned applications, see e.g. [3],[6] to cite only a few of these.

<footnote n="2010">2010 Mathematics Subject Classification. 65N30,35K15,47A07.</footnote>

<footnote>Key words and phrases. linear elliptic problem with minimal regularity, optimal error estimate, finite volume method, linear parabolic problem.</footnote>
Important features of the TPFA scheme are $L^\infty$ stability and monotony [18]. Thanks to these properties, stable and convergent numerical schemes based on TPFA were designed for nonlinear problems [5] or linear problems with singular source terms [11].

A mathematical proof of convergence of the TPFA scheme for the Laplace equation was first given for rectangular grids in [20], and the first comprehensive mathematical analysis of this scheme fore more general grids (the so called admissible grids) dates back to 2000 [18]. Therein, under suitable assumptions on the physical domain and the right hand side, it is proven that if the exact solution of the Poisson equation posed on a open polytopal bounded subset of $\mathbb{R}^d$ with $d = 2$ or 3 is assumed to be in $H^2$, then an error estimate of order 1 can be obtained. Superconvergence has been observed and proven for a class of triangular meshes, see [13] and Remark 3.13 below.

Our aim in the present paper is to give an error estimate in the case where the right-hand side belongs to the dual space of the energy space associated to the equation, namely $H^{-1}$ and without any supplementary regularity assumption on the exact solution other than the natural regularity given by the problem, namely $H_0^1$. The case of a right-hand side in $H^{-1}$ had already been studied in [12], where the convergence of the scheme is proven for a non coercive elliptic operator. In the present paper, we study the same scheme for the Laplace equation, and prove an optimal bound on the error; more precisely we show that the discretization error is bounded by below and by above, in each case up to a constant, by the sum of two errors, namely the interpolation error and the conformity error, whose precise definitions are given in the sequel. We also obtain a bound of the error between an approximate gradient of the solution and the gradient of the exact solution. As far as we know, both results are original.

In order to introduce these results, let us first give some definitions. Consider $f \in L^2(\Omega)$ and $\Omega$ an open polytopal bounded subset of $\mathbb{R}^d$ with $d = 2$ or 3; denote by $\bar{u} \in H_0^1(\Omega)$ the solution of the problem

$$ \forall v \in H_0^1(\Omega), \quad \int_{\Omega} \nabla \bar{u}(x) \cdot \nabla v(x) dx = \int_{\Omega} f(x) v(x) dx. \qquad (1.1) $$

As above mentioned, if $\bar{u} \in H^2(\Omega)$, there exists $C > 0$ only depending on a regularity factor of the finite volume mesh $\mathcal{T}$ such that the solution $u$ of the finite volume scheme satisfies

$$ \|u - \bar{u}\|_{L^2(\Omega)} \le C h_{\mathcal{T}} \|\bar{u}\|_{H^2(\Omega)}, \qquad (1.2) $$

where $h_{\mathcal{T}}$ is the size of the mesh.

Now when seeking an approximate solution $u_h \in V_h$ of the same Problem (1.1) by a conforming finite element method, which means that

$$ \forall v_h \in V_h, \quad \int_{\Omega} \nabla v_h(x) \cdot \nabla v_h(x) dx = \int_{\Omega} f(x) v_h(x) dx, \qquad (1.3) $$

we get that

$$ \|\bar{u} - u_h\|_{H_0^1(\Omega)} = \inf_{v_h \in V_h} \|\bar{u} - v_h\|_{H_0^1(\Omega)}, $$

where:

* $\| \cdot \|_{H_0^1(\Omega)}$ is the energy norm associated to Problem (1.1), that is to say $\|v\|_{H_0^1(\Omega)} = \int_{\Omega} |\nabla v(x)|^2 dx$,

* $V_h \subset H_0^1(\Omega)$ is the (finite dimensional) finite element space.

This result, which is a particular case of Céa's lemma for the problem at hand, states that the energy norm of the discretization error is equal to the energy norm of the interpolation error; it
Important features of the TPFA scheme are $L^\infty$ stability and monotony [18]. Thanks to these properties, stable and convergent numerical schemes based on TPFA were designed for nonlinear problems [5] or linear problems with singular source terms [11].

A mathematical proof of convergence of the TPFA scheme for the Laplace equation was first given for rectangular grids in [20], and the first comprehensive mathematical analysis of this scheme fore more general grids (the so called admissible grids) dates back to 2000 [18]. Therein, under suitable assumptions on the physical domain and the right hand side, it is proven that if the exact solution of the Poisson equation posed on a open polytopal bounded subset of $\mathbb{R}^d$ with $d = 2$ or 3 is assumed to be in $H^2$, then an error estimate of order 1 can be obtained. Superconvergence has been observed and proven for a class of triangular meshes, see [13] and Remark 3.13 below.

Our aim in the present paper is to give an error estimate in the case where the right-hand side belongs to the dual space of the energy space associated to the equation, namely $H^{-1}$ and without any supplementary regularity assumption on the exact solution other than the natural regularity given by the problem, namely $H_0^1$. The case of a right-hand side in $H^{-1}$ had already been studied in [12], where the convergence of the scheme is proven for a non coercive elliptic operator. In the present paper, we study the same scheme for the Laplace equation, and prove an optimal bound on the error; more precisely we show that the discretization error is bounded by below and by above, in each case up to a constant, by the sum of two errors, namely the interpolation error and the conformity error, whose precise definitions are given in the sequel. We also obtain a bound of the error between an approximate gradient of the solution and the gradient of the exact solution. As far as we know, both results are original.

In order to introduce these results, let us first give some definitions. Consider $f \in L^2(\Omega)$ and $\Omega$ an open polytopal bounded subset of $\mathbb{R}^d$ with $d = 2$ or 3; denote by $\bar{u} \in H_0^1(\Omega)$ the solution of the problem

$$ \forall v \in H_0^1(\Omega), \quad \int_{\Omega} \nabla \bar{u}(x) \cdot \nabla v(x)dx = \int_{\Omega} f(x)v(x)dx. \quad (1.1) $$

As above mentioned, if $\bar{u} \in H^2(\Omega)$, there exists $C > 0$ only depending on a regularity factor of the finite volume mesh $\mathcal{T}$ such that the solution $u$ of the finite volume scheme satisfies

$$ \|u - \bar{u}\|_{L^2(\Omega)} \le C h_{\mathcal{T}} \|\bar{u}\|_{H^2(\Omega)}, \quad (1.2) $$

where $h_{\mathcal{T}}$ is the size of the mesh.

Now when seeking an approximate solution $u_h \in V_h$ of the same Problem (1.1) by a conforming finite element method, which means that

$$ \forall v_h \in V_h, \quad \int_{\Omega} \nabla v_h(x) \cdot \nabla v_h(x)dx = \int_{\Omega} f(x)v_h(x)dx, \quad (1.3) $$

we get that

$$ \|\bar{u} - u_h\|_{H_0^1(\Omega)} = \inf_{v_h \in V_h} \|\bar{u} - v_h\|_{H_0^1(\Omega)}, $$

where:

* $\| \cdot \|_{H_0^1(\Omega)}$ is the energy norm associated to Problem (1.1), that is to say $\|v\|_{H_0^1(\Omega)} = \int_{\Omega} |\nabla v(x)|^2 dx$,

* $V_h \subset H_0^1(\Omega)$ is the (finite dimensional) finite element space.

This result, which is a particular case of Céa's lemma for the problem at hand, states that the energy norm of the discretization error is equal to the energy norm of the interpolation error; it
holds without any further regularity assumption than $\bar{u} \in H_0^1(\Omega)$. It is optimal in the sense that the order of convergence is given by the order of the interpolation error $\inf_{v \in V_h} \|\bar{u} - v\|_{H_0^1(\Omega)}$.

Such an optimal error bound is extended in [7] to some nonconforming schemes that fall in the framework of the Gradient Discretization Method (GDM). The GDM was invented and analysed some time ago in order to construct a framework for a number of schemes that have been devised in the past twenty years in order to deal with anisotropic diffusion problems and/or distorted meshes. Let us briefly recall the basics of the method. Let $X_\tau$ be a finite dimensional vector space of the degrees of freedom generated by a mesh $\mathcal{T}$, $\Pi_\tau v$ denote a function that is reconstructed from any $v \in X_\tau$ and defined a.e. in $\Omega$, and $\nabla_\tau v$ be an approximate gradient, also reconstructed from $v$. Then Problem (1.1) is approximated by $u \in X_\tau$ such that

$$ \forall v \in X_{\tau}, \quad \int_{\Omega} \nabla_{\tau} u(x) \cdot \nabla_{\tau} v(x) dx = \int_{\Omega} f(x) \Pi_{\tau} v(x) dx. \quad (1.4) $$

Note that for non conforming schemes, $\Pi_\tau(X_\tau) \not\subset H_0^1(\Omega)$, the approximate gradient $\nabla_\tau v$ cannot be defined as the continuous gradient of $\Pi_\tau v$. This difficulty is overcome by defining the conformity error of the scheme by

$$ \zeta_{\tau}(\nabla \bar{u}) = \max_{v \in X_{\tau} \setminus \{0\}} \frac{\langle \nabla \bar{u}, \nabla_{\tau} v \rangle_{L^2} + \langle \operatorname{div}(\nabla \bar{u}), \Pi_{\tau} v \rangle_{L^2}}{\| \nabla_{\tau} v \|_{L^2}}. \quad (1.5) $$

An optimal error bound (see [7, Theorem 2.28]) may then be obtained in the spirit of second Strang's lemma [26]; introducing

$$ \delta_{\tau}(\bar{u}, v)^2 = \|\bar{u} - \Pi_{\tau} v\|_{L^2}^2 + \|\nabla \bar{u} - \nabla_{\tau} v\|_{L^2}^2, \quad (1.6) $$

this bound states that the error $\delta_{\tau}(\bar{u}, u)$ committed on the solution and its gradient satisfies

$$ \frac{1}{2} \left( \zeta_{\tau}(\nabla \bar{u}) + \inf_{v \in X_{\tau}} \delta_{\tau}(\bar{u}, v) \right) \le \delta_{\tau}(\bar{u}, u) \le C \left( \zeta_{\tau}(\nabla \bar{u}) + \inf_{v \in X_{\tau}} \delta_{\tau}(\bar{u}, v) \right). \quad (1.7) $$

without requiring any more regularity than $\bar{u} \in H_0^1(\Omega)$ on the exact solution. This bound is said to be optimal in the sense that the order of the approximation error $\delta_{\tau}(\bar{u}, u)$ is the same as that of the sum of the interpolation error $\inf_{v \in X_\tau} \delta_{\tau}(\bar{u}, v)$ and the conformity error $\zeta_{\tau}(\nabla \bar{u})$.

In some particular situations (rectangular or acute triangular cells), the TPFA FV scheme can be shown to be a GDM (see [7, Lemma 13.20]), and therefore the optimal bound (1.7) holds. However, in the case of general admissible meshes of Definition 2.1 below, this is no more the case; for instance the TPFA FV scheme on Voronoï meshes [16] cannot be seen as a GDM.

The aim of the present paper is to prove an error estimate for the TPFA finite volume scheme which

- is optimal in the sense of the bound (1.7),

- still holds for meshes which only satisfy the assumptions of Definition 2.1,

- applies to general right-hand sides in $H^{-1}(\Omega)$ so that the exact solution belongs only to $H_0^1(\Omega)$,

- and is formulated through a stronger norm than the $L^2(\Omega)$ norm.

To this purpose, we first propose to use an approximate gradient in the right-hand side (the so-called “inflated approximate gradient”) of the scheme in order to deal with for right-hand sides in $H^{-1}(\Omega)$. Unfortunately this approximate gradient cannot be used for a GDM scheme, due to the fact that it can only weakly converge. A variational formulation of the TPFA FV scheme is then obtained by considering only the normal component of the approximate gradient. We may then
holds without any further regularity assumption than $\bar{u} \in H_0^1(\Omega)$. It is optimal in the sense that the order of convergence is given by the order of the interpolation error $\inf_{v \in V_h} \|\bar{u} - v\|_{H_0^1(\Omega)}$.

Such an optimal error bound is extended in [7] to some nonconforming schemes that fall in the framework of the Gradient Discretization Method (GDM). The GDM was invented and analysed some time ago in order to construct a framework for a number of schemes that have been devised in the past twenty years in order to deal with anisotropic diffusion problems and/or distorted meshes. Let us briefly recall the basics of the method. Let $X_\tau$ be a finite dimensional vector space of the degrees of freedom generated by a mesh $\mathcal{T}$, $\Pi_\tau v$ denote a function that is reconstructed from any $v \in X_\tau$ and defined a.e. in $\Omega$, and $\nabla_\tau v$ be an approximate gradient, also reconstructed from $v$. Then Problem (1.1) is approximated by $u \in X_\tau$ such that

$$ \forall v \in X_{\tau}, \quad \int_{\Omega} \nabla_{\tau} u(x) \cdot \nabla_{\tau} v(x) dx = \int_{\Omega} f(x) \Pi_{\tau} v(x) dx. \quad (1.4) $$

Note that for non conforming schemes, $\Pi_\tau(X_\tau) \not\subset H_0^1(\Omega)$, the approximate gradient $\nabla_\tau v$ cannot be defined as the continuous gradient of $\Pi_\tau v$. This difficulty is overcome by defining the conformity error of the scheme by

$$ \zeta_{\tau}(\nabla \bar{u}) = \max_{v \in X_{\tau} \setminus \{0\}} \frac{\langle \nabla \bar{u}, \nabla_{\tau} v \rangle_{L^2} + \langle \operatorname{div}(\nabla \bar{u}), \Pi_{\tau} v \rangle_{L^2}}{\| \nabla_{\tau} v \|_{L^2}}. \quad (1.5) $$

An optimal error bound (see [7, Theorem 2.28]) may then be obtained in the spirit of second Strang's lemma [26]; introducing

$$ \delta_{\tau}(\bar{u}, v)^2 = \|\bar{u} - \Pi_{\tau} v\|_{L^2}^2 + \|\nabla \bar{u} - \nabla_{\tau} v\|_{L^2}^2, \quad (1.6) $$

this bound states that the error $\delta_{\tau}(\bar{u}, u)$ committed on the solution and its gradient satisfies

$$ \frac{1}{2} \left( \zeta_{\tau}(\nabla \bar{u}) + \inf_{v \in X_{\tau}} \delta_{\tau}(\bar{u}, v) \right) \le \delta_{\tau}(\bar{u}, u) \le C \left( \zeta_{\tau}(\nabla \bar{u}) + \inf_{v \in X_{\tau}} \delta_{\tau}(\bar{u}, v) \right). \quad (1.7) $$

without requiring any more regularity than $\bar{u} \in H_0^1(\Omega)$ on the exact solution. This bound is said to be optimal in the sense that the order of the approximation error $\delta_{\tau}(\bar{u}, u)$ is the same as that of the sum of the interpolation error $\inf_{v \in X_\tau} \delta_{\tau}(\bar{u}, v)$ and the conformity error $\zeta_{\tau}(\nabla \bar{u})$.

In some particular situations (rectangular or acute triangular cells), the TPFA FV scheme can be shown to be a GDM (see [7, Lemma 13.20]), and therefore the optimal bound (1.7) holds. However, in the case of general admissible meshes of Definition 2.1 below, this is no more the case; for instance the TPFA FV scheme on Voronoï meshes [16] cannot be seen as a GDM.

The aim of the present paper is to prove an error estimate for the TPFA finite volume scheme which

- is optimal in the sense of the bound (1.7),

- still holds for meshes which only satisfy the assumptions of Definition 2.1,

- applies to general right-hand sides in $H^{-1}(\Omega)$ so that the exact solution belongs only to $H_0^1(\Omega)$,

- and is formulated through a stronger norm than the $L^2(\Omega)$ norm.

To this purpose, we first propose to use an approximate gradient in the right-hand side (the so-called “inflated approximate gradient”) of the scheme in order to deal with for right-hand sides in $H^{-1}(\Omega)$. Unfortunately this approximate gradient cannot be used for a GDM scheme, due to the fact that it can only weakly converge. A variational formulation of the TPFA FV scheme is then obtained by considering only the normal component of the approximate gradient. We may then
follow the line of thought of [7, Theorem 2.28] and obtain an optimal error bound involving an interpolation error between the approximation of the normal gradient and the normal component of the gradient of the continuous solution. This result is, to our knowledge, original for at least these two reasons: it does not require more regularity on the exact solution than the natural regularity obtained from the continuous problem, and it leads to an error estimate between a consistent reconstruction of an approximate gradient (necessarily different from the inflated approximate gradient) and the full gradient of the continuous solution.

We then apply these two results in the case where the continuous solution is in $H^2$, extending the results of [18] to the convergence of the consistent approximate gradient.

A numerical example in the case of a solution with minimal regularity illustrates in Section 4 the optimality of the error estimate; indeed the considered solution is unbounded and does not belong to any $W_0^{1,p}(\Omega)$, for $p > 2$.

A short conclusion follows.

In Appendix A, following the methods of proof provided by [2] and [8], we extend the optimal error bound obtained for the Laplace problem to the transient heat equation, using an implicit Euler time discretization.

## 2. THE FINITE VOLUME SCHEME FOR THE LAPLACE PROBLEM

Let $\Omega \subset \mathbb{R}^d$ be an open polytopal domain (with $d = 2$ or $d = 3$), with boundary $\partial\Omega$, and $L \in H^{-1}(\Omega)$. We consider the Dirichlet problem

$$ - \Delta \bar{u} = L \qquad (2.1) $$

$$ \bar{u} = 0 \text{ on } \partial\Omega. \qquad (2.2) $$

It is wellknown that there exists a unique solution to this problem, which satisfies

$$ \bar{u} \in H_0^1(\Omega) \\ \int_{\Omega} \nabla \bar{u} \cdot \nabla \varphi \, dx = \langle L, \varphi \rangle_{H^{-1}(\Omega), H_0^1(\Omega)}, \quad \forall \varphi \in H_0^1(\Omega). $$

It is also wellknown that any linear form $L \in H^{-1}(\Omega)$ may be decomposed as $L = f + \operatorname{div} \mathbf{F}$ with $f \in L^2(\Omega)$ and $\mathbf{F} \in L^2(\Omega)^d$. Since $\langle \operatorname{div} \mathbf{F}, \varphi \rangle_{H^{-1}(\Omega), H_0^1(\Omega)} = - \int_{\Omega} \mathbf{F} \cdot \nabla \varphi \, dx$, the above weak formulation may be recast as

$$ \bar{u} \in H_0^1(\Omega), \quad \forall \varphi \in H_0^1(\Omega), \qquad (2.3) $$

$$ \int_{\Omega} \nabla \bar{u}(x) \cdot \nabla \varphi(x) dx = \int_{\Omega} (f(x)\varphi(x) - \mathbf{F}(x) \cdot \nabla \varphi(x)) dx. \qquad (2.4) $$

Observe that

$$ \nabla \bar{u} + \mathbf{F} \in \mathbf{H}_{\operatorname{div}}(\Omega) \text{ with } -\operatorname{div}(\nabla \bar{u} + \mathbf{F}) = f \text{ a.e in } \Omega. \qquad (2.5) $$

The remaining part of this section is dedicated to the definition of a finite volume for the approximation of this problem.

Let $\mathcal{T}$ be an admissible mesh of $\Omega$ in the following sense, close to that given in [18, Definition 3.1].
follow the line of thought of [7, Theorem 2.28] and obtain an optimal error bound involving an interpolation error between the approximation of the normal gradient and the normal component of the gradient of the continuous solution. This result is, to our knowledge, original for at least these two reasons: it does not require more regularity on the exact solution than the natural regularity obtained from the continuous problem, and it leads to an error estimate between a consistent reconstruction of an approximate gradient (necessarily different from the inflated approximate gradient) and the full gradient of the continuous solution.

We then apply these two results in the case where the continuous solution is in $H^2$, extending the results of [18] to the convergence of the consistent approximate gradient.

A numerical example in the case of a solution with minimal regularity illustrates in Section 4 the optimality of the error estimate; indeed the considered solution is unbounded and does not belong to any $W_0^{1,p}(\Omega)$, for $p > 2$.

A short conclusion follows.

In Appendix A, following the methods of proof provided by [2] and [8], we extend the optimal error bound obtained for the Laplace problem to the transient heat equation, using an implicit Euler time discretization.

## 2. THE FINITE VOLUME SCHEME FOR THE LAPLACE PROBLEM

Let $\Omega \subset \mathbb{R}^d$ be an open polytopal domain (with $d = 2$ or $d = 3$), with boundary $\partial\Omega$, and $L \in H^{-1}(\Omega)$. We consider the Dirichlet problem

$$ - \Delta \bar{u} = L \qquad (2.1) $$

$$ \bar{u} = 0 \text{ on } \partial\Omega. \qquad (2.2) $$

It is wellknown that there exists a unique solution to this problem, which satisfies

$$ \bar{u} \in H_0^1(\Omega) \\ \int_{\Omega} \nabla \bar{u} \cdot \nabla \varphi \, dx = \langle L, \varphi \rangle_{H^{-1}(\Omega), H_0^1(\Omega)}, \quad \forall \varphi \in H_0^1(\Omega). $$

It is also wellknown that any linear form $L \in H^{-1}(\Omega)$ may be decomposed as $L = f + \operatorname{div} \mathbf{F}$ with $f \in L^2(\Omega)$ and $\mathbf{F} \in L^2(\Omega)^d$. Since $\langle \operatorname{div} \mathbf{F}, \varphi \rangle_{H^{-1}(\Omega), H_0^1(\Omega)} = - \int_{\Omega} \mathbf{F} \cdot \nabla \varphi \, dx$, the above weak formulation may be recast as

$$ \bar{u} \in H_0^1(\Omega), \quad \forall \varphi \in H_0^1(\Omega), \qquad (2.3) $$

$$ \int_{\Omega} \nabla \bar{u}(x) \cdot \nabla \varphi(x) dx = \int_{\Omega} (f(x)\varphi(x) - \mathbf{F}(x) \cdot \nabla \varphi(x)) dx. \qquad (2.4) $$

Observe that

$$ \nabla \bar{u} + \mathbf{F} \in \mathbf{H}_{\operatorname{div}}(\Omega) \text{ with } -\operatorname{div}(\nabla \bar{u} + \mathbf{F}) = f \text{ a.e in } \Omega. \qquad (2.5) $$

The remaining part of this section is dedicated to the definition of a finite volume for the approximation of this problem.

Let $\mathcal{T}$ be an admissible mesh of $\Omega$ in the following sense, close to that given in [18, Definition 3.1].
**Definition 2.1** (Admissible meshes). An admissible finite volume mesh of $\Omega$, denoted by $\mathcal{T}$, is given by a finite family of “control volumes”, which are disjoint open polytopal convex subsets of $\Omega$, a finite family of disjoint subsets of $\bar{\Omega}$ contained in hyperplanes of $\mathbb{R}^d$, denoted by $\mathcal{E}$ (these are the edges (two-dimensional) or sides (three-dimensional) of the control volumes), with strictly positive $(d-1)$-dimensional measure, and a family of points of $\Omega$ denoted by $\mathcal{P}$ satisfying the following properties (in fact, we shall denote, somewhat incorrectly, by $\mathcal{T}$ the family of control volumes):

(i) The closure of the union of all the control volumes is $\bar{\Omega}$;

(ii) For any $K \in \mathcal{T}$, there exists a subset $\mathcal{E}_K$ of $\mathcal{E}$ such that $\partial K = \bar{K} \setminus K = \cup_{\sigma \in \mathcal{E}_K} \bar{\sigma}$. Furthermore, $\mathcal{E} = \cup_{K \in \mathcal{T}} \mathcal{E}_K$.

(iii) The family $\mathcal{P} = (x_K)_{K \in \mathcal{T}}$ is such that $x_K \in K$ for all $K \in \mathcal{T}$.

(iv) The set $\mathcal{E}$ is partitioned into $\mathcal{E} = \mathcal{E}_{\text{ext}} \cup \mathcal{E}_{\text{int}}$. For any $\sigma \in \mathcal{E}$:

- either $\sigma \in \mathcal{E}_{\text{ext}}$; then there exists exactly one $K \in \mathcal{T}$ such that $\sigma \in \mathcal{E}_K$, $\sigma \subset \partial\Omega$ and the straight line $D_{K,\sigma}$ going through $x_K$ and orthogonal to $\sigma$ is such that $D_{K,\sigma} \cap \sigma \neq \emptyset$;

- or $\sigma \in \mathcal{E}_{\text{int}}$; then there exist exactly two elements of $\mathcal{T}$ denoted $K$ and $L$ such that $\sigma \in \mathcal{E}_K \cap \mathcal{E}_L$, $\partial K \cap \partial L = \bar{\sigma}$ and the straight line $D_{K,L}$ going through $x_K$ and $x_L$ is orthogonal to $\sigma$ ($D_{K,L} \cap \bar{\sigma} = \emptyset$ is not excluded); in this case $\sigma$ is denoted by $K|L$.

Figure 1 shows an example of mesh with a few notations.

<figure><img src="image_2.png" /><figcaption>FIGURE 1. Two neighbouring control volumes of an admissible mesh.</figcaption></figure>

The following notations are used

• $|A|$ the $d$ (resp. $(d-1)$)-dimensional measure of any subset $A$ of $\mathbb{R}^d$ (resp. $\mathbb{R}^{d-1}$).

• $\bar{x}_{\sigma}$: center of gravity of $\sigma \in \mathcal{E}$,

• $d_{K,\sigma} > 0$: orthogonal distance from $x_K$ to $\sigma$, for $K \in \mathcal{T}$ and $\sigma \in \mathcal{E}_K$,

• $D_{K,\sigma}$: cone with vertex $x_K$ and basis $\sigma$, for $K \in \mathcal{T}$ and $\sigma \in \mathcal{E}_K$, so that $|D_{K,\sigma}| = \frac{|\sigma|d_{K,\sigma}}{d}$,

• $\mathbf{n}_{K,\sigma}$: unit vector, normal to $\sigma$ and outward to $K$

• $h_{\mathcal{T}} = \sup\{\text{diam}(K), K \in \mathcal{T}\}$: mesh size,

• $\theta_{\mathcal{T}} = \inf\{\frac{d_{K,\sigma}}{\text{diam}(K)}, K \in \mathcal{T}, \sigma \in \mathcal{E}_K\}$: mesh regularity parameter,

• $X_{\mathcal{T}}$: set of all real families $u = ((u_K)_{K \in \mathcal{T}}, (u_{\sigma})_{\sigma \in \mathcal{E}})$ such that $u_{\sigma} = 0$ if $\sigma \in \mathcal{E}_{\text{ext}}$.
**Definition 2.1** (Admissible meshes). An admissible finite volume mesh of $\Omega$, denoted by $\mathcal{T}$, is given by a finite family of “control volumes”, which are disjoint open polytopal convex subsets of $\Omega$, a finite family of disjoint subsets of $\bar{\Omega}$ contained in hyperplanes of $\mathbb{R}^d$, denoted by $\mathcal{E}$ (these are the edges (two-dimensional) or sides (three-dimensional) of the control volumes), with strictly positive $(d-1)$-dimensional measure, and a family of points of $\Omega$ denoted by $\mathcal{P}$ satisfying the following properties (in fact, we shall denote, somewhat incorrectly, by $\mathcal{T}$ the family of control volumes):

(i) The closure of the union of all the control volumes is $\bar{\Omega}$;

(ii) For any $K \in \mathcal{T}$, there exists a subset $\mathcal{E}_K$ of $\mathcal{E}$ such that $\partial K = \bar{K} \setminus K = \cup_{\sigma \in \mathcal{E}_K} \bar{\sigma}$. Furthermore, $\mathcal{E} = \cup_{K \in \mathcal{T}} \mathcal{E}_K$.

(iii) The family $\mathcal{P} = (x_K)_{K \in \mathcal{T}}$ is such that $x_K \in K$ for all $K \in \mathcal{T}$.

(iv) The set $\mathcal{E}$ is partitioned into $\mathcal{E} = \mathcal{E}_{\text{ext}} \cup \mathcal{E}_{\text{int}}$. For any $\sigma \in \mathcal{E}$:

- either $\sigma \in \mathcal{E}_{\text{ext}}$; then there exists exactly one $K \in \mathcal{T}$ such that $\sigma \in \mathcal{E}_K$, $\sigma \subset \partial\Omega$ and the straight line $D_{K,\sigma}$ going through $x_K$ and orthogonal to $\sigma$ is such that $D_{K,\sigma} \cap \sigma \neq \emptyset$;

- or $\sigma \in \mathcal{E}_{\text{int}}$; then there exist exactly two elements of $\mathcal{T}$ denoted $K$ and $L$ such that $\sigma \in \mathcal{E}_K \cap \mathcal{E}_L$, $\partial K \cap \partial L = \bar{\sigma}$ and the straight line $D_{K,L}$ going through $x_K$ and $x_L$ is orthogonal to $\sigma$ ($D_{K,L} \cap \bar{\sigma} = \emptyset$ is not excluded); in this case $\sigma$ is denoted by $K|L$.

Figure 1 shows an example of mesh with a few notations.

<figure><img src="image_2.png" /><figcaption>FIGURE 1. Two neighbouring control volumes of an admissible mesh.</figcaption></figure>

The following notations are used

• $|A|$ the $d$ (resp. $(d-1)$)-dimensional measure of any subset $A$ of $\mathbb{R}^d$ (resp. $\mathbb{R}^{d-1}$).

• $\bar{x}_{\sigma}$: center of gravity of $\sigma \in \mathcal{E}$,

• $d_{K,\sigma} > 0$: orthogonal distance from $x_K$ to $\sigma$, for $K \in \mathcal{T}$ and $\sigma \in \mathcal{E}_K$,

• $D_{K,\sigma}$: cone with vertex $x_K$ and basis $\sigma$, for $K \in \mathcal{T}$ and $\sigma \in \mathcal{E}_K$, so that $|D_{K,\sigma}| = \frac{|\sigma|d_{K,\sigma}}{d}$,

• $\mathbf{n}_{K,\sigma}$: unit vector, normal to $\sigma$ and outward to $K$

• $h_{\mathcal{T}} = \sup\{\text{diam}(K), K \in \mathcal{T}\}$: mesh size,

• $\theta_{\mathcal{T}} = \inf\{\frac{d_{K,\sigma}}{\text{diam}(K)}, K \in \mathcal{T}, \sigma \in \mathcal{E}_K\}$: mesh regularity parameter,

• $X_{\mathcal{T}}$: set of all real families $u = ((u_K)_{K \in \mathcal{T}}, (u_{\sigma})_{\sigma \in \mathcal{E}})$ such that $u_{\sigma} = 0$ if $\sigma \in \mathcal{E}_{\text{ext}}$.
With a slight abuse of notation, for any $u \in X_{\mathcal{T}}$, we also denote by $u$ the element of $L^2(\Omega)$ which is a.e. equal to the value $u_K$ in any $K \in \mathcal{T}$.

Let us observe that this definition leads to admissible meshes in the sense of [18, Definition 3.1], but is slightly more restrictive since, for the sake of simplicity and as in [12], we do not allow $x_K \in \bar{\sigma}$ if $\sigma \in \mathcal{E}_K$. There are ways to relax this assumption by eliminating the face unknowns in the scheme, but this leads to additional technical difficulties.

As mentioned in the introduction, our analysis requires several discrete derivation operators which we now precisely define.

**Definition 2.2** (Discrete derivative and gradients).

• $G_{\mathcal{T}}$: normal discrete derivative

$$ G_{\mathcal{T}}:X_{\mathcal{T}} \to L^2(\Omega) \\ u \mapsto G_{\mathcal{T}}u, \text{ with } G_{\mathcal{T}}u(x) = \frac{u_{\sigma} - u_K}{d_{K,\sigma}} \text{ for a.e. } x \in D_{K,\sigma}. \quad (2.6) $$

• $\nabla_{\mathcal{T}}$: inflated discrete gradient

$$ \nabla_{\mathcal{T}}:X_{\mathcal{T}} \to (L^2(\Omega))^d \\ u \mapsto \nabla_{\mathcal{T}}u, \text{ with } \nabla_{\mathcal{T}}u(x) = d \frac{u_{\sigma} - u_K}{d_{K,\sigma}} \mathbf{n}_{K,\sigma} \text{ for a.e. } x \in D_{K,\sigma}. \quad (2.7) $$

• $\hat{\nabla}_{\mathcal{T}}$: consistent discrete gradient

$$ \hat{\nabla}_{\mathcal{T}}:X_{\mathcal{T}} \to (L^2(\Omega))^d \\ u \mapsto \hat{\nabla}_{\mathcal{T}}u, \text{ with } \hat{\nabla}_{\mathcal{T}}u(x) = d \frac{1}{|K|} \sum_{\sigma \in \mathcal{E}_K} |\sigma|(\bar{x}_{\sigma} - x_K) \frac{u_{\sigma} - u_K}{d_{K,\sigma}} \text{ for a.e. } x \in D_{K,\sigma}. \quad (2.8) $$

**Remark 2.3** (On the inflated and consistent gradients). The inflated discrete gradient, first introduced in [17, Definition 2] only involves the normal discrete gradient, but with a factor d, hence the term inflated; it also appears, but somewhat hidden, in the weak formulation (2.6) of the FV scheme in [12]. The consistent discrete gradient $\hat{\nabla}_{\mathcal{T}}u$ was first introduced in [19, Definition 2.3] in one of the first attempts to generalize finite volume schemes to anisotropic diffusion problems.

Following [12], integrating (2.1) on a cell $K$ and (formally) integrating by parts yields the following balance equation on each cell $K$:

$$ - \sum_{\sigma \in \mathcal{E}_K} \int_{\sigma} \nabla u \cdot \mathbf{n}_{K,\sigma} \mathrm{d}s = \int_{K} f \mathrm{d}x + \sum_{\sigma \in \mathcal{E}_K} \int_{\sigma} \mathbf{F} \cdot \mathbf{n}_{K,\sigma} \mathrm{d}s. $$
With a slight abuse of notation, for any $u \in X_T$, we also denote by $u$ the element of $L^2(\Omega)$ which is a.e. equal to the value $u_K$ in any $K \in \mathcal{T}$.

Let us observe that this definition leads to admissible meshes in the sense of [18, Definition 3.1], but is slightly more restrictive since, for the sake of simplicity and as in [12], we do not allow $x_K \in \bar{\sigma}$ if $\sigma \in \mathcal{E}_K$. There are ways to relax this assumption by eliminating the face unknowns in the scheme, but this leads to additional technical difficulties.

As mentioned in the introduction, our analysis requires several discrete derivation operators which we now precisely define.

**Definition 2.2** (Discrete derivative and gradients).

• $G_T$: normal discrete derivative

$$ G_T: X_T \to L^2(\Omega) \\ u \mapsto G_T u, \text{ with } G_T u(x) = \frac{u_\sigma - u_K}{d_{K,\sigma}} \text{ for a.e. } x \in D_{K,\sigma}. \qquad (2.6) $$

• $\nabla_T$: inflated discrete gradient

$$ \nabla_T: X_T \to (L^2(\Omega))^d \\ u \mapsto \nabla_T u, \text{ with } \nabla_T u(x) = d \frac{u_\sigma - u_K}{d_{K,\sigma}} \mathbf{n}_{K,\sigma} \text{ for a.e. } x \in D_{K,\sigma}. \qquad (2.7) $$

• $\hat{\nabla}_T$: consistent discrete gradient

$$ \hat{\nabla}_T: X_T \to (L^2(\Omega))^d \\ u \mapsto \hat{\nabla}_T u, \text{ with } \hat{\nabla}_T u(x) = d \frac{1}{|K|} \sum_{\sigma \in \mathcal{E}_K} |\sigma|(\bar{x}_\sigma - x_K) \frac{u_\sigma - u_K}{d_{K,\sigma}} \text{ for a.e. } x \in D_{K,\sigma}. \qquad (2.8) $$

**Remark 2.3** (On the inflated and consistent gradients). The inflated discrete gradient, first introduced in [17, Definition 2] only involves the normal discrete gradient, but with a factor d, hence the term inflated; it also appears, but somewhat hidden, in the weak formulation (2.6) of the FV scheme in [12]. The consistent discrete gradient $\hat{\nabla}_T u$ was first introduced in [19, Definition 2.3] in one of the first attempts to generalize finite volume schemes to anisotropic diffusion problems.

Following [12], integrating (2.1) on a cell $K$ and (formally) integrating by parts yields the following balance equation on each cell $K$:

$$ - \sum_{\sigma \in \mathcal{E}_K} \int_{\sigma} \nabla u \cdot \mathbf{n}_{K,\sigma} \mathrm{d}s = \int_{K} f \mathrm{d}x + \sum_{\sigma \in \mathcal{E}_K} \int_{\sigma} \mathbf{F} \cdot \mathbf{n}_{K,\sigma} \mathrm{d}s. $$
Introducing the discrete unknowns $(u_K)_{K \in \mathcal{T}}$ and $(u_\sigma)_{\sigma \in \mathcal{E}}$, and following [12], we propose the following scheme:

$$ \forall K \in \mathcal{T}, -\sum_{\sigma \in \mathcal{E}_K} |\sigma| \frac{u_\sigma - u_K}{d_{K,\sigma}} = \int_K f \, dx + \sum_{\sigma \in \mathcal{E}_K} |\sigma| \frac{1}{|D_{K,\sigma}|} \int_{D_{K,\sigma}} \mathbf{F} \cdot \mathbf{n}_{K,\sigma}, \quad (2.9a) $$

$$ \forall \sigma = K |L \in \mathcal{E}_{\text{int}}, $$

$$ \frac{u_\sigma - u_K}{d_{K,\sigma}} + \frac{1}{|D_{K,\sigma}|} \int_{D_{K,\sigma}} \mathbf{F} \cdot \mathbf{n}_{K,\sigma} \, dx = - \frac{u_\sigma - u_L}{d_{L,\sigma}} - \frac{1}{|D_{L,\sigma}|} \int_{D_{L,\sigma}} \mathbf{F} \cdot \mathbf{n}_{L,\sigma} \, dx. \quad (2.9b) $$

$$ \forall \sigma \in \mathcal{E}_{\text{ext}}, u_{\sigma} = 0. \quad (2.9c) $$

Equation (2.9a) is the discretization of the local mass balance on the cell $K$, while (2.9b) expresses the conservativity of the discrete fluxes. By [12, Theorem 2.1], there exists a unique solution to the scheme (2.9). Moreover, the scheme admits a weak formulation [12, Lemma 2.1].

**Lemma 2.4** (Weak formulation of the scheme). The scheme (2.9) is equivalent to the following weak formulation:

$$ \text{Find } u \in X_{\mathcal{T}} \text{ such that, for any } v = ((v_K)_{K \in \mathcal{T}}, (v_\sigma)_{\sigma \in \mathcal{E}}) \in X_{\mathcal{T}}, \\ d \int_{\Omega} G_{\mathcal{T}} u(x) G_{\mathcal{T}} v(x) \, dx = \int_{\Omega} (f(x)v(x) - \mathbf{F}(x) \cdot \nabla_{\mathcal{T}} v(x)) \, dx, \quad (2.10) $$

where $G_{\mathcal{T}}v$ and $\nabla_{\mathcal{T}}v$ are the discrete derivative and gradient given in Definition 2.2.

*Proof.* Owing to the definitions (2.6) and (2.7) of $G_{\mathcal{T}}v$ and $\nabla_{\mathcal{T}}v$ and noting that $|D_{K,\sigma}| = \frac{|\sigma| d_{K,\sigma}}{d}$, the scheme (2.10) also reads

$$ \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K} \frac{|\sigma|}{d_{K,\sigma}} (u_\sigma - u_K)(v_\sigma - v_K) = \sum_{K \in \mathcal{T}} \int_K f v_K \, dx - \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K} \int_{D_{K,\sigma}} \mathbf{F} d \frac{v_\sigma - v_K}{d_{K,\sigma}} \cdot \mathbf{n}_{K,\sigma}. \quad (2.11) $$

The proof that if $u$ is a solution to (2.9) then $u$ satisfies (2.11) may be found in [12, Lemma 2.1]. Conversely, letting $v_K = 1$ and $v_\sigma = 0$ in (2.11) leads to (2.9a) and letting $v_K = 0$ and $v_\sigma = 1$ in (2.11) leads to (2.9b). The condition (2.9c) is ensured by the fact that $u \in X_{\mathcal{T}}$. $\square$

Note that the inflated approximate gradient is expected to weakly converge in $L^2$ toward the gradient of the exact solution but can never converge in $L^2$ toward the gradient of the exact solution except if the exact solution is 0 (see [17, Lemma 2 and Remark 2]). The convergence in $L^2$ of the consistent approximate gradient is proved in [19], in the case where $\mathbf{F} = 0$, and using a modified finite volume scheme in order to handle anisotropic diffusion problems, under the condition that $\theta_{\mathcal{T}}$ is uniformly bounded by below.

### 3. CONVERGENCE ANALYSIS

3.1. **An optimal error estimate.** A norm on $X_{\mathcal{T}}$ is defined by

$$ \|u\|_{\mathcal{T}}^2 = d \|G_{\mathcal{T}}u\|_{L^2(\Omega)}^2 = \frac{1}{d} \|\nabla_{\mathcal{T}} u\|_{L^2(\Omega)^d}^2 = \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K} \frac{|\sigma|}{d_{K,\sigma}} (u_\sigma - u_K)^2. \quad (3.1) $$
Introducing the discrete unknowns $(u_K)_{K \in \mathcal{T}}$ and $(u_\sigma)_{\sigma \in \mathcal{E}}$, and following [12], we propose the following scheme:

$$ \forall K \in \mathcal{T}, - \sum_{\sigma \in \mathcal{E}_K} |\sigma| \frac{u_\sigma - u_K}{d_{K,\sigma}} = \int_K f \, dx + \sum_{\sigma \in \mathcal{E}_K} |\sigma| \frac{1}{|D_{K,\sigma}|} \int_{D_{K,\sigma}} \mathbf{F} \cdot \mathbf{n}_{K,\sigma}, \quad (2.9a) $$

$$ \forall \sigma = K |L \in \mathcal{E}_{\text{int}}, $$

$$ \frac{u_\sigma - u_K}{d_{K,\sigma}} + \frac{1}{|D_{K,\sigma}|} \int_{D_{K,\sigma}} \mathbf{F} \cdot \mathbf{n}_{K,\sigma} \, dx = - \frac{u_\sigma - u_L}{d_{L,\sigma}} - \frac{1}{|D_{L,\sigma}|} \int_{D_{L,\sigma}} \mathbf{F} \cdot \mathbf{n}_{L,\sigma} \, dx. \quad (2.9b) $$

$$ \forall \sigma \in \mathcal{E}_{\text{ext}}, u_{\sigma} = 0. \quad (2.9c) $$

Equation (2.9a) is the discretization of the local mass balance on the cell $K$, while (2.9b) expresses the conservativity of the discrete fluxes. By [12, Theorem 2.1], there exists a unique solution to the scheme (2.9). Moreover, the scheme admits a weak formulation [12, Lemma 2.1].

**Lemma 2.4** (Weak formulation of the scheme). The scheme (2.9) is equivalent to the following weak formulation:

$$ \text{Find } u \in X_{\mathcal{T}} \text{ such that, for any } v = ((v_K)_{K \in \mathcal{T}}, (v_\sigma)_{\sigma \in \mathcal{E}}) \in X_{\mathcal{T}}, \\ d \int_{\Omega} G_{\mathcal{T}} u(x) G_{\mathcal{T}} v(x) \, dx = \int_{\Omega} (f(x)v(x) - \mathbf{F}(x) \cdot \nabla_{\mathcal{T}} v(x)) \, dx, \quad (2.10) $$

where $G_{\mathcal{T}}v$ and $\nabla_{\mathcal{T}}v$ are the discrete derivative and gradient given in Definition 2.2.

*Proof.* Owing to the definitions (2.6) and (2.7) of $G_{\mathcal{T}}v$ and $\nabla_{\mathcal{T}}v$ and noting that $|D_{K,\sigma}| = \frac{|\sigma| d_{K,\sigma}}{d}$, the scheme (2.10) also reads

$$ \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K} \frac{|\sigma|}{d_{K,\sigma}} (u_\sigma - u_K)(v_\sigma - v_K) = \sum_{K \in \mathcal{T}} \int_K f v_K \, dx - \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K} \int_{D_{K,\sigma}} \mathbf{F} d \frac{v_\sigma - v_K}{d_{K,\sigma}} \cdot \mathbf{n}_{K,\sigma}. \quad (2.11) $$

The proof that if $u$ is a solution to (2.9) then $u$ satisfies (2.11) may be found in [12, Lemma 2.1]. Conversely, letting $v_K = 1$ and $v_\sigma = 0$ in (2.11) leads to (2.9a) and letting $v_K = 0$ and $v_\sigma = 1$ in (2.11) leads to (2.9b). The condition (2.9c) is ensured by the fact that $u \in X_{\mathcal{T}}$. $\square$

Note that the inflated approximate gradient is expected to weakly converge in $L^2$ toward the gradient of the exact solution but can never converge in $L^2$ toward the gradient of the exact solution except if the exact solution is 0 (see [17, Lemma 2 and Remark 2]). The convergence in $L^2$ of the consistent approximate gradient is proved in [19], in the case where $\mathbf{F} = 0$, and using a modified finite volume scheme in order to handle anisotropic diffusion problems, under the condition that $\theta_{\mathcal{T}}$ is uniformly bounded by below.

### 3. CONVERGENCE ANALYSIS

**3.1. An optimal error estimate.** A norm on $X_{\mathcal{T}}$ is defined by

$$ \|u\|_{\mathcal{T}}^2 = d \|G_{\mathcal{T}}u\|_{L^2(\Omega)}^2 = \frac{1}{d} \|\nabla_{\mathcal{T}} u\|_{L^2(\Omega)^d}^2 = \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K} \frac{|\sigma|}{d_{K,\sigma}} (u_\sigma - u_K)^2. \quad (3.1) $$
The discrete Poincaré inequality [18, Lemma 9.1] states that for any piecewise function $u$ equal to $u_K$ on the cell $K$,

$$ \|u\|_{L^2(\Omega)} \le \mathrm{diam}(\Omega) \left( \sum_{\sigma=K|L \in \mathcal{E}_{\mathrm{int}}} \frac{|\sigma|}{d_{K,\sigma} + d_{L,\sigma}} (u_K - u_L)^2 + \sum_{\sigma \in \mathcal{E}_{\mathrm{ext}} \cap \mathcal{E}_K} \frac{|\sigma|}{d_{K,\sigma}} u_K^2 \right). \quad (3.2) $$

Now consider $u \in X_\mathcal{T}$ satisfying

$$ \frac{u_\sigma - u_K}{d_{K,\sigma}} + \frac{u_\sigma - u_L}{d_{L,\sigma}} = 0, \text{ for all } \sigma = K|L \in \mathcal{E}_{\mathrm{int}}; \quad (3.3) $$

this yields $\frac{(u_\sigma - u_K)^2}{d_{K,\sigma}} + \frac{(u_\sigma - u_L)^2}{d_{L,\sigma}} = \frac{(u_L - u_K)^2}{d_{K,\sigma} + d_{L,\sigma}}$ and therefore, recalling that $u_\sigma = 0$ for any $\sigma \in \mathcal{E}_{\mathrm{ext}}$,

$$ \sum_{\sigma=K|L \in \mathcal{E}_{\mathrm{int}}} \frac{|\sigma|}{d_{K,\sigma} + d_{L,\sigma}} (u_K - u_L)^2 + \sum_{\sigma \in \mathcal{E}_{\mathrm{ext}} \cap \mathcal{E}_K} \frac{|\sigma|}{d_{K,\sigma}} u_K^2 = \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K} \frac{|\sigma|}{d_{K,\sigma}} (u_\sigma - u_K)^2. $$

Observe then that the minimal value of the function $s \mapsto \frac{(s-u_K)^2}{d_{K,\sigma}} + \frac{(s-u_L)^2}{d_{L,\sigma}}$ is obtained for $s = u_\sigma$ such that (3.3) holds, so that for any $u \in X_\mathcal{T}$, whatever the values $u_\sigma$ for any $\sigma = K|L \in \mathcal{E}_{\mathrm{int}}$,

$$ \|u\|_{\mathcal{T}}^2 \ge \sum_{\sigma=K|L \in \mathcal{E}_{\mathrm{int}}} \frac{|\sigma|}{d_{K,\sigma} + d_{L,\sigma}} (u_K - u_L)^2 + \sum_{\sigma \in \mathcal{E}_{\mathrm{ext}} \cap \mathcal{E}_K} \frac{|\sigma|}{d_{K,\sigma}} u_K^2. $$

Hence

$$ \|u\|_{L^2(\Omega)} \le \mathrm{diam}(\Omega) \|u\|_{\mathcal{T}} = \mathrm{diam}(\Omega) \sqrt{d} \|G_{\mathcal{T}} u\|_{L^2(\Omega)}, \quad \forall u \in X_{\mathcal{T}}. \quad (3.4) $$

Note that (3.3) expresses the conservativity of the discrete fluxes (2.9b) in the case $F = 0$.

Next, we recall that the space $H_{\mathrm{div}}(\Omega) = \{\varphi \in L^2(\Omega)^d : \mathrm{div}\varphi \in L^2(\Omega)\}$ is a Hilbert space when equipped with the norm

$$ \|\varphi\|_{H_{\mathrm{div}}}^2 = \|\varphi\|_{L^2(\Omega)^d}^2 + \|\mathrm{div}\varphi\|_{L^2(\Omega)}^2. $$

We introduce the conformity error: $\zeta_\mathcal{T}: H_{\mathrm{div}}(\Omega) \to \mathbb{R}$, defined by

$$ \zeta_{\mathcal{T}}(\varphi) = \sup \left\{ \int_{\Omega} (\mathrm{div}\varphi(x)v(x) + \varphi(x) \cdot \nabla_{\mathcal{T}} v(x)) dx, v \in X_{\mathcal{T}} \text{ with } \|v\|_{\mathcal{T}} = 1 \right\}. \quad (3.5) $$

Let $G_\mathcal{T}: H_0^1(\Omega) \to L^2(\Omega)$, called in the following the mean normal gradient of a function belonging to $H_0^1(\Omega)$, defined for $\varphi \in H_0^1(\Omega)$ by

$$ G_{\mathcal{T}}\varphi(x) = \frac{1}{|D_{K,\sigma}|} \int_{D_{K,\sigma}} \nabla\varphi(x) \cdot n_{K,\sigma} dx \quad \text{for a.e. } x \in D_{K,\sigma}, \quad (3.6) $$

Observe that Definition (3.6) implies the following equality, for $\bar{u} \in H_0^1(\Omega)$ and $v \in H_\mathcal{T}$,

$$ \int_{\Omega} \nabla \bar{u} \cdot \nabla_{\mathcal{T}} v(x) dx = d \int_{\Omega} G_{\mathcal{T}} \bar{u}(x) G_{\mathcal{T}} v(x) dx. \quad (3.7) $$

Next, we define a distance between any function $v \in X_\mathcal{T}$ and $\varphi \in H_0^1(\Omega)$ in the following way:

$$ \delta_{\mathcal{T}}(\varphi, v) = \frac{1}{\mathrm{diam}(\Omega)} \|\varphi - v\|_{L^2(\Omega)} + \sqrt{d} \|G_{\mathcal{T}}\varphi - G_{\mathcal{T}}v\|_{L^2(\Omega)}. \quad (3.8) $$

We may then define a (generalized) interpolation error by:

$$ \forall \varphi \in H_0^1(\Omega), \quad \mathcal{I}_{\mathcal{T}}(\varphi) = \inf_{v \in X_{\mathcal{T}}} \delta_{\mathcal{T}}(\varphi, v). \quad (3.9) $$
The discrete Poincaré inequality [18, Lemma 9.1] states that for any piecewise function $u$ equal to $u_K$ on the cell $K$,

$$ \|u\|_{L^2(\Omega)} \le \mathrm{diam}(\Omega) \left( \sum_{\sigma=K|L \in \mathcal{E}_{\mathrm{int}}} \frac{|\sigma|}{d_{K,\sigma} + d_{L,\sigma}} (u_K - u_L)^2 + \sum_{\sigma \in \mathcal{E}_{\mathrm{ext}} \cap \mathcal{E}_K} \frac{|\sigma|}{d_{K,\sigma}} u_K^2 \right). \quad (3.2) $$

Now consider $u \in X_\mathcal{T}$ satisfying

$$ \frac{u_\sigma - u_K}{d_{K,\sigma}} + \frac{u_\sigma - u_L}{d_{L,\sigma}} = 0, \quad \text{for all } \sigma = K|L \in \mathcal{E}_{\mathrm{int}}; \quad (3.3) $$

this yields $\frac{(u_\sigma - u_K)^2}{d_{K,\sigma}} + \frac{(u_\sigma - u_L)^2}{d_{L,\sigma}} = \frac{(u_L - u_K)^2}{d_{K,\sigma} + d_{L,\sigma}}$ and therefore, recalling that $u_\sigma = 0$ for any $\sigma \in \mathcal{E}_{\mathrm{ext}}$,

$$ \sum_{\sigma=K|L \in \mathcal{E}_{\mathrm{int}}} \frac{|\sigma|}{d_{K,\sigma} + d_{L,\sigma}} (u_K - u_L)^2 + \sum_{\sigma \in \mathcal{E}_{\mathrm{ext}} \cap \mathcal{E}_K} \frac{|\sigma|}{d_{K,\sigma}} u_K^2 = \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K} \frac{|\sigma|}{d_{K,\sigma}} (u_\sigma - u_K)^2. $$

Observe then that the minimal value of the function $s \mapsto \frac{(s-u_K)^2}{d_{K,\sigma}} + \frac{(s-u_L)^2}{d_{L,\sigma}}$ is obtained for $s = u_\sigma$ such that (3.3) holds, so that for any $u \in X_\mathcal{T}$, whatever the values $u_\sigma$ for any $\sigma = K|L \in \mathcal{E}_{\mathrm{int}}$,

$$ \|u\|_{\mathcal{T}}^2 \ge \sum_{\sigma=K|L \in \mathcal{E}_{\mathrm{int}}} \frac{|\sigma|}{d_{K,\sigma} + d_{L,\sigma}} (u_K - u_L)^2 + \sum_{\sigma \in \mathcal{E}_{\mathrm{ext}} \cap \mathcal{E}_K} \frac{|\sigma|}{d_{K,\sigma}} u_K^2. $$

Hence

$$ \|u\|_{L^2(\Omega)} \le \mathrm{diam}(\Omega) \|u\|_{\mathcal{T}} = \mathrm{diam}(\Omega) \sqrt{d} \|G_{\mathcal{T}} u\|_{L^2(\Omega)}, \quad \forall u \in X_{\mathcal{T}}. \quad (3.4) $$

Note that (3.3) expresses the conservativity of the discrete fluxes (2.9b) in the case $F = 0$.

Next, we recall that the space $H_{\mathrm{div}}(\Omega) = \{\varphi \in L^2(\Omega)^d : \mathrm{div}\varphi \in L^2(\Omega)\}$ is a Hilbert space when equipped with the norm

$$ \|\varphi\|_{H_{\mathrm{div}}}^2 = \|\varphi\|_{L^2(\Omega)^d}^2 + \|\mathrm{div}\varphi\|_{L^2(\Omega)}^2. $$

We introduce the conformity error: $\zeta_\mathcal{T} : H_{\mathrm{div}}(\Omega) \to \mathbb{R}$, defined by

$$ \zeta_{\mathcal{T}}(\varphi) = \sup \left\{ \int_{\Omega} (\mathrm{div}\varphi(x)v(x) + \varphi(x) \cdot \nabla_{\mathcal{T}} v(x)) dx, v \in X_{\mathcal{T}} \text{ with } \|v\|_{\mathcal{T}} = 1 \right\}. \quad (3.5) $$

Let $G_\mathcal{T} : H_0^1(\Omega) \to L^2(\Omega)$, called in the following the mean normal gradient of a function belonging to $H_0^1(\Omega)$, defined for $\varphi \in H_0^1(\Omega)$ by

$$ G_{\mathcal{T}}\varphi(x) = \frac{1}{|D_{K,\sigma}|} \int_{D_{K,\sigma}} \nabla\varphi(x) \cdot n_{K,\sigma} dx \quad \text{for a.e. } x \in D_{K,\sigma}, \quad (3.6) $$

Observe that Definition (3.6) implies the following equality, for $\bar{u} \in H_0^1(\Omega)$ and $v \in H_\mathcal{T}$,

$$ \int_{\Omega} \nabla \bar{u} \cdot \nabla_{\mathcal{T}} v(x) dx = d \int_{\Omega} G_{\mathcal{T}} \bar{u}(x) G_{\mathcal{T}} v(x) dx. \quad (3.7) $$

Next, we define a distance between any function $v \in X_\mathcal{T}$ and $\varphi \in H_0^1(\Omega)$ in the following way:

$$ \delta_{\mathcal{T}}(\varphi, v) = \frac{1}{\mathrm{diam}(\Omega)} \|\varphi - v\|_{L^2(\Omega)} + \sqrt{d} \|G_{\mathcal{T}}\varphi - G_{\mathcal{T}}v\|_{L^2(\Omega)}. \quad (3.8) $$

We may then define a (generalized) interpolation error by:

$$ \forall \varphi \in H_0^1(\Omega), \quad \mathcal{I}_{\mathcal{T}}(\varphi) = \inf_{v \in X_{\mathcal{T}}} \delta_{\mathcal{T}}(\varphi, v). \quad (3.9) $$
The following theorem yields an optimal error bound, in the sense that the distance between the solution $\bar{u}$ to (2.1) and the solution $u$ to (2.10) is bounded by above and below by the sum of the conformity error and the interpolation error.

**Theorem 3.1** (Optimal error bound for the approximation of the elliptic problem (2.3)). Let $\bar{u}$ be the solution to the elliptic problem (2.3). Then the solution $u$ to the numerical scheme (2.10) satisfies

$$ \frac{1}{2} (\zeta_{\mathcal{T}}(\nabla \bar{u} + \mathbf{F}) + \mathcal{I}_{\mathcal{T}}(\bar{u})) \leq \delta_{\mathcal{T}}(\bar{u}, u) \leq 3 (\zeta_{\mathcal{T}}(\nabla \bar{u} + \mathbf{F}) + \mathcal{I}_{\mathcal{T}}(\bar{u})). \quad (3.10) $$

**Remark 3.2** (On the definition of $\delta_{\mathcal{T}}(\varphi, v)$). If we replace the definition (3.8) of $\delta_{\mathcal{T}}$ by

$$ \delta_{\mathcal{T},w}(\varphi, v) = \frac{1}{\mathrm{diam}(\Omega)} \|w - v\|_{L^2(\Omega)} + \sqrt{d} \|\mathcal{G}_{\mathcal{T}}\varphi - G_{\mathcal{T}}v\|_{L^2(\Omega)}, $$

with $w \in L^2(\Omega)$ possibly depending on $\mathcal{T}$, then Theorem 3.1 still holds. Nevertheless, in order for the bound (3.10) to yield an error estimate, the interpolation error $\mathcal{I}_{\mathcal{T}}(\bar{u})$ must tend to 0 with the size of the mesh. This implies that $w$ tends to $\bar{u}$ in $L^2(\Omega)$ with the size of the mesh; this is indeed the case if, for instance, $w$ is the piecewise constant function defined by the mean value over the cells, or, for regular enough functions, by the value of $\bar{u}$ at point $x_K$ (for the cell $K$).

**Remark 3.3** (On the definition of $\mathcal{G}_{\mathcal{T}}$). The results of Theorem 3.1 remain true if we define $\mathcal{G}_{\mathcal{T}}: H_0^1(\Omega) \to L^2(\Omega)$ by

$$ \mathcal{G}_{\mathcal{T}}\bar{u}(x) = \nabla\bar{u}(x) \cdot \mathbf{n}_{K,\sigma} \quad \text{for a.e. } x \in D_{K,\sigma}, \quad (3.11) $$

that is by the value of the normal gradient instead of its mean value in $D_{K,\sigma}$. However, this choice leads to larger values for $\delta_{\mathcal{T}}(\bar{u}, u)$. Moreover, it leads to expressions which can be more difficult to evaluate in the numerical implementation.

*Proof of Theorem 3.1.* From the definition (3.5) of the conformity error $\zeta_{\mathcal{T}}$ and owing to (2.5), we obtain that

$$ \int_{\Omega} ((\nabla \bar{u} + \mathbf{F})(x) \cdot \nabla_{\mathcal{T}} v(x) - f(x)v(x)) dx \leq \zeta_{\mathcal{T}}(\nabla \bar{u} + \mathbf{F}) \|v\|_{\mathcal{T}}, \forall v \in X_{\mathcal{T}}. $$

Since $u$ is the solution to (2.10), we get

$$ \int_{\Omega} (\nabla \bar{u}(x) \cdot \nabla_{\mathcal{T}} v(x) - d G_{\mathcal{T}} u(x) G_{\mathcal{T}} v(x)) dx \leq \zeta_{\mathcal{T}} (\nabla \bar{u} + \mathbf{F}) \|v\|_{\mathcal{T}}, \forall v \in X_{\mathcal{T}}. $$

Owing to (3.7), we have

$$ \int_{\Omega} d(\mathcal{G}_{\mathcal{T}}\bar{u}(x) - G_{\mathcal{T}}u(x))G_{\mathcal{T}}v(x)dx = \int_{\Omega} d(\mathcal{G}_{\mathcal{T}}\bar{u}(x) - G_{\mathcal{T}}u(x))G_{\mathcal{T}}v(x)dx \leq \zeta_{\mathcal{T}}(\nabla \bar{u} + \mathbf{F})\|v\|_{\mathcal{T}}. $$

This inequality, together with the triangle and the Cauchy-Schwarz inequalities yield that for any $w \in X_{\mathcal{T}}$,

$$ \int_{\Omega} d(G_{\mathcal{T}}w(x) - G_{\mathcal{T}}u(x))G_{\mathcal{T}}v(x)dx \\ \leq \zeta_{\mathcal{T}}(\nabla \bar{u} + \mathbf{F})\|v\|_{\mathcal{T}} + d\|\mathcal{G}_{\mathcal{T}}\bar{u} - G_{\mathcal{T}}w\|_{L^2(\Omega)}\|G_{\mathcal{T}}v\|_{L^2(\Omega)}, \forall v \in X_{\mathcal{T}}, $$

Choosing $v = w - u$ and simplifying by $\sqrt{d}\|\mathcal{G}_{\mathcal{T}}v\|_{L^2(\Omega)}$, we get

$$ \sqrt{d} \|G_{\mathcal{T}}u - G_{\mathcal{T}}w\|_{L^2(\Omega)} \leq \zeta_{\mathcal{T}}(\nabla \bar{u} + \mathbf{F}) + \sqrt{d} \|\mathcal{G}_{\mathcal{T}}\bar{u} - G_{\mathcal{T}}w\|_{L^2(\Omega)}, \forall w \in X_{\mathcal{T}}. \quad (3.12) $$
The following theorem yields an optimal error bound, in the sense that the distance between the solution $\bar{u}$ to (2.1) and the solution $u$ to (2.10) is bounded by above and below by the sum of the conformity error and the interpolation error.

**Theorem 3.1** (Optimal error bound for the approximation of the elliptic problem (2.3)). Let $\bar{u}$ be the solution to the elliptic problem (2.3). Then the solution $u$ to the numerical scheme (2.10) satisfies

$$ \frac{1}{2} (\zeta_{\mathcal{T}}(\nabla \bar{u} + \mathbf{F}) + \mathcal{I}_{\mathcal{T}}(\bar{u})) \leq \delta_{\mathcal{T}}(\bar{u}, u) \leq 3 (\zeta_{\mathcal{T}}(\nabla \bar{u} + \mathbf{F}) + \mathcal{I}_{\mathcal{T}}(\bar{u})). \quad (3.10) $$

**Remark 3.2** (On the definition of $\delta_{\mathcal{T}}(\varphi, v)$). If we replace the definition (3.8) of $\delta_{\mathcal{T}}$ by

$$ \delta_{\mathcal{T},w}(\varphi, v) = \frac{1}{\mathrm{diam}(\Omega)} \|w - v\|_{L^2(\Omega)} + \sqrt{d} \|\mathcal{G}_{\mathcal{T}}\varphi - G_{\mathcal{T}}v\|_{L^2(\Omega)}, $$

with $w \in L^2(\Omega)$ possibly depending on $\mathcal{T}$, then Theorem 3.1 still holds. Nevertheless, in order for the bound (3.10) to yield an error estimate, the interpolation error $\mathcal{I}_{\mathcal{T}}(\bar{u})$ must tend to 0 with the size of the mesh. This implies that $w$ tends to $\bar{u}$ in $L^2(\Omega)$ with the size of the mesh; this is indeed the case if, for instance, $w$ is the piecewise constant function defined by the mean value over the cells, or, for regular enough functions, by the value of $\bar{u}$ at point $x_K$ (for the cell $K$).

**Remark 3.3** (On the definition of $\mathcal{G}_{\mathcal{T}}$). The results of Theorem 3.1 remain true if we define $\mathcal{G}_{\mathcal{T}}: H_0^1(\Omega) \to L^2(\Omega)$ by

$$ \mathcal{G}_{\mathcal{T}}\bar{u}(x) = \nabla\bar{u}(x) \cdot \mathbf{n}_{K,\sigma} \quad \text{for a.e. } x \in D_{K,\sigma}, \quad (3.11) $$

that is by the value of the normal gradient instead of its mean value in $D_{K,\sigma}$. However, this choice leads to larger values for $\delta_{\mathcal{T}}(\bar{u}, u)$. Moreover, it leads to expressions which can be more difficult to evaluate in the numerical implementation.

*Proof of Theorem 3.1.* From the definition (3.5) of the conformity error $\zeta_{\mathcal{T}}$ and owing to (2.5), we obtain that

$$ \int_{\Omega} ((\nabla \bar{u} + \mathbf{F})(x) \cdot \nabla_{\mathcal{T}} v(x) - f(x)v(x)) dx \leq \zeta_{\mathcal{T}} (\nabla \bar{u} + \mathbf{F}) \|v\|_{\mathcal{T}}, \forall v \in X_{\mathcal{T}}. $$

Since $u$ is the solution to (2.10), we get

$$ \int_{\Omega} (\nabla \bar{u}(x) \cdot \nabla_{\mathcal{T}} v(x) - d G_{\mathcal{T}} u(x) G_{\mathcal{T}} v(x)) dx \leq \zeta_{\mathcal{T}} (\nabla \bar{u} + \mathbf{F}) \|v\|_{\mathcal{T}}, \forall v \in X_{\mathcal{T}}. $$

Owing to (3.7), we have

$$ \int_{\Omega} d(\mathcal{G}_{\mathcal{T}}\bar{u}(x) - G_{\mathcal{T}}u(x))G_{\mathcal{T}}v(x)dx = \int_{\Omega} d(\mathcal{G}_{\mathcal{T}}\bar{u}(x) - G_{\mathcal{T}}u(x))G_{\mathcal{T}}v(x)dx \leq \zeta_{\mathcal{T}}(\nabla\bar{u} + \mathbf{F})\|v\|_{\mathcal{T}}. $$

This inequality, together with the triangle and the Cauchy-Schwarz inequalities yield that for any $w \in X_{\mathcal{T}}$,

$$ \int_{\Omega} d(G_{\mathcal{T}}w(x) - G_{\mathcal{T}}u(x))G_{\mathcal{T}}v(x)dx \\ \leq \zeta_{\mathcal{T}}(\nabla \bar{u} + \mathbf{F})\|v\|_{\mathcal{T}} + d\|\mathcal{G}_{\mathcal{T}}\bar{u} - G_{\mathcal{T}}w\|_{L^2(\Omega)}\|G_{\mathcal{T}}v\|_{L^2(\Omega)}, \forall v \in X_{\mathcal{T}}, $$

Choosing $v = w - u$ and simplifying by $\sqrt{d}\|\mathcal{G}_{\mathcal{T}}v\|_{L^2(\Omega)}$, we get

$$ \sqrt{d} \|G_{\mathcal{T}}u - G_{\mathcal{T}}w\|_{L^2(\Omega)} \leq \zeta_{\mathcal{T}}(\nabla \bar{u} + \mathbf{F}) + \sqrt{d} \|\mathcal{G}_{\mathcal{T}}\bar{u} - G_{\mathcal{T}}w\|_{L^2(\Omega)}, \forall w \in X_{\mathcal{T}}. \quad (3.12) $$
Owing to the Poincaré inequality (3.4), this latter inequality implies that

$$ \frac{1}{\mathrm{diam}(\Omega)} \|u - w\|_{L^2(\Omega)} \le \zeta_{\mathcal{T}}(\nabla \bar{u} + \mathbf{F}) + \sqrt{d} \|\mathcal{G}_{\mathcal{T}}\bar{u} - G_{\mathcal{T}}w\|_{L^2(\Omega)}, $$

and thanks to the triangle inequality,

$$ \frac{1}{\mathrm{diam}(\Omega)} \|u - \bar{u}\|_{L^2(\Omega)} \le \zeta_{\mathcal{T}}(\nabla \bar{u} + \mathbf{F}) + \sqrt{d} \|\mathcal{G}_{\mathcal{T}}\bar{u} - G_{\mathcal{T}}w\|_{L^2(\Omega)} + \frac{1}{\mathrm{diam}(\Omega)} \|\bar{u} - w\|_{L^2(\Omega)}. $$

Again by the triangle inequality, we get from (3.12) that

$$ \sqrt{d} \|G_{\mathcal{T}}u - G_{\mathcal{T}}\bar{u}\|_{L^2(\Omega)} \le \zeta_{\mathcal{T}}(\nabla \bar{u} + \mathbf{F}) + 2\sqrt{d} \|G_{\mathcal{T}}\bar{u} - G_{\mathcal{T}}w\|_{L^2(\Omega)}. $$

Adding the two previous inequalities and taking the infimum on $w \in X_{\mathcal{T}}$ yields the inequality on the right of (3.10).

Let us now prove the left inequality of (3.10). Let $v \in X_{\mathcal{T}}$; from (3.7) and the FV scheme (2.10), and owing to the definitions (3.8) of and (3.1) of the norm, we get that,

$$ \int_{\Omega} ((\nabla \bar{u} + \mathbf{F})(x) \cdot \nabla_{\mathcal{T}} v(x) - f(x)v(x)) dx = \int_{\Omega} d(\mathcal{G}_{\mathcal{T}}\bar{u}(x) - G_{\mathcal{T}}u(x)) G_{\mathcal{T}}v(x) dx \\ \le \delta_{\mathcal{T}}(\bar{u}, u) \|v\|_{\mathcal{T}}. $$

Passing to the supremum on the functions $v \in H_{\mathcal{T}}$ which are such that $\|v\|_{\mathcal{T}} = 1$ yields

$$ \zeta_{\mathcal{T}}(\nabla \bar{u} + \mathbf{F}) \le \delta_{\mathcal{T}}(\bar{u}, u), $$

and since $\mathcal{I}_{\mathcal{T}}\bar{u} = \inf_{v \in X_{\mathcal{T}}} \delta_{\mathcal{T}}(\bar{u}, v) \le \delta_{\mathcal{T}}(\bar{u}, u)$, we get the left part of (3.10).

□

**Remark 3.4** (Existence and uniqueness). Note that the existence and uniqueness result of $u$, solution to the numerical scheme (2.10) which was proven in [12, Theorem 2.1], may also be seen as a consequence of (3.10). Indeed, the components of $u$ are solution to the square linear system given by the numerical scheme. A null right hand side to this linear system is obtained by setting $\mathbf{F} = 0$ and $f = 0$, which leads to $\bar{u} = 0$. The error estimate result yields a 0 value in the right hand side of (3.10) and therefore $\delta_{\mathcal{T}}(0, u) = 0$, implying $u = 0$. Hence the linear system is invertible so that there exists a unique solution to (2.10).

**Lemma 3.5** (Error bound for the consistent approximate gradient). Let $\bar{u} \in H_0^1(\Omega)$ and let $u \in X_{\mathcal{T}}$, then the following bound holds:

$$ \|\hat{\nabla}_{\mathcal{T}} u - \nabla \bar{u}\|_{L^2(\Omega)^d} \le \frac{d}{\theta_{\mathcal{T}}} (\|G_{\mathcal{T}} u - \mathcal{G}_{\mathcal{T}} \bar{u}\|_{L^2(\Omega)} + \Theta_{\mathcal{T}}(\nabla \bar{u})), $$

with, for any $\varphi \in L^2(\Omega)^d$,

$$ \Theta_{\mathcal{T}}(\varphi)^2 = \sum_{K \in \mathcal{T}} \int_K \int_K \frac{1}{|K|} |\varphi(y) - \varphi(x)|^2 dx dy. \quad (3.13) $$

As a consequence, if $\bar{u}$ is the solution to (2.3) and $u$ is the solution of (2.10), then

$$ \|\hat{\nabla}_{\mathcal{T}} u - \nabla \bar{u}\|_{L^2(\Omega)^d} \le \frac{d}{\theta_{\mathcal{T}}} \left( 2(1 + \mathrm{diam}(\Omega)) \left( \zeta_{\mathcal{T}}(\nabla \bar{u} + \mathbf{F}) + \inf_{v \in X_{\mathcal{T}}} \delta_{\mathcal{T}}(\bar{u}, v) \right) + \Theta_{\mathcal{T}}(\nabla \bar{u}) \right). \quad (3.14) $$
Owing to the Poincaré inequality (3.4), this latter inequality implies that

$$ \frac{1}{\text{diam}(\Omega)} \|u - w\|_{L^2(\Omega)} \le \zeta_{\mathcal{T}}(\nabla \bar{u} + \mathbf{F}) + \sqrt{d} \|\mathcal{G}_{\mathcal{T}}\bar{u} - G_{\mathcal{T}}w\|_{L^2(\Omega)}, $$

and thanks to the triangle inequality,

$$ \frac{1}{\text{diam}(\Omega)} \|u - \bar{u}\|_{L^2(\Omega)} \le \zeta_{\mathcal{T}}(\nabla \bar{u} + \mathbf{F}) + \sqrt{d} \|\mathcal{G}_{\mathcal{T}}\bar{u} - G_{\mathcal{T}}w\|_{L^2(\Omega)} + \frac{1}{\text{diam}(\Omega)} \|\bar{u} - w\|_{L^2(\Omega)}. $$

Again by the triangle inequality, we get from (3.12) that

$$ \sqrt{d} \|G_{\mathcal{T}}u - G_{\mathcal{T}}\bar{u}\|_{L^2(\Omega)} \le \zeta_{\mathcal{T}}(\nabla \bar{u} + \mathbf{F}) + 2\sqrt{d} \|G_{\mathcal{T}}\bar{u} - G_{\mathcal{T}}w\|_{L^2(\Omega)}. $$

Adding the two previous inequalities and taking the infimum on $w \in X_{\mathcal{T}}$ yields the inequality on the right of (3.10).

Let us now prove the left inequality of (3.10). Let $v \in X_{\mathcal{T}}$; from (3.7) and the FV scheme (2.10), and owing to the definitions (3.8) of and (3.1) of the norm, we get that,

$$ \int_{\Omega} ((\nabla \bar{u} + \mathbf{F})(x) \cdot \nabla_{\mathcal{T}} v(x) - f(x)v(x)) dx = \int_{\Omega} d(\mathcal{G}_{\mathcal{T}}\bar{u}(x) - G_{\mathcal{T}}u(x)) G_{\mathcal{T}}v(x) dx \\ \le \delta_{\mathcal{T}}(\bar{u}, u) \|v\|_{\mathcal{T}}. $$

Passing to the supremum on the functions $v \in H_{\mathcal{T}}$ which are such that $\|v\|_{\mathcal{T}} = 1$ yields

$$ \zeta_{\mathcal{T}}(\nabla \bar{u} + \mathbf{F}) \le \delta_{\mathcal{T}}(\bar{u}, u), $$

and since $\mathcal{I}_{\mathcal{T}}\bar{u} = \inf_{v \in X_{\mathcal{T}}} \delta_{\mathcal{T}}(\bar{u}, v) \le \delta_{\mathcal{T}}(\bar{u}, u)$, we get the left part of (3.10).

□

**Remark 3.4** (Existence and uniqueness). Note that the existence and uniqueness result of $u$, solution to the numerical scheme (2.10) which was proven in [12, Theorem 2.1], may also be seen as a consequence of (3.10). Indeed, the components of $u$ are solution to the square linear system given by the numerical scheme. A null right hand side to this linear system is obtained by setting $\mathbf{F} = 0$ and $f = 0$, which leads to $\bar{u} = 0$. The error estimate result yields a 0 value in the right hand side of (3.10) and therefore $\delta_{\mathcal{T}}(0, u) = 0$, implying $u = 0$. Hence the linear system is invertible so that there exists a unique solution to (2.10).

**Lemma 3.5** (Error bound for the consistent approximate gradient). Let $\bar{u} \in H_0^1(\Omega)$ and let $u \in X_{\mathcal{T}}$, then the following bound holds:

$$ \|\hat{\nabla}_{\mathcal{T}} u - \nabla \bar{u}\|_{L^2(\Omega)^d} \le \frac{d}{\theta_{\mathcal{T}}} (\|G_{\mathcal{T}} u - \mathcal{G}_{\mathcal{T}} \bar{u}\|_{L^2(\Omega)} + \Theta_{\mathcal{T}}(\nabla \bar{u})), $$

with, for any $\varphi \in L^2(\Omega)^d$,

$$ \Theta_{\mathcal{T}}(\varphi)^2 = \sum_{K \in \mathcal{T}} \int_K \int_K \frac{1}{|K|} |\varphi(y) - \varphi(x)|^2 dx dy. \quad (3.13) $$

As a consequence, if $\bar{u}$ is the solution to (2.3) and $u$ is the solution of (2.10), then

$$ \|\hat{\nabla}_{\mathcal{T}} u - \nabla \bar{u}\|_{L^2(\Omega)^d} \le \frac{d}{\theta_{\mathcal{T}}} \left( 2(1 + \text{diam}(\Omega)) \left( \zeta_{\mathcal{T}}(\nabla \bar{u} + \mathbf{F}) + \inf_{v \in X_{\mathcal{T}}} \delta_{\mathcal{T}}(\bar{u}, v) \right) + \Theta_{\mathcal{T}}(\nabla \bar{u}) \right). \quad (3.14) $$
Proof. We define $\overline{\nabla}_{K,\sigma}\bar{u} = \frac{1}{|D_{K,\sigma}|} \int_{D_{K,\sigma}} \nabla\bar{u}(x)dx$. Note that $G_T\bar{u}$ is a.e. equal to $n_{K,\sigma} \cdot \overline{\nabla}_{K,\sigma}\bar{u}$ in $D_{K,\sigma}$. Owing to [19, Lemme 2.4],

$$ \frac{1}{|K|} \sum_{\sigma \in \mathcal{E}_K} |\sigma|(\bar{x}_{\sigma} - x_K)\mathbf{n}_{K,\sigma}^t = \mathrm{Id}_d, $$

where $\mathrm{Id}_d$ is the $d \times d$ identity matrix. This is true for any choice of $x_K$, however in the sequel we use it for the choice of $x_K$ satisfying the definition 2.1 of admissible meshes. Let $\overline{\nabla}_T\bar{u}$ be defined a.e. on $\Omega$ by:

$$ \forall x \in K, \overline{\nabla}_T \bar{u}(x) \overline{\nabla}_K \bar{u} = \frac{1}{|K|} \sum_{\sigma \in \mathcal{E}_K} |\sigma|(\bar{x}_{\sigma} - x_K)\mathbf{n}_{K,\sigma} \cdot \overline{\nabla}_{K,\sigma} \bar{u}, $$

By definition of the mesh regularity parameter $\theta_{mesh}$, we have $|\bar{x}_{\sigma} - x_K| \le \frac{d_{K,\sigma}}{\theta_T}$,

$$ \begin{aligned} \|\overline{\nabla}_T \bar{u} - \nabla \bar{u}\|_{L^2(\Omega)^d}^2 &= \sum_{K \in \mathcal{T}} \int_K \left( \frac{1}{|K|} \sum_{\sigma \in \mathcal{E}_K} |\sigma|(\bar{x}_{\sigma} - x_K) \mathbf{n}_{K,\sigma} \cdot (\overline{\nabla}_{K,\sigma} \bar{u} - \nabla \bar{u}(x)) \right)^2 dx \\ &\le \frac{d^2}{\theta_T^2} \sum_{K \in \mathcal{T}} \int_K \left( \frac{1}{|K|} \sum_{\sigma \in \mathcal{E}_K} |D_{K,\sigma}| |\overline{\nabla}_{K,\sigma} \bar{u} - \nabla \bar{u}(x)| \right)^2 dx \\ &\le \frac{d^2}{\theta_T^2} \sum_{K \in \mathcal{T}} \int_K \frac{1}{|K|} \sum_{\sigma \in \mathcal{E}_K} |D_{K,\sigma}| \left| \overline{\nabla}_{K,\sigma} \bar{u} - \nabla \bar{u}(x) \right|^2 dx. \end{aligned} $$

On the one hand, by the Cauchy-Schwarz inequality,

$$ \left| \overline{\nabla}_{K,\sigma} \bar{u} - \nabla \bar{u}(x) \right|^2 = \left| \frac{1}{|D_{K,\sigma}|} \int_{D_{K,\sigma}} (\nabla \bar{u}(y) - \nabla \bar{u}(x)) dy \right|^2 \le \frac{1}{|D_{K,\sigma}|} \int_{D_{K,\sigma}} |\nabla \bar{u}(y) - \nabla \bar{u}(x)|^2 dy, $$

so that

$$ \|\overline{\nabla}_T \bar{u} - \nabla \bar{u}\|_{L^2(\Omega)^d}^2 \le \frac{d^2}{\theta_T^2} \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K} \int_{D_{K,\sigma}} \int_K \frac{1}{|K|} |\nabla \bar{u}(y) - \nabla \bar{u}(x)|^2 dx dy = \frac{d^2}{\theta_T^2} \Theta_T (\nabla \bar{u})^2. $$

On the other hand,

$$ \begin{aligned} \|\hat{\nabla}_T u - \overline{\nabla}_T \bar{u}\|_{L^2(\Omega)^d}^2 &= \sum_{K \in \mathcal{T}} \frac{1}{|K|} \left| \sum_{\sigma \in \mathcal{E}_K} |\sigma|(\bar{x}_{\sigma} - x_K) \left( \frac{u_{\sigma} - u_K}{d_{K,\sigma}} - \mathbf{n}_{K,\sigma} \cdot \overline{\nabla}_{K,\sigma} \bar{u} \right) \right|^2 \\ &\le \sum_{K \in \mathcal{T}} \frac{1}{|K|} \left( \sum_{\sigma \in \mathcal{E}_K} |\sigma||\bar{x}_{\sigma} - x_K| \left( \frac{u_{\sigma} - u_K}{d_{K,\sigma}} - \mathbf{n}_{K,\sigma} \cdot \overline{\nabla}_{K,\sigma} \bar{u} \right) \right)^2 \\ &\le \frac{d^2}{\theta_T^2} \sum_{K \in \mathcal{T}} \frac{1}{|K|} \left( \sum_{\sigma \in \mathcal{E}_K} |D_{K,\sigma}| \left( \frac{u_{\sigma} - u_K}{d_{K,\sigma}} - \mathbf{n}_{K,\sigma} \cdot \overline{\nabla}_{K,\sigma} \bar{u} \right) \right)^2 \end{aligned} $$
Proof. We define $\overline{\nabla}_{K,\sigma}\bar{u} = \frac{1}{|D_{K,\sigma}|} \int_{D_{K,\sigma}} \nabla\bar{u}(x)dx$. Note that $G_T\bar{u}$ is a.e. equal to $n_{K,\sigma} \cdot \overline{\nabla}_{K,\sigma}\bar{u}$ in $D_{K,\sigma}$. Owing to [19, Lemme 2.4],

$$ \frac{1}{|K|} \sum_{\sigma \in \mathcal{E}_K} |\sigma|(\bar{x}_{\sigma} - x_K)\mathbf{n}_{K,\sigma}^t = \mathrm{Id}_d, $$

where $\mathrm{Id}_d$ is the $d \times d$ identity matrix. This is true for any choice of $x_K$, however in the sequel we use it for the choice of $x_K$ satisfying the definition 2.1 of admissible meshes. Let $\overline{\nabla}_T\bar{u}$ be defined a.e. on $\Omega$ by:

$$ \forall x \in K, \quad \overline{\nabla}_T \bar{u}(x) \overline{\nabla}_K \bar{u} = \frac{1}{|K|} \sum_{\sigma \in \mathcal{E}_K} |\sigma|(\bar{x}_{\sigma} - x_K)\mathbf{n}_{K,\sigma} \cdot \overline{\nabla}_{K,\sigma} \bar{u}, $$

By definition of the mesh regularity parameter $\theta_{mesh}$, we have $|\bar{x}_{\sigma} - x_K| \le \frac{d_{K,\sigma}}{\theta_T}$,

$$ \begin{aligned} \|\overline{\nabla}_T \bar{u} - \nabla \bar{u}\|_{L^2(\Omega)^d}^2 &= \sum_{K \in \mathcal{T}} \int_K \left( \frac{1}{|K|} \sum_{\sigma \in \mathcal{E}_K} |\sigma|(\bar{x}_{\sigma} - x_K) \mathbf{n}_{K,\sigma} \cdot (\overline{\nabla}_{K,\sigma} \bar{u} - \nabla \bar{u}(x)) \right)^2 dx \\ &\le \frac{d^2}{\theta_T^2} \sum_{K \in \mathcal{T}} \int_K \left( \frac{1}{|K|} \sum_{\sigma \in \mathcal{E}_K} |D_{K,\sigma}| |\overline{\nabla}_{K,\sigma} \bar{u} - \nabla \bar{u}(x)| \right)^2 dx \\ &\le \frac{d^2}{\theta_T^2} \sum_{K \in \mathcal{T}} \int_K \frac{1}{|K|} \sum_{\sigma \in \mathcal{E}_K} |D_{K,\sigma}| \left| \overline{\nabla}_{K,\sigma} \bar{u} - \nabla \bar{u}(x) \right|^2 dx. \end{aligned} $$

On the one hand, by the Cauchy-Schwarz inequality,

$$ \left| \overline{\nabla}_{K,\sigma} \bar{u} - \nabla \bar{u}(x) \right|^2 = \left| \frac{1}{|D_{K,\sigma}|} \int_{D_{K,\sigma}} (\nabla \bar{u}(y) - \nabla \bar{u}(x)) dy \right|^2 \le \frac{1}{|D_{K,\sigma}|} \int_{D_{K,\sigma}} |\nabla \bar{u}(y) - \nabla \bar{u}(x)|^2 dy, $$

so that

$$ \|\overline{\nabla}_T \bar{u} - \nabla \bar{u}\|_{L^2(\Omega)^d}^2 \le \frac{d^2}{\theta_T^2} \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K} \int_{D_{K,\sigma}} \int_K \frac{1}{|K|} |\nabla \bar{u}(y) - \nabla \bar{u}(x)|^2 dx dy = \frac{d^2}{\theta_T^2} \Theta_T (\nabla \bar{u})^2. $$

On the other hand,

$$ \begin{aligned} \|\hat{\nabla}_T u - \overline{\nabla}_T \bar{u}\|_{L^2(\Omega)^d}^2 &= \sum_{K \in \mathcal{T}} \frac{1}{|K|} \left| \sum_{\sigma \in \mathcal{E}_K} |\sigma|(\bar{x}_{\sigma} - x_K) \left( \frac{u_{\sigma} - u_K}{d_{K,\sigma}} - \mathbf{n}_{K,\sigma} \cdot \overline{\nabla}_{K,\sigma} \bar{u} \right) \right|^2 \\ &\le \sum_{K \in \mathcal{T}} \frac{1}{|K|} \left( \sum_{\sigma \in \mathcal{E}_K} |\sigma||\bar{x}_{\sigma} - x_K| \left( \frac{u_{\sigma} - u_K}{d_{K,\sigma}} - \mathbf{n}_{K,\sigma} \cdot \overline{\nabla}_{K,\sigma} \bar{u} \right) \right)^2 \\ &\le \frac{d^2}{\theta_T^2} \sum_{K \in \mathcal{T}} \frac{1}{|K|} \left( \sum_{\sigma \in \mathcal{E}_K} |D_{K,\sigma}| \left( \frac{u_{\sigma} - u_K}{d_{K,\sigma}} - \mathbf{n}_{K,\sigma} \cdot \overline{\nabla}_{K,\sigma} \bar{u} \right) \right)^2 \end{aligned} $$
Applying the Cauchy-Schwarz inequality, we obtain

$$
\begin{align*}
\| \hat{\nabla}_{\mathcal{T}} u - \bar{\nabla}_{\mathcal{T}} \bar{u} \|_{L^2(\Omega)^d}^2 &\le \frac{d^2}{\theta_{\mathcal{T}}^2} \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K} |D_{K,\sigma}| (\frac{u_\sigma - u_K}{d_{K,\sigma}} - \mathbf{n}_{K,\sigma} \cdot \bar{\nabla}_{K,\sigma} \bar{u})^2 \\
&\le \frac{d^2}{\theta_{\mathcal{T}}^2} \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K} |D_{K,\sigma}| (\frac{u_\sigma - u_K}{d_{K,\sigma}} - \mathbf{n}_{K,\sigma} \cdot \bar{\nabla}_{K,\sigma} \bar{u})^2 \\
&= \left( \frac{d}{\theta_{\mathcal{T}}} \| G_{\mathcal{T}} u - G_{\mathcal{T}} \bar{u} \|_{L^2(\Omega)} \right)^2.
\end{align*}
$$

□

3.2. **Convergence of the interpolation and conformity errors.** Owing to Theorem 3.1, the convergence of the scheme (3.16) relies on the convergence of the conformity error and the interpolation error.

**Lemma 3.6** (Convergence of the conformity error). Let $(\mathcal{T}_n)_{n \in \mathbb{N}}$ be a sequence of admissible meshes such that $h_{\mathcal{T}_n} = \max_{K \in \mathcal{T}_n} \text{diam}(K)$ tends to 0 as $n \to \infty$. For any $\varphi \in \mathbf{H}_{\text{div}}(\Omega)$,

$$
\lim_{n \to \infty} \zeta_{\mathcal{T}_n}(\varphi) = 0.
$$

*Proof.* We first consider the case where $\varphi \in C^1(\bar{\Omega})^d$. Let us momentarily drop the index $n$ for simplicity, and let $v \in X_{\mathcal{T}}$ with $\|v\|_{\mathcal{T}} = 1$. Owing to the definition (1.5) of $\zeta_{\mathcal{T}}$,

$$
\zeta_{\mathcal{T}}(\varphi) = \sup_{\substack{v \in X_{\mathcal{T}} \\ \|v\|_{\mathcal{T}} = 1}} z_{\mathcal{T}}(\varphi, v) \text{ where } z_{\mathcal{T}}(\varphi, v) = \int_{\Omega} (\text{div}\varphi(x)v(x) + \varphi(x) \cdot \nabla_{\mathcal{T}} v(x))dx.
$$

Let $\varphi_{K,\sigma} = \frac{1}{|D_{K,\sigma}|} \int_{D_{K,\sigma}} \varphi(x)dx$ and $\varphi_{\sigma} = \frac{1}{|\sigma|} \int_{\sigma} \varphi(x)ds$; by the definition 2.7 of $\nabla_{\mathcal{T}}$, we have

$$
z_{\mathcal{T}}(\varphi, v) = \sum_{K \in \mathcal{T}} v_K \sum_{\sigma \in \mathcal{E}_K} |\sigma| \varphi_{\sigma} \cdot \mathbf{n}_{K,\sigma} + \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K} |D_{K,\sigma}| (v_{\sigma} - v_K) \varphi_{K,\sigma} \cdot \mathbf{n}_{K,\sigma}.
$$

Observing that $|D_{K,\sigma}| = \frac{|\sigma|d_{K,\sigma}}{d}$ and that $\sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K} |\sigma|v_\sigma \varphi_{K,\sigma} \cdot \mathbf{n}_{K,\sigma} = 0$, we get that

$$
\begin{align*}
z_{\mathcal{T}}(\varphi, v) &= \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K} |\sigma|(v_\sigma - v_K)(\varphi_{K,\sigma} \cdot \mathbf{n}_{K,\sigma} - \varphi_\sigma \cdot \mathbf{n}_{K,\sigma}) \\
&\le \|v\|_{\mathcal{T}} \left( \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K} |\sigma| d_{K,\sigma} |\varphi_{K,\sigma} - \varphi_\sigma|^2 \right)^{1/2},
\end{align*}
$$

owing to the Cauchy-Schwarz inequality. Hence

$$
\zeta_{\mathcal{T}}(\varphi) \leq C_{\varphi} h_{\mathcal{T}} (d |\Omega|)^{1/2}.
$$

We now consider the general case $\varphi \in \mathbf{H}_{\text{div}}(\Omega)$. Let $\psi \in C^1(\bar{\Omega})^d$ be given. Again dropping the subscript $n$, we have for $v \in X_{\mathcal{T}}$ with $\|v\|_{\mathcal{T}} = 1$

$$
\begin{align*}
\int_{\Omega} (\operatorname{div}\varphi(x)v(x) + \varphi(x) \cdot \nabla_{\mathcal{T}} v(x))dx &= \int_{\Omega} (\operatorname{div}\psi(x)v(x) + \psi(x) \cdot \nabla_{\mathcal{T}} v(x))dx + \delta \\
&\le C_{\psi} h_{\mathcal{T}} (d |\Omega|)^{1/2} + \delta,
\end{align*}
$$
Applying the Cauchy-Schwarz inequality, we obtain

$$ \begin{aligned} \|\hat{\nabla}_{\mathcal{T}} u - \bar{\nabla}_{\mathcal{T}} \bar{u}\|_{L^2(\Omega)^d}^2 &\le \frac{d^2}{\theta_{\mathcal{T}}^2} \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K} |D_{K,\sigma}| (\frac{u_\sigma - u_K}{d_{K,\sigma}} - \mathbf{n}_{K,\sigma} \cdot \bar{\nabla}_{K,\sigma} \bar{u})^2 \\ &\le \frac{d^2}{\theta_{\mathcal{T}}^2} \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K} |D_{K,\sigma}| (\frac{u_\sigma - u_K}{d_{K,\sigma}} - \mathbf{n}_{K,\sigma} \cdot \bar{\nabla}_{K,\sigma} \bar{u})^2 \\ &= \left( \frac{d}{\theta_{\mathcal{T}}} \|G_{\mathcal{T}} u - G_{\mathcal{T}} \bar{u}\|_{L^2(\Omega)} \right)^2. \end{aligned} $$

□

## 3.2. Convergence of the interpolation and conformity errors.
Owing to Theorem 3.1, the convergence of the scheme (3.16) relies on the convergence of the conformity error and the interpolation error.

**Lemma 3.6** (Convergence of the conformity error). Let $(\mathcal{T}_n)_{n \in \mathbb{N}}$ be a sequence of admissible meshes such that $h_{\mathcal{T}_n} = \max_{K \in \mathcal{T}_n} \text{diam}(K)$ tends to 0 as $n \to \infty$. For any $\varphi \in \mathbf{H}_{\text{div}}(\Omega)$,

$$ \lim_{n \to \infty} \zeta_{\mathcal{T}_n}(\varphi) = 0. $$

*Proof.* We first consider the case where $\varphi \in C^1(\bar{\Omega})^d$. Let us momentarily drop the index $n$ for simplicity, and let $v \in X_{\mathcal{T}}$ with $\|v\|_{\mathcal{T}} = 1$. Owing to the definition (1.5) of $\zeta_{\mathcal{T}}$,

$$ \zeta_{\mathcal{T}}(\varphi) = \sup_{\substack{v \in X_{\mathcal{T}} \\ \|v\|_{\mathcal{T}} = 1}} z_{\mathcal{T}}(\varphi, v) \text{ where } z_{\mathcal{T}}(\varphi, v) = \int_{\Omega} (\text{div}\varphi(x)v(x) + \varphi(x) \cdot \nabla_{\mathcal{T}} v(x))dx. $$

Let $\varphi_{K,\sigma} = \frac{1}{|D_{K,\sigma}|} \int_{D_{K,\sigma}} \varphi(x)dx$ and $\varphi_{\sigma} = \frac{1}{|\sigma|} \int_{\sigma} \varphi(x)ds$; by the definition 2.7 of $\nabla_{\mathcal{T}}$, we have

$$ z_{\mathcal{T}}(\varphi, v) = \sum_{K \in \mathcal{T}} v_K \sum_{\sigma \in \mathcal{E}_K} |\sigma| \varphi_{\sigma} \cdot \mathbf{n}_{K,\sigma} + \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K} |D_{K,\sigma}| (v_{\sigma} - v_K) \varphi_{K,\sigma} \cdot \mathbf{n}_{K,\sigma}. $$

Observing that $|D_{K,\sigma}| = \frac{|\sigma|d_{K,\sigma}}{d}$ and that $\sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K} |\sigma|v_{\sigma}\varphi_{K,\sigma} \cdot \mathbf{n}_{K,\sigma} = 0$, we get that

$$ \begin{aligned} z_{\mathcal{T}}(\varphi, v) &= \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K} |\sigma|(v_{\sigma} - v_K)(\varphi_{K,\sigma} \cdot \mathbf{n}_{K,\sigma} - \varphi_{\sigma} \cdot \mathbf{n}_{K,\sigma}) \\ &\le \|v\|_{\mathcal{T}} \left( \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K} |\sigma|d_{K,\sigma}|\varphi_{K,\sigma} - \varphi_{\sigma}|^2 \right)^{1/2}, \end{aligned} $$

owing to the Cauchy-Schwarz inequality. Hence

$$ \zeta_{\mathcal{T}}(\varphi) \le C_{\varphi} h_{\mathcal{T}} (d |\Omega|)^{1/2}. $$

We now consider the general case $\varphi \in \mathbf{H}_{\text{div}}(\Omega)$. Let $\psi \in C^1(\bar{\Omega})^d$ be given. Again dropping the subscript $n$, we have for $v \in X_{\mathcal{T}}$ with $\|v\|_{\mathcal{T}} = 1$

$$ \begin{aligned} \int_{\Omega} (\text{div}\varphi(x)v(x) + \varphi(x) \cdot \nabla_{\mathcal{T}} v(x))dx &= \int_{\Omega} (\text{div}\psi(x)v(x) + \psi(x) \cdot \nabla_{\mathcal{T}} v(x))dx + \delta \\ &\le C_{\psi} h_{\mathcal{T}} (d |\Omega|)^{1/2} + \delta, \end{aligned} $$
where thanks to (3.4) and to the Cauchy-Schwarz inequality,

$$ \begin{aligned} |\delta| &\le \|\varphi - \psi\|_{L^2(\Omega)^d} \|\nabla \varphi\|_{L^2(\Omega)^d} + \|\operatorname{div}(\varphi - \psi)\|_{L^2(\Omega)} \|v\|_{L^2(\Omega)} \\ &\le (d + \operatorname{diam}(\Omega)^2)^{1/2} \|\varphi - \psi\|_{H_{\operatorname{div}}(\Omega)}. \end{aligned} $$

Hence

$$ \zeta_{\mathcal{T}}(\varphi) \le C_{\psi} h_{\mathcal{T}}(d |\Omega|)^{1/2} + (d + \operatorname{diam}(\Omega)^2)^{1/2} \|\varphi - \psi\|_{H_{\operatorname{div}}(\Omega)}. $$

We conclude thanks to [27, Theorem 1.1] which states the density of $C^1(\bar{\Omega})^d$ in $H_{\operatorname{div}}(\Omega)$ for any Lipschitz open set $\Omega$ of $\mathbb{R}^d$. $\square$

**Lemma 3.7** (Convergence of the interpolation error). Let $(T_n)_{n \in \mathbb{N}}$ be a sequence of admissible meshes such that $h_{T_n} = \max_{K \in T_n} \operatorname{diam}(K)$ tends to 0 as $n \to \infty$. Then, for any $\varphi \in H_0^1(\Omega)$,

$$ \lim_{n \to \infty} \mathcal{I}_{T_n}(\varphi). $$

*Proof.* We first consider $\varphi \in C_c^\infty(\Omega)$. Let $v^{(n)} \in H_{T_n}$ be defined by $v_K^{(n)} = \varphi(x_K)$ and $v_\sigma^{(n)}$ be such that

$$ |\sigma| \frac{v_{\sigma}^{(n)} - v_{K}^{(n)}}{d_{K,\sigma}} + |\sigma| \frac{v_{\sigma}^{(n)} - v_{L}^{(n)}}{d_{L,\sigma}} = 0, \quad \text{for all } \sigma \in \mathcal{E}_K \cap \mathcal{E}_L, $$

By the definition (3.8),

$$ \delta_{T_n}(\varphi, v^{(n)}) = \frac{1}{\operatorname{diam}(\Omega)} \|\varphi - v^{(n)}\|_{L^2(\Omega)} + \sqrt{d} \|\mathcal{G}_{T_n}\varphi - G_{T_n}v^{(n)}\|_{L^2(\Omega)}. $$

It is clear that $\|\varphi - v^{(n)}\|_{L^2(\Omega)} \to 0$ as $n \to +\infty$. Now observe that

$$ \frac{v_{\sigma}^{(n)} - v_{K}^{(n)}}{d_{K,\sigma}} = \frac{v_{L}^{(n)} - v_{K}^{(n)}}{d_{K,\sigma} + d_{L,\sigma}}; $$

therefore, since $\varphi \in C_c^\infty(\Omega)$, we have for any $x \in D_{K,\sigma}$,

$$ |\mathcal{G}_{T_n}\varphi(x) - G_{T_n}v^{(n)}(x)| \le C_{\varphi}h_{T_n}, $$

and the result follows.

We now consider the case $\varphi \in H_0^1(\Omega)$; let $\psi \in C_c^\infty(\Omega)$ be given. By the triangle inequality,

$$ \delta_{T_n}(\varphi, v) \le \delta_{T_n}(\psi, v) + \frac{1}{\operatorname{diam}(\Omega)} \|\varphi - \psi\|_{L^2} + \sqrt{d} \|\mathcal{G}_{T_n}\varphi - \mathcal{G}_{T_n}\psi\|_{L^2}. $$

Observing that for a.e. $x \in \Omega$,

$$ |\mathcal{G}_{T_n}(\varphi - \psi)(x)| \le |\nabla(\varphi - \psi)(x)|, $$

we get

$$ \delta_{T_n}(\varphi, v) \le \delta_{T_n}(\psi, v) + \max(\frac{1}{\operatorname{diam}(\Omega)}, \sqrt{d}) \|\varphi - \psi\|_{H^1}, $$

and the result follows by density. $\square$

**Lemma 3.8.** Let $\varphi \in L^2(\Omega)^d$, and let $\Theta_{T_n}(\varphi)$ be defined by (3.13) for any $n \in \mathbb{N}$. Then

$$ \lim_{n \to \infty} \Theta_{T_n}(\varphi) = 0. \qquad (3.15) $$
where thanks to (3.4) and to the Cauchy-Schwarz inequality,

$$
\begin{align*}
|\delta| &\le \|\varphi - \psi\|_{L^2(\Omega)^d} \|\nabla \tau v\|_{L^2(\Omega)^d} + \|\mathrm{div}(\varphi - \psi)\|_{L^2(\Omega)} \|v\|_{L^2(\Omega)} \\
&\le (d + \mathrm{diam}(\Omega)^2)^{1/2} \|\varphi - \psi\|_{H_{\mathrm{div}}(\Omega)}.
\end{align*}
$$

Hence

$$
\zeta_{\mathcal{T}}(\varphi) \le C_{\psi} h_{\mathcal{T}}(d |\Omega|)^{1/2} + (d + \mathrm{diam}(\Omega)^2)^{1/2} \|\varphi - \psi\|_{\mathcal{H}_{\mathrm{div}}(\Omega)}.
$$

We conclude thanks to [27, Theorem 1.1] which states the density of $C^1(\bar{\Omega})^d$ in $\mathbf{H}_{\text{div}}(\Omega)$ for any Lipschitz open set $\Omega$ of $\mathbb{R}^d$. $\square$

**Lemma 3.7** (Convergence of the interpolation error). Let ($\mathcal{T}_n$)$_{n \in \mathbb{N}}$ be a sequence of admissible meshes such that $h_{\mathcal{T}_n} = \max_{K \in \mathcal{T}_n} \text{diam}(K)$ tends to 0 as $n \to \infty$. Then, for any $\varphi \in H_0^1(\Omega)$,

$$
\lim_{n \to \infty} \mathcal{I}_{\mathcal{T}_n}(\varphi).
$$

Proof. We first consider $\varphi \in C_c^\infty(\Omega)$. Let $v^{(n)} \in H_{\mathcal{T}_n}$ be defined by $v_K^{(n)} = \varphi(x_K)$ and $v_\sigma^{(n)}$ be such that

$$
|\sigma| \frac{v_{\sigma}^{(n)} - v_{K}^{(n)}}{d_{K,\sigma}} + |\sigma| \frac{v_{\sigma}^{(n)} - v_{L}^{(n)}}{d_{L,\sigma}} = 0, \text{ for all } \sigma \in \mathcal{E}_K \cap \mathcal{E}_L,
$$

By the definition (3.8),

$$
\delta_{\mathcal{T}_n}(\varphi, v^{(n)}) = \frac{1}{\text{diam}(\Omega)} \|\varphi - v^{(n)}\|_{L^2(\Omega)} + \sqrt{d} \|\mathcal{G}_{\mathcal{T}_n} \varphi - G_{\mathcal{T}_n} v^{(n)}\|_{L^2(\Omega)}.
$$

It is clear that $\|\varphi - v^{(n)}\|_{L^2(\Omega)} \to 0$ as $n \to +\infty$. Now observe that

$$
\frac{v_{\sigma}^{(n)} - v_{K}^{(n)}}{d_{K,\sigma}} = \frac{v_{L}^{(n)} - v_{K}^{(n)}}{d_{K,\sigma} + d_{L,\sigma}};
$$

therefore, since $\varphi \in C_c^\infty(\Omega)$, we have for any $x \in D_{K,\sigma}$,

$$
|\mathcal{G}_{\mathcal{T}_n} \varphi(x) - G_{\mathcal{T}_n} v^{(n)}(x)| \leq C_{\varphi} h_{\mathcal{T}_n},
$$

and the result follows.

We now consider the case $\varphi \in H_0^1(\Omega)$; let $\psi \in C_c^\infty(\Omega)$ be given. By the triangle inequality,

$$
\delta_{\mathcal{T}_n}(\varphi, v) \le \delta_{\mathcal{T}_n}(\psi, v) + \frac{1}{\mathrm{diam}(\Omega)} \|\varphi - \psi\|_{L^2} + \sqrt{d} \|\mathcal{G}_{\mathcal{T}_n}\varphi - \mathcal{G}_{\mathcal{T}_n}\psi\|_{L^2}.
$$

Observing that for a.e. $x \in \Omega$,

$$
|\mathcal{G}_{\mathcal{T}_n}(\varphi - \psi)(x)| \leq |\nabla(\varphi - \psi)(x)|,
$$

we get

$$
\delta_{\mathcal{T}_n}(\varphi, v) \le \delta_{\mathcal{T}_n}(\psi, v) + \max\left(\frac{1}{\mathrm{diam}(\Omega)}, \sqrt{d}\right) \|\varphi - \psi\|_{H^1},
$$

and the result follows by density.

**Lemma 3.8.** Let $\varphi \in L^2(\Omega)^d$, and let $\Theta_{\mathcal{T}_n}(\varphi)$ be defined by (3.13) for any $n \in \mathbb{N}$. Then

$$
\lim_{n \to \infty} \Theta_{\mathcal{T}_n}(\varphi) = 0. \qquad (3.15)
$$
Proof. We observe that, for a given admissible mesh $\mathcal{T}$ and for $\psi \in C^1(\bar{\Omega})^d$, using

$$|\psi(y) - \psi(x)| \le h_{\mathcal{T}} C_{\psi},$$

we get

$$|\Theta_{\mathcal{T}}(\psi)| \le h_{\mathcal{T}} C_{\psi} \sqrt{|\bar{\Omega}|}.$$

By the triangle inequality

$$|\varphi(y) - \varphi(x)| \le |\varphi(y) - \psi(y)| + |\psi(y) - \psi(x)| + |\psi(x) - \varphi(x)|,$$

we obtain

$$\Theta_{\mathcal{T}}(\varphi)^2 \le 3(2\|\varphi - \psi\|_{L^2(\Omega)^d}^2 + \Theta_{\mathcal{T}}(\psi)^2)$$

which concludes the proof by density of $C^1(\bar{\Omega})^d$ in $L^2(\Omega)^d$. $\square$

We can then conclude the convergence of the finite volume scheme.

**Theorem 3.9** (Convergence of the approximate solution and gradient to the solution of (2.3)).
Let $\bar{u} \in H_0^1(\Omega)$ be the solution of (2.3) and let $u_n \in H_{T_n}$ be the solution of (2.10) for $\mathcal{T} = \mathcal{T}_n$.
Then the following holds:

$$\lim_{n \to \infty} (\|\bar{u} - u_n\|_{L^2(\Omega)}^2 + d\|G_{T_n}\bar{u} - G_{T_n}u_n\|_{L^2(\Omega)}^2) = 0. \quad (3.16)$$

Moreover, if the sequence $(\theta_{T_n})_{n \in \mathbb{N}}$ is bounded by below by $\theta_0 > 0$, then

$$\lim_{n \to \infty} \| \hat{\nabla}_{T_n} u_n - \nabla \bar{u} \|_{L^2(\Omega)} = 0. \quad (3.17)$$

Proof. Applying Lemmas 3.6 and 3.7, (3.16) is a consequence of Theorem 3.1. We conclude (3.17) using Lemmas 3.8 and 3.5. $\square$

**3.3. Error estimate in the $H^2(\Omega)$ case.** We now suppose that the exact solution belongs to $H^2(\Omega)$.

**Lemma 3.10** (Estimate of difference with average value in the $H^1$ case). Let $V \subset \mathbb{R}^d$ be a bounded convex open set and let $h_V$ be the diameter of $V$. Let $\varphi \in H^1(V)$. Then

$$\int_V \left( \frac{1}{|V|} \int_V \varphi(y) dy - \varphi(x) \right)^2 dx \le h_V^2 \frac{C_d h_V^d}{|V|} \int_V |\nabla \varphi(z)|^2 dz, \quad (3.18)$$

where $C_d > 0$ is the $d$-dimensional measure of the unit ball of $\mathbb{R}^d$.

Proof. The proof of (3.18) is inspired by that of [18, Lemma 10.2], which is expressed in a discrete setting and is in fact trickier than the present proof. By density of $C_c^1(\mathbb{R}^d)$ in $H^1(V)$, we only have to prove (3.18) for $\varphi \in C_c^1(\mathbb{R}^d)$. For $x, y \in V$,

$$[\varphi(y) - \varphi(x)]^2 = \left[ \int_0^1 \nabla\varphi(tx + (1-t)y) \cdot (x-y)dt \right]^2 \le h_V^2 \int_0^1 |\nabla\varphi(tx + (1-t)y)|^2 dt,$$

so that owing to the Fubini-Tonelli theorem,

$$\int_V \int_V [\varphi(y) - \varphi(x)]^2 dx dy \le h_V^2 \int_V \int_V \int_0^1 |\nabla \varphi(tx + (1-t)y)|^2 dt dx dy \\ \le h_V^2 \int_V \int_0^1 \int_V |\nabla \varphi(tx + (1-t)y)|^2 dx dt dy.$$
**Proof.** We observe that, for a given admissible mesh $\mathcal{T}$ and for $\psi \in C^1(\bar{\Omega})^d$, using

$$|\psi(y) - \psi(x)| \le h_{\mathcal{T}} C_{\psi},$$

we get

$$|\Theta_{\mathcal{T}}(\psi)| \le h_{\mathcal{T}} C_{\psi} \sqrt{|\bar{\Omega}|}.$$

By the triangle inequality

$$|\varphi(y) - \varphi(x)| \le |\varphi(y) - \psi(y)| + |\psi(y) - \psi(x)| + |\psi(x) - \varphi(x)|,$$

we obtain

$$\Theta_{\mathcal{T}}(\varphi)^2 \le 3(2\|\varphi - \psi\|_{L^2(\Omega)^d}^2 + \Theta_{\mathcal{T}}(\psi)^2)$$

which concludes the proof by density of $C^1(\bar{\Omega})^d$ in $L^2(\Omega)^d$. $\square$

We can then conclude the convergence of the finite volume scheme.

**Theorem 3.9** (Convergence of the approximate solution and gradient to the solution of (2.3)).
Let $\bar{u} \in H_0^1(\Omega)$ be the solution of (2.3) and let $u_n \in H_{T_n}$ be the solution of (2.10) for $\mathcal{T} = \mathcal{T}_n$.
Then the following holds:

$$\lim_{n \to \infty} (\|\bar{u} - u_n\|_{L^2(\Omega)}^2 + d\|G_{T_n}\bar{u} - G_{T_n}u_n\|_{L^2(\Omega)}^2) = 0. \qquad (3.16)$$

Moreover, if the sequence $(\theta_{T_n})_{n \in \mathbb{N}}$ is bounded by below by $\theta_0 > 0$, then

$$\lim_{n \to \infty} \| \hat{\nabla}_{T_n} u_n - \nabla \bar{u} \|_{L^2(\Omega)} = 0. \qquad (3.17)$$

**Proof.** Applying Lemmas 3.6 and 3.7, (3.16) is a consequence of Theorem 3.1. We conclude (3.17) using Lemmas 3.8 and 3.5. $\square$

**3.3. Error estimate in the $H^2(\Omega)$ case.** We now suppose that the exact solution belongs to $H^2(\Omega)$.

**Lemma 3.10** (Estimate of difference with average value in the $H^1$ case). Let $V \subset \mathbb{R}^d$ be a bounded convex open set and let $h_V$ be the diameter of $V$. Let $\varphi \in H^1(V)$. Then

$$\int_V \left( \frac{1}{|V|} \int_V \varphi(y) dy - \varphi(x) \right)^2 dx \le h_V^2 \frac{C_d h_V^d}{|V|} \int_V |\nabla \varphi(z)|^2 dz, \qquad (3.18)$$

where $C_d > 0$ is the $d$-dimensional measure of the unit ball of $\mathbb{R}^d$.

**Proof.** The proof of (3.18) is inspired by that of [18, Lemma 10.2], which is expressed in a discrete setting and is in fact trickier than the present proof. By density of $C_c^1(\mathbb{R}^d)$ in $H^1(V)$, we only have to prove (3.18) for $\varphi \in C_c^1(\mathbb{R}^d)$. For $x, y \in V$,

$$[\varphi(y) - \varphi(x)]^2 = \left[ \int_0^1 \nabla\varphi(tx + (1-t)y) \cdot (x-y)dt \right]^2 \le h_V^2 \int_0^1 |\nabla\varphi(tx + (1-t)y)|^2 dt,$$

so that owing to the Fubini-Tonelli theorem,

$$\int_V \int_V [\varphi(y) - \varphi(x)]^2 dx dy \le h_V^2 \int_V \int_V \int_0^1 |\nabla \varphi(tx + (1-t)y)|^2 dt dx dy \\ \le h_V^2 \int_V \int_0^1 \int_V |\nabla \varphi(tx + (1-t)y)|^2 dx dt dy.$$
Applying the change of variable (with $t$ and $y$ fixed) $z = tx + (1-t)y$, noticing that $z \in V \cap B(y, th_V)$, and again applying the Fubini-Tonelli theorem, we thus get that

$$ 
\begin{aligned}
\int_V \int_V (\varphi(y) - \varphi(x))^2 dx dy &\le h_V^2 \int_V \left[ \int_0^1 \left( \int_V 1_{B(z,th_V)}(y) dy \right) t^{-d} dt \right] |\nabla \varphi(z)|^2 dz \\
&\le h_V^2 C_d h_V^d \int_V |\nabla \varphi(z)|^2 dz.
\end{aligned}
$$

The proof of (3.18) is thus complete.

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

**Lemma 3.11** (Estimate of $\Theta_{\mathcal{T}}(\nabla \bar{u})$ in the $H^2$ case). Let $\bar{u} \in H^2(\Omega)$. Then

$$ \Theta_{\mathcal{T}}(\nabla \bar{u}) \le \frac{h_{\mathcal{T}}}{\theta_{\mathcal{T}}^{d/2}} \|\bar{u}\|_{H^2}. \quad (3.19) $$

*Proof.* Let $K \in \mathcal{T}$ be a given control volume, and let $h_K$ be the diameter of $K$. By definition of $\theta_{\mathcal{T}}$, for any edge or face $\sigma$ of $K$, $d_{K,\sigma} \ge \theta_{\mathcal{T}}h_K$; thus, the ball with center $x_K$ and radius $\theta_{\mathcal{T}}h_K$ has $d$-dimensional measure $C_d \times (\theta_{\mathcal{T}}h_K)^d$ and is included in $K$. Now $|K| \ge C_d \inf_{\sigma \in \mathcal{E}_K} d_{K,\sigma}^d$, and by definition of $\theta_{\mathcal{T}}$, $\inf_{\sigma \in \mathcal{E}_K} d_{K,\sigma} \ge \theta_{\mathcal{T}}h_K$; therefore, $h_K^d \le \frac{|K|}{C_d \theta_{\mathcal{T}}^d}$; applying (3.18) to each of the components of $\nabla \bar{u}$ provides (3.19). □

**Lemma 3.12** (Error estimate in the $H^2$ case). Let $d \le 3$, $\bar{u} \in H^2(\Omega)$ and $\mathbf{F} = 0$. Letting $u \in X_{\mathcal{T}}$ be defined by $u_K = \bar{u}(x_K)$, Then there exists $C$, only depending on $d$, $\Omega$ and continuously depending on $\theta_{\mathcal{T}}$, such that

$$ \|u - \bar{u}\|_{L^2(\Omega)} + \|\hat{\nabla}_{\mathcal{T}}u - \nabla \bar{u}\|_{L^2(\Omega)^d} \le C h_{\mathcal{T}} \|\bar{u}\|_{H^2}. $$

*Proof.* In this proof, we denote by $C_i$, for $i \in \mathbb{N}$, various real functions only depending on $d$, $\Omega$ and continuously depending on $\theta_{\mathcal{T}}$. Let us first observe that

$$ \inf_{v \in X_{\mathcal{T}}} \delta_{\mathcal{T}}(\bar{u}, v) \le C_1 h_{\mathcal{T}} \|\bar{u}\|_{H^2}. \quad (3.20) $$

Indeed, following the computations of [18, Theorem 9.4] (which involves the regularity factor $\theta_{\mathcal{T}}$), and denoting $\tilde{u} \in X_{\mathcal{T}}$ defined by $\tilde{u}_K = \bar{u}(x_K)$ and $\tilde{u}_{\sigma}$ such that

$$ |\sigma| \frac{\tilde{u}_{\sigma} - \tilde{u}_K}{d_{K,\sigma}} + |\sigma| \frac{\tilde{u}_{\sigma} - \tilde{u}_L}{d_{L,\sigma}} = 0, \text{ for all } \sigma \in \mathcal{E}_K \cap \mathcal{E}_L, $$

we get the existence of $C_1$ such that

$$ \delta_{\mathcal{T}}(\bar{u}, \tilde{u}) \le C_1 h_{\mathcal{T}} \|\bar{u}\|_{H^2}. $$

Since $\nabla \bar{u} \in H^1(\Omega)^d \subset \mathbf{H}_{\text{div}}(\Omega)$, we have, for some $v \in X_{\mathcal{T}}$ with $\|v\|_{\mathcal{T}} = 1$, using the Cauchy-Schwarz inequality,

$$ 
\begin{aligned}
\int_{\Omega} (\operatorname{div} \nabla \bar{u}(x)v(x) + \nabla \bar{u}(x) \cdot \nabla_{\mathcal{T}} v(x)) dx &= \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K} |\sigma|(v_{\sigma} - v_K)(\overline{\nabla} \bar{u}_{K,\sigma} \cdot \mathbf{n}_{K,\sigma} - \overline{\nabla} \bar{u}_{\sigma} \cdot \mathbf{n}_{K,\sigma}) \\
&\le \|v\|_{\mathcal{T}} \left( \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K} |\sigma| d_{K,\sigma} |\overline{\nabla} \bar{u}_{K,\sigma} - \overline{\nabla} \bar{u}_{\sigma}|^2 \right)^{1/2},
\end{aligned}
$$

where

$$ \overline{\nabla} \bar{u}_{K,\sigma} = \frac{1}{|D_{K,\sigma}|} \int_{D_{K,\sigma}} \nabla \bar{u}(x) dx \text{ and } \overline{\nabla} \bar{u}_{\sigma} = \frac{1}{|\sigma|} \int_{\sigma} \nabla \bar{u}(x) ds. $$
Applying the change of variable (with $t$ and $y$ fixed) $z = tx + (1-t)y$, noticing that $z \in V \cap B(y, th_V)$, and again applying the Fubini-Tonelli theorem, we thus get that

$$ 
\begin{aligned}
\int_V \int_V (\varphi(y) - \varphi(x))^2 dx dy &\le h_V^2 \int_V \left[ \int_0^1 \left( \int_V 1_{B(z,th_V)}(y) dy \right) t^{-d} dt \right] |\nabla \varphi(z)|^2 dz \\
&\le h_V^2 C_d h_V^d \int_V |\nabla \varphi(z)|^2 dz.
\end{aligned}
$$

The proof of (3.18) is thus complete.

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

**Lemma 3.11** (Estimate of $\Theta_{\mathcal{T}}(\nabla \bar{u})$ in the $H^2$ case). Let $\bar{u} \in H^2(\Omega)$. Then

$$ \Theta_{\mathcal{T}}(\nabla \bar{u}) \le \frac{h_{\mathcal{T}}}{\theta_{\mathcal{T}}^{d/2}} \|\bar{u}\|_{H^2}. \quad (3.19) $$

*Proof.* Let $K \in \mathcal{T}$ be a given control volume, and let $h_K$ be the diameter of $K$. By definition of $\theta_{\mathcal{T}}$, for any edge or face $\sigma$ of $K$, $d_{K,\sigma} \ge \theta_{\mathcal{T}}h_K$; thus, the ball with center $x_K$ and radius $\theta_{\mathcal{T}}h_K$ has $d$-dimensional measure $C_d \times (\theta_{\mathcal{T}}h_K)^d$ and is included in $K$. Now $|K| \ge C_d \inf_{\sigma \in \mathcal{E}_K} d_{K,\sigma}^d$, and by definition of $\theta_{\mathcal{T}}$, $\inf_{\sigma \in \mathcal{E}_K} d_{K,\sigma} \ge \theta_{\mathcal{T}}h_K$; therefore, $h_K^d \le \frac{|K|}{C_d \theta_{\mathcal{T}}^d}$; applying (3.18) to each of the components of $\nabla \bar{u}$ provides (3.19).

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

**Lemma 3.12** (Error estimate in the $H^2$ case). Let $d \le 3$, $\bar{u} \in H^2(\Omega)$ and $\mathbf{F} = 0$. Letting $u \in X_{\mathcal{T}}$ be defined by $u_K = \bar{u}(x_K)$, Then there exists $C$, only depending on $d$, $\Omega$ and continuously depending on $\theta_{\mathcal{T}}$, such that

$$ \|u - \bar{u}\|_{L^2(\Omega)} + \|\hat{\nabla}_{\mathcal{T}}u - \nabla \bar{u}\|_{L^2(\Omega)^d} \le C h_{\mathcal{T}} \|\bar{u}\|_{H^2}. $$

*Proof.* In this proof, we denote by $C_i$, for $i \in \mathbb{N}$, various real functions only depending on $d$, $\Omega$ and continuously depending on $\theta_{\mathcal{T}}$. Let us first observe that

$$ \inf_{v \in X_{\mathcal{T}}} \delta_{\mathcal{T}}(\bar{u}, v) \le C_1 h_{\mathcal{T}} \|\bar{u}\|_{H^2}. \quad (3.20) $$

Indeed, following the computations of [18, Theorem 9.4] (which involves the regularity factor $\theta_{\mathcal{T}}$), and denoting $\tilde{u} \in X_{\mathcal{T}}$ defined by $\tilde{u}_K = \bar{u}(x_K)$ and $\tilde{u}_{\sigma}$ such that

$$ |\sigma| \frac{\tilde{u}_{\sigma} - \tilde{u}_K}{d_{K,\sigma}} + |\sigma| \frac{\tilde{u}_{\sigma} - \tilde{u}_L}{d_{L,\sigma}} = 0, \text{ for all } \sigma \in \mathcal{E}_K \cap \mathcal{E}_L, $$

we get the existence of $C_1$ such that

$$ \delta_{\mathcal{T}}(\bar{u}, \tilde{u}) \le C_1 h_{\mathcal{T}} \|\bar{u}\|_{H^2}. $$

Since $\nabla \bar{u} \in H^1(\Omega)^d \subset \mathbf{H}_{\text{div}}(\Omega)$, we have, for some $v \in X_{\mathcal{T}}$ with $\|v\|_{\mathcal{T}} = 1$, using the Cauchy-Schwarz inequality,

$$ 
\begin{aligned}
\int_{\Omega} (\operatorname{div} \nabla \bar{u}(x)v(x) + \nabla \bar{u}(x) \cdot \nabla_{\mathcal{T}} v(x)) dx &= \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K} |\sigma|(v_{\sigma} - v_K)(\overline{\nabla} \bar{u}_{K,\sigma} \cdot \mathbf{n}_{K,\sigma} - \overline{\nabla} \bar{u}_{\sigma} \cdot \mathbf{n}_{K,\sigma}) \\
&\le \|v\|_{\mathcal{T}} \left( \sum_{K \in \mathcal{T}} \sum_{\sigma \in \mathcal{E}_K} |\sigma| d_{K,\sigma} |\overline{\nabla} \bar{u}_{K,\sigma} - \overline{\nabla} \bar{u}_{\sigma}|^2 \right)^{1/2},
\end{aligned}
$$

where

$$ \overline{\nabla} \bar{u}_{K,\sigma} = \frac{1}{|D_{K,\sigma}|} \int_{D_{K,\sigma}} \nabla \bar{u}(x) dx \text{ and } \overline{\nabla} \bar{u}_{\sigma} = \frac{1}{|\sigma|} \int_{\sigma} \nabla \bar{u}(x) ds. $$
Hence

$$ \zeta_{\mathcal{T}}(\nabla \bar{u}) \leq C_2 h_{\mathcal{T}} \|\bar{u}\|_{H^2}. $$

Using (3.19), (3.10), (3.20) and (3.14), the conclusion follows.

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

**Remark 3.13 (Superconvergence).** The superconvergence of the TPFA scheme on a class of acute triangles in 2D was observed numerically long ago [22, 4, 21], and was proved theoretically in [13] more recently; the proof of [13] relies on the fact that on a large class of 2D triangular grids, the TPFA scheme is a hybrid mixed method as defined in [10].

## 4. NUMERICAL EXAMPLE OF A MINIMAL REGULARITY PROBLEM

We illustrate here the error estimate of Theorem 3.1 in the case where the regularity of the solution of Problem (2.3) is minimal. Let $\Omega = (0, 1)^2$, and let $\bar{u} \in H_0^1(\Omega)$ be defined by

$$ \bar{u}(x) = \left( -\log \left( \max \left( |x_1 - \frac{1}{2}|, |x_2 - \frac{1}{2}| \right) \right) \right)^{\gamma} - \left( -\log(r_0) \right)^{\gamma}, \quad x = (x_1, x_2) \in \Omega, \quad (4.1) $$

with $\gamma = \frac{1}{4}$ and $r_0 = 1/2$. Then $\bar{u}$ is solution to (2.3), letting $f = 0$ and $F(x) = -\nabla \bar{u}(x)$. We have, for any $p > 1$,

$$ \int_{\Omega} |\nabla \bar{u}(x)|^p dx = 4 \int_{0}^{r_0} \int_{-x}^{x} \frac{\gamma^p}{x^p} (-\log(x))^{p(\gamma-1)} dy dx = 8 \int_{0}^{r_0} \frac{\gamma^p}{x^{p-1}} (-\log(x))^{p(\gamma-1)} dy dx. $$

This integral is infinite for any $p > 2$, which shows that $\bar{u} \notin W_0^{1,p}(\Omega)$ if $p > 2$; this was expected since $\bar{u}$ is discontinuous and $W_0^{1,p}(\Omega) \subset C^0(\bar{\Omega})$ for any $p > 2$ in 2D. Now if $p = 2$, we get, since $2\gamma - 1 < 0$,

$$ \int_{\Omega} |\nabla \bar{u}(x)|^2 dx = 8 \int_{0}^{r_0} \frac{\gamma^2}{x} (-\log(x))^{2\gamma-2} dx = \frac{8\gamma^2}{1-2\gamma} (-\log(r_0))^{2\gamma-1}, $$

which shows that $\nabla \bar{u}(x) \in L^2(\Omega)^2$ with $\|\nabla \bar{u}\|_{L^2} \simeq 1.0959573$. Let us now compute $\|\bar{u}\|_{L^2}$. We have

$$ \int_{\Omega} \bar{u}(x)^2 dx = 4 \int_{0}^{r_0} \int_{-x}^{x} (-\log(x))^{2\gamma} dy dx - 2 (-\log(r_0))^{\gamma} 4 \int_{0}^{r_0} \int_{-x}^{x} (-\log(x))^{\gamma} dy dx + (-\log(r_0))^{2\gamma}, $$

which gives

$$ \int_{\Omega} \bar{u}(x)^2 dx = 8 \int_{0}^{r_0} x (-\log(x))^{2\gamma} dx - 2 (-\log(r_0))^{\gamma} 8 \int_{0}^{r_0} x (-\log(x))^{\gamma} dx + (2r_0)^2 (-\log(r_0))^{2\gamma}. $$

Owing to the change of variable $-\log(x) = y/2$, we get

$$ \int_{0}^{r_0} x (-\log(x))^{2\gamma} dx = \frac{1}{2^{2\gamma+1}} \int_{-2\log(r_0)}^{\infty} y^{2\gamma} e^{-y} dy. $$

Introducing the upper incomplete gamma function $\Gamma$ defined by

$$ \Gamma(a, x) = \int_{x}^{+\infty} t^{a-1} e^{-t} dt, $$

we get

$$ \int_{0}^{r_0} x (-\log(x))^{2\gamma} dx = \frac{1}{2^{2\gamma+1}} \Gamma(2\gamma+1, -2\log(r_0)). $$
Hence

$$ \zeta_{\mathcal{T}}(\nabla \bar{u}) \leq C_2 h_{\mathcal{T}} \|\bar{u}\|_{H^2}. $$

Using (3.19), (3.10), (3.20) and (3.14), the conclusion follows.

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

**Remark 3.13 (Superconvergence).** The superconvergence of the TPFA scheme on a class of acute triangles in 2D was observed numerically long ago [22, 4, 21], and was proved theoretically in [13] more recently; the proof of [13] relies on the fact that on a large class of 2D triangular grids, the TPFA scheme is a hybrid mixed method as defined in [10].

## 4. NUMERICAL EXAMPLE OF A MINIMAL REGULARITY PROBLEM

We illustrate here the error estimate of Theorem 3.1 in the case where the regularity of the solution of Problem (2.3) is minimal. Let $\Omega = (0, 1)^2$, and let $\bar{u} \in H_0^1(\Omega)$ be defined by

$$ \bar{u}(x) = \left( -\log\left(\max\left(|x_1 - \frac{1}{2}|, |x_2 - \frac{1}{2}|\right)\right) \right)^{\gamma} - \left(-\log(r_0)\right)^{\gamma}, \quad x = (x_1, x_2) \in \Omega, \quad (4.1) $$

with $\gamma = \frac{1}{4}$ and $r_0 = 1/2$. Then $\bar{u}$ is solution to (2.3), letting $f = 0$ and $F(x) = -\nabla \bar{u}(x)$. We have, for any $p > 1$,

$$ \int_{\Omega} |\nabla \bar{u}(x)|^p dx = 4 \int_{0}^{r_0} \int_{-x}^{x} \frac{\gamma^p}{x^p} (-\log(x))^{p(\gamma-1)} dy dx = 8 \int_{0}^{r_0} \frac{\gamma^p}{x^{p-1}} (-\log(x))^{p(\gamma-1)} dy dx. $$

This integral is infinite for any $p > 2$, which shows that $\bar{u} \notin W_0^{1,p}(\Omega)$ if $p > 2$; this was expected since $\bar{u}$ is discontinuous and $W_0^{1,p}(\Omega) \subset C^0(\bar{\Omega})$ for any $p > 2$ in 2D. Now if $p = 2$, we get, since $2\gamma - 1 < 0$,

$$ \int_{\Omega} |\nabla \bar{u}(x)|^2 dx = 8 \int_{0}^{r_0} \frac{\gamma^2}{x} (-\log(x))^{2\gamma-2} dx = \frac{8\gamma^2}{1-2\gamma} (-\log(r_0))^{2\gamma-1}, $$

which shows that $\nabla \bar{u}(x) \in L^2(\Omega)^2$ with $\|\nabla \bar{u}\|_{L^2} \simeq 1.0959573$. Let us now compute $\|\bar{u}\|_{L^2}$. We have

$$ \int_{\Omega} \bar{u}(x)^2 dx = 4 \int_{0}^{r_0} \int_{-x}^{x} (-\log(x))^{2\gamma} dy dx - 2 (-\log(r_0))^{\gamma} 4 \int_{0}^{r_0} \int_{-x}^{x} (-\log(x))^{\gamma} dy dx + (-\log(r_0))^{2\gamma}, $$

which gives

$$ \int_{\Omega} \bar{u}(x)^2 dx = 8 \int_{0}^{r_0} x (-\log(x))^{2\gamma} dx - 2 (-\log(r_0))^{\gamma} 8 \int_{0}^{r_0} x (-\log(x))^{\gamma} dx + (2r_0)^2 (-\log(r_0))^{2\gamma}. $$

Owing to the change of variable $-\log(x) = y/2$, we get

$$ \int_{0}^{r_0} x (-\log(x))^{2\gamma} dx = \frac{1}{2^{2\gamma+1}} \int_{-2\log(r_0)}^{\infty} y^{2\gamma} e^{-y} dy. $$

Introducing the upper incomplete gamma function $\Gamma$ defined by

$$ \Gamma(a, x) = \int_{x}^{+\infty} t^{a-1} e^{-t} dt, $$

we get

$$ \int_{0}^{r_0} x (-\log(x))^{2\gamma} dx = \frac{1}{2^{2\gamma+1}} \Gamma(2\gamma+1, -2\log(r_0)). $$
Again using the change of variable $-\log(x) = y/2$, we have

$$ \int_{0}^{r_0} x(-\log(x))^{\gamma} dx = \frac{1}{2^{\gamma+1}} \int_{-2\log(r_0)}^{\infty} y^{\gamma} e^{-y} dy = \frac{1}{2^{\gamma+1}} \Gamma(\gamma+1, -2\log(r_0)) $$

Finally we obtain

$$ \int_{\Omega} \bar{u}(x)^2 dx = 2^{2-2\gamma} \Gamma(2\gamma + 1, -2\log(r_0)) \\ - 2^{3-\gamma} (-\log(r_0))^{\gamma} \Gamma(\gamma + 1, -2\log(r_0)) + (2r_0)^2 (-\log(r_0))^{2\gamma}, $$

leading to $\|\bar{u}\|_{L^2} \simeq 0.1519926$.

Let us then set $f = 0$ and $F(x) = -\nabla \bar{u}(x)$ in the finite volume scheme (2.10). Since $\nabla \bar{u} + F = 0$, we get that the conformity error $\zeta_{\mathcal{T}}(\nabla \bar{u} + F)$ vanishes. Hence the numerical validation of the error estimate of Theorem 3.1 follows by a comparison of the approximation error $\delta_{\mathcal{T}}(\bar{u}, u)$ and the upper bound of the interpolation error $\delta_{\mathcal{T}}(\bar{u}, \bar{u}_{\mathcal{T}})$ for the interpolation $\bar{u}_{\mathcal{T}} \in X_{\mathcal{T}}$ of the exact solution $\bar{u}$, defined by the set of values $((\bar{u}(x_K))_{K \in \mathcal{T}}, (\bar{u}(x_{\sigma}))_{\sigma \in \mathcal{E}})$, where $x_{\sigma}$ is the middle of the edge $\sigma$ (it is then the intersection of the segment $[x_K, x_L]$ and of $\sigma$ when $\sigma$ is the common edge between $K$ and $L$).

The right-hand side $\int_{\Omega} \nabla \bar{u}(x) \cdot \nabla_{\mathcal{T}} v(x) dx$ can be computed by using the relation

$$ \int_{D_{K,\sigma}} \nabla \bar{u}(x) \cdot \mathbf{n}_{K,\sigma} dx = \int_{\partial D_{K,\sigma}} \bar{u}(x) \mathbf{n}(y) \cdot \mathbf{n}_{K,\sigma} ds(y), $$

where $\mathbf{n}(y)$ is the outward unit vector normal to $\partial D_{K,\sigma}$ at the point $y$. This latter integral is then a linear combination of terms under the form

$$ \int_0^1 (-\log(\alpha + (\beta - \alpha)s))^\gamma ds = -\frac{1}{\beta - \alpha} \int_{-\log(\alpha)}^{-\log(\beta)} t^\gamma e^{-t} dt \\ = \frac{1}{\beta - \alpha} \left( \Gamma(\gamma + 1, -\log(\beta)) - \Gamma(\gamma + 1, -\log(\alpha)) \right), $$

Hence we can accurately compute the right-hand side for the computation of the approximate solution, even in the case where the elements have the point $(1/2, 1/2)$ as vertex since we can extend by continuity the value of the function $g(\alpha) = \Gamma(a, -\log(\alpha))$ to $0$ by setting $g(0) = 0$.

Note that the main part of $\|\nabla \bar{u}\|_{L^2(\Omega)^2}$ is concentrated inside a small ball. Indeed, denoting by $B_{\infty}(r) = \{x \in \Omega, \max(|x_1 - \frac{1}{2}|, |x_2 - \frac{1}{2}|) < r\}$ (hence $\Omega = B_{\infty}(r_0)$), we get the following values for $\|\nabla \bar{u}\|_{L^2(B_{\infty}(r))^2}$ and $\|\bar{u}\|_{L^2(B_{\infty}(r))}$:

<table><thead><tr><th>r</th><th>0.5</th><th>exp(-10)</th><th>exp(-1000)</th><th>exp(-1000000)</th><th>exp(-1000000000)</th></tr></thead><tbody><tr><td>||∇u||<sub>L<sup>2</sup>(B<sub>∞</sub>(r))<sup>2</sup></sub></td><td>1.096</td><td>0.562</td><td>0.178</td><td>0.032</td><td>0.006</td></tr><tr><td>||u||<sub>L<sup>2</sup>(B<sub>∞</sub>(r))</sub></td><td>0.152</td><td>2.71e-6</td><td>&lt; 1.e-200</td><td>&lt; 1.e-200000</td><td>&lt;1.e-200000000</td></tr></tbody></table>

It is clear from this table that $\lim_{r \to 0} \|\nabla \bar{u}\|_{L^2(B_{\infty}(r))^2} = 0$; however, this convergence is so slow that a numerical scheme based on a standard mesh cannot be expected to show a numerical convergence of the approximation of $\nabla \bar{u}$ for the $L^2$ norm, even with a refined mesh at the point $(0.5, 0.5)$ (see Remark 3.2). Nevertheless, we can numerically assess the result proven in Theorem 3.1 by considering the family of benchmark meshes [21, Mesh1_k] with $k = 1, 2, 3, 4, 5, 6, 7$, which can be respectively characterized by the values $h_{\mathcal{T}} = \frac{1}{4}, \frac{1}{8}, \frac{1}{16}, \frac{1}{32}, \frac{1}{64}, \frac{1}{128}, \frac{1}{256}$ for the size of the
Again using the change of variable $-\log(x) = y/2$, we have

$$ \int_{0}^{r_0} x(-\log(x))^{\gamma} dx = \frac{1}{2^{\gamma+1}} \int_{-2\log(r_0)}^{\infty} y^{\gamma} e^{-y} dy = \frac{1}{2^{\gamma+1}} \Gamma(\gamma+1, -2\log(r_0)) $$

Finally we obtain

$$ \int_{\Omega} \bar{u}(x)^2 dx = 2^{2-2\gamma} \Gamma(2\gamma+1, -2\log(r_0)) \\ - 2^{3-\gamma} (-\log(r_0))^{\gamma} \Gamma(\gamma+1, -2\log(r_0)) + (2r_0)^2 (-\log(r_0))^{2\gamma}, $$

leading to $\|\bar{u}\|_{L^2} \simeq 0.1519926$.

Let us then set $f = 0$ and $F(x) = -\nabla \bar{u}(x)$ in the finite volume scheme (2.10). Since $\nabla \bar{u} + F = 0$, we get that the conformity error $\zeta_{\mathcal{T}}(\nabla \bar{u} + F)$ vanishes. Hence the numerical validation of the error estimate of Theorem 3.1 follows by a comparison of the approximation error $\delta_{\mathcal{T}}(\bar{u}, u)$ and the upper bound of the interpolation error $\delta_{\mathcal{T}}(\bar{u}, \bar{u}_{\mathcal{T}})$ for the interpolation $\bar{u}_{\mathcal{T}} \in X_{\mathcal{T}}$ of the exact solution $\bar{u}$, defined by the set of values $((\bar{u}(x_K))_{K \in \mathcal{T}}, (\bar{u}(x_{\sigma}))_{\sigma \in \mathcal{E}})$, where $x_{\sigma}$ is the middle of the edge $\sigma$ (it is then the intersection of the segment $[x_K, x_L]$ and of $\sigma$ when $\sigma$ is the common edge between $K$ and $L$).

The right-hand side $\int_{\Omega} \nabla \bar{u}(x) \cdot \nabla_{\mathcal{T}} v(x) dx$ can be computed by using the relation

$$ \int_{D_{K,\sigma}} \nabla \bar{u}(x) \cdot \mathbf{n}_{K,\sigma} dx = \int_{\partial D_{K,\sigma}} \bar{u}(x) \mathbf{n}(y) \cdot \mathbf{n}_{K,\sigma} ds(y), $$

where $\mathbf{n}(y)$ is the outward unit vector normal to $\partial D_{K,\sigma}$ at the point $y$. This latter integral is then a linear combination of terms under the form

$$ \int_0^1 (-\log(\alpha + (\beta - \alpha)s))^\gamma ds = -\frac{1}{\beta - \alpha} \int_{-\log(\alpha)}^{-\log(\beta)} t^\gamma e^{-t} dt \\ = \frac{1}{\beta - \alpha} \left( \Gamma(\gamma + 1, -\log(\beta)) - \Gamma(\gamma + 1, -\log(\alpha)) \right), $$

Hence we can accurately compute the right-hand side for the computation of the approximate solution, even in the case where the elements have the point $(1/2, 1/2)$ as vertex since we can extend by continuity the value of the function $g(\alpha) = \Gamma(a, -\log(\alpha))$ to $0$ by setting $g(0) = 0$.

Note that the main part of $\|\nabla \bar{u}\|_{L^2(\Omega)^2}$ is concentrated inside a small ball. Indeed, denoting by $B_{\infty}(r) = \{x \in \Omega, \max(|x_1 - \frac{1}{2}|, |x_2 - \frac{1}{2}|) < r\}$ (hence $\Omega = B_{\infty}(r_0)$), we get the following values for $\|\nabla \bar{u}\|_{L^2(B_{\infty}(r))^2}$ and $\|\bar{u}\|_{L^2(B_{\infty}(r))}$:

<table><thead><tr><th>r</th><th>0.5</th><th>exp(-10)</th><th>exp(-1000)</th><th>exp(-1000000)</th><th>exp(-1000000000)</th></tr></thead><tbody><tr><td>||∇u||<sub>L<sup>2</sup>(B<sub>∞</sub>(r))<sup>2</sup></sub></td><td>1.096</td><td>0.562</td><td>0.178</td><td>0.032</td><td>0.006</td></tr><tr><td>||u||<sub>L<sup>2</sup>(B<sub>∞</sub>(r))</sub></td><td>0.152</td><td>2.71e-6</td><td>&lt; 1.e-200</td><td>&lt; 1.e-200000</td><td>&lt;1.e-200000000</td></tr></tbody></table>

It is clear from this table that $\lim_{r \to 0} \|\nabla \bar{u}\|_{L^2(B_{\infty}(r))^2} = 0$; however, this convergence is so slow that a numerical scheme based on a standard mesh cannot be expected to show a numerical convergence of the approximation of $\nabla \bar{u}$ for the $L^2$ norm, even with a refined mesh at the point $(0.5, 0.5)$ (see Remark 3.2). Nevertheless, we can numerically assess the result proven in Theorem 3.1 by considering the family of benchmark meshes [21, Mesh1_k] with $k = 1, 2, 3, 4, 5, 6, 7$, which can be respectively characterized by the values $h_{\mathcal{T}} = \frac{1}{4}, \frac{1}{8}, \frac{1}{16}, \frac{1}{32}, \frac{1}{64}, \frac{1}{128}, \frac{1}{256}$ for the size of the
<figure><img src="image_5.png" /><figcaption>FIGURE 2. The approximate solution of Problem (2.3) computed on Mesh1_3.</figcaption></figure>

mesh. The following table gives the number of control volumes (cv), vertices and edges for each mesh.

<table><thead><tr><th>mesh</th><th>h<sub>τ</sub></th><th># cv</th><th># vertices</th><th># edges</th></tr></thead><tbody><tr><td>Mesh1_1</td><td><sup>1</sup>&frasl;<sub>4</sub></td><td>56</td><td>37</td><td>92</td></tr><tr><td>Mesh1_2</td><td><sup>1</sup>&frasl;<sub>8</sub></td><td>224</td><td>129</td><td>352</td></tr><tr><td>Mesh1_3</td><td><sup>1</sup>&frasl;<sub>16</sub></td><td>896</td><td>481</td><td>1376</td></tr><tr><td>Mesh1_4</td><td><sup>1</sup>&frasl;<sub>32</sub></td><td>3584</td><td>1857</td><td>5440</td></tr><tr><td>Mesh1_5</td><td><sup>1</sup>&frasl;<sub>64</sub></td><td>14336</td><td>7297</td><td>21632</td></tr><tr><td>Mesh1_6</td><td><sup>1</sup>&frasl;<sub>128</sub></td><td>57344</td><td>28929</td><td>86272</td></tr><tr><td>Mesh1_7</td><td><sup>1</sup>&frasl;<sub>256</sub></td><td>229376</td><td>115201</td><td>344576</td></tr></tbody></table>

We show in Figure 2 a numerical solution obtained on Mesh1_3.

Note that these meshes share the same regularity factor $\theta_{\tau}$ but do not feature any symmetries that might increase the order of convergence.

For each mesh, we consider the interpolation $\bar{u}_{\tau} \in X_{\tau}$ of $\bar{u}$ defined as specified above. Figure 3 provides the approximation error and the interpolation error obtained on each mesh.

We observe that, as in the regular case, the error $\|\bar{u}_{\tau} - u\|_{L^2}$ is much smaller than $\|\bar{u} - u\|_{L^2}$ (recall that superconvergence is observed on regular solutions), and that we observe an order 1
<figure><img src="image_6.png" /><figcaption>FIGURE 2. The approximate solution of Problem (2.3) computed on Mesh1_3.</figcaption></figure>

mesh. The following table gives the number of control volumes (cv), vertices and edges for each mesh.

<table><thead><tr><th>mesh</th><th>h<sub>τ</sub></th><th># cv</th><th># vertices</th><th># edges</th></tr></thead><tbody><tr><td>Mesh1_1</td><td><sup>1</sup>&frasl;<sub>4</sub></td><td>56</td><td>37</td><td>92</td></tr><tr><td>Mesh1_2</td><td><sup>1</sup>&frasl;<sub>8</sub></td><td>224</td><td>129</td><td>352</td></tr><tr><td>Mesh1_3</td><td><sup>1</sup>&frasl;<sub>16</sub></td><td>896</td><td>481</td><td>1376</td></tr><tr><td>Mesh1_4</td><td><sup>1</sup>&frasl;<sub>32</sub></td><td>3584</td><td>1857</td><td>5440</td></tr><tr><td>Mesh1_5</td><td><sup>1</sup>&frasl;<sub>64</sub></td><td>14336</td><td>7297</td><td>21632</td></tr><tr><td>Mesh1_6</td><td><sup>1</sup>&frasl;<sub>128</sub></td><td>57344</td><td>28929</td><td>86272</td></tr><tr><td>Mesh1_7</td><td><sup>1</sup>&frasl;<sub>256</sub></td><td>229376</td><td>115201</td><td>344576</td></tr></tbody></table>

We show in Figure 2 a numerical solution obtained on Mesh1_3.

Note that these meshes share the same regularity factor $\theta_{\tau}$ but do not feature any symmetries that might increase the order of convergence.

For each mesh, we consider the interpolation $\bar{u}_{\tau} \in X_{\tau}$ of $\bar{u}$ defined as specified above. Figure 3 provides the approximation error and the interpolation error obtained on each mesh.

We observe that, as in the regular case, the error $\|\bar{u}_{\tau} - u\|_{L^2}$ is much smaller than $\|\bar{u} - u\|_{L^2}$ (recall that superconvergence is observed on regular solutions), and that we observe an order 1
in the convergence of $||\bar{u} - u||_{L^2}$ to 0. When accounting for the error on the gradient, we get that the approximation error behaves similarly for the small values of $h_{\tau}$ as the interpolation error. Nevertheless, Figure 3 does not show that the quantities $\delta_{\tau}(\bar{u}, \bar{u}_{\tau})$ and $||G\bar{u} - G_{\tau}u||_{L^2}$ behave as $(h_{\tau})^{\alpha}$ for some $\alpha > 0$. Our conjecture is that these terms behave as $(-\log(h_{\tau}))^{-\alpha}$ but the range of mesh sizes that can be tested is not sufficient to assess it.

<figure><img src="image_6.png" /><figcaption>FIGURE 3. Approximation and interpolation errors vs. the size of the mesh.</figcaption></figure>

5. CONCLUSION AND PERSPECTIVES

In the present study, we were able to show an optimal error estimate for the discretization of a Poisson equation with minimal regularity, namely when the solution possesses no greater regularity than belonging to $H_0^1(\Omega)$, even though the two-point flux approximation (TPFA) cannot be placed in the framework of the gradient discretization method for which such an estimate already exists [9]. Error bounds are derived for both the solution and the approximation of the gradient component orthogonal to the mesh faces. Our estimates are optimal, in the sense that the approximation error is shown to be of the same order as the sum of the interpolation error and the conformity error. Furthermore, numerical experiments show the validity of our error for a solution with minimal regularity. The results are extended to evolution problems discretized via the implicit Euler scheme, as detailed in the accompanying appendix.

While our study sheds some light on error estimates for finite volume schemes on admissible meshes for irregular and regular right-hand sides, the superconvergence of the method for
in the convergence of $||\bar{u} - u||_{L^2}$ to 0. When accounting for the error on the gradient, we get that the approximation error behaves similarly for the small values of $h_{\tau}$ as the interpolation error. Nevertheless, Figure 3 does not show that the quantities $\delta_{\tau}(\bar{u}, u_{\tau})$ and $||G\bar{u} - G_{\tau}u||_{L^2}$ behave as $(h_{\tau})^{\alpha}$ for some $\alpha > 0$. Our conjecture is that these terms behave as $(-\log(h_{\tau}))^{-\alpha}$ but the range of mesh sizes that can be tested is not sufficient to assess it.

<figure><img src="image_7.png" /><figcaption>FIGURE 3. Approximation and interpolation errors vs. the size of the mesh.</figcaption></figure>

## 5. CONCLUSION AND PERSPECTIVES

In the present study, we were able to show an optimal error estimate for the discretization of a Poisson equation with minimal regularity, namely when the solution possesses no greater regularity than belonging to $H_0^1(\Omega)$, even though the two-point flux approximation (TPFA) cannot be placed in the framework of the gradient discretization method for which such an estimate already exists [9]. Error bounds are derived for both the solution and the approximation of the gradient component orthogonal to the mesh faces. Our estimates are optimal, in the sense that the approximation error is shown to be of the same order as the sum of the interpolation error and the conformity error. Furthermore, numerical experiments show the validity of our error for a solution with minimal regularity. The results are extended to evolution problems discretized via the implicit Euler scheme, as detailed in the accompanying appendix.

While our study sheds some light on error estimates for finite volume schemes on admissible meshes for irregular and regular right-hand sides, the superconvergence of the method for
regular right-hand sides remains an open problem; indeed, even though numerical experiments demonstrate that the two point scheme superconverges on triangular [4, 22] and Voronoï meshes [24, 14, 15], the theoretical proof remains an open problem except for the case of a (rather large) class of triangular meshes [13] or for a double finite volume scheme known as “DDFV” which involves two TPFA schemes, one on the Delaunay triangulation and one on the associated Voronoï mesh [25].

## APPENDIX A. THE TRANSIENT CASE

This appendix is focused on the study of an optimal error bound in the transient case. For a given $T > 0$, the solution of the continuous problem is $u: \Omega \times [0, T] \to \mathbb{R}$ such that

$$ \partial_t u - \Delta u = L \qquad (A.1) $$

$$ u(\cdot, t) = 0 \text{ on } \partial\Omega \text{ for } t \in [0, T], \qquad (A.2) $$

$$ u(0) - \Phi u(T) = \xi_0, \qquad (A.3) $$

where $L \in L^2(0, T; H^{-1}(\Omega))$, $\xi_0 \in L^2(\Omega)$ and $\Phi: L^2(\Omega) \to L^2(\Omega)$ is a linear contractive mapping, i.e. such that $\|\Phi v\|_{L^2(\Omega)} \le \|v\|_{L^2(\Omega)}$ for all $v \in L^2(\Omega)$. Note that if $\Phi = 0$, the problem (A.1)-(A.3) is the standard Cauchy problem, and if $\Phi = \text{Id}$, (A.3) becomes a periodic condition.

Identifying $L^2(\Omega)$ with its dual space thanks to the Riesz-Fréchet representation theorem, we have

$$ H_0^1(\Omega) \subset L^2(\Omega) = (L^2(\Omega))' \subset (H_0^1(\Omega))' = H^{-1}(\Omega), $$

so that the space

$$ W(0, T; \Omega) = \{u \in L^2(0, T; H_0^1(\Omega)) : u' \in L^2(0, T; H^{-1}(\Omega))\}, $$

is well defined, with $u'$ the weak derivative of the function $t \mapsto u(\cdot, t)$. Recall that if $u \in W(0, T; \Omega)$, then $u \in C([0, T], L^2(\Omega))$. As in the elliptic case, we decompose $L = f + \text{div } \mathbf{F}$ with $f \in L^2(0, T; L^2(\Omega))$ and $\mathbf{F} \in L^2(0, T; L^2(\Omega)^d)$. A weak formulation of the problem (A.1)-(A.3) then reads

$$ \begin{gathered} \bar{u} \in W(0, T; \Omega) \text{ s.t. } \forall \bar{v} \in L^2(0, T; H_0^1(\Omega)), \\ \langle \bar{u}'(t), \bar{v}(t) \rangle_{H^{-1}(\Omega), H_0^1(\Omega)} + \langle \nabla \bar{u}(t) + \mathbf{F}(t), \nabla \bar{v}(t) \rangle_{L^2(\Omega)^d} = \langle f(t), \bar{v}(t) \rangle_{L^2(\Omega)} \quad \text{for a.e. } t \in [0, T] \\ \text{and } \bar{u}(0) - \Phi \bar{u}(T) = \xi_0, \end{gathered} $$

which also reads

$$ \bar{u} \in W(0, T; \Omega) \text{ s.t. } \bar{u}' - \text{div}(\nabla \bar{u} + \mathbf{F}) = f \text{ and } \bar{u}(0) - \Phi \bar{u}(T) = \xi_0. \qquad (A.4) $$

Recall that $-\Delta$ is a linear continuous bijective operator from $H_0^1(\Omega)$ to $H^{-1}(\Omega)$ and therefore also from $L^2(0, T; H_0^1(\Omega))$ to $L^2(0, T; H^{-1}(\Omega))$. We denote by $R$ its inverse which is linear continuous from $L^2(0, T; H^{-1}(\Omega))$ to $L^2(0, T; H_0^1(\Omega))$. More precisely, noting that

$$ \forall (u, v) \in H_0^1(\Omega)^2, \langle -\Delta u, v \rangle_{H^{-1}(\Omega), H_0^1(\Omega)} = \langle \nabla u, \nabla v \rangle_{L^2(\Omega)^d}, $$

and that $u = R\xi \iff \xi = -\Delta u$, the operator $R$ is characterized by

$$ \forall (\xi, v) \in H^{-1}(\Omega) \times H_0^1(\Omega), \quad \langle \nabla R \xi, \nabla v \rangle_{L^2(\Omega)^d} = \langle \xi, v \rangle_{H^{-1}(\Omega), H_0^1(\Omega)}. \qquad (A.5) $$

The problem (A.4) is then equivalent to

$$ \text{find } \bar{u} \in W(0, T; \Omega) \text{ s.t. } -\text{div}(\nabla R \bar{u}' + \nabla \bar{u} + \mathbf{F}) = f \text{ and } \bar{u}(0) - \Phi \bar{u}(T) = \xi_0, \qquad (A.6) $$
regular right-hand sides remains an open problem; indeed, even though numerical experiments demonstrate that the two point scheme superconverges on triangular [4, 22] and Voronoï meshes [24, 14, 15], the theoretical proof remains an open problem except for the case of a (rather large) class of triangular meshes [13] or for a double finite volume scheme known as “DDFV” which involves two TPFA schemes, one on the Delaunay triangulation and one on the associated Voronoï mesh [25].

## APPENDIX A. THE TRANSIENT CASE

This appendix is focused on the study of an optimal error bound in the transient case. For a given $T > 0$, the solution of the continuous problem is $u: \Omega \times [0, T] \to \mathbb{R}$ such that

$$ \partial_t u - \Delta u = L \qquad (A.1) $$

$$ u(\cdot, t) = 0 \text{ on } \partial\Omega \text{ for } t \in [0, T], \qquad (A.2) $$

$$ u(0) - \Phi u(T) = \xi_0, \qquad (A.3) $$

where $L \in L^2(0, T; H^{-1}(\Omega))$, $\xi_0 \in L^2(\Omega)$ and $\Phi: L^2(\Omega) \to L^2(\Omega)$ is a linear contractive mapping, i.e. such that $\|\Phi v\|_{L^2(\Omega)} \le \|v\|_{L^2(\Omega)}$ for all $v \in L^2(\Omega)$. Note that if $\Phi = 0$, the problem (A.1)-(A.3) is the standard Cauchy problem, and if $\Phi = \text{Id}$, (A.3) becomes a periodic condition.

Identifying $L^2(\Omega)$ with its dual space thanks to the Riesz-Fréchet representation theorem, we have

$$ H_0^1(\Omega) \subset L^2(\Omega) = (L^2(\Omega))' \subset (H_0^1(\Omega))' = H^{-1}(\Omega), $$

so that the space

$$ W(0, T; \Omega) = \{u \in L^2(0, T; H_0^1(\Omega)) : u' \in L^2(0, T; H^{-1}(\Omega))\}, $$

is well defined, with $u'$ the weak derivative of the function $t \mapsto u(\cdot, t)$. Recall that if $u \in W(0, T; \Omega)$, then $u \in C([0, T], L^2(\Omega))$. As in the elliptic case, we decompose $L = f + \text{div } \mathbf{F}$ with $f \in L^2(0, T; L^2(\Omega))$ and $\mathbf{F} \in L^2(0, T; L^2(\Omega)^d)$. A weak formulation of the problem (A.1)-(A.3) then reads

$$ \bar{u} \in W(0, T; \Omega) \text{ s.t. } \forall \bar{v} \in L^2(0, T; H_0^1(\Omega)), $$

$$ \langle \bar{u}'(t), \bar{v}(t) \rangle_{H^{-1}(\Omega), H_0^1(\Omega)} + \langle \nabla \bar{u}(t) + \mathbf{F}(t), \nabla \bar{v}(t) \rangle_{L^2(\Omega)^d} = \langle f(t), \bar{v}(t) \rangle_{L^2(\Omega)} \text{ for a.e. } t \in [0, T] $$

$$ \text{and } \bar{u}(0) - \Phi \bar{u}(T) = \xi_0, $$

which also reads

$$ \bar{u} \in W(0, T; \Omega) \text{ s.t. } \bar{u}' - \text{div}(\nabla \bar{u} + \mathbf{F}) = f \text{ and } \bar{u}(0) - \Phi \bar{u}(T) = \xi_0. \qquad (A.4) $$

Recall that $-\Delta$ is a linear continuous bijective operator from $H_0^1(\Omega)$ to $H^{-1}(\Omega)$ and therefore also from $L^2(0, T; H_0^1(\Omega))$ to $L^2(0, T; H^{-1}(\Omega))$. We denote by $R$ its inverse which is linear continuous from $L^2(0, T; H^{-1}(\Omega))$ to $L^2(0, T; H_0^1(\Omega))$. More precisely, noting that

$$ \forall (u, v) \in H_0^1(\Omega)^2, \langle -\Delta u, v \rangle_{H^{-1}(\Omega), H_0^1(\Omega)} = \langle \nabla u, \nabla v \rangle_{L^2(\Omega)^d}, $$

and that $u = R\xi \iff \xi = -\Delta u$, the operator $R$ is characterized by

$$ \forall (\xi, v) \in H^{-1}(\Omega) \times H_0^1(\Omega), \quad \langle \nabla R \xi, \nabla v \rangle_{L^2(\Omega)^d} = \langle \xi, v \rangle_{H^{-1}(\Omega), H_0^1(\Omega)}. \qquad (A.5) $$

The problem (A.4) is then equivalent to

$$ \text{find } \bar{u} \in W(0, T; \Omega) \text{ s.t. } -\text{div}(\nabla R \bar{u}' + \nabla \bar{u} + \mathbf{F}) = f \text{ and } \bar{u}(0) - \Phi \bar{u}(T) = \xi_0, \qquad (A.6) $$
which in particular implies that $\nabla R\bar{u}' + \nabla \bar{u} + \mathbf{F} \in L^2(0, T; \mathbf{H}_{\text{div}}(\Omega))$. Let us first state a known result on the well-posedness of the problem (A.4).

**Theorem A.1 (Existence and uniqueness of $\bar{u}$).** For all $f \in L^2(0, T; L^2(\Omega))$, $\mathbf{F} \in L^2(0, T; L^2(\Omega)^d)$ and $\xi_0 \in L^2(\Omega)$, Problem (A.4) has a unique solution.

The proof of this theorem may be found in [23] in the case that $\Phi = 0$ and in [1] for any general contractive mapping $\Phi$.

A.1. **Description of the implicit Euler scheme.** Let $N \in \mathbb{N} \setminus \{0\}$ and define the time step (taken to be uniform for simplicity of presentation) $k = \frac{T}{N}$. For all $m \in [1, N]$, $f^{(m)} \in L^2(\Omega)$, $\mathbf{F}^{(m)} \in L^2(\Omega)^d$ are defined by

$$f^{(m)} = \frac{1}{k} \int_{(m-1)k}^{mk} f(t)dt \quad \text{and} \quad \mathbf{F}^{(m)} = \frac{1}{k} \int_{(m-1)k}^{mk} \mathbf{F}(t)dt.$$

The implicit Euler scheme consists in seeking $N+1$ elements of the space $X_{\mathcal{T}}$ given in Definition 2.1, denoted by $(u^{(m)})_{m=0,...,N}$, such that

$$\text{Initialization: } \langle u^{(0)} - \Phi u^{(N)}, v \rangle_{L^2(\Omega)} = \langle \xi_0, v \rangle_{L^2(\Omega)} \quad \text{for all } v \in X_{\mathcal{T}} \qquad (A.7a)$$

$$\text{Step } m: \quad \begin{aligned} & \langle \frac{u^{(m)} - u^{(m-1)}}{k}, v \rangle_{L^2(\Omega)} + d\langle G_{\mathcal{T}}u^{(m)}, G_{\mathcal{T}}v \rangle_{L^2(\Omega)} = \langle f^{(m)}, v \rangle_{L^2(\Omega)} \\ & \quad - \langle \mathbf{F}^{(m)}, \nabla_{\mathcal{T}}v \rangle_{L^2(\Omega)^d}, \quad \text{for all } m \in [1, N] \text{ and } v \in X_{\mathcal{T}}. \end{aligned} \qquad (A.7b)$$

Note that, if $\Phi = 0$, the scheme is the usual implicit scheme, and the existence and uniqueness of a solution $(u^{(0)}, (u^{(m)})_{m \in [1, N]})$ to (A.7) is standard. In the general case, a linear system involving $u^{(0)}$ must be solved, and its invertibility is proved by Theorem 3.1. Let us now introduce the discrete functional space $W_{\mathcal{T}}$, which consists of all functions $w: [0, T] \to X_{\mathcal{T}}$ in the sense that there exist $N+1$ elements of $X_{\mathcal{T}}$, denoted by $(w^{(m)})_{m=0,...,N}$, such that

$$\begin{aligned} w(0) &= w^{(0)}, \quad \text{and} \\ w(t) &= w^{(m)} \quad \text{for all } t \in ((m-1)k, mk], \quad \text{for all } m \in [1, N]. \end{aligned} \qquad (A.8)$$

Note that such a function $w \in W_{\mathcal{T}}$ is a function of time with values in a (discrete) functional space; however, observe that there is a bijection from the space $W_{\mathcal{T}}$ to the space $X_{\mathcal{T}}^{N+1}$, through the mapping $w \mapsto (w(mk))_{m \in [0, N]}$, so that $w$ may also be seen as a function of space and time, defined, with an abuse of notation, by $w: (x, t) \in \Omega \times [0, T] \mapsto w(x, t) = w^{(0)}(x)$ if $t=0$ and $w(x, t) = w^{(m)}(x)$ if $t \in ((m-1)k, mk]$. Therefore, to avoid any confusion, we denote the discrete time derivative of $w$ by $\partial_t w$, defined as follows:

$$\partial_t w(t) = \frac{w^{(m)} - w^{(m-1)}}{k}, \quad \text{for a.e. } t \in ((m-1)k, mk), \quad \text{for all } m \in [1, N], \qquad (A.9)$$

and its discrete gradient $G_{\mathcal{T}}v: [0, T] \to X_{\mathcal{T}}$ by

$$G_{\mathcal{T}}w(x, t) = \begin{cases} G_{\mathcal{T}}w^{(0)}(x) & \text{if } t=0, \\ G_{\mathcal{T}}w^{(m)}(x) & \text{if } t \in ((m-1)k, mk], \quad \text{for } m \in [1, N], \end{cases}$$

where $G_{\mathcal{T}}: X_{\mathcal{T}} \to L^2(\Omega)$ be the discrete normal gradient defined by (2.6). We also need to define the space $V_{\mathcal{T}}$ of all functions $v \in L^2(0, T; X_{\mathcal{T}})$ for which there exist $N$ elements of $X_{\mathcal{T}}$,
which in particular implies that $\nabla R\bar{u}' + \nabla \bar{u} + \mathbf{F} \in L^2(0, T; \mathbf{H}_{\text{div}}(\Omega))$. Let us first state a known result on the well-posedness of the problem (A.4).

**Theorem A.1 (Existence and uniqueness of $\bar{u}$).** For all $f \in L^2(0, T; L^2(\Omega))$, $\mathbf{F} \in L^2(0, T; L^2(\Omega)^d)$ and $\xi_0 \in L^2(\Omega)$, Problem (A.4) has a unique solution.

The proof of this theorem may be found in [23] in the case that $\Phi = 0$ and in [1] for any general contractive mapping $\Phi$.

**A.1. Description of the implicit Euler scheme.** Let $N \in \mathbb{N} \setminus \{0\}$ and define the time step (taken to be uniform for simplicity of presentation) $k = \frac{T}{N}$. For all $m \in [1, N]$, $f^{(m)} \in L^2(\Omega)$, $\mathbf{F}^{(m)} \in L^2(\Omega)^d$ are defined by

$$f^{(m)} = \frac{1}{k} \int_{(m-1)k}^{mk} f(t)dt \quad \text{and} \quad \mathbf{F}^{(m)} = \frac{1}{k} \int_{(m-1)k}^{mk} \mathbf{F}(t)dt.$$

The implicit Euler scheme consists in seeking $N+1$ elements of the space $X_{\mathcal{T}}$ given in Definition 2.1, denoted by $(u^{(m)})_{m=0,...,N}$, such that

$$\text{Initialization: } \langle u^{(0)} - \Phi u^{(N)}, v \rangle_{L^2(\Omega)} = \langle \xi_0, v \rangle_{L^2(\Omega)} \quad \text{for all } v \in X_{\mathcal{T}} \qquad (A.7a)$$

$$\text{Step } m: \quad \begin{aligned} \langle \frac{u^{(m)} - u^{(m-1)}}{k}, v \rangle_{L^2(\Omega)} + d\langle G_{\mathcal{T}}u^{(m)}, G_{\mathcal{T}}v \rangle_{L^2(\Omega)} = \langle f^{(m)}, v \rangle_{L^2(\Omega)} \\ - \langle \mathbf{F}^{(m)}, \nabla_{\mathcal{T}}v \rangle_{L^2(\Omega)^d}, \quad \text{for all } m \in [1, N] \text{ and } v \in X_{\mathcal{T}}. \end{aligned} \qquad (A.7b)$$

Note that, if $\Phi = 0$, the scheme is the usual implicit scheme, and the existence and uniqueness of a solution $(u^{(0)}, (u^{(m)})_{m \in [1, N]})$ to (A.7) is standard. In the general case, a linear system involving $u^{(0)}$ must be solved, and its invertibility is proved by Theorem 3.1. Let us now introduce the discrete functional space $W_{\mathcal{T}}$, which consists of all functions $w: [0, T] \to X_{\mathcal{T}}$ in the sense that there exist $N+1$ elements of $X_{\mathcal{T}}$, denoted by $(w^{(m)})_{m=0,...,N}$, such that

$$\begin{aligned} w(0) &= w^{(0)}, \quad \text{and} \\ w(t) &= w^{(m)} \quad \text{for all } t \in ((m-1)k, mk], \quad \text{for all } m \in [1, N]. \end{aligned} \qquad (A.8)$$

Note that such a function $w \in W_{\mathcal{T}}$ is a function of time with values in a (discrete) functional space; however, observe that there is a bijection from the space $W_{\mathcal{T}}$ to the space $X_{\mathcal{T}}^{N+1}$, through the mapping $w \mapsto (w(mk))_{m \in [0, N]}$, so that $w$ may also be seen as a function of space and time, defined, with an abuse of notation, by $w: (x, t) \in \Omega \times [0, T] \mapsto w(x, t) = w^{(0)}(x)$ if $t = 0$ and $w(x, t) = w^{(m)}(x)$ if $t \in ((m-1)k, mk]$. Therefore, to avoid any confusion, we denote the discrete time derivative of $w$ by $\partial_t w$, defined as follows:

$$\partial_t w(t) = \frac{w^{(m)} - w^{(m-1)}}{k}, \quad \text{for a.e. } t \in ((m-1)k, mk), \quad \text{for all } m \in [1, N], \qquad (A.9)$$

and its discrete gradient $G_{\mathcal{T}}v: [0, T] \to X_{\mathcal{T}}$ by

$$G_{\mathcal{T}}w(x, t) = \begin{cases} G_{\mathcal{T}}w^{(0)}(x) & \text{if } t = 0, \\ G_{\mathcal{T}}w^{(m)}(x) & \text{if } t \in ((m-1)k, mk], \quad \text{for } m \in [1, N], \end{cases}$$

where $G_{\mathcal{T}}: X_{\mathcal{T}} \to L^2(\Omega)$ be the discrete normal gradient defined by (2.6). We also need to define the space $V_{\mathcal{T}}$ of all functions $v \in L^2(0, T; X_{\mathcal{T}})$ for which there exist $N$ elements of $X_{\mathcal{T}}$,
denoted by $(v^{(m)})_{m \in [1, N]}$, such that

$$v(t) = v^{(m)} \text{ for a.e. } t \in ((m-1)k, mk), \text{ for all } m \in [1, N]. \quad (A.10)$$

Observe that there is now a bijection from the space $V_T$ to $X_T^N$, through the mapping $w \mapsto (w(mk))_{m \in [1, N]}$.

**Remark A.2** (Difference between $W_T$ and $V_T$). The functions belonging to $W_T$ are defined pointwise and everywhere, whereas those belonging to $V_T$ are only defined almost everywhere on $(0, T)$.

The scheme (A.7a)-(A.7b) can then be written under the following form (for the sake of brevity, here and in the sequel, we omit $0, T$ in the expression of some spaces of time dependent functions: for example, we write $L^2(L^2(\Omega))$ instead of $L^2(0, T; L^2(\Omega))$):

$$\text{Find } u \in W_T \text{ such that } \forall (v, z) \in V_T \times X_T, \quad b(u, (v, z)) = \mathcal{L}((v, z)), \text{ with} \quad (A.11a)$$

$$b(u, (v, z)) = \langle \partial_t u, v \rangle_{L^2(L^2(\Omega))} + d\langle G_T u, G_T v \rangle_{L^2(L^2(\Omega))} + \langle u(0) - \Phi u(T), z \rangle_{L^2(\Omega)}, \quad (A.11b)$$

$$\mathcal{L}((v, z)) = \langle f, v \rangle_{L^2(L^2(\Omega))} - \langle \mathbf{F}, \nabla_T v \rangle_{L^2(L^2(\Omega)^d)} + \langle \xi_0, z \rangle_{L^2(\Omega)}.$$

**A.2. An optimal error estimate.** Let $R_T: v \in X_T \mapsto R_T v \in X_T$, with $R_T v$ defined by

$$d\langle G_T R_T v, G_T w \rangle_{L^2(\Omega)^d} = \langle v, w \rangle_{L^2(\Omega)} \quad \text{for all } w \in X_T. \quad (A.12)$$

With this definition, the scheme (A.11) can be recast as: for all $(v, z) \in V_T \times X_T$,

$$d\langle G_T R_T \partial_t u + G_T u, G_T v \rangle_{L^2(L^2(\Omega)^d)} + \langle \mathbf{F}, \nabla_T v \rangle_{L^2(L^2(\Omega)^d)} \\ + \langle u(0) - \Phi u(T), z \rangle_{L^2(\Omega)} = \langle f, v \rangle_{L^2(L^2(\Omega))} + \langle \xi_0, z \rangle_{L^2(\Omega)}. \quad (A.13)$$

We define a distance between any functions $\varphi \in W(0, T; \Omega)$ and $v \in W_T$ by

$$\delta_T^{(T)}(\varphi, v) = d\|G_T R \varphi' - G_T R_T \partial_t v\|_{L^2(L^2(\Omega))} \\ + d\|G_T \varphi - G_T v\|_{L^2(L^2(\Omega))} + \max_{t \in [0,T]} \|\varphi(t) - v(t)\|_{L^2(\Omega)}. \quad (A.14)$$

(Recall that $G_T$ is the mean normal gradient defined by (3.6).) The interpolation error for a given function $\varphi \in W(0, T; \Omega)$ is then defined by

$$\mathcal{I}_T^{(T)}(\varphi) = \inf_{v \in W_T} \delta_T^{(T)}(\varphi, v) \quad (A.15)$$

We also define the time-space conformity error $\zeta_T^{(T)}: L^2(0, T; H_{\text{div}}(\Omega)) \to [0, +\infty)$ by

$$\forall \varphi \in L^2(0, T; H_{\text{div}}(\Omega)), \quad \zeta_T^{(T)}(\varphi) = \sup_{w \in V_T \setminus \{0\}} \frac{\left| \langle \varphi, \nabla_T w \rangle_{L^2(L^2(\Omega)^d)} + \langle \text{div}\varphi, w \rangle_{L^2(L^2(\Omega))} \right|}{\|G_T w\|_{L^2(L^2(\Omega)^d)}}. \quad (A.16)$$

Before proving the optimal error result, we need the following preliminary lemma.

**Lemma A.3.** For $w \in X_T$, the following inequalities hold

$$\max_{t \in [0, T]} \|w(t)\|_{L^2(\Omega)} \le \|dG_T R_T \partial_t w\|_{L^2(L^2(\Omega))} + \|dG_T w\|_{L^2(L^2(\Omega))} + \|w(0)\|_{L^2(\Omega)}, \quad (A.17)$$

$$\langle G_T R_T \partial_t w, dG_T w \rangle_{L^2(L^2(\Omega))} \ge \frac{1}{2} \|w(T)\|_{L^2(\Omega)}^2 - \frac{1}{2} \|w(0)\|_{L^2(\Omega)}^2, \quad (A.18)$$
denoted by $(v^{(m)})_{m \in \llbracket 1, N \rrbracket}$, such that

$$v(t) = v^{(m)} \text{ for a.e. } t \in ((m-1)k, mk), \text{ for all } m \in \llbracket 1, N \rrbracket. \quad (A.10)$$

Observe that there is now a bijection from the space $V_T$ to $X_T^N$, through the mapping $w \mapsto (w(mk))_{m \in \llbracket 1, N \rrbracket}$.

**Remark A.2** (Difference between $W_T$ and $V_T$). The functions belonging to $W_T$ are defined pointwise and everywhere, whereas those belonging to $V_T$ are only defined almost everywhere on $(0, T)$.

The scheme (A.7a)–(A.7b) can then be written under the following form (for the sake of brevity, here and in the sequel, we omit 0, T in the expression of some spaces of time dependent functions: for example, we write $L^2(L^2(\Omega))$ instead of $L^2(0, T; L^2(\Omega))$):

$$\text{Find } u \in W_T \text{ such that } \forall (v, z) \in V_T \times X_T, \quad b(u, (v, z)) = \mathcal{L}((v, z)), \text{ with} \quad (A.11a)$$

$$b(u, (v, z)) = \langle \partial_t u, v \rangle_{L^2(L^2(\Omega))} + d\langle G_T u, G_T v \rangle_{L^2(L^2(\Omega))} + \langle u(0) - \Phi u(T), z \rangle_{L^2(\Omega)}, \quad (A.11b)$$

$$\mathcal{L}((v, z)) = \langle f, v \rangle_{L^2(L^2(\Omega))} - \langle \mathbf{F}, \nabla_T v \rangle_{L^2(L^2(\Omega)^d)} + \langle \xi_0, z \rangle_{L^2(\Omega)}.$$

**A.2. An optimal error estimate.** Let $R_T: v \in X_T \mapsto R_T v \in X_T$, with $R_T v$ defined by

$$d\langle G_T R_T v, G_T w \rangle_{L^2(\Omega)^d} = \langle v, w \rangle_{L^2(\Omega)} \quad \text{for all } w \in X_T. \quad (A.12)$$

With this definition, the scheme (A.11) can be recast as: for all $(v, z) \in V_T \times X_T$,

$$d\langle G_T R_T \partial_t u + G_T u, G_T v \rangle_{L^2(L^2(\Omega)^d)} + \langle \mathbf{F}, \nabla_T v \rangle_{L^2(L^2(\Omega)^d)} \\ + \langle u(0) - \Phi u(T), z \rangle_{L^2(\Omega)} = \langle f, v \rangle_{L^2(L^2(\Omega))} + \langle \xi_0, z \rangle_{L^2(\Omega)}. \quad (A.13)$$

We define a distance between any functions $\varphi \in W(0, T; \Omega)$ and $v \in W_T$ by

$$\delta_T^{(T)}(\varphi, v) = d\|G_T R \varphi' - G_T R_T \partial_t v\|_{L^2(L^2(\Omega))} \\ + d\|G_T \varphi - G_T v\|_{L^2(L^2(\Omega))} + \max_{t \in [0,T]} \|\varphi(t) - v(t)\|_{L^2(\Omega)}. \quad (A.14)$$

(Recall that $G_T$ is the mean normal gradient defined by (3.6).) The interpolation error for a given function $\varphi \in W(0, T; \Omega)$ is then defined by

$$\mathcal{I}_T^{(T)}(\varphi) = \inf_{v \in W_T} \delta_T^{(T)}(\varphi, v) \quad (A.15)$$

We also define the time-space conformity error $\zeta_T^{(T)}: L^2(0, T; H_{\text{div}}(\Omega)) \to [0, +\infty)$ by

$$\forall \varphi \in L^2(0, T; H_{\text{div}}(\Omega)), \quad \zeta_T^{(T)}(\varphi) = \sup_{w \in V_T \setminus \{0\}} \frac{\left| \langle \varphi, \nabla_T w \rangle_{L^2(L^2(\Omega)^d)} + \langle \text{div} \varphi, w \rangle_{L^2(L^2(\Omega))} \right|}{\|G_T w\|_{L^2(L^2(\Omega)^d)}}. \quad (A.16)$$

Before proving the optimal error result, we need the following preliminary lemma.

**Lemma A.3.** For $w \in X_T$, the following inequalities hold

$$\max_{t \in [0, T]} \|w(t)\|_{L^2(\Omega)} \le \|dG_T R_T \partial_t w\|_{L^2(L^2(\Omega))} + \|dG_T w\|_{L^2(L^2(\Omega))} + \|w(0)\|_{L^2(\Omega)}, \quad (A.17)$$

$$\langle G_T R_T \partial_t w, dG_T w \rangle_{L^2(L^2(\Omega))} \ge \frac{1}{2} \|w(T)\|_{L^2(\Omega)}^2 - \frac{1}{2} \|w(0)\|_{L^2(\Omega)}^2, \quad (A.18)$$
and, recalling that $\mathrm{diam}(\Omega)$ is defined by (3.4),

$$ \|d G_{\mathcal{T}} R_{\mathcal{T}} \partial_t w\|_{L^2(L^2(\Omega))}^2 + \left(1 + \frac{\mathrm{diam}(\Omega)^2}{T}\right) \|d G_{\mathcal{T}} w\|_{L^2(L^2(\Omega))}^2 \geq \|w(T)\|_{L^2(\Omega)}^2. \quad (A.19) $$

*Proof.* Let $w \in X_{\mathcal{T}}$; thanks to the equality $(a-b)a = \frac{1}{2}a^2 + \frac{1}{2}(a-b)^2 - \frac{1}{2}b^2$, the definition (A.12) of $R_{\mathcal{T}}$ yields, for $0 \le m \le n \le N$,

$$ \begin{aligned} \int_{mk}^{nk} \langle G_{\mathcal{T}} R_{\mathcal{T}} \partial_t w(t), dG_{\mathcal{T}} w(t) \rangle_{L^2(\Omega)} dt &= \int_{mk}^{nk} \langle \partial_t w(t), w(t) \rangle_{L^2(\Omega)} dt \\ &= \sum_{p=m}^{n-1} k \langle \frac{w^{(p+1)} - w^{(p)}}{k}, w^{(p+1)} \rangle_{L^2(\Omega)} \\ &= \frac{1}{2} \|w^{(n)}\|_{L^2(\Omega)}^2 + \frac{1}{2} \sum_{p=m}^{n-1} \|(w^{(p+1)} - w^{(p)})\|_{L^2(\Omega)}^2 - \frac{1}{2} \|w^{(m)}\|_{L^2(\Omega)}^2. \end{aligned} \quad (A.20) $$

Applying the Cauchy–Schwarz inequality to the left-hand side, we get

$$ \begin{aligned} \frac{1}{2} \|w^{(n)}\|_{L^2(\Omega)}^2 &\le \|d G_{\mathcal{T}} R_{\mathcal{T}} \partial_t w\|_{L^2(L^2(\Omega)^d)} \|d G_{\mathcal{T}} w\|_{L^2(L^2(\Omega)^d)} + \frac{1}{2} \|w^{(m)}\|_{L^2(\Omega)}^2 \\ &\le \frac{1}{2} \|d G_{\mathcal{T}} R_{\mathcal{T}} \partial_t w\|_{L^2(L^2(\Omega)^d)}^2 + \frac{1}{2} \|d G_{\mathcal{T}} w\|_{L^2(L^2(\Omega)^d)}^2 + \frac{1}{2} \|w^{(m)}\|_{L^2(\Omega)}^2, \end{aligned} \quad (A.21) $$

where the second line follows from the Young inequality. Setting $m = 0$ allows us to take any $n = 0, \dots, N$. Taking the square root of the above inequality and using $(a^2+b^2+c^2)^{1/2} \le a+b+c$ then concludes the proof of (A.17).

The inequality (A.18) is obtained letting $m = 0$ and $n = N$ in (A.20). To prove (A.19), we come back to (A.21) and set $n = N$ to get, after multiplication by $2k$, for all $m = 0, \dots, N$,

$$ k \|w(T)\|_{L^2(\Omega)}^2 \le k \|d G_{\mathcal{T}} R_{\mathcal{T}} \partial_t w\|_{L^2(L^2(\Omega)^d)}^2 + k \|d G_{\mathcal{T}} w\|_{L^2(L^2(\Omega)^d)}^2 + k \|w^{(m)}\|_{L^2(\Omega)}^2. $$

Summing over $m \in [1, N]$ yields

$$ \begin{aligned} T \|w(T)\|_{L^2(\Omega)}^2 &\le T \|d G_{\mathcal{T}} R_{\mathcal{T}} \partial_t w\|_{L^2(L^2(\Omega)^d)}^2 + T \|d G_{\mathcal{T}} w\|_{L^2(L^2(\Omega)^d)}^2 + \|w\|_{L^2(L^2(\Omega))}^2 \\ &\le T \|d G_{\mathcal{T}} R_{\mathcal{T}} \partial_t w\|_{L^2(L^2(\Omega)^d)}^2 + (T + \mathrm{diam}(\Omega)^2) \|d G_{\mathcal{T}} w\|_{L^2(L^2(\Omega)^d)}^2, \end{aligned} $$

which proves (A.19).

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

We now turn to the error estimate result.

**Theorem A.4.** There exists one and only one solution $u$ to (A.7). Moreover, letting $\bar{u}$ be the solution to (A.4) and $v$ be defined by

$$ v := \nabla R \bar{u}' + \nabla \bar{u} + \mathbf{F} \in L^2(0, T; \mathbf{H}_{\mathrm{div}}(\Omega)), \quad (A.22) $$

there exists $C_{\Omega,T} \ge 0$, depending only on $\Omega$ and $T$, such that:

$$ \frac{1}{2} \left[ \zeta_{\mathcal{T}}^{(T)}(\mathbf{v}) + \mathcal{I}_{\mathcal{T}}^{(T)}(\bar{u}, v) \right] \le \delta_{\mathcal{T}}^{(T)}(\bar{u}, u) \le C_{\Omega,T} \left[ \zeta_{\mathcal{T}}^{(T)}(\mathbf{v}) + \mathcal{I}_{\mathcal{T}}^{(T)}(\bar{u}, v) \right]. \quad (A.23) $$

**Remark A.5** (Optimality of the error estimate (A.23)). As in the steady case, the second inequality in (A.23) gives an error estimate for the scheme, while the first inequality shows its optimality.
and, recalling that $\mathrm{diam}(\Omega)$ is defined by (3.4),

$$ \|dG_{\mathcal{T}}R_{\mathcal{T}}\partial_t w\|_{L^2(L^2(\Omega))}^2 + \left(1 + \frac{\mathrm{diam}(\Omega)^2}{T}\right) \|dG_{\mathcal{T}}w\|_{L^2(L^2(\Omega))}^2 \ge \|w(T)\|_{L^2(\Omega)}^2. \quad (A.19) $$

*Proof.* Let $w \in X_{\mathcal{T}}$; thanks to the equality $(a-b)a = \frac{1}{2}a^2 + \frac{1}{2}(a-b)^2 - \frac{1}{2}b^2$, the definition (A.12) of $R_{\mathcal{T}}$ yields, for $0 \le m \le n \le N$,

$$ \begin{aligned} \int_{mk}^{nk} \langle G_{\mathcal{T}} R_{\mathcal{T}} \partial_t w(t), dG_{\mathcal{T}} w(t) \rangle_{L^2(\Omega)} dt &= \int_{mk}^{nk} \langle \partial_t w(t), w(t) \rangle_{L^2(\Omega)} dt \\ &= \sum_{p=m}^{n-1} k \langle \frac{w^{(p+1)} - w^{(p)}}{k}, w^{(p+1)} \rangle_{L^2(\Omega)} \\ &= \frac{1}{2} \|w^{(n)}\|_{L^2(\Omega)}^2 + \frac{1}{2} \sum_{p=m}^{n-1} \|(w^{(p+1)} - w^{(p)})\|_{L^2(\Omega)}^2 - \frac{1}{2} \|w^{(m)}\|_{L^2(\Omega)}^2. \end{aligned} \quad (A.20) $$

Applying the Cauchy–Schwarz inequality to the left-hand side, we get

$$ \begin{aligned} \frac{1}{2} \|w^{(n)}\|_{L^2(\Omega)}^2 &\le \|dG_{\mathcal{T}}R_{\mathcal{T}}\partial_t w\|_{L^2(L^2(\Omega)^d)} \|dG_{\mathcal{T}}w\|_{L^2(L^2(\Omega)^d)} + \frac{1}{2} \|w^{(m)}\|_{L^2(\Omega)}^2 \\ &\le \frac{1}{2} \|dG_{\mathcal{T}}R_{\mathcal{T}}\partial_t w\|_{L^2(L^2(\Omega)^d)}^2 + \frac{1}{2} \|dG_{\mathcal{T}}w\|_{L^2(L^2(\Omega)^d)}^2 + \frac{1}{2} \|w^{(m)}\|_{L^2(\Omega)}^2, \end{aligned} \quad (A.21) $$

where the second line follows from the Young inequality. Setting $m = 0$ allows us to take any $n = 0, \dots, N$. Taking the square root of the above inequality and using $(a^2+b^2+c^2)^{1/2} \le a+b+c$ then concludes the proof of (A.17).

The inequality (A.18) is obtained letting $m = 0$ and $n = N$ in (A.20). To prove (A.19), we come back to (A.21) and set $n = N$ to get, after multiplication by $2k$, for all $m = 0, \dots, N$,

$$ k \|w(T)\|_{L^2(\Omega)}^2 \le k \|dG_{\mathcal{T}}R_{\mathcal{T}}\partial_t w\|_{L^2(L^2(\Omega)^d)}^2 + k \|dG_{\mathcal{T}}w\|_{L^2(L^2(\Omega)^d)}^2 + k \|w^{(m)}\|_{L^2(\Omega)}^2. $$

Summing over $m \in [1, N]$ yields

$$ \begin{aligned} T \|w(T)\|_{L^2(\Omega)}^2 &\le T \|dG_{\mathcal{T}}R_{\mathcal{T}}\partial_t w\|_{L^2(L^2(\Omega)^d)}^2 + T \|dG_{\mathcal{T}}w\|_{L^2(L^2(\Omega)^d)}^2 + \|w\|_{L^2(L^2(\Omega))}^2 \\ &\le T \|dG_{\mathcal{T}}R_{\mathcal{T}}\partial_t w\|_{L^2(L^2(\Omega)^d)}^2 + (T + \mathrm{diam}(\Omega)^2) \|dG_{\mathcal{T}}w\|_{L^2(L^2(\Omega)^d)}^2, \end{aligned} $$

which proves (A.19).

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

We now turn to the error estimate result.

**Theorem A.4.** There exists one and only one solution $u$ to (A.7). Moreover, letting $\bar{u}$ be the solution to (A.4) and $v$ be defined by

$$ v := \nabla R \bar{u}' + \nabla \bar{u} + \mathbf{F} \in L^2(0, T; \mathbf{H}_{\mathrm{div}}(\Omega)), \quad (A.22) $$

there exists $C_{\Omega,T} \ge 0$, depending only on $\Omega$ and $T$, such that:

$$ \frac{1}{2} \left[ \zeta_{\mathcal{T}}^{(T)}(\mathbf{v}) + \mathcal{I}_{\mathcal{T}}^{(T)}(\bar{u}, v) \right] \le \delta_{\mathcal{T}}^{(T)}(\bar{u}, u) \le C_{\Omega,T} \left[ \zeta_{\mathcal{T}}^{(T)}(\mathbf{v}) + \mathcal{I}_{\mathcal{T}}^{(T)}(\bar{u}, v) \right]. \quad (A.23) $$

**Remark A.5** (Optimality of the error estimate (A.23)). As in the steady case, the second inequality in (A.23) gives an error estimate for the scheme, while the first inequality shows its optimality.
Proof of Theorem A.4. Let $v \in V_T$ and $z \in X_T$ be given. Definition (A.16) of $\zeta_T^{(T)}(v)$ gives

$$ \int_0^T \left( \langle \mathbf{v}(t), \nabla_T v(t) \rangle_{L^2(\Omega)^d} + \langle \operatorname{div} \mathbf{v}(t), v(t) \rangle_{L^2(\Omega)} \right) dt \le \zeta_T^{(T)}(\mathbf{v}) \sqrt{d} \|G_T v\|_{L^2(L^2(\Omega))}. $$

This yields, using (A.6) (which reads $\operatorname{div} \mathbf{v} = -f$), and (A.13),

$$ \int_0^T d \left( \langle G_T R \bar{u}'(t) + G_T \bar{u}(t) - (G_T R \partial_t u(t) + G_T u(t)), G_T v(t) \rangle_{L^2(\Omega)} \right) dt \\ + \langle \xi_0 - (u(0) - \Phi u(T)), z \rangle_{L^2(\Omega)} \le \zeta_T^{(T)}(\mathbf{v}) \sqrt{d} \|G_T v\|_{L^2(L^2(\Omega))}. $$

Using $\bar{u}(0) - \Phi \bar{u}(T) = \xi_0$, we get

$$ \int_0^T d \left( \langle G_T R \bar{u}'(t) + G_T \bar{u}(t) - (G_T R \partial_t u(t) + G_T u(t)), G_T v(t) \rangle_{L^2(\Omega)} \right) dt \\ + \langle \bar{u}(0) - \Phi \bar{u}(T) - (u(0) - \Phi u(T)), z \rangle_{L^2(\Omega)} \\ \le \zeta_T^{(T)}(\mathbf{v}) \sqrt{d} \|G_T v\|_{L^2(L^2(\Omega))}. \quad (A.24) $$

We then take an arbitrary element $\tilde{v} \in W_T$ and notice that, by definition (A.14) of $\delta_T^{(T)}(\bar{u}, \tilde{v})$ and since $\Phi$ is a contraction,

$$ \int_0^T d \left[ \langle [G_T R_T \partial \tilde{v} - G_T R \bar{u}'](t), G_T v(t) \rangle_{L^2(\Omega)} + \langle [G_T \tilde{v} - G_T \bar{u}](t), G_T v(t) \rangle_{L^2(\Omega)} dt \right. \\ \left. + \langle \tilde{v}(0) - \bar{u}(0) - \Phi(\tilde{v}(T) - \bar{u}(T)), z \rangle_{L^2(\Omega)} \right. \\ \left. \le \delta_T^{(T)}(\bar{u}, \tilde{v}) \|G_T v\|_{L^2(L^2(\Omega))} + 2\delta_T^{(T)}(\bar{u}, \tilde{v}) \|z\|_{L^2(\Omega)} \right]. $$

Adding this inequality to (A.24) yields

$$ \int_0^T d \langle G_T R_T \partial(\tilde{v} - u)(t) + G_T(\tilde{v}(t) - u(t)), G_T v(t) \rangle_{L^2(\Omega)} dt \\ \langle (\tilde{v}(0) - u(0)) - \Phi(\tilde{v}(T) - u(T)), z \rangle_{L^2(\Omega)} \\ \le (\delta_T^{(T)}(\bar{u}, \tilde{v}) + \zeta_T^{(T)}(\mathbf{v})) \|G_T v\|_{L^2(L^2(\Omega)^d)} + 2\delta_T^{(T)}(\bar{u}, \tilde{v}) \|z\|_{L^2(\Omega)}. \quad (A.25) $$

Let us now introduce the notations used in Lemma A.6, and prove that the hypotheses of the lemma hold. We denote by $V = L^2(0, T; L^2(\Omega))$ the Hilbert space endowed with the inner product $d\langle \cdot, \cdot \rangle_{L^2(L^2(\Omega))}$, and by $L = L^2(\Omega)$. Let $Z$ and $Y$ be the Hilbert spaces defined by $Z = V \times V \times L \times L$ and $Y = V \times L$. We define the bilinear form $\hat{b}: Z \times Y \to \mathbb{R}$ by

$$ \hat{b}((z_1, z_2, z_3, z_4), (y_1, y_2)) = d\langle z_1 + Az_2, y_1 \rangle_V + \langle z_3 - \Phi z_4, y_2 \rangle_L, \quad (A.26) $$

where $A = \operatorname{Id}$. We then define the subspace $X \subset Z$ by

$$ X = G_T(V_T) \times G_T(V_T) \times X_T \times X_T. $$

The operator $A = \operatorname{Id}$ satisfies (A.30) with $\alpha = M = 1$. Condition (A.32) is satisfied since $X_1 = X_2 = G_T(V_T)$. Let us prove that we can find $\zeta > 0$ and $\delta > 0$ such that (A.33) holds.
Proof of Theorem A.4. Let $v \in V_T$ and $z \in X_T$ be given. Definition (A.16) of $\zeta_T^{(T)}(v)$ gives

$$ \int_0^T \left( \langle \mathbf{v}(t), \nabla_T v(t) \rangle_{L^2(\Omega)^d} + \langle \operatorname{div} \mathbf{v}(t), v(t) \rangle_{L^2(\Omega)} \right) dt \le \zeta_T^{(T)}(\mathbf{v}) \sqrt{d} \|G_T v\|_{L^2(L^2(\Omega))}. $$

This yields, using (A.6) (which reads $\operatorname{div} \mathbf{v} = -f$), and (A.13),

$$ \int_0^T d\left( \langle G_T R \bar{u}'(t) + G_T \bar{u}(t) - (G_T R \partial_t u(t) + G_T u(t)), G_T v(t) \rangle_{L^2(\Omega)} \right) dt \\ + \langle \xi_0 - (u(0) - \Phi u(T)), z \rangle_{L^2(\Omega)} \le \zeta_T^{(T)}(\mathbf{v}) \sqrt{d} \|G_T v\|_{L^2(L^2(\Omega))}. $$

Using $\bar{u}(0) - \Phi \bar{u}(T) = \xi_0$, we get

$$ \int_0^T d\left( \langle G_T R \bar{u}'(t) + G_T \bar{u}(t) - (G_T R \partial_t u(t) + G_T u(t)), G_T v(t) \rangle_{L^2(\Omega)} \right) dt \\ + \langle \bar{u}(0) - \Phi \bar{u}(T) - (u(0) - \Phi u(T)), z \rangle_{L^2(\Omega)} \\ \le \zeta_T^{(T)}(\mathbf{v}) \sqrt{d} \|G_T v\|_{L^2(L^2(\Omega))}. \quad (A.24) $$

We then take an arbitrary element $\tilde{v} \in W_T$ and notice that, by definition (A.14) of $\delta_T^{(T)}(\bar{u}, \tilde{v})$ and since $\Phi$ is a contraction,

$$ \int_0^T d\left[ \langle G_T R_T \partial \tilde{v} - G_T R \bar{u}'(t), G_T v(t) \rangle_{L^2(\Omega)} + \langle [G_T \tilde{v} - G_T \bar{u}](t), G_T v(t) \rangle_{L^2(\Omega)} dt \right. \\ \left. + \langle \tilde{v}(0) - \bar{u}(0) - \Phi(\tilde{v}(T) - \bar{u}(T)), z \rangle_{L^2(\Omega)} \right. \\ \left. \le \delta_T^{(T)}(\bar{u}, \tilde{v}) \|G_T v\|_{L^2(L^2(\Omega))} + 2\delta_T^{(T)}(\bar{u}, \tilde{v}) \|z\|_{L^2(\Omega)} \right]. $$

Adding this inequality to (A.24) yields

$$ \int_0^T d\langle G_T R_T \partial(\tilde{v} - u)(t) + G_T(\tilde{v}(t) - u(t)), G_T v(t) \rangle_{L^2(\Omega)} dt \\ \langle (\tilde{v}(0) - u(0)) - \Phi(\tilde{v}(T) - u(T)), z \rangle_{L^2(\Omega)} \\ \le (\delta_T^{(T)}(\bar{u}, \tilde{v}) + \zeta_T^{(T)}(\mathbf{v})) \|G_T v\|_{L^2(L^2(\Omega)^d)} + 2\delta_T^{(T)}(\bar{u}, \tilde{v}) \|z\|_{L^2(\Omega)}. \quad (A.25) $$

Let us now introduce the notations used in Lemma A.6, and prove that the hypotheses of the lemma hold. We denote by $V = L^2(0, T; L^2(\Omega))$ the Hilbert space endowed with the inner product $d\langle \cdot, \cdot \rangle_{L^2(L^2(\Omega))}$, and by $L = L^2(\Omega)$. Let $Z$ and $Y$ be the Hilbert spaces defined by $Z = V \times V \times L \times L$ and $Y = V \times L$. We define the bilinear form $\hat{b}: Z \times Y \to \mathbb{R}$ by

$$ \hat{b}((z_1, z_2, z_3, z_4), (y_1, y_2)) = d\langle z_1 + Az_2, y_1 \rangle_V + \langle z_3 - \Phi z_4, y_2 \rangle_L, \quad (A.26) $$

where $A = \operatorname{Id}$. We then define the subspace $X \subset Z$ by

$$ X = G_T(V_T) \times G_T(V_T) \times X_T \times X_T. $$

The operator $A = \operatorname{Id}$ satisfies (A.30) with $\alpha = M = 1$. Condition (A.32) is satisfied since $X_1 = X_2 = G_T(V_T)$. Let us prove that we can find $\zeta > 0$ and $\delta > 0$ such that (A.33) holds.
Applying Lemma A.3, we add (A.18) to $(1 + \frac{\mathrm{diam}(\Omega)^2}{T})^{-1} \frac{1}{12} \times (\mathrm{A}.19)$. This shows that the hypothesis (A.33) is satisfied with

$$ \zeta := \mu = \frac{1}{2} + \left( 1 + \frac{\mathrm{diam}(\Omega)^2}{T} \right)^{-1} \frac{1}{12} \quad \text{and} \quad \nu = \frac{1}{2}. $$

We note that $\mu - \nu \|\Phi\|^2 \ge \mu - \nu = \left(1 + \frac{\mathrm{diam}(\Omega)^2}{T}\right)^{-1} \frac{1}{12} =: \delta$.

We can therefore apply Lemma A.6, which yields the existence of $\hat{\beta} > 0$ depending only on $\mathrm{diam}(\Omega)$ and $T$ such that

$$ \sup_{(y_1,y_2) \in X_2 \times X_3, \|(y_1,y_2)\|_Y = 1} \hat{b}((z_1, z_2, z_3, z_4), (y_1, y_2)) \ge \hat{\beta} \|(z_1, z_2, z_3, z_4)\|_Z \\ \forall (z_1, z_2, z_3, z_4) \in X. \quad (\mathrm{A}.27) $$

We now remark that Inequality (A.25) can be expressed by

$$ \hat{b}((z_1, z_2, z_3, z_4), (y_1, y_2)) \le \hat{c}_1 \|y_1\|_V + \hat{c}_2 \|y_2\|_{L^2(\Omega)}, \quad (\mathrm{A}.28) $$

with

$$ \begin{align*} z_1 &= G_{\mathcal{T}} R_{\mathcal{T}} \partial(\tilde{v} - u), & z_2 &= G_{\mathcal{T}}(\tilde{v} - u), \\ z_3 &= (\tilde{v}(0) - u(0)), & z_4 &= (\tilde{v}(T) - u(T)), \\ y_1 &= G_{\mathcal{T}} v, & y_2 &= z, \\ \hat{c}_1 &= \delta_{\mathcal{T}}^{(T)}(\bar{u}, \tilde{v}) + \zeta_{\mathcal{T}}^{(T)}(\mathbf{v}), & \hat{c}_2 &= 2\delta_{\mathcal{T}}^{(T)}(\bar{u}, \tilde{v}). \end{align*} $$

Taking in the right-hand-side of (A.28) the maximum over all $(y_1, y_2) \in G_{\mathcal{T}}(V_{\mathcal{T}}) \times X_{\mathcal{T}}$ with norm in $V \times L^2(\Omega)$ equal to 1 and using (A.27), we deduce that

$$ \begin{align*} \hat{\beta} \left( \|G_{\mathcal{T}} R_{\mathcal{T}} \partial(\tilde{v} - u)\|_{L^2(L^2(\Omega)^d)}^2 + \|G_{\mathcal{T}}(\tilde{v} - u)\|_{L^2(L^2(\Omega)^d)}^2 \right. \\ \left. + \|(\tilde{v} - u)(0)\|_{L^2(\Omega)}^2 + \|(\tilde{v} - u)(T)\|_{L^2(\Omega)}^2 \right)^{1/2} \\ \le \left[ \left( \delta_{\mathcal{T}}^{(T)}(\bar{u}, \tilde{v}) + \zeta_{\mathcal{T}}^{(T)}(\mathbf{v}) \right)^2 + 4\delta_{\mathcal{T}}^{(T)}(\bar{u}, \tilde{v})^2 \right]^{1/2} \\ \le \delta_{\mathcal{T}}^{(T)}(\bar{u}, \tilde{v}) + \zeta_{\mathcal{T}}^{(T)}(\mathbf{v}) + 2\delta_{\mathcal{T}}^{(T)}(\bar{u}, \tilde{v}), \quad (\mathrm{A}.29) \end{align*} $$

where we use $(a^2 + b^2)^{1/2} \le a + b$ for positive $a$ and $b$ in the last inequality. Plugging this into (A.29) and using (A.17) in Lemma A.3 together with $a + b + c \le \sqrt{3}(a^2 + b^2 + c^2)^{1/2}$, we infer

$$ \hat{\beta} \left( \|G_{\mathcal{T}} R_{\mathcal{T}} \partial(\tilde{v} - u)\|_{L^2(L^2(\Omega)^d)} + \|G_{\mathcal{T}}(\tilde{v} - u)\|_{L^2(L^2(\Omega)^d)} + \max_{t \in [0,T]} \|(\tilde{v} - u)(t)\|_{L^2(\Omega)} \right) \\ \le \sqrt{3} \left( 3\delta_{\mathcal{T}}^{(T)}(\bar{u}, \tilde{v}) + \zeta_{\mathcal{T}}^{(T)}(\mathbf{v}) \right). $$

Using the triangle inequality in the definition (A.14) of $\delta_{\mathcal{T}}^{(T)}$, we infer

$$ \hat{\beta} \delta_{\mathcal{T}}^{(T)}(\bar{u}, u) \le \sqrt{3} \left( 3 \delta_{\mathcal{T}}^{(T)}(\bar{u}, \tilde{v}) + \zeta_{\mathcal{T}}^{(T)}(\mathbf{v}) \right) + \hat{\beta} \delta_{\mathcal{T}}^{(T)}(\bar{u}, \tilde{v}). $$

Since $\tilde{v}$ is arbitrary in $W_{\mathcal{T}}$, this concludes the proof of the second inequality in (A.23).
Applying Lemma A.3, we add (A.18) to $(1 + \frac{\mathrm{diam}(\Omega)^2}{T})^{-1} \frac{1}{12} \times (\mathrm{A}.19)$. This shows that the hypothesis (A.33) is satisfied with

$$ \zeta := \mu = \frac{1}{2} + \left( 1 + \frac{\mathrm{diam}(\Omega)^2}{T} \right)^{-1} \frac{1}{12} \quad \text{and} \quad \nu = \frac{1}{2}. $$

We note that $\mu - \nu \|\Phi\|^2 \ge \mu - \nu = \left(1 + \frac{\mathrm{diam}(\Omega)^2}{T}\right)^{-1} \frac{1}{12} =: \delta$.

We can therefore apply Lemma A.6, which yields the existence of $\hat{\beta} > 0$ depending only on $\mathrm{diam}(\Omega)$ and $T$ such that

$$ \sup_{(y_1,y_2) \in X_2 \times X_3, \|(y_1,y_2)\|_Y = 1} \hat{b}((z_1, z_2, z_3, z_4), (y_1, y_2)) \ge \hat{\beta} \|(z_1, z_2, z_3, z_4)\|_Z \\ \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \qquad \forall (z_1, z_2, z_3, z_4) \in X. \quad (\mathrm{A}.27) $$

We now remark that Inequality (A.25) can be expressed by

$$ \hat{b}((z_1, z_2, z_3, z_4), (y_1, y_2)) \le \hat{c}_1 \|y_1\|_V + \hat{c}_2 \|y_2\|_{L^2(\Omega)}, \quad (\mathrm{A}.28) $$

with

$$ \begin{aligned} z_1 &= G_T R_T \partial(\tilde{v} - u), & z_2 &= G_T(\tilde{v} - u), \\ z_3 &= (\tilde{v}(0) - u(0)), & z_4 &= (\tilde{v}(T) - u(T)), \\ y_1 &= G_T v, & y_2 &= z, \\ \hat{c}_1 &= \delta_T^{(T)}(\bar{u}, \tilde{v}) + \zeta_T^{(T)}(\mathbf{v}), & \hat{c}_2 &= 2\delta_T^{(T)}(\bar{u}, \tilde{v}). \end{aligned} $$

Taking in the right-hand-side of (A.28) the maximum over all $(y_1, y_2) \in G_T(V_T) \times X_T$ with norm in $V \times L^2(\Omega)$ equal to 1 and using (A.27), we deduce that

$$ \begin{aligned} & \hat{\beta} \left( \|G_T R_T \partial(\tilde{v} - u)\|_{L^2(L^2(\Omega)^d)}^2 + \|G_T(\tilde{v} - u)\|_{L^2(L^2(\Omega)^d)}^2 \right. \\ & \qquad \left. + \|(\tilde{v} - u)(0)\|_{L^2(\Omega)}^2 + \|(\tilde{v} - u)(T)\|_{L^2(\Omega)}^2 \right)^{1/2} \\ & \le \left[ \left( \delta_T^{(T)}(\bar{u}, \tilde{v}) + \zeta_T^{(T)}(\mathbf{v}) \right)^2 + 4\delta_T^{(T)}(\bar{u}, \tilde{v})^2 \right]^{1/2} \\ & \le \delta_T^{(T)}(\bar{u}, \tilde{v}) + \zeta_T^{(T)}(\mathbf{v}) + 2\delta_T^{(T)}(\bar{u}, \tilde{v}), \quad (\mathrm{A}.29) \end{aligned} $$

where we use $(a^2 + b^2)^{1/2} \le a + b$ for positive $a$ and $b$ in the last inequality. Plugging this into (A.29) and using (A.17) in Lemma A.3 together with $a + b + c \le \sqrt{3}(a^2 + b^2 + c^2)^{1/2}$, we infer

$$ \begin{aligned} & \hat{\beta} \left( \|G_T R_T \partial(\tilde{v} - u)\|_{L^2(L^2(\Omega)^d)} + \|G_T(\tilde{v} - u)\|_{L^2(L^2(\Omega)^d)} + \max_{t \in [0,T]} \|(\tilde{v} - u)(t)\|_{L^2(\Omega)} \right) \\ & \le \sqrt{3} \left( 3\delta_T^{(T)}(\bar{u}, \tilde{v}) + \zeta_T^{(T)}(\mathbf{v}) \right). \end{aligned} $$

Using the triangle inequality in the definition (A.14) of $\delta_T^{(T)}$, we infer

$$ \hat{\beta} \delta_T^{(T)}(\bar{u}, u) \le \sqrt{3} \left( 3 \delta_T^{(T)}(\bar{u}, \tilde{v}) + \zeta_T^{(T)}(\mathbf{v}) \right) + \hat{\beta} \delta_T^{(T)}(\bar{u}, \tilde{v}). $$

Since $\tilde{v}$ is arbitrary in $W_T$, this concludes the proof of the second inequality in (A.23).
Let us now turn to the first inequality in (A.23); first note that

$$ \inf_{v \in W_{\tau}} \delta_{\tau}^{(T)}(\bar{u}, v) \leq \delta_{\tau}^{(T)}(\bar{u}, u). $$

To bound $\zeta_{\tau}^{(T)}(\mathbf{v})$ we recall that $\mathbf{v} = \nabla R \bar{u}' + \nabla \bar{u} + \mathbf{F}$ satisfies $-\mathrm{div} \mathbf{v} = f$ (see (A.6)), and use the scheme (A.13) (with $z=0$) to write, for any $v \in V_{\tau} \setminus \{0\}$,

$$ \begin{aligned} & \langle \mathbf{v}, G_{\tau} v \rangle_{L^2(\Omega)^d} + \langle \mathrm{div} \mathbf{v}, v \rangle_{L^2(\Omega)} \\ &= \int_0^T \left( \langle \nabla R \bar{u}'(t) + \nabla \bar{u}(t) \right. \\ & \qquad \left. - (G_{\tau} R \tau \partial_t u(t) + G_{\tau} u(t)), G_{\tau} v(t) \rangle_{L^2(\Omega)^d} \right) dt \\ & \leq \left( \| (\nabla R \bar{u}' - G_{\tau} R \tau \partial_t u) \|_{L^2(L^2(\Omega)^d)} + \| (\nabla \bar{u} - G_{\tau} u) \|_{L^2(L^2(\Omega)^d)} \right) \\ & \qquad \times \| G_{\tau} v \|_{L^2(L^2(\Omega)^d)}. \end{aligned} $$

Dividing by $\|G_{\tau}v\|_{L^2(L^2(\Omega)^d)}$ and taking the supremum over $v \in V_{\tau} \setminus \{0\}$ shows that $\zeta_{\tau}^{(T)}(\mathbf{v}) \le \delta_{\tau}^{(T)}(\bar{u}, u)$, which concludes the proof.

Finally, as in the steady case, we notice that $(u^{(0)}, (u^{(m)})_{m \in [1, N]})$ is solution to a square linear system. Therefore the error estimate Theorem 3.1 proves that, for a null right-hand-side, the solution is null, which shows that the system is invertible. $\square$

The following lemma is proven in [9].

**Lemma A.6** ([9, Lemma 4.7]). Let $V$ and $L$ be Hilbert spaces. Let $Z$ and $Y$ be the Hilbert spaces defined by $Z = V \times V \times L \times L$ and $Y = V \times L$. Let $A: V \to V$ be an $M$-continuous and $\alpha$-coercive linear operator (with $M \ge 1$ and $\alpha > 0$), which means that

$$ \|Av\|_V \leq M \|v\|_V \quad \text{and} \quad \langle Av, v \rangle_V \geq \alpha \|v\|_V^2 \quad \forall v \in V. \quad (A.30) $$

Let $\Phi: L \to L$ be a linear operator such that $\|\Phi\| \le 1$, and let $\hat{b}: Z \times Y \to \mathbb{R}$ be defined by

$$ \hat{b}((z_1, z_2, z_3, z_4), (y_1, y_2)) = \langle z_1 + Az_2, y_1 \rangle_V + \langle z_3 - \Phi z_4, y_2 \rangle_L, \quad (A.31) $$

for all $(z_1, z_2, z_3, z_4) \in Z$ and for all $(y_1, y_2) \in Y$.

Let $X \subset Z$ be a subspace of $Z$. We define the Hilbert spaces $X_1 \subset V$, $X_2 \subset V$, $X_3 \subset L$ and $X_4 \subset L$ by: for $i = 1, \dots, 4$,

$$ X_i = \overline{\{x_i : x \in X\}}, \text{ where } x_i \text{ is the } i\text{-th component of } x \in Z. $$

Assume that

$$ X_1 \subset X_2, \quad (A.32) $$

and that there exist $\zeta > 0$ and $\delta > 0$ such that, for all $x \in X$,

$$ \langle x_1, x_2 \rangle_V + \frac{\alpha^2}{12 M^3} (\|x_2\|_V^2 + \|x_1\|_V^2) \geq \mu \|x_4\|_L^2 - \nu \|x_3\|_L^2, \quad (A.33) $$

for some $\mu \in (0, \zeta]$ and $\nu \in [0, \mu]$ with $\mu - \nu \|\Phi\|^2 \ge \delta$. Then, there exists $\hat{\beta} > 0$, only depending on $\alpha, M, \zeta$ and $\delta$ (and not on $\mu, \nu$ and $\|\Phi\|$) such that

$$ \sup_{y \in X_2 \times X_3, \|y\|_Y = 1} \hat{b}(x, y) \geq \hat{\beta} \|x\|_Z \quad \forall x \in X. \quad (A.34) $$
Let us now turn to the first inequality in (A.23); first note that

$$ \inf_{v \in W_{\tau}} \delta_{\tau}^{(T)}(\bar{u}, v) \leq \delta_{\tau}^{(T)}(\bar{u}, u). $$

To bound $\zeta_{\tau}^{(T)}(\mathbf{v})$ we recall that $\mathbf{v} = \nabla R \bar{u}' + \nabla \bar{u} + \mathbf{F}$ satisfies $-\mathrm{div} \mathbf{v} = f$ (see (A.6)), and use the scheme (A.13) (with $z=0$) to write, for any $v \in V_{\tau} \setminus \{0\}$,

$$ \begin{aligned} & \langle \mathbf{v}, G_{\tau} v \rangle_{L^2(\Omega)^d} + \langle \mathrm{div} \mathbf{v}, v \rangle_{L^2(\Omega)} \\ &= \int_0^T \left( \langle \nabla R \bar{u}'(t) + \nabla \bar{u}(t) \right. \\ & \qquad \left. - (G_{\tau} R \tau \partial_t u(t) + G_{\tau} u(t)), G_{\tau} v(t) \rangle_{L^2(\Omega)^d} \right) dt \\ & \leq \left( \| (\nabla R \bar{u}' - G_{\tau} R \tau \partial_t u) \|_{L^2(L^2(\Omega)^d)} + \| (\nabla \bar{u} - G_{\tau} u) \|_{L^2(L^2(\Omega)^d)} \right) \\ & \qquad \times \| G_{\tau} v \|_{L^2(L^2(\Omega)^d)}. \end{aligned} $$

Dividing by $\|G_{\tau}v\|_{L^2(L^2(\Omega)^d)}$ and taking the supremum over $v \in V_{\tau} \setminus \{0\}$ shows that $\zeta_{\tau}^{(T)}(\mathbf{v}) \le \delta_{\tau}^{(T)}(\bar{u}, u)$, which concludes the proof.

Finally, as in the steady case, we notice that $(u^{(0)}, (u^{(m)})_{m \in [1, N]})$ is solution to a square linear system. Therefore the error estimate Theorem 3.1 proves that, for a null right-hand-side, the solution is null, which shows that the system is invertible. $\square$

The following lemma is proven in [9].

**Lemma A.6** ([9, Lemma 4.7]). Let $V$ and $L$ be Hilbert spaces. Let $Z$ and $Y$ be the Hilbert spaces defined by $Z = V \times V \times L \times L$ and $Y = V \times L$. Let $A: V \to V$ be an $M$-continuous and $\alpha$-coercive linear operator (with $M \ge 1$ and $\alpha > 0$), which means that

$$ \|Av\|_V \le M \|v\|_V \quad \text{and} \quad \langle Av, v \rangle_V \ge \alpha \|v\|_V^2 \quad \forall v \in V. \quad (A.30) $$

Let $\Phi: L \to L$ be a linear operator such that $\|\Phi\| \le 1$, and let $\hat{b}: Z \times Y \to \mathbb{R}$ be defined by

$$ \hat{b}((z_1, z_2, z_3, z_4), (y_1, y_2)) = \langle z_1 + Az_2, y_1 \rangle_V + \langle z_3 - \Phi z_4, y_2 \rangle_L, \quad (A.31) $$

for all $(z_1, z_2, z_3, z_4) \in Z$ and for all $(y_1, y_2) \in Y$.

Let $X \subset Z$ be a subspace of $Z$. We define the Hilbert spaces $X_1 \subset V$, $X_2 \subset V$, $X_3 \subset L$ and $X_4 \subset L$ by: for $i = 1, \dots, 4$,

$$ X_i = \overline{\{x_i : x \in X\}}, \text{ where } x_i \text{ is the } i\text{-th component of } x \in Z. $$

Assume that

$$ X_1 \subset X_2, \quad (A.32) $$

and that there exist $\zeta > 0$ and $\delta > 0$ such that, for all $x \in X$,

$$ \langle x_1, x_2 \rangle_V + \frac{\alpha^2}{12 M^3} (\|x_2\|_V^2 + \|x_1\|_V^2) \ge \mu \|x_4\|_L^2 - \nu \|x_3\|_L^2, \quad (A.33) $$

for some $\mu \in (0, \zeta]$ and $\nu \in [0, \mu]$ with $\mu - \nu \|\Phi\|^2 \ge \delta$. Then, there exists $\hat{\beta} > 0$, only depending on $\alpha, M, \zeta$ and $\delta$ (and not on $\mu, \nu$ and $\|\Phi\|$) such that

$$ \sup_{y \in X_2 \times X_3, \|y\|_Y = 1} \hat{b}(x, y) \ge \hat{\beta} \|x\|_Z \quad \forall x \in X. \quad (A.34) $$
REFERENCES

[1] W. Arendt, I. Chalendar, and R. Eymard. Lions' representation theorem and applications. *J. Math. Anal. Appl.*, 522(2):Paper No. 126946, 23, 2023.

[2] W. Arendt, I. Chalendar, and R. Eymard. Space-time error estimates for approximations of linear parabolic problems with generalized time boundary conditions. accepted for publication in IMAJNA, Apr. 2023.

[3] K. Aziz and A. Settari. *Petroleum Reservoir Simulation*. Society of Petroleum Engineers.

[4] S. Boivin, F. Cayré, and J.-M. Hérard. A finite volume method to solve the navier–stokes equations for incompressible flows on unstructured meshes. *International Journal of Thermal Sciences*, 39(8):806–825, 2000.

[5] C. Cancès, C. Chainais-Hillairet, M. Herda, and S. Krell. Large time behavior of nonlinear finite volume schemes for convection-diffusion equations. *SIAM J. Numer. Anal.*, 58(5):2544–2571, 2020.

[6] Y. Coudière, M. Saad, and A. Uzureau. Analysis of a finite volume method for a bone growth system in vivo. *Computers & Mathematics with Applications*, 66(9):1581–1594, 2013. *BioMath 2012*.

[7] J. Droniou, R. Eymard, T. Gallouët, C. Guichard, and R. Herbin. *The gradient discretisation method*, volume 82 of Mathematics & Applications. Springer, 2018.

[8] J. Droniou, R. Eymard, T. Gallouët, C. Guichard, and R. Herbin. Optimal error estimates for non-conforming approximations of linear parabolic problems with minimal regularity. working paper or preprint, Aug. 2023.

[9] J. Droniou, R. Eymard, T. Gallouët, C. Guichard, and R. Herbin. Optimal error estimates for non-conforming approximations of linear parabolic problems with minimal regularity. accepted for publication in *SeMA, Aug. 2023*.

[10] J. Droniou, R. Eymard, T. Gallouët, and R. Herbin. A unified approach to mimetic finite difference, hybrid finite volume and mixed finite volume methods. *Math. Models Methods Appl. Sci.*, 20(2):265–295, 2010.

[11] J. Droniou, T. Gallouët, and R. Herbin. A finite volume scheme for a noncoercive elliptic equation with measure data. *SIAM J. Numer. Anal.*, 41(6):1997–2031, 2003.

[12] J. Droniou and T. Gallouët. Finite volume methods for convection-diffusion equations with right-hand side in $h^{-1}$. *ESAIM: Mathematical Modelling and Numerical Analysis*, 36(4):705–724, 2002.

[13] J. Droniou and N. Nataraj. Improved $L^2$ estimate for gradient schemes and super-convergence of the TPFA finite volume scheme. *IMA J. Numer. Anal.*, 38(3):1254–1293, 2018.

[14] Q. Du, M. D. Gunzburger, and L. Ju. Voronoï-based finite volume methods, optimal Voronoï meshes, and pdes on the sphere. *Computer Methods in Applied Mechanics and Engineering*, 192(35):3933–3957, 2003.

[15] Q. Du and L. Ju. Finite volume methods on spheres and spherical centroidal Voronoï meshes. *SIAM Journal on Numerical Analysis*, 43(4):1673–1692, 2005.

[16] R. Eymard, J. Fuhrmann, and K. Gärtner. A finite volume scheme for nonlinear parabolic equations derived from one-dimensional local Dirichlet problems. *Numer. Math.*, 102(3):463–495, 2006.

[17] R. Eymard and T. Gallouët. H-convergence and numerical schemes for elliptic problems. *SIAM J. Numer. Anal.*, 41(2):539–562, 2003.

[18] R. Eymard, T. Gallouët, and R. Herbin. Finite volume methods. In P. Ciarlet and J. Lions, editors, *Handbook of Numerical Analysis, Volume VII*, pages 713–1020. North Holland, 2000.

[19] R. Eymard, T. Gallouët, and R. Herbin. A cell-centred finite-volume approximation for anisotropic diffusion operators on unstructured meshes in any space dimension. *IMA Journal of Numerical Analysis*, 26(2):326–353, 04 2006.

[20] P. Forsyth and P. Sammon. Quadratic convergence for cell-centered grids. *Applied Numerical Mathematics*, 4(5):377–394, 1988.

[21] R. Herbin and F. Hubert. Benchmark on discretization schemes for anisotropic diffusion problems on general grids. In *Finite volumes for complex applications V*, pages 659–692. ISTE, London, 2008.

[22] R. Herbin and O. Labergerie. Finite volume schemes for elliptic and elliptic-hyperbolic problems on triangular meshes. *Comput. Methods Appl. Mech. Engrg.*, 147(1-2):85–103, 1997.

[23] J. L. Lions. *Équations différentielles opérationnelles et problèmes aux limites*, volume 111. Springer, Cham, 1961.

[24] I. D. Mishev. Finite volume methods on Voronoï meshes. *Numerical Methods for Partial Differential Equations*, 14:193–212, 1998.

[25] P. Omnes. On the second-order convergence of a function reconstructed from finite volume approximations of the laplace equation on delaunay-Voronoï meshes. *ESAIM: M2AN*, 45(4):627–650, 2011.
REFERENCES

[1] W. Arendt, I. Chalendar, and R. Eymard. Lions' representation theorem and applications. *J. Math. Anal. Appl.*, 522(2):Paper No. 126946, 23, 2023.

[2] W. Arendt, I. Chalendar, and R. Eymard. Space-time error estimates for approximations of linear parabolic problems with generalized time boundary conditions. accepted for publication in IMAJNA, Apr. 2023.

[3] K. Aziz and A. Settari. *Petroleum Reservoir Simulation*. Society of Petroleum Engineers.

[4] S. Boivin, F. Cayré, and J.-M. Hérard. A finite volume method to solve the navier–stokes equations for incompressible flows on unstructured meshes. *International Journal of Thermal Sciences*, 39(8):806–825, 2000.

[5] C. Cancès, C. Chainais-Hillairet, M. Herda, and S. Krell. Large time behavior of nonlinear finite volume schemes for convection-diffusion equations. *SIAM J. Numer. Anal.*, 58(5):2544–2571, 2020.

[6] Y. Coudière, M. Saad, and A. Uzureau. Analysis of a finite volume method for a bone growth system in vivo. *Computers & Mathematics with Applications*, 66(9):1581–1594, 2013. *BioMath 2012*.

[7] J. Droniou, R. Eymard, T. Gallouët, C. Guichard, and R. Herbin. *The gradient discretisation method*, volume 82 of Mathematics & Applications. Springer, 2018.

[8] J. Droniou, R. Eymard, T. Gallouët, C. Guichard, and R. Herbin. Optimal error estimates for non-conforming approximations of linear parabolic problems with minimal regularity. working paper or preprint, Aug. 2023.

[9] J. Droniou, R. Eymard, T. Gallouët, C. Guichard, and R. Herbin. Optimal error estimates for non-conforming approximations of linear parabolic problems with minimal regularity. accepted for publication in SeMA, Aug. 2023.

[10] J. Droniou, R. Eymard, T. Gallouët, and R. Herbin. A unified approach to mimetic finite difference, hybrid finite volume and mixed finite volume methods. *Math. Models Methods Appl. Sci.*, 20(2):265–295, 2010.

[11] J. Droniou, T. Gallouët, and R. Herbin. A finite volume scheme for a noncoercive elliptic equation with measure data. *SIAM J. Numer. Anal.*, 41(6):1997–2031, 2003.

[12] J. Droniou and T. Gallouët. Finite volume methods for convection-diffusion equations with right-hand side in $h^{-1}$. *ESAIM: Mathematical Modelling and Numerical Analysis*, 36(4):705–724, 2002.

[13] J. Droniou and N. Nataraj. Improved $L^2$ estimate for gradient schemes and super-convergence of the TPFA finite volume scheme. *IMA J. Numer. Anal.*, 38(3):1254–1293, 2018.

[14] Q. Du, M. D. Gunzburger, and L. Ju. Voronoï-based finite volume methods, optimal Voronoï meshes, and pdes on the sphere. *Computer Methods in Applied Mechanics and Engineering*, 192(35):3933–3957, 2003.

[15] Q. Du and L. Ju. Finite volume methods on spheres and spherical centroidal Voronoï meshes. *SIAM Journal on Numerical Analysis*, 43(4):1673–1692, 2005.

[16] R. Eymard, J. Fuhrmann, and K. Gärtner. A finite volume scheme for nonlinear parabolic equations derived from one-dimensional local Dirichlet problems. *Numer. Math.*, 102(3):463–495, 2006.

[17] R. Eymard and T. Gallouët. H-convergence and numerical schemes for elliptic problems. *SIAM J. Numer. Anal.*, 41(2):539–562, 2003.

[18] R. Eymard, T. Gallouët, and R. Herbin. Finite volume methods. In P. Ciarlet and J. Lions, editors, *Handbook of Numerical Analysis, Volume VII*, pages 713–1020. North Holland, 2000.

[19] R. Eymard, T. Gallouët, and R. Herbin. A cell-centred finite-volume approximation for anisotropic diffusion operators on unstructured meshes in any space dimension. *IMA Journal of Numerical Analysis*, 26(2):326–353, 04 2006.

[20] P. Forsyth and P. Sammon. Quadratic convergence for cell-centered grids. *Applied Numerical Mathematics*, 4(5):377–394, 1988.

[21] R. Herbin and F. Hubert. Benchmark on discretization schemes for anisotropic diffusion problems on general grids. In *Finite volumes for complex applications V*, pages 659–692. ISTE, London, 2008.

[22] R. Herbin and O. Labergerie. Finite volume schemes for elliptic and elliptic-hyperbolic problems on triangular meshes. *Comput. Methods Appl. Mech. Engrg.*, 147(1-2):85–103, 1997.

[23] J. L. Lions. *Équations différentielles opérationnelles et problèmes aux limites*, volume 111. Springer, Cham, 1961.

[24] I. D. Mishev. Finite volume methods on Voronoï meshes. *Numerical Methods for Partial Differential Equations*, 14:193–212, 1998.

[25] P. Omnes. On the second-order convergence of a function reconstructed from finite volume approximations of the laplace equation on delaunay-Voronoï meshes. *ESAIM: M2AN*, 45(4):627–650, 2011.
[26] G. Strang. Variational crimes in the finite element method. In *The mathematical foundations of the finite element method with applications to partial differential equations (Proc. Sympos., Univ. Maryland, Baltimore, Md., 1972)*, pages 689–710. Academic Press, New York, 1972.

[27] R. Temam. *Navier-Stokes equations*, volume 2 of *Studies in Mathematics and its Applications*. North-Holland Publishing Co., Amsterdam, third edition, 1984. Theory and numerical analysis, With an appendix by F. Thomasset.

ROBERT EYMARD: **UNIVERSITÉ GUSTAVE EIFFEL, LAMA, (UMR 8050), UPEM, UPEC, CNRS, F-77454, MARNE-LA-VALLEÉ (FRANCE)**

Email address: robert.eymard@univ-eiffel.fr

THIERRY GALLOUËT: I2M UMR 7373, AIX-MARSEILLE UNIVERSITÉ, CNRS, ECOLE CENTRALE DE MARSEILLE, F-13453 MARSEILLE, FRANCE

Email address: thierry.gallouet@univ-amu.fr

RAPHAÈLE HERBIN: I2M UMR 7373, AIX-MARSEILLE UNIVERSITÉ, CNRS, ECOLE CENTRALE DE MARSEILLE, F-13453 MARSEILLE, FRANCE

Email address: raphaele.herbin@univ-amu.fr
[26] G. Strang. Variational crimes in the finite element method. In *The mathematical foundations of the finite element method with applications to partial differential equations (Proc. Sympos., Univ. Maryland, Baltimore, Md., 1972)*, pages 689–710. Academic Press, New York, 1972.

[27] R. Temam. *Navier-Stokes equations*, volume 2 of *Studies in Mathematics and its Applications*. North-Holland Publishing Co., Amsterdam, third edition, 1984. Theory and numerical analysis, With an appendix by F. Thomasset.

ROBERT EYMARD: **UNIVERSITÉ GUSTAVE EIFFEL, LAMA, (UMR 8050), UPEM, UPEC, CNRS, F-77454, MARNE-LA-VALLEE (FRANCE)**

Email address: robert.eymard@univ-eiffel.fr

THIERRY GALLOUËT: I2M UMR 7373, AIX-MARSEILLE UNIVERSITÉ, CNRS, ECOLE CENTRALE DE MARSEILLE, F-13453 MARSEILLE, FRANCE

Email address: thierry.gallouet@univ-amu.fr

RAPHAÈLE HERBIN: I2M UMR 7373, AIX-MARSEILLE UNIVERSITÉ, CNRS, ECOLE CENTRALE DE MARSEILLE, F-13453 MARSEILLE, FRANCE

Email address: raphaele.herbin@univ-amu.fr