**Remark 6.6 (Error structure).** From (6.26), (6.28), and (6.32a), we see that

$$ 
\begin{aligned}
c_{\tilde{K}}^{1/2} \|\mathbf{u} - \mathbf{u}_h^{k,i}\| & \le c_{\tilde{K}}^{-1/2} C_{\tilde{K}} \max_{\varphi \in H_0^1(\Omega), \|\nabla \varphi\|=1} \{(f, \varphi) + (\mathbf{u}_h^{k,i}, \nabla \varphi)\} + c_{\tilde{K}}^{-1/2} \min_{v \in H_0^1(\Omega)} \|\tilde{K}(|\mathbf{u}_h^{k,i}|)\mathbf{u}_h^{k,i} + \nabla v\| \\
& \le 2c_{\tilde{K}}^{-1/2} C_{\tilde{K}} \|\mathbf{u} - \mathbf{u}_h^{k,i}\|;
\end{aligned}
\quad (6.33) $$

in the second inequality, we have used the Cauchy-Schwarz inequality and (6.8) together with (6.3b) for the second term. This means that the energy error $c_{\tilde{K}}^{1/2} \|\mathbf{u}-\mathbf{u}_h^{k,i}\|$ is equivalent to the sum of the dual norm of the residual and nonconformity evaluated as the distance to the $H_0^1(\Omega)$ space. This is an immediate extension of Theorem 4.9 from the linear case. Compared to Corollary 4.6 for the linear case with an exact algebraic solver giving $\nabla \cdot \mathbf{u}_h = f$, this in particular enables to treat the nonlinear case with inexact linearization and algebraic solvers, so that $\nabla \cdot \mathbf{u}_h^{k,i} \neq f$. Our motivation here was to evaluate the error as the $L^2(\Omega)$-norm, which gives rise to the weights $c_{\tilde{K}}$ and $C_{\tilde{K}}$ and the factor 2 in (6.33). Alternatively, we could have stick to the intrinsic error measure of the form

$$ \left\{ \max_{\varphi \in H_0^1(\Omega), \|\nabla \varphi\|=1} \{(f, \varphi) + (\mathbf{u}_h^{k,i}, \nabla \varphi)\}^2 + \min_{v \in H_0^1(\Omega)} \|\tilde{K}(|\mathbf{u}_h^{k,i}|)\mathbf{u}_h^{k,i} + \nabla v\|^2 \right\}^{\frac{1}{2}}, $$

cf. (4.9), which would lead to constant-free estimates without the weights $c_{\tilde{K}}$ and $C_{\tilde{K}}$. In this way, we define the error measure $\mathcal{N}^{n,k,i}$ in (7.17a) in the complex multiphase flow below.

## 6.8 Numerical experiments

In this section, we numerically illustrate the efficiency of our theoretical results of Section 6.7 on two different examples. Our main goals are to assess the sharpness of the guaranteed bound (6.24) and to examine the robustness of our estimates with respect to the ratio $C_{\tilde{K}}/c_{\tilde{K}}$ where we consider $C_{\tilde{K}}/c_{\tilde{K}} = 10^i, i \in \{1, ..., 6\}$. As in previous sections, the effectivity index is defined by the ratio of the estimator from Theorem 6.2 to the exact error,

$$ I_{\text{eff}} := \frac{\eta^{k,i}}{c_{\tilde{K}}^{1/2} \|\mathbf{u} - \mathbf{u}_h^{k,i}\|}. \quad (6.34) $$

Additionally, we study in the following examples a stopping criterion for the linearization algorithm based on balancing the error components of Theorem 6.2. The criterion is that the linearization iteration is pursued until step $k$ such that

$$ \eta_{\text{lin}}^{k,i} \le \Gamma_{\text{lin}} \eta_{\text{sp}}^{k,i}, \quad (6.35) $$

with $\Gamma_{\text{lin}} \in (0, 1)$ (we focus on the linearization only, the algebraic solver is used with a very small tolerance yielding $\eta_{\text{alg}}^{k,i}, \eta_{\text{rem}}^{k,i} \approx 0$). We compare this criterion with a classic one based on a fixed threshold on the relative linearization residual

$$ e_{\text{lin}}^{k} \le 10^{-8}, \quad (6.36) $$

where the relative linearization residual, for a resolution of a system of nonlinear algebraic equations $F(X) = 0$ by the Newton method, is given by

$$ e_{\text{lin}}^{k} := \frac{\|F(X^{k})\|}{\|F(X^{0})\|}. $$

Note that the criterion (6.35) can be used in more general adaptive algorithm balancing all the error components such as in Algorithm 7.1.

### 6.8.1 Regular solution

In this test, we consider on the domain $\Omega = (0, 1) \times (0, 1)$ the exact solution

$$ p(x, y) = 16x(1-x)y(1-y). \quad (6.37) $$

The source term $f$ is then constructed through (6.1). We perform the numerical tests with the nonlinearity example given by (6.5) on a uniformly refined rectangular mesh of $10^4$ elements. Figure 24 shows the
<figure><img src="image_27.png" /><figcaption>Figure 24: [Section 6.8.1] The approximate solution given by the values $p_K^{k,i}$, $K \in \mathcal{T}_h$, (left), the energy errors $c_{\tilde{K}}^{1/2} \|\mathbf{u}-\mathbf{u}_h^{k,i}\|_K$ (middle), and the error estimators $\eta_{\text{sp},K}^{k,i}+\eta_{\text{lin},K}^{k,i}+\eta_{\text{alg},K}^{k,i}+\eta_{\text{rem},K}^{k,i}$ from Theorem 6.2 (right).</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Figure 25: [Section 6.8.1], Relative errors $c_{\tilde{K}}^{1/2} \|\mathbf{u}-\mathbf{u}_h^{k,i}\|/\|\mathbf{u}_h^{k,i}\|$ and relative estimates $\eta^{k,i}/\|\mathbf{u}_h^{k,i}\|$ from Theorem 6.2 (left), effectivity indices $I_{\text{eff}}$ from (6.34) (right).</figcaption></figure>

approximate solution given by the values $p_K^{k,i}$, $K \in \mathcal{T}_h$, the elementwise errors $c_{\tilde{K}}^{1/2} \|\mathbf{u}-\mathbf{u}_h^{k,i}\|_K$, and the corresponding a posteriori error estimators $\eta_K^{k,i}$ from Theorem 6.2.

The robustness of our estimates with respect to the ratio $C_{\tilde{K}}/c_{\tilde{K}}$ is examined in Figure 25, where we fix the value of $c_{\tilde{K}} = 1$ and plot, for different values of $C_{\tilde{K}} \in \{10, 100, \dots, 10^6\}$, the values of the relative energy error $c_{\tilde{K}}^{1/2} \|\mathbf{u}-\mathbf{u}_h^{k,i}\|/\|\mathbf{u}_h^{k,i}\|$ and the corresponding relative total error estimate $\eta^{k,i}/\|\mathbf{u}_h^{k,i}\|$ from Theorem 6.2. The effectivity index plotted in the right part of Figure 25 reveal that, in this case and for the choice of nonlinearity example (6.5), the estimate is robust and not affected by the values of the ratio $C_{\tilde{K}}/c_{\tilde{K}}$.

Figure 26 is dedicated to study the stopping criterion (6.35) and its comparison with the usual stopping criterion (6.36). The left part of Figure 26 depicts the evolution of the total estimator $\eta^{k,i}$, the linearization estimator $\eta_{\text{lin}}^{k,i}$, and the relative linearization residual $e_{\text{lin}}^k$ as a function of the number of Newton iterations for a uniform fixed mesh of $10^4$ elements and with the ratio $C_{\tilde{K}}/c_{\tilde{K}} = 10$. We observe that the linearization estimator and the relative linearization residual steadily decrease, while the total estimator stagnate starting from the third iteration. The stopping criterion (6.35) with $\Gamma_{\text{lin}} = 0.1$ indicates that the resolution can be stopped at the sixth iteration, which already yields a sufficiently accurate approximate solution and helps avoid unnecessary additional iterations. The right part of Figure 26 shows a comparison of the required number of Newton iteration in order to satisfy the standard criterion $e_{\text{lin}}^k \le 10^{-8}$ and the adaptive one $\eta_{\text{lin}}^{k,i} \le 0.1\eta_{\text{sp},K}^{k,i}$ for different values of the ratio $C_{\tilde{K}}/c_{\tilde{K}}$. We see that we can save at least the half of Newton iterations on each resolution without altering the precision.
<figure><img src="image_29.png" /><figcaption>Figure 26: [Section 6.8.1], Standard linearization vs. adaptive linearization</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>Figure 27: [Section 6.8.2] The approximate solution given by the values $p_K^{k,i}$, $K \in \mathcal{T}_h$, (left), the energy errors $c_{\bar{K}}^{1/2} \|u - u_h^{k,i}\|_K$ (middle), and the error estimators $\eta_{sp,K}^{k,i} + \eta_{lin,K}^{k,i} + \eta_{alg,K}^{k,i} + \eta_{rem,K}^{k,i}$ from Theorem 6.2 (right).</figcaption></figure>

## 6.8.2 Singular solution

In this test, we consider the L-shaped domain

$$ \Omega = (-1, 1) \times (-1, 1) \setminus (-1, 0] \times (-1, 0] $$

with the exact solution

$$ p(r, \theta) = r^{\frac{2}{3}} \sin \left( \frac{2}{3}\theta + \frac{3}{2}\pi \right), $$

where $(r, \theta)$ are the polar coordinates. The Dirichlet boundary condition is partly inhomogeneous and given by the value of the exact solution on $\partial\Omega$. The following results are carried out with the choice of the nonlinearity example given by (6.5) and on a uniformly refined rectangular mesh of 10800 elements.

In Figure 27, we show the distribution of the approximate solution and we compare the distributions of the errors $c_{\bar{K}}^{1/2} \|u - u_h^{k,i}\|_K$ and the corresponding a posteriori error estimators $\eta_K^{k,i}$ from Theorem 6.2. We can observe that the estimators identify perfectly the error in the zone around the singularity at (0,0).

In Figure 28, we study the effect of the ratio $C_{\bar{K}}/c_{\bar{K}}$ on our estimators. The left part of Figure 28 shows the relative exact errors $c_{\bar{K}}^{1/2} \|u - u_h^{k,i}\|/\|u_h^{k,i}\|$ together with their relative a posteriori error estimates $\eta^{k,i}/\|u_h^{k,i}\|$ from Theorem 6.2. The right part of Figure 28 displays the corresponding effectivity indices,
<figure><img src="image_31.png" /><figcaption>Figure 28: [Section 6.8.2], Relative errors $c_{\bar{K}}^{1/2} \|\mathbf{u} - \mathbf{u}_h^{k,i}\| / \|\mathbf{u}_h^{k,i}\|$ and relative estimates $\eta^{k,i} / \|\mathbf{u}_h^{k,i}\|$ from Theorem 6.2 (left), effectivity indices $I_{\text{eff}}$ from (6.34) (right)</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>Figure 29: [Section 6.8.2], Standard linearization vs. adaptive linearization</figcaption></figure>

given by the ratio of the estimate over the error (6.34). We see that the effectivity indices are not affected by the values of the ratio $C_{\bar{K}}/c_{\bar{K}}$, which numerically gives the robustness of the estimate.

In order to evaluate the potential of a stopping criterion for the linearization algorithm, we illustrate in the left part of Figure 29 the evolution of the total estimator $\eta^{k,i}$, the linearization estimator $\eta_{\text{lin}}^{k,i}$, and the relative linearization residual $e_{\text{lin}}^k$ as a function of the number of Newton iterations for a uniform fixed mesh of 10800 elements and with the ratio $C_{\bar{K}}/c_{\bar{K}} = 10$. These results are given with a standard resolution stopped when $e_{\text{lin}}^k \le 10^{-8}$. The linearization estimator and the relative linearization residual decrease as expected while the total estimator start to stagnate from the fifth iteration. Considering the stopping criterion (6.35) with $\Gamma_{\text{lin}} = 0.1$ suggests that the resolution process can be halted after the seventh iteration, as it already provides a sufficiently accurate approximate solution while avoiding a significant number of unnecessary linearization iterations. In the right part of Figure 29, computational savings in terms of linearization iterations is observed for different values of the ratio $C_{\bar{K}}/c_{\bar{K}}$. These results make it apparent that using criteria based on balancing the error components gives a powerful tool for adaptivity algorithms. This is shown in the next section by considering Algorithm 7.1 applied to a realistic problem and achieving a significant gain in terms of the computation effort.

## 6.9 Bibliographic resources

One of the first a posteriori error estimates for finite volume methods for nonlinear problems is presented in Bergam et al. [48], though it concerns vertex-centered finite volumes, yielding $H_0^1(\Omega)$-conforming approximations as in the finite element case. Kim [163] was one of the first to consider locally conservative methods in a Prager-Synge-type framework of the form (6.32).

The distinction of different error components (and the associated adaptivity) has been initiated in
Han [143], Becker et al. [39], Kaltenbacher [156], Axelsson and Kaporin [27], Patera and Rønquist [211], Arioli [24], Bakushinsky and Smirnova [30], Chaillou and Suri [75, 74], Picasso [215], Ahusborde et al. [11], Bernardi et al. [49, 50], Jiránek et al. [154], El Alaoui et al. [114], Nordbotten et al. [200], Arioli et al. [26, 25], Rannacher and Vihharev [217], Ern and Vohralík [124], Keilegavlen and Nordbotten [158, 159], Rey et al. [222, 221, 220], and Becker et al. [38], see also the references therein. It has been reconsidered more recently, namely in the framework of conforming finite elements, in Heid and Wihler [147], Heid et al. [146], Gantner et al. [136], Haberl et al. [142], Kumar et al. [171], Chaudhry et al. [77], Dolgov and Vejchodský [106], Liu and Keyes [182], Févotte et al. [134], Bringmann et al. [66], Harnist et al. [145], Mitra and Vohralík [189], Stokke et al. [233], and Papež [209]. The classical books of Kelley [160], Saad [229], or Olshanskii and Tyrtyshnikov [205] do not discuss total a posteriori error estimates during iterative linearization and iterative linear algebraic solvers iteration but rather only treat the nonlinear algebra level (for iterative linearization) and linear algebra level (for iterative linear algebraic solvers), not addressing the discretization/partial differential equation level.

# 7 Unsteady nonlinear coupled degenerate advection-diffusion-reaction problems. Temporal discretization error

In this section, the methodology derived in Sections 4–6 is applied to polytopal discretizations of a model real-life multiphase compositional Darcy flow in porous media. In the following, we describe briefly the multiphase model, discuss its discretization by an implicit finite volume scheme, the linearization by the Newton or combined Newton and fixed point method, and the algebraic solution of the arising linear system. We explain that the required fictitious flux and potential reconstructions stay unchanged with respect to Section 6 and finally develop our a posteriori error estimate with inexpensive implementation and evaluation distinguishing the different error components.

By decomposing the estimators into the space, time, linearization, and algebraic error components, we can design an algorithm involving adaptivity of both solvers of meshes. We namely formulate criteria for stopping the iterative algebraic solver and the iterative linearization solver when the corresponding error components do not affect significantly the total error. Moreover, the spatial and temporal error components are balanced respectively by time step and space mesh adaptation, materializing the abstract considerations of Section 2.6. Denote by $\eta_{sp}^{n,k,i}$, $\eta_{tm}^{n,k,i}$, $\eta_{lin}^{n,k,i}$, and $\eta_{alg}^{n,k,i}$ respectively the total spatial, temporal, linearization, and algebraic estimators, on time step $n$, linearization step $k$, and algebraic solver step $i$. Let $\gamma_{lin}, \gamma_{alg} \in (0, 1)$ and $\Gamma_{tm} > \gamma_{tm} > 0$ be user-given fixed balancing parameters and fix the fractions (percentages) of cells to refine $\delta_{ref}$ and to derefine $\delta_{deref}$, $0 < \delta_{deref} < \delta_{ref} < 1$. Following [124, 253, 70, 103, 104] and the references therein, we in particular derive in this section the following fully adaptive algorithm:

**Algorithm 7.1 (Adaptive stopping criteria and adaptive time and space mesh refinement).**

Set $n := 0$.
while $t^n \le t_F$ do {time loop}
    Set $n := n + 1$.
    loop {spatial and temporal errors balancing loop}
        Set $k := 0$.
        loop {linearization loop}
            Set $k := k + 1$.
            Set up the linear system.
            Set $i := 0$.
            loop {algebraic solver loop}
                Perform a step of the iterative algebraic solver and set $i := i + 1$.
                Evaluate the different a posteriori error estimators.
                Terminate if: $\eta_{alg}^{n,k,i} \le \gamma_{alg}\eta_{sp}^{n,k,i}$.
                end loop {algebraic solver loop}
                Terminate if: $\eta_{lin}^{n,k,i} \le \gamma_{lin}\eta_{sp}^{n,k,i}$.
                end loop {linearization loop}
                Terminate if

$$ \eta_{sp,K}^{n,k,i} \ge \delta_{ref} \max_{L \in \mathcal{T}_H^n} \{\eta_{sp,L}^{n,k,i}\} \quad \forall K \in \mathcal{T}_H^n \\ \gamma_{tm} \eta_{sp}^{n,k,i} \le \eta_{tm}^{n,k,i} \le \Gamma_{tm} \eta_{sp}^{n,k,i}; $$
**else**

Refine the cells $K \in \mathcal{T}_H^n$ such that $\eta_{sp,K}^{n,k,i} \ge \delta_{\text{ref}} \max_{L \in \mathcal{T}_H^n} \{\eta_{sp,L}^{n,k,i}\}$.

Derefine the cells $K \in \mathcal{T}_H^n$ such that $\eta_{sp,K}^{n,k,i} \le \delta_{\text{deref}} \max_{L \in \mathcal{T}_H^n} \{\eta_{sp,L}^{n,k,i}\}$.

Refine the time step if $\eta_{\text{tm}}^{n,k,i} > \Gamma_{\text{tm}} \eta_{\text{sp}}^{n,k,i}$, derefine the time step if $\gamma_{\text{tm}} \eta_{\text{sp}}^{n,k,i} > \eta_{\text{tm}}^{n,k,i}$.

**end loop** {spatial and temporal errors balancing loop}

Update data.

**end while** {time loop}

## 7.1 Multiphase compositional unsteady Darcy flow

We consider the multiphase compositional Darcy flow model in porous media [4, 28, 93, 83, 261]. We use the generalization of Coats’ formulation [94] to an arbitrary number of phases by [129].

Let $\mathcal{P} = \{p\}$ be the set of phases, $\mathcal{C} = \{c\}$ the set of components, and, for a given phase $p \in \mathcal{P}$, let $C_p \subset \mathcal{C}$ be the set of its components. $S_p$ then denotes the saturation of the phase $p$ and $C_{p,c}$ the molar fraction of the component $c$ in the phase $p$. For a given component $c \in \mathcal{C}$, denote by $\mathcal{P}_c$ the set of the phases which contain $c$.

We denote by $\mathcal{P}$ the reference pressure such that the phase pressures $P_p, p \in \mathcal{P}$, are expressed as

$$P_p := P + P_{c_p}(S_p), \qquad (7.1)$$

where $P_{c_p}: \mathbb{R} \to \mathbb{R}$ is a given generalized capillary pressure function. We collect the unknowns of the model in the vector

$$\mathbf{X} := (P, (S_p)_{p \in \mathcal{P}}, (C_{p,c})_{p \in \mathcal{P}, c \in \mathcal{C}_p}). \qquad (7.2)$$

For a phase $p \in \mathcal{P}$, let $k_{r,p}: \mathbb{R} \to \mathbb{R}$ be the relative permeability function of the unknown saturation $S_p$. We denote by $\phi$ the porosity of the medium, by $\mathbf{K}$ the permeability tensor, by $\mu_p$ the dynamic viscosity, by $\zeta_p$ the molar density, and by $\rho_p$ the mass density. For simplicity, we suppose that all these porous media parameters are given as constants on each mesh element and on each time step; $\phi$ is actually supposed to be constant in time. We also define the mobility by $\nu_p := \zeta_p \frac{k_{r,p}}{\mu_p}$. Finally, let $q_c$ be the source or sink term, again piecewise constant on the space-time mesh.

Let $t_F > 0$ be the final simulation time. The system of governing equations is posed in $\Omega \times (0, t_F)$ by

$$\partial_t l_c + \nabla \cdot \boldsymbol{\theta}_c = q_c, \quad \forall c \in \mathcal{C}, \qquad (7.3)$$

where $l_c$ is the amount (in moles) of component $c$ per unit volume,

$$l_c := \phi \sum_{p \in \mathcal{P}_c} \zeta_p S_p C_{p,c}. \qquad (7.4)$$

In (7.3), for each component $c \in \mathcal{C}$, the component flux $\boldsymbol{\theta}_c$ has the following expression:

$$\boldsymbol{\theta}_c := \sum_{p \in \mathcal{P}_c} \boldsymbol{\theta}_{p,c}, \quad \boldsymbol{\theta}_{p,c} := \boldsymbol{\theta}_{p,c}(\mathbf{X}) = \nu_p C_{p,c} \mathbf{v}_p(P_p), \qquad (7.5)$$

where for all phases $p \in \mathcal{P}$, $\mathbf{v}_p(P_p)$ represents the average phase velocity given by Darcy’s law,

$$\mathbf{v}_p(P_p) := -\underline{\mathbf{K}}(\nabla P_p + \rho_p g \nabla z). \qquad (7.6)$$

Here, $z$ is the vertical space coordinate and $g$ is the gravitation acceleration constant. Note at this occasion that only the first part of the Darcy velocity,

$$\mathbf{u}_p(P_p) := -\underline{\mathbf{K}} \nabla P_p, \qquad (7.7)$$

is a gradient of the phase pressure. We will use this observation in our error measure below. We assume that no-flow boundary conditions are prescribed for all the component fluxes,

$$\boldsymbol{\theta}_c \cdot \mathbf{n}_{\Omega} = 0 \quad \text{on } \partial\Omega \times (0, t_F) \quad \forall c \in \mathcal{C}, \qquad (7.8)$$

where $\partial\Omega$ denotes the boundary of $\Omega$ and $\mathbf{n}_{\Omega}$ its unit outward normal. At $t = 0$, we prescribe the initial amount of each component,

$$l_c(\cdot, 0) = l_c^0 \quad \forall c \in \mathcal{C}. \qquad (7.9)$$
The previous PDEs is supplemented by a system of algebraic equations imposing the volume conservation

$$ \sum_{p \in \mathcal{P}} S_p = 1, \qquad (7.10) $$

the conservation of the quantity of matter

$$ \sum_{c \in \mathcal{C}_p} C_{p,c} = 1 \quad \forall p \in \mathcal{P}, \qquad (7.11) $$

and local thermodynamic equilibrium expressed by

$$ \sum_{c \in \mathcal{C}} (|\mathcal{P}_c| - 1) = \sum_{p \in \mathcal{P}} |\mathcal{C}_p| - |\mathcal{C}| \qquad (7.12) $$

equalities of fugacities, where we refer to [103] for details.

## 7.2 Space-time mesh

We suppose that the time simulation interval $[0, t_F]$ is partitioned into time intervals $I^n := [t^{n-1}, t^n]$, $1 \le n \le N$, given by the discrete times $0 = t^0 < \dots < t^n < \dots < t^N = t_F$. The discrete time steps are then $\tau^n := t^n - t^{n-1}$. Then, we let $(\mathcal{T}_H^n)_{0 \le n \le N}$ be a family of polytopal meshes of the space domain $\Omega$ defined in the sense of Section 3.6; we suppose that $\mathcal{T}_H^n$ has been obtained from $\mathcal{T}_H^{n-1}$ by refinement of some elements and by coarsening of some other. For all $0 \le n \le N$, we define $\mathcal{F}_H^n, \mathcal{F}_h^n, \mathcal{F}_{H,h}^n, \mathcal{F}_{K,h}^{n,\text{ext}}$, and $\mathcal{F}_{K,h}^{n,\text{int}}$ as in Section 3.6.

## 7.3 Weak solution and its properties

At this stage, we need to characterize a weak solution for the multiphase compositional model (7.1)–(7.12). We define

$$ X := L^2((0, t_F); H^1(\Omega)), \qquad (7.13a) $$

$$ Y := H^1((0, t_F); L^2(\Omega)), \qquad (7.13b) $$

We equip the space $X$ with the norm

$$ \|\varphi\|_X := \left\{ \sum_{n=1}^{N} \|\varphi\|_{X^n}^2 \right\}^{\frac{1}{2}}, \quad \|\varphi\|_{X^n}^2 := \int_{I^n} \sum_{K \in \mathcal{T}_H^n} \left\{ h_K^{-2} c_{\underline{K},K} \|\varphi\|_K^2 + \|\underline{K}^{\frac{1}{2}} \nabla \varphi\|_K^2 \right\} dt \quad \varphi \in X, \quad (7.14) $$

where $c_{\underline{K},K}$ is the smallest eigenvalue of $\underline{K}$ on the mesh element $K \in \mathcal{T}_H^n$ and $h_K$ is the diameter of $K$. Below, it will be convenient to use $X^n := L^2(I^n; H^1(\Omega))$. We suppose that there exists a unique weak solution of (7.1)–(7.12) in the following sense:

**Assumption 7.2 (Weak solution).** There exists a weak solution $X$ of (7.1)–(7.12) which can be characterized as follows:

$$ l_c \in Y \quad \forall c \in \mathcal{C}, \qquad (7.15a) $$

$$ P_p \in X \quad \forall p \in \mathcal{P}, \qquad (7.15b) $$

$$ \theta_c \in [L^2((0, t_F); L^2(\Omega))]^d \quad \forall c \in \mathcal{C}, \qquad (7.15c) $$

$$ \int_0^{t_F} \{ (\partial_t l_c, \varphi)(t) - (\theta_c, \nabla \varphi)(t) \} dt = \int_0^{t_F} (q_c, \varphi)(t) dt \quad \forall \varphi \in X, \forall c \in \mathcal{C}, \qquad (7.15d) $$

the initial condition (7.9) holds, \hfill (7.15e)

and the algebraic closure equations (7.10)–(7.12) hold, \hfill (7.15f)

where $P_p, l_c$, and $\theta_c$ are defined, respectively, by (7.1), (7.4), and (7.5).

We have formulated the existence and uniqueness of a weak solution as an assumption, since, in contrast to the model problems of Sections 4–6, the problem (7.1)–(7.12) is too abstract and complicated to prove it. In simplified settings, with possibly only two phases present and each phase composed of a single component, such results can be found in [169, 78, 23, 81, 69, 161, 19] and the references therein.
**Remark 7.3** (Component fluxes). Since we assume $q_c$ piecewise constant and thus $q_c \in L^2((0, t_F); L^2(\Omega))$, it follows from (7.15) that actually

$$ \theta_c \in L^2((0, t_F); \mathbf{H}(\mathrm{div}, \Omega)) \quad \forall c \in \mathcal{C}, \qquad (7.16a) $$

$$ \nabla \cdot \theta_c = q_c - \partial_t l_c \quad \forall c \in \mathcal{C}, \qquad (7.16b) $$

$$ \theta_c \cdot n_{\Omega} = 0 \quad \text{on } \partial\Omega \times (0, t_F) \quad \forall c \in \mathcal{C}, \qquad (7.16c) $$

so that the component fluxes $\theta_c$ have the normal trace continuous in a proper sense, the governing equation (7.3) is satisfied with a weak divergence, and the boundary conditions (7.8) hold in the normal trace sense. This extends Proposition 4.1, (5.4), and (6.9) to the multiphase compositional model (7.1)–(7.12).

## 7.4 Total error and its spatial discretization, temporal discretization, linearization, and algebraic components

Suppose that Assumption 7.2 holds, so that we dispose of a weak solution. Then, we describe the distance of the current approximate solution to the exact solution by the following intrinsic error measure:

$$ \mathcal{N}^{n,k,i} := \left\{ \sum_{c \in \mathcal{C}} (\mathcal{N}_c^{n,k,i})^2 \right\}^{\frac{1}{2}} + \left\{ \sum_{p \in \mathcal{P}} (\mathcal{N}_p^{n,k,i})^2 \right\}^{\frac{1}{2}}, \qquad (7.17a) $$

where

$$ \mathcal{N}_c^{n,k,i} := \sup_{\varphi \in X^n, \|\varphi\|_{X^n} = 1} \int_{I^n} \left\{ (q_c, \varphi) - (\partial_t l_{c,h\tau}^{n,k,i}, \varphi) + (\theta_{c,h\tau}^{n,k,i}, \nabla \varphi) \right\} dt \qquad (7.17b) $$

is the dual norm of the residual of the weak formulation for each component flux $\theta_{c,h\tau}^{n,k,i}$ on the time interval $I^n$ and where

$$ \mathcal{N}_p^{n,k,i} := \inf_{\zeta_p \in X^n} \left\{ \sum_{c \in \mathcal{C}_p} \int_{I^n} \left\{ \sum_{K \in \mathcal{T}_H^n} (\nu_{p,K}^{n,k,i} C_{p,c,K}^{n,k,i})^2 \|\mathbf{u}_{p,h\tau}^{n,k,i} + \underline{\mathbf{K}} \nabla \zeta_p\|_{\underline{\mathbf{K}}^{-\frac{1}{2}, K}}^2 \right\} dt \right\}^{\frac{1}{2}} \qquad (7.17c) $$

evaluates the nonconformity (the distance of $\mathbf{u}_{p,h\tau}^{n,k,i}$ to $\underline{\mathbf{K}}\nabla H_0^1(\Omega)$) for each Darcy phase flux $\mathbf{u}_{p,h\tau}^{n,k,i}$. Indeed, from (7.16b), the first term should be zero and from (7.7), $\mathbf{u}_p(P_p)$ should be minus a gradient of a scalar-valued function, which is not necessarily the case for the discrete $\mathbf{u}_{p,h\tau}^{n,k,i}$ given at time $t^n$.

**Remark 7.4** (Error measure $\mathcal{N}^{n,k,i}$). The error measure $\mathcal{N}^{n,k,i}$ is a natural extension of the energy error from the simpler cases. Indeed, component by component and phase by phase, it extends Theorem 4.9 from the Poisson case as well as Remark 6.6 from the steady nonlinear case. The sum $\sum_{n=1}^N (\mathcal{N}^{n,k,i})^2$ (if the initial condition is satisfied exactly) then corresponds to the square of the energy error for the heat equation, see [120, Theorem 2.1 and equation (2.7)].

In (7.17), we have denoted the approximate component fluxes by $\theta_{c,h\tau}^{n,k,i}$ and the approximate Darcy phase fluxes by $\mathbf{u}_{p,h\tau}^{n,k,i}$, anticipating that our numerical approximation will here consist in temporal discretization (with characteristic time step $\tau$ and discrete times $t^n$), spatial discretization (with characteristic mesh size $h$), iterative linearization (with iteration index $k$), and iterative algebraic resolution (with iteration index $i$). Our goal will be to distinguish the error components as in (2.14) in Section 2.4.6, to design an a posteriori error estimate of the form (2.23b), and to use them in the fully adaptive Algorithm 7.1.

## 7.5 Backward Euler time discretization and locally conservative space discretization on polytopal meshes

To discretize problem (7.1)–(7.12), we consider a fully implicit backward Euler time stepping and an abstract locally conservative space discretization on a polytopal mesh of the form of Assumption 5.1, using phase-upwinding. For all $1 \le n \le N$, we let $\mathbf{X}_{\mathcal{T}_H}^n := (\mathbf{X}_K^n)_{K \in \mathcal{T}_H^n}$, with, following (7.2),

$$ \mathbf{X}_K^n := (P_K^n, (S_{p,K}^n)_{p \in \mathcal{P}}, (C_{p,c,K}^n)_{p \in \mathcal{P}, c \in \mathcal{C}_p}) \qquad (7.18) $$

be the algebraic vector of discrete unknowns on an element $K$. System (7.1)–(7.12) is then discretized as follows: for all time steps $1 \le n \le N$, all polytopal cells $K \in \mathcal{T}_H^n$, and each component $c \in \mathcal{C}$, we require

$$ \frac{|K|}{\tau^n} (l_{c,K}(\mathbf{X}_K^n) - l_{c,K}(\mathbf{X}_K^{n-1})) + \sum_{\sigma \in \mathcal{F}_K \cap \mathcal{F}_H^{\mathrm{int}}} \theta_{c,K,\sigma,\uparrow}(\mathbf{X}_{\mathcal{T}_H}^n) = |K|q_{c,K}^n. \qquad (7.19) $$
This equation expresses the mass balance for the element K. Here $q_{c,K}^n$ is the value of the source or sink $q_c$ on the element K and time interval $I^n$ – recall that we consider $q$ piecewise constant. Recall also that $|K|$ stands for the volume of the element K.

In a discrete version of (7.4), we set

$$l_{c,K}(\mathbf{X}_K^n) := \phi_K \sum_{p \in \mathcal{P}_c} \zeta_{p,K}^n S_{p,K}^n C_{p,c,K}^n. \qquad (7.20)$$

Similarly, for each component $c \in \mathcal{C}$, following (7.5), the total flux across the face $\sigma$ results from the sum of the corresponding fluxes for each phase $p \in \mathcal{P}_c$, i.e., for all elements $K \in \mathcal{T}_H^n$ and all faces $\sigma \in \mathcal{F}_K \cap \mathcal{F}_H^{\text{int}}$ with $\sigma \subset \partial K \cap \partial L$,

$$\theta_{c,K,\sigma,\uparrow}(\mathbf{X}_{\mathcal{T}_H}^n) := \sum_{p \in \mathcal{P}_c} \nu_{p,K_p^\uparrow}^n C_{p,c,K_p^\uparrow}^n v_{p,K,\sigma}(\mathbf{X}_{\mathcal{T}_H}^n), \quad K_p^\uparrow := \begin{cases} K & \text{if } v_{p,K,\sigma}(\mathbf{X}_{\mathcal{T}_H}^n) \ge 0, \\ L & \text{otherwise,} \end{cases} \qquad (7.21)$$

with $C_{p,c,K_p^\uparrow}^n$ and $\nu_{p,K_p^\uparrow}^n$ denoting, respectively, the upstream molar fraction and the upstream mobility. Let us also introduce

$$\theta_{c,K,\sigma}(\mathbf{X}_{\mathcal{T}_H}^n) := \sum_{p \in \mathcal{P}_c} \nu_{p,K}^n C_{p,c,K}^n v_{p,K,\sigma}(\mathbf{X}_{\mathcal{T}_H}^n), \qquad (7.22)$$

where no upwinding is used. In (7.21)–(7.22), we have employed a finite volume approximation of the normal component of the average phase velocity on face $\sigma$ given by

$$v_{p,K,\sigma}(\mathbf{X}_{\mathcal{T}_H}^n) := v_{p,K,\sigma}(\{\mathbf{X}_{K'}^n\}_{K' \in \mathcal{S}_\sigma}) := \sum_{K' \in \mathcal{S}_\sigma} \tau_{K'}^\sigma (P_{p,K'}^n + \rho_{p,\sigma}^n g z_{K'}), \qquad (7.23)$$

where, following (7.1),

$$P_{p,K}^n := P_K^n + P_{c_p}(S_{p,K}^n), \qquad (7.24)$$

and, for all $\sigma \in \mathcal{F}_K \cap \mathcal{F}_H^{\text{int}}$, $S_\sigma$ is the flux stencil collecting the elements in $\mathcal{T}_H^n$ with nonzero flux contribution $(S_\sigma = \{K, L\}$ for $\sigma$ such that $\sigma = \sigma_{K,L} = \partial K \cap \partial L \in \mathcal{F}_H^{\text{int}})$ in the two-point scheme (4.19b). Moreover, for all $L \in S_\sigma$, $\tau_L^\sigma \in \mathbb{R}$ is the transmissibility coefficient of the face $\sigma$ (if $\underline{K} = \underline{I}$, $\tau_K^\sigma = \frac{|\sigma_{K,L}|}{d_{K,L}}$ and $\tau_L^\sigma = -\frac{|\sigma_{K,L}|}{d_{K,L}}$ for the two-point scheme (4.19b)). Please remark that the boundary fluxes are set to zero to account for the homogeneous no-flow boundary conditions (7.8).

Finally, the initial condition comes from (7.9), and the algebraic constraints (7.10)–(7.12) are trivially discretized separately on all $K \in \mathcal{T}_H^n$, namely

$$\sum_{p \in \mathcal{P}} S_{p,K}^n = 1, \qquad (7.25)$$

$$\sum_{c \in \mathcal{C}_p} C_{p,c,K}^n = 1 \quad \forall p \in \mathcal{P}. \qquad (7.26)$$

## 7.6 Iterative linearization

At this stage, we need to solve, at each time step, the system of nonlinear algebraic equations resulting from the discretization (7.19)–(7.26). To this purpose, for all times $1 \le n \le N$, we apply an iterative linearization, in extension of (6.12). We consider the Newton linearization, as in (6.13). This generates, for an initial guess $\mathbf{X}_{\mathcal{T}_H}^{n,0}$ (typically given by the previous time step/initial condition), a sequence $(\mathbf{X}_{\mathcal{T}_H}^{n,k})_{k \ge 1}$ with $\mathbf{X}_{\mathcal{T}_H}^{n,k}$ solution to the following system of linear algebraic equations: for all components $c \in \mathcal{C}$ and all mesh elements $K \in \mathcal{T}_H^n$,

$$\sum_{K' \in \mathcal{T}_H^n} \frac{\partial R_{c,K}}{\partial \mathbf{X}_{K'}^n} (\mathbf{X}_{\mathcal{T}_H}^{n,k-1}) \cdot (\mathbf{X}_{K'}^{n,k} - \mathbf{X}_{K'}^{n,k-1}) + R_{c,K} (\mathbf{X}_{\mathcal{T}_H}^{n,k-1}) = 0, \qquad (7.27)$$

with, for all $c \in \mathcal{C}$ and all $K \in \mathcal{T}_H^n$,

$$R_{c,K}(\mathbf{X}_{\mathcal{T}_H}^n) := \frac{|K|}{\tau^n} (l_{c,K}(\mathbf{X}_K^n) - l_{c,K}(\mathbf{X}_K^{n-1})) + \sum_{\sigma \in \mathcal{F}_K \cap \mathcal{F}_H^{\text{int}}} \theta_{c,K,\sigma,\uparrow}(\mathbf{X}_{\mathcal{T}_H}^n) - |K|q_{c,K}^n.$$
Developing the derivatives in (7.27), we can write

$$ \frac{|K|}{\tau^n} (l_{c,K}^{n,k} - l_{c,K}(\mathbf{X}_K^{n-1})) + \sum_{\sigma \in \mathcal{F}_K \cap \mathcal{F}_H^{\text{int}}} \theta_{c,K,\sigma,\uparrow}^{n,k} - |K|q_{c,K}^n = 0, \quad (7.28) $$

where the linearization of the amount of component $l_{c,K}(\mathbf{X}_K^{n,k})$ of (7.20) is given by

$$ l_{c,K}^{n,k} := l_{c,K}(\mathbf{X}_K^{n,k-1}) + \frac{\partial l_{c,K}}{\partial \mathbf{X}_K^n}(\mathbf{X}_K^{n,k-1}) \cdot (\mathbf{X}_K^{n,k} - \mathbf{X}_K^{n,k-1}), $$

and where $\theta_{c,K,\sigma,\uparrow}^{n,k}$ are appropriate linearizations of the component fluxes $\theta_{c,K,\sigma,\uparrow}(\mathbf{X}_K^{n,k})$ of (7.21).

**Remark 7.5** (Blending Newton and fixed-point linearizations). In numerical experiments below in Section 7.11, we will actually define $\theta_{c,K,\sigma,\uparrow}^{n,k}$ by evaluation of the mobility $\nu_{p,K_p^\uparrow}^n$ and the molar fraction $C_{p,c,K_p^\uparrow}^n$ at the previous linearization iterate $k-1$ as

$$ \theta_{c,K,\sigma,\uparrow}^{n,k} := \sum_{p \in \mathcal{P}_c} \left\{ \nu_{p,K_p^\uparrow}^{n,k-1} C_{p,c,K_p^\uparrow}^{n,k-1} \left[ v_{p,K,\sigma}(\mathbf{X}_{T_H}^{n,k-1}) + \sum_{K' \in \mathcal{T}_H^n} \frac{\partial v_{p,K,\sigma}}{\partial \mathbf{X}_{K'}^n}(\mathbf{X}_{T_H}^{n,k-1}) \cdot (\mathbf{X}_{K'}^{n,k} - \mathbf{X}_{K'}^{n,k-1}) \right] \right\}, $$

blending thus the Newton and the fixed-point linearizations.

## 7.7 Iterative algebraic system solution

Now, for a given time step $1 \le n \le N$ and a given linearization iteration $k \ge 1$, (7.28) is a system of linear algebraic equations. To find its approximate solution, we use an iterative algebraic solver. Let $\mathbf{X}_{T_H}^{n,k,0}$ be an initial guess. Typically, we take $\mathbf{X}_{T_H}^{n,k-1} = \mathbf{X}_{T_H}^{n,k-1}$. This gives a sequence $(\mathbf{X}_{T_H}^{n,k,i})_{i \ge 1}$ solving (7.28) up to the residuals $R_{c,K}^{n,k,i}$ for all $c \in \mathcal{C}$ and all $K \in \mathcal{T}_H^n$, i.e.,

$$ \frac{|K|}{\tau^n} (l_{c,K}^{n,k,i} - l_{c,K}(\mathbf{X}_K^{n-1})) + \sum_{\sigma \in \mathcal{F}_K \cap \mathcal{F}_H^{\text{int}}} \theta_{c,K,\sigma,\uparrow}^{n,k,i} - |K|q_{c,K}^n =: R_{c,K}^{n,k,i}, \quad (7.29) $$

cf. (6.14)–(6.15). Here

$$ l_{c,K}^{n,k,i} := l_{c,K}(\mathbf{X}_K^{n,k-1}) + \frac{\partial l_{c,K}}{\partial \mathbf{X}_K^n}(\mathbf{X}_K^{n,k-1}) \cdot (\mathbf{X}_K^{n,k,i} - \mathbf{X}_K^{n,k-1}), \quad (7.30) $$

and, in the case of Remark 7.5,

$$ \theta_{c,K,\sigma,\uparrow}^{n,k,i} := \sum_{p \in \mathcal{P}_c} \left\{ \nu_{p,K_p^\uparrow}^{n,k-1} C_{p,c,K_p^\uparrow}^{n,k-1} \left[ v_{p,K,\sigma}(\mathbf{X}_{T_H}^{n,k-1}) + \sum_{K' \in \mathcal{T}_H^n} \frac{\partial v_{p,K,\sigma}}{\partial \mathbf{X}_{K'}^n}(\mathbf{X}_{T_H}^{n,k-1}) \cdot (\mathbf{X}_{K'}^{n,k,i} - \mathbf{X}_{K'}^{n,k-1}) \right] \right\}. \quad (7.31) $$

We will also consider $j \ge 1$ additional algebraic solver steps in (7.29), similarly as we have done it in (6.16).

## 7.8 Face normal fluxes and potential point values

We proceed here as in Section 6.5. In particular, in extension of (6.19), we rely on the finite volume face normal fluxes. These employ the approximate phase normal velocities (7.23), the total normal flux nonlinear expressions (7.22) and their upwindings (7.21), as well as the linearizations (7.31) of the latter. For all phases $p \in \mathcal{P}$, components $c \in \mathcal{C}$, time levels $n \ge 1$, mesh elements $K \in \mathcal{T}_H^n$ and their faces not lying on the boundary of $\Omega$, $\sigma \in \mathcal{F}_K \cap \mathcal{F}_H^{\text{int}}$, linearization steps $k \ge 1$, and algebraic iterations $i \ge 1$, at times $t \in I^n$ where necessary, we define

$$ (U_{K,p}^{t,n,k,i})_{\sigma} := \frac{t-t^{n-1}}{\tau^n} \sum_{K' \in \mathcal{S}_{\sigma}} \tau_{K'}^{\sigma} P_{p,K'}^{n,k,i} + \frac{t^n-t}{\tau^n} \sum_{K' \in \mathcal{S}_{\sigma}} \tau_{K'}^{\sigma} P_{p,K'}^{n-1}, \quad (7.32a) $$

$$ (\Theta_{\text{tm},K,c}^{t,n,k,i})_{\sigma} := \frac{t^n-t}{\tau^n} (\theta_{c,K,\sigma}(\mathbf{X}_{T_H}^{n,k,i}) - \theta_{c,K,\sigma}(\mathbf{X}_{T_H}^{n-1})), \quad (7.32b) $$

$$ (\Theta_{\text{upw},K,c}^{n,k,i})_{\sigma} := \theta_{c,K,\sigma,\uparrow}(\mathbf{X}_{T_H}^{n,k,i}) - \theta_{c,K,\sigma}(\mathbf{X}_{T_H}^{n,k,i}), \quad (7.32c) $$

$$ (\Theta_{\text{lin},K,c}^{n,k,i})_{\sigma} := \theta_{c,K,\sigma,\uparrow}^{n,k,i} - \theta_{c,K,\sigma,\uparrow}(\mathbf{X}_{T_H}^{n,k,i}), \quad (7.32d) $$

$$ (\Theta_{\text{alg},K,c}^{n,k,i})_{\sigma} := \theta_{c,K,\sigma,\uparrow}^{n,k,i+j} - \theta_{c,K,\sigma,\uparrow}^{n,k,i}, \quad (7.32e) $$
All these face normal fluxes are set to zero on faces $\sigma$ located on the boundary of $\Omega$, in accordance with the no-flow boundary condition (7.8), with the exception of $(U_{K,p}^{t,n,k,i})_{\sigma}$ that is set to $(\underline{K}\rho_g g \nabla z \cdot \mathbf{n}_{\Omega}, 1)_{\sigma}$, in accordance with (7.6)–(7.7). Note that (7.32a) is an affine-in-time interpolation of the previous time step $P_{p,K}^{n-1}$ and the current $P_{p,K}^{n,k,i}$ values.

We also need to define the phase pressure point values $Z_{p,K}^{t,n,k,i}$ and $Z_{p,K}^{\text{ext},t,n,k,i}$. We again proceed as in Section 6.5, using directly the cells pressure values in Definition 5.11 or 5.13. More precisely, the starting values are the phase pressure values obtained from the finite volume scheme, the current $P_{p,K}^{n,k,i}$ and the previous time step $P_{p,K}^{n-1}$, that we combine for each time $t \in I^n$ to get an affine-in-time interpolation

$$ \frac{t - t^{n-1}}{\tau^n} P_{p,K}^{n,k,i} + \frac{t^n - t}{\tau^n} P_{p,K}^{n-1}. \quad (7.33) $$

Note that all quantities of (7.32) and (7.33) are readily available from the finite volume discretization of Section 7.5, more precisely on each linearization step $k \ge 1$ and each algebraic solver step $i \ge 1$ as described in Sections 7.6–7.7.

## 7.9 Fictitious flux and potential reconstructions

We now identify the fictitious flux and potential reconstructions, following Section 6.6.

Let a phase $p \in \mathcal{P}$, a component $c \in \mathcal{C}$, a time step $1 \le n \le N$, a mesh element $K \in \mathcal{T}_H^n$, a linearization iteration $k \ge 1$, and an algebraic solver iteration $i \ge 1$ be fixed. For the fictitious flux reconstructions, we lift the face normal fluxes from (7.32), using them to define the spaces $\mathbf{V}_{h,N}^K$ in (5.16) of Definition 5.5. We define:

* The Darcy phase fluxes reconstructions $u_{p,h}^{n,k,i}|_K \in \mathcal{R}\mathcal{T}_0(\mathcal{T}_K) \cap \mathbf{H}(\text{div}, K)$, as a discrete counterpart of (7.7) at the discrete time $t^n$, using the face fluxes $U_{K,p}^{t,n,k,i}$ of (7.32a) at $t^n$.

* The component fluxes reconstructions $\theta_{c,h}^{n,k,i}|_K \in \mathcal{R}\mathcal{T}_0(\mathcal{T}_K) \cap \mathbf{H}(\text{div}, K)$, as a discrete counterpart of (7.5) at the discrete time $t^n$, prescribed by $\theta_{c,K,\sigma}(\mathbf{X}_{\mathcal{T}_H}^{n,k,i})$ appearing in (7.32b).

* The upwinding error component fluxes reconstructions $\theta_{\text{upw},c,h}^{n,k,i}|_K \in \mathcal{R}\mathcal{T}_0(\mathcal{T}_K) \cap \mathbf{H}(\text{div}, K)$ prescribed by $\Theta_{\text{upw},K,c}^{n,k,i}$ of (7.32c).

* The linearization error flux reconstruction $\theta_{\text{lin},c,h}^{n,k,i}|_K \in \mathcal{R}\mathcal{T}_0(\mathcal{T}_K) \cap \mathbf{H}(\text{div}, K)$ prescribed by $\Theta_{\text{lin},K,c}^{n,k,i}$ of (7.32d).

* The algebraic error flux reconstruction $\theta_{\text{alg},c,h}^{n,k,i}|_K \in \mathcal{R}\mathcal{T}_0(\mathcal{T}_K) \cap \mathbf{H}(\text{div}, K)$ prescribed by $\Theta_{\text{alg},K,c}^{n,k,i}$ of (7.32e).

The error measure $\mathcal{N}^{n,k,i}$ of (7.17) works with approximate Darcy phase fluxes $u_{p,h\tau}^{n,k,i}$ together with approximate component fluxes $\theta_{c,h\tau}^{n,k,i}$. We define them as affine-in-time on the time interval $I^n$, given respectively by $u_{p,h}^{n,k,i}$ and $\theta_{c,h}^{n,k,i}$ at the current time $t^n$, and by $u_{p,h}^{n-1}$ and $\theta_{c,h}^{n-1}$ obtained in the same way but on time previous time $t^{n-1}$, once the linearization and algebraic iterations have been stopped. All $u_{p,h\tau}^{n,k,i}$, $\theta_{c,h\tau}^{n,k,i}$, and the other fluxes $\theta_{\text{upw},c,h}^{n,k,i}$, $\theta_{\text{lin},c,h}^{n,k,i}$, and $\theta_{\text{alg},c,h}^{n,k,i}$ are only used for the proof of Theorem 7.6 below and need not be constructed in practice (whence, recall, the naming fictitious). It namely follows from (7.29), on algebraic step $i + j$, the above flux reconstructions, and the Green theorem that there holds

$$ q_{c,K}^n - \frac{l_{c,K}^{n,k,i+j} - l_{c,K}(\mathbf{X}_K^{n-1})}{\tau^n} - (\nabla \cdot (\theta_{c,h}^{n,k,i} + \theta_{\text{upw},c,h}^{n,k,i} + \theta_{\text{lin},c,h}^{n,k,i} + \theta_{\text{alg},c,h}^{n,k,i}))|_K = -|K|^{-1} R_{c,K}^{n,k,i+j} \quad \forall K \in \mathcal{T}_H^n, \quad (7.34) $$

similarly as in (6.22).

The fictitious reconstruction of the phase pressures is again achieved as in Definition 5.14, i.e., setting

$$ \zeta_{p,h\tau}^{n,k,i}(\mathbf{a}) := Z_{p,\mathbf{a}}^{t,n,k,i} \quad \forall \mathbf{a} \in \mathcal{V}_h $$

from the potential point values identified in Section 7.8. It is piecewise affine with respect to the simplicial submeshes $\mathcal{T}_h^n$ and $H^1(\Omega)$-conforming and piecewise affine in time. In the error measure $\mathcal{N}^{n,k,i}$ of (7.17),
we also need $l_{c,h\tau}^{n,k,i}$. It is defined as affine-in-time on each time interval $I^n$, given respectively by the values $l_{c,K}(\mathbf{X}_K^{n-1})$ and $l_{c,K}(\mathbf{X}_K^{n,k,i})$ at times $t^{n-1}$ and $t^n$, for all $K \in \mathcal{T}_H^n$. This in particular gives

$$ \partial_t l_{c,h\tau}^{n,k,i}|_K = \frac{l_{c,K}(\mathbf{X}_K^{n,k,i}) - l_{c,K}(\mathbf{X}_K^{n-1})}{\tau^n} \quad (7.35) $$

on each time interval $I^n$.

## 7.10 A guaranteed a posteriori error estimate distinguishing the error components with inexpensive implementation and evaluation

Proceeding as for the steady linear problem in Section 5.9 and the steady nonlinear problem in Section 6.7 (under Assumption 5.1), we obtain a posteriori error estimates for the problem (7.1)-(7.12). We first treat the simpler case where there is no coarsening of the space meshes between the time steps and postpone the general case to Remark 7.8 below.

**Theorem 7.6** (A guaranteed a posteriori error estimate distinguishing the error components with inexpensive implementation and evaluation).

For the multiphase compositional Darcy flow (7.1)-(7.12), let the weak solution $\mathbf{X}$ satisfy Assumption 7.2. Consider any polytopal discretization of the form (7.19), an iterative linearization (7.27) on step $k \ge 1$, and any iterative algebraic solver (7.29) on step $i \ge 1$, for a given time step $1 \le n \le N$. Consider $j \ge 1$ additional algebraic solver steps. Let, for each polytopal element $K \in \mathcal{T}_H^n$, the element vectors of the face normal fluxes $\Theta_{\bullet,K,c}^{n,k,i}$, $\bullet$ = upw, lin, alg, together with $\Theta_{\mathrm{tm},K,c}^{t,n,k,i}$ and $\mathcal{U}_{K,p}^{t,n,k,i}$ be given by (7.32). Let the element vectors of potential point values $Z_{p,K}^{t,n,k,i}$ and $Z_{p,K}^{\mathrm{ext},t,n,k,i}$ be constructed from the time-interpolated values (7.33) of $P_{p,K}^{n,k,i}$ and $P_{p,K}^{n-1}$ following Section 7.8. Consider the approximate solution $\mathbf{u}_{p,h\tau}^{n,k,i}$, $\theta_{c,h\tau}^{n,k,i}$ as defined in Section 7.9 and the intrinsic error measure $\mathcal{N}^{n,k,i}$ of (7.17). Let finally the element matrices $\hat{\mathbf{A}}_{\mathrm{MFE},K}$, $\hat{\mathbf{S}}_{\mathrm{FE},K}$, and $\hat{\mathbf{M}}_{\mathrm{FE},K}$ be respectively defined by (5.20), (5.32), and (5.33). Let $\mathcal{T}_H^n \subset \mathcal{T}_H^{n-1}$, i.e., there is no coarsening. Then

$$ \mathcal{N}^{n,k,i} \le \left\{ \sum_{c \in \mathcal{C}} \left( \eta_{\mathrm{sp},c}^{n,k,i} + \eta_{\mathrm{tm},c}^{n,k,i} + \eta_{\mathrm{lin},c}^{n,k,i} + \eta_{\mathrm{alg},c}^{n,k,i} + \eta_{\mathrm{rem},c}^{n,k,i} \right)^2 \right\}^{\frac{1}{2}} \quad (7.36) $$

with

$$ \eta_{\bullet,c}^{n,k,i} := \left\{ \delta_{\bullet} \int_{I^n} \sum_{K \in \mathcal{T}_H^n} \left( \eta_{\bullet,K,c}^{n,k,i} \right)^2 dt \right\}^{\frac{1}{2}}, \quad (7.37) $$

where $\bullet$ = sp, tm, lin, alg, rem and $\delta_{\bullet} := 2$, except for $\delta_{\mathrm{sp}} := 4$. Here, for $c \in \mathcal{C}$, we prescribe the elementwise spatial estimators

$$ \eta_{\mathrm{sp},K,c}^{n,k,i} := \eta_{\mathrm{upw},K,c}^{n,k,i} + \left\{ \sum_{p \in \mathcal{P}_c} \left( \eta_{\mathrm{NC},K,c,p}^{t,n,k,i} \right)^2 \right\}^{\frac{1}{2}}, \quad (7.38a) $$

with

$$ \left( \eta_{\mathrm{NC},K,c,p}^{t,n,k,i} \right)^2 := \left( \nu_{p,K}^{n,k,i} C_{p,c,K}^{n,k,i} \right)^2 \left[ \left( \mathcal{U}_{K,p}^{t,n,k,i} \right)^{\mathrm{t}} \hat{\mathbf{A}}_{\mathrm{MFE},K} \mathcal{U}_{K,p}^{t,n,k,i} + \left( \mathcal{Z}_{p,K}^{t,n,k,i} \right)^{\mathrm{t}} \hat{\mathbf{S}}_{\mathrm{FE},K} \mathcal{Z}_{p,K}^{t,n,k,i} \right. \\ \left. + 2 \left( \mathcal{U}_{K,p}^{t,n,k,i} \right)^{\mathrm{t}} \mathcal{Z}_{p,K}^{\mathrm{ext},t,n,k,i} - 2 \sum_{\sigma \in \mathcal{F}_K} \left( \mathcal{U}_{K,p}^{t,n,k,i} \right)_{\sigma} |K|^{-1} \mathcal{M}_{\mathrm{FE},K} \mathcal{Z}_{p,K}^{t,n,k,i} \right], \quad (7.38b) $$

and the upwinding estimators

$$ \left( \eta_{\mathrm{upw},K,c}^{n,k,i} \right)^2 := \left( \Theta_{\mathrm{upw},K,c}^{n,k,i} \right)^{\mathrm{t}} \hat{\mathbf{A}}_{\mathrm{MFE},K} \left( \Theta_{\mathrm{upw},K,c}^{n,k,i} \right) \quad (7.38c) $$

together with the temporal estimators

$$ \left( \eta_{\mathrm{tm},K,c}^{n,k,i} \right)^2 := \left( \Theta_{\mathrm{tm},K,c}^{t,n,k,i} \right)^{\mathrm{t}} \hat{\mathbf{A}}_{\mathrm{MFE},K} \Theta_{\mathrm{tm},K,c}^{t,n,k,i} \quad (7.38d) $$
the linearization estimators

$$ \eta_{\text{lin},K,c}^{n,k,i} := \left( (\Theta_{\text{lin},K,c}^{n,k,i})^{\text{t}} \hat{\Lambda}_{\text{MFE},K} \Theta_{\text{lin},K,c}^{n,k,i} \right)^{1/2} + h_K c_{\underline{K},K}^{-1/2} (\tau^n)^{-1} \|l_{c,K}^{n,k,i} - l_{c,K}(\mathbf{X}_K^{n,k,i})\|_K, \quad (7.38e) $$

the algebraic estimators

$$ \eta_{\text{alg},K,c}^{n,k,i} := \left( (\Theta_{\text{alg},K,c}^{n,k,i})^{\text{t}} \hat{\Lambda}_{\text{MFE},K} \Theta_{\text{alg},K,c}^{n,k,i} \right)^{1/2} + h_K c_{\underline{K},K}^{-1/2} (\tau^n)^{-1} \|l_{c,K}^{n,k,i+j} - l_{c,K}^{n,k,i}\|_K, \quad (7.38f) $$

and the algebraic remainder estimators by

$$ \eta_{\text{rem},K,c}^{n,k,i} := h_K c_{\underline{K},K}^{-1/2} |K|^{-1/2} |R_{c,K}^{n,k,i+j}|. \quad (7.38g) $$

*Proof.* Let the time step $1 \le n \le N$, the iterative linearization step $k \ge 1$, and the iterative algebraic solver $i \ge 1$ be fixed. We bound the two building terms $N_c^{n,k,i}$ and $N_p^{n,k,i}$ of the error measure $N^{n,k,i}$ from (7.17) separately, for each component $c \in C$ and each phase $p \in \mathcal{P}$. We rely on the fictitious flux and potential reconstructions of Section 7.9 and namely on (7.34).

We first bound $N_c^{n,k,i}$. Let $\varphi \in X^n$ with $\|\varphi\|_{X^n} = 1$ be fixed. We have, using (7.35), (7.34), and the Green theorem

$$ \begin{aligned} & (q_c, \varphi) - (\partial_t l_{c,h\tau}^{n,k,i}, \varphi) + (\theta_{c,h\tau}^{n,k,i}, \nabla\varphi) \\ &= \sum_{K \in \mathcal{T}_H^n} \left( q_c - \frac{l_{c,K}^{n,k,i+j} - l_{c,K}(\mathbf{X}_K^{n-1})}{\tau^n} - \nabla \cdot (\theta_{c,h}^{n,k,i} + \theta_{\text{upw},c,h}^{n,k,i} + \theta_{\text{lin},c,h}^{n,k,i} + \theta_{\text{alg},c,h}^{n,k,i}), \varphi \right)_K \\ &\quad + \sum_{K \in \mathcal{T}_H^n} \left( \frac{l_{c,K}^{n,k,i+j} - l_{c,K}(\mathbf{X}_K^{n-1})}{\tau^n} - \frac{l_{c,K}(\mathbf{X}_K^{n,k,i}) - l_{c,K}(\mathbf{X}_K^{n-1})}{\tau^n}, \varphi \right)_K \\ &\quad + (\nabla \cdot (\theta_{c,h}^{n,k,i} + \theta_{\text{upw},c,h}^{n,k,i} + \theta_{\text{lin},c,h}^{n,k,i} + \theta_{\text{alg},c,h}^{n,k,i}), \varphi) + (\theta_{c,h\tau}^{n,k,i}, \nabla\varphi) \\ &= - \sum_{K \in \mathcal{T}_H^n} (|K|^{-1} R_{c,K}^{n,k,i+j}, \varphi)_K + \sum_{K \in \mathcal{T}_H^n} \left( \frac{l_{c,K}^{n,k,i+j} - l_{c,K}(\mathbf{X}_K^{n,k,i})}{\tau^n}, \varphi \right)_K \\ &\quad + (\theta_{c,h\tau}^{n,k,i} - \theta_{c,h}^{n,k,i} - \theta_{\text{upw},c,h}^{n,k,i} - \theta_{\text{lin},c,h}^{n,k,i} - \theta_{\text{alg},c,h}^{n,k,i}, \nabla\varphi) \end{aligned} $$

We now treat the terms herein separately. First, the triangle inequality gives

$$ \begin{aligned} & (\theta_{c,h\tau}^{n,k,i} - \theta_{c,h}^{n,k,i} - \theta_{\text{upw},c,h}^{n,k,i} - \theta_{\text{lin},c,h}^{n,k,i} - \theta_{\text{alg},c,h}^{n,k,i}, \nabla\varphi) \\ &\le \sum_{K \in \mathcal{T}_H^n} \left\{ \|\theta_{c,h\tau}^{n,k,i} - \theta_{c,h}^{n,k,i}\|_{\underline{K}^{-1/2},K} + \|\theta_{\text{upw},c,h}^{n,k,i}\|_{\underline{K}^{-1/2},K} \right. \\ &\qquad \left. + \|\theta_{\text{lin},c,h}^{n,k,i}\|_{\underline{K}^{-1/2},K} + \|\theta_{\text{alg},c,h}^{n,k,i}\|_{\underline{K}^{-1/2},K} \right\} \|\underline{K}^{1/2} \nabla\varphi\|_K. \end{aligned} $$

Second, the Cauchy–Schwarz inequality gives

$$ \sum_{K \in \mathcal{T}_H^n} \left( \frac{l_{c,K}^{n,k,i+j} - l_{c,K}(\mathbf{X}_K^{n,k,i})}{\tau^n}, \varphi \right)_K \le \sum_{K \in \mathcal{T}_H^n} h_K c_{\underline{K},K}^{-1/2} (\tau^n)^{-1} \|l_{c,K}^{n,k,i+j} - l_{c,K}(\mathbf{X}_K^{n,k,i})\|_K h_K^{-1} c_{\underline{K},K}^{1/2} \|\varphi\|_K. $$

And third,

$$ - \sum_{K \in \mathcal{T}_H^n} (|K|^{-1} R_{c,K}^{n,k,i+j}, \varphi)_K \le \sum_{K \in \mathcal{T}_H^n} h_K c_{\underline{K},K}^{-1/2} |K|^{-1/2} |R_{c,K}^{n,k,i+j}| h_K^{-1} c_{\underline{K},K}^{1/2} \|\varphi\|_K. $$
Thus, regrouping these estimates and using the Cauchy–Schwarz inequality together with (7.14), we have

$$
\begin{align}
\mathcal{N}_c^{n,k,i} &\le \sup_{\varphi \in X^n, \|\varphi\|_{X^n} = 1} \int_{I^n} \sum_{K \in \mathcal{T}_H^n} \left( \eta_{\text{tm},K,c}^{n,k,i} + \eta_{\text{upw},K,c}^{n,k,i} + \eta_{\text{lin},K,c}^{n,k,i} + \eta_{\text{alg},K,c}^{n,k,i} + \eta_{\text{rem},K,c}^{n,k,i} \right. \nonumber \\
&\qquad \left. \{h_K^{-2} c_{\underline{K},K} \|\varphi\|_K^2 + \|\underline{K}^{\frac{1}{2}} \nabla \varphi\|_K^2\}^{\frac{1}{2}} dt \right. \nonumber \\
&\le \sup_{\varphi \in X^n, \|\varphi\|_{X^n} = 1} \left\{ \int_{I^n} \sum_{K \in \mathcal{T}_H^n} \left( \eta_{\text{tm},K,c}^{n,k,i} + \eta_{\text{upw},K,c}^{n,k,i} + \eta_{\text{lin},K,c}^{n,k,i} + \eta_{\text{alg},K,c}^{n,k,i} + \eta_{\text{rem},K,c}^{n,k,i} \right)^2 dt \right\}^{\frac{1}{2}} \underbrace{\|\varphi\|_{X^n}}_{=1} \tag{7.39}
\end{align}
$$

where, for all components $c \in \mathcal{C}$ and all $K \in \mathcal{T}_H^n$,

$$
\eta_{\text{tm},K,c}^{n,k,i} := \| \theta_{c,h\tau}^{n,k,i} - \theta_{c,h}^{n,k,i} \|_{\underline{K}^{-\frac{1}{2}},K}, \qquad (7.40a)
$$

$$
\eta_{\text{upw},K,c}^{n,k,i} := \| \theta_{\text{upw},c,h}^{n,k,i} \|_{\underline{K}^{-\frac{1}{2}},K}, \qquad (7.40\text{b})
$$

$$
\eta_{\text{lin},K,c}^{n,k,i} := \| \theta_{\text{lin},c,h}^{n,k,i} \|_{\underline{K}^{-\frac{1}{2}},K} + h_K c_{\underline{K},K}^{-1/2} (\tau^n)^{-1} \| l_{c,K}^{n,k,i} - l_{c,K}(\mathbf{X}_K^{n,k,i}) \|_K, \quad (7.40c)
$$

$$
\eta_{\text{alg},K,c}^{n,k,i} := \| \theta_{\text{alg},c,h}^{n,k,i} \|_{\underline{K}^{-\frac{1}{2}},K} + h_K c_{\underline{K},K}^{-1/2} (\tau^n)^{-1} \| l_{c,K}^{n,k,i+j} - l_{c,K}^{n,k,i} \|_K, \quad (7.40\text{d})
$$

$$
\eta_{\text{rem},K,c}^{n,k,i} := h_K c_{\underline{K},K}^{-1/2} |K|^{-1/2} |R_{c,K}^{n,k,i+j}|. \quad (7.40e)
$$

We now bound $\mathcal{N}_p^{n,k,i}$. Let

$$
\eta_{\text{NC},K,c,p}^{t,n,k,i} := \nu_{p,K}^{n,k,i} C_{p,c,K}^{n,k,i} \|\mathbf{u}_{p,h\tau}^{n,k,i} + \underline{\mathbf{K}} \nabla \zeta_{p,h\tau}^{n,k,i}\|_{\underline{K}^{-1/2},K'} \quad (7.40\text{f})
$$

We have

$$
\mathcal{N}_p^{n,k,i} \le \left\{ \sum_{c \in C_p} \int_{I^n} \left\{ \sum_{K \in \mathcal{T}_H^n} (\eta_{\text{NC},K,c,p}^{t,n,k,i})^2 \right\} dt \right\}^{\frac{1}{2}} . \quad (7.41)
$$

We now remark that the estimators in (7.40a)–(7.40f) and (7.38b)–(7.38g) are indeed the same, pro-
ceeding as in Theorem 5.18 to evaluate the nonconformity estimator given by (7.40f) and by applying
Lemma 5.8 to evaluate the upwinding, temporal, linearization, and algebraic estimators. Regrouping the
spatial discretization estimators

$$
\eta_{\text{sp},K,c}^{n,k,i} := \eta_{\text{upw},K,c}^{n,k,i} + \left\{ \sum_{p \in \mathcal{P}_c} \left( \eta_{\text{NC},K,c,p}^{t,n,k,i} \right)^2 \right\}^{\frac{1}{2}} \quad (7.42)
$$

and all the estimators per component $c \in \mathcal{C}$, using the Cauchy–Schwarz inequality, and recalling (7.17),
we obtain (7.36). $\square$

**Remark 7.7 (Sharper estimates).** *Sharper estimates, namely avoiding the factors* $\delta_\bullet$, **(but not identifying/regrouping the error component estimators) immediately follow from (7.39) and (7.41).*

**Remark 7.8 (Coarsening).** If a mesh element $K \in \mathcal{T}_H^n$ has been obtained by coarsening of elements of the previous mesh $\mathcal{T}_H^{n-1}$ collected in the set $\mathcal{T}_{H,K}^{n-1}$, then the estimators $\eta_{\text{NC},K,c,p}^{t,n,k,i}$ of (7.38b) and $\eta_{\text{tm},K,c}^{n,k,i}$ of (7.38d) need to be evaluated using the matrices $\hat{A}_{\text{MFE},K}$, $\hat{S}_{\text{FE},K}$, and $\hat{M}_{\text{FE},K}$ stemming from a common refinement of the element K by the simplicial meshes $\mathcal{T}_h^{n-1}|_{\mathcal{T}_{H,K}^{n-1}}$ and $\mathcal{T}_K$ and not $\mathcal{T}_K$ solely. In this sense, Theorem 7.6 also holds when $\mathcal{T}_H^n \not\subset \mathcal{T}_H^{n-1}$.

7.11 Numerical experiments

The proposed a posteriori error estimate framework summarized in Theorem 7.6 and Algorithm 7.1 has been implemented in a reservoir prototype simulator [223], a thermal multi-purpose simulator written in C++, which is a part of the next generation IFPEn research simulators based on the Arcane framework [141]. The execution platform is a public computer Intel Core i7, 8 cores, 3.7Ghz with 16GB of memory.
### 7.11.1 Setup

We study two different test cases. The first one is taken from [89], relying on the tenth SPE comparative solution project model. It is an incompressible water-oil two-phase flow problem built on a Cartesian regular geometry. This test corresponds to the layer 85 of SPE10. We choose here first the Cartesian regular mesh so as to compare our new approach with approaches already validated on this type of meshes. The second case is a simulation of a black-oil model. For this test, the problem is built on a three-dimensional corner-point geometry (distorted grids), well-known and most often used in reservoir simulation due to the flexibility allowing a good representation of reservoir description, see [105] and the references therein.

In our computational experiments presented below, we perform a nonconforming-type mesh refinement and allow one hanging node per cell face for the space mesh adaptation. The resulting mesh is a polytopal mesh in the sense of Section 3.6. Coarsening is then obtained by retrogressive cell agglomeration. During the coarsening and refinement phases, we use some upscaling and interpolating mechanisms to ensure the mass balance. Geological data (rock properties like porosity and permeability) are stored at the finest mesh level as data on a cloud of space points. These data are used during mesh adaptation by standard upscaling algorithms to compute the values of the rock model properties on the newly created cells following [185].

### 7.11.2 Two-phase Darcy flow

We consider a two-dimensional spatial domain discretized by a grid of 60 × 220 rectangular cells of size 6.096m in the x direction and 3.048m in the y direction. We choose the initial time step as $\tau^0 = 4.32 \times 10^4$s, which equals to 0.5 days, and the process is simulated to $t_F = 2000$ days. The reservoir is initially saturated with hydrocarbons and we consider the injection of water by a well located at the center of the grid. Four production wells are placed at the four corners of the domain. Therefore, we have a water component W and an oil component O collected in the set of components $C = \{W, O\}$ and two phases $P = \{w, o\}$ corresponding to water and oil. The model is actually simplified in that the components can be identified with the phases, so that $P_W = \{w\}$ and $P_O = \{o\}$, $C_{w,W} = C_{o,O} = 1$, $C_{o,W} = C_{w,O} = 0$, and the vector of unknowns $X$ reduces to $(P, (S_p)_{p \in P})$. Note that then the accumulation term $l_c$ becomes linear in the only unknown $S_p$, and one particular consequence is that the second term in the definition (7.38e) of $\eta_{\text{lin},K,c}^{n,k,i}$ vanishes. The porosity $\phi$ and the permeability field $K$ (scalar coefficient times an identity matrix) are shown in Figure 30. The other parameters of Section 7.1 are given by (see [89] for more details):

* $\mu_o = 10^{-3}$ Pa·s and $\mu_w = 0.3 \times 10^{-3}$ Pa·s,

* $\zeta_o = \zeta_w = 1$ mole·m$^{-3}$,

* $S_{wr} = S_{or} = 0.2$,

* $k_{r,w}(S_w) = \left( \frac{S_w - S_{wr}}{1 - S_{wr} - S_{or}} \right)^2$ and $k_{r,o}(S_o) = \left( \frac{S_o - S_{or}}{1 - S_{wr} - S_{or}} \right)^2$,

* $P_{c_p}(S_p) = 0$,

* there is no gravitational force, $z = 0$, so that the mass densities $\rho_p$ need not be specified.

We consider a discretization by the implicit multi-point finite volume scheme of Section 7.5 with the linearization detailed in Section 7.6, using the blending of Newton and fixed-point of Remark 7.5. For the linear solver needed in Section 7.7, we use the Bi-Conjugated Gradient Stabilized (BiCGStab) [239] with an ILU{0} preconditioner.

Figure 31 shows, at 500 days of the simulation, the evolution of the approximate water saturation, the spatial estimator computed by the formula on Cartesian grid proposed in [262, Chapter 4], and the spatial estimator given by (7.38a). Similar behavior is observed for the two a posteriori estimators. Furthermore, we see that they both detect well the error following the saturation front despite the strong heterogeneity of the domain.

The previous result motivates an adaptive mesh refinement/coarsening (AMR) strategy based on these estimators. We now verify that the simple a posteriori estimators on polyhedral mesh of Theorem 7.6 give the expected results, while comparing it with the already validated case in [262, Chapter 4]. We
<figure><img src="image_33.png" /><figcaption>Figure 30: [Section 7.11.2] Porosity (**left**) and permeability (**right**), layer 85 of the 10th SPE case</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>Figure 31: [Section 7.11.2] Results on the fine mesh at 500 days</figcaption></figure>

focus on the mesh space adaptivity and we do not act on the time step. Actually, in our simulator, the time step size is controlled by the behavior of the Newton/fixed-point algorithm (increased systematically in general, by multiplying by 2, and divided by 2 if the Newton/fixed-point algorithm diverges), so the space-mesh adaptation only will enable an easy comparison with the non-adaptive code. We first apply Algorithm 7.1 with $\delta_{\text{ref}} = 0.7$, $\delta_{\text{deref}} = 0.2$, and “exact” (i.e., negligible-error) algebraic and linearization solvers. On the coarse scale, the domain is discretized by a grid of 30 × 110 cells and we allow one refinement level. Figure 32 shows the evolution of the approximate water saturation and of the meshes at two different simulation times. We remark that the refinement follows the saturation front as times evolves. Additionally, since we have a model with highly heterogeneous permeability, we are lead to perform a slow derefinement process in the zone abandoned by the front of water saturation. Figure 33
<figure><img src="image_35.png" /><figcaption>Figure 32: [Section 7.11.2] Results on the adaptive mesh at 400 days and 1100 days</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>Figure 33: [Section 7.11.2] Fine mesh, coarse mesh, and adaptive meshes: cumulated oil rate (**left**) and water cut (**right**)</figcaption></figure>

depicts the cumulated oil rate (**left**) and the water-cut<fnref n="2" /> (**right**) during the simulation. We compare there the results on the fine grid, the results on the coarse grid, and the results of the two AMR strategies. We remark that the accuracy of the results on the fine grid is almost recovered by the AMR strategy, and appears much better than the coarse-grid result.

The details of the efficiency of the AMR strategy based on the a posteriori error estimator can finally be appreciated in Table 1. We compare the global CPU time of the simulation for the different strategies. We detail the CPU time spent on the evaluation of the estimators and the mesh adaptation. We remark a low cost of the estimators evaluation compared with the total computation CPU time, thanks to the use of the simple and fast-to-evaluate form of a posteriori error estimates on polygonal meshes of Theorem 7.6. In Table 1, when applying the mesh adaptation, the overall CPU time is the sum of the resolution time, AMR time, and estimators evaluation time. We remark that applying the AMR strategy on this two-dimensional test case leads to a gain factor in the overall CPU time at around 2, for both the polygonal estimate proposed in this paper and the Cartesian estimate already validated on this type of meshes. Note that, however, the simulation is slightly faster when adapting using the Cartesian estimate, due to

<footnote n="2">²The ratio of water produced compared to the volume of total liquids produced.</footnote>
<table><caption>Table 1: [Section 7.11.2] Fine grid vs. adaptive mesh refinements</caption><thead><tr><th>-</th><th>Resolution</th><th>AMR</th><th>Estimators evaluation</th><th>Gain factor</th></tr></thead><tbody><tr><td>Fine</td><td>603s</td><td>-</td><td>-</td><td>-</td></tr><tr><td>AMR Cartesian est.</td><td>229s</td><td>39s</td><td>19s</td><td>2.1</td></tr><tr><td>AMR polygonal est.</td><td>242s</td><td>46s</td><td>27s</td><td>1.9</td></tr></tbody></table>

<figure><img src="image_37.png" /><figcaption>Figure 34: [Section 7.11.3] Permeability</figcaption></figure>

the fact that we can here directly compute the $RT$ basis functions on the rectangular cells.

### 7.11.3 Three-phases, three-components Darcy flow

In this section, we present a simulation of a black-oil model. Here, we have three phases constituted by water, oil, and gas, represented by lowercase letters w, o, g as indices, respectively. The oil phase contains two types of components: nonvolatile oil and volatile oil, which we call here oil component and gas component, respectively. This is due to the fact that in this model, the hydrocarbon components are divided into light and heavy components. The light component can dissolve into the liquid oil phase or volatilize in the gas phase according to the pressure and temperature. The gas phase only contains the gas components and the water phase only contains the water component. The components are represented by uppercase letters W for the water component, O for the oil component, and G for the gas component. Therefore, we have a problem with three phases $P = \{w, o, g\}$ and three components $C = \{W, O, G\}$.

#### Model setting

The reservoir considered in this test case is a 3-dimensional domain $\Omega := 4750\text{m} \times 3000\text{m} \times 114\text{m}$ discretized by a corner-point geometry grid. We consider a heterogeneous anisotropic reservoir with porosity 0.3 and permeability $K$ in a form of a diagonal matrix with the $x$ and $y$ components identical and forming three horizontal layers, and the $z$ component forming two horizontal layers, see Figure 34. We consider a gas injection in a reservoir initially unsaturated. A vertical gas injection well perforates a corner of the reservoir in the $z$ direction and a production well is located in the opposite corner. On the fine scale, the domain is discretized by a grid of $76 \times 48 \times 10$ elements and on the coarse scale by a grid of $38 \times 24 \times 5$ cells, leading to one refinement level. The process is simulated to $t_F = 2000$ days with initial time step $\tau^0 = 4.32 \times 10^4$ s, which equals to 0.5 days. Data, constraints, and pressure-volume properties are adapted from the first SPE comparative solution project model (SPE1) designed to simulate a three-dimensional black-oil reservoir, given in [201, Tables 1,2, and 3]. We consider a discretization by the implicit multi-point finite volume scheme of Section 7.5 with the blended Newton and fixed-point linearization detailed in Section 7.6. For the linear solver of Section 7.7, we again use the BiCGStab [239] with an ILU{0} preconditioner. Figure 35 shows the evolution of the gas saturation and of the spatial estimator at 1000 days. Note that, for the spatial estimator, the data are normalized by max value in order to have a [0, 1] range. We observe that the spatial estimator follows the saturation front though the heterogeneous anisotropic medium with time evolution.
<figure><img src="image_38.png" /><figcaption>Figure 35: [Section 7.11.3] Results on the fine mesh at 1000 days, gas saturation (left) and normalized polygonal estimate (right)</figcaption></figure>

<figure><img src="image_39.png" /><figcaption>Figure 36: [Section 7.11.3] Gas saturation: results on the adaptive mesh at 500 days (left) and 1500 days (right)</figcaption></figure>

## Adaptive space mesh refinement and stopping criteria for the linear solver

In the standard resolution of our reservoir prototype simulator, the chosen grid is the fine-scale one and the initial time step, as mentioned before, is chosen as $\tau^0 = 4.32 \times 10^4$ s. As in the test of Section 7.11.2, the time step is increased systematically, by multiplying by 2, but also controlled by the convergence of the linearization loop in such a way that we divide it by 2 if a divergence of the linearization algorithm occurs. We thus stick to this setting and only focus on the stopping criteria for the linear solver (and not for the linearization one) and on mesh adaptivity in space (and not in time). For the adaptive resolution, we start on the coarse-level grid allowing to one refinement level and we fix $\delta_{\text{ref}} = 0.7$, $\delta_{\text{deref}} = 0.2$ in Algorithm 7.1.

Figure 36 illustrates the evolution of the approximate gas saturation at two different time steps. We remark that the refinement follows the saturation front as times evolves. Additionally, the fact that the light component G (gas component) can dissolve into the liquid oil phase or volatilize in the gas phase, we are lead to perform some localized refinement in zones abandoned by the front of gas saturation.

We show in the left part of Figure 37, at a fixed time of 500 days and for the first linearization iteration, the evolution of the total estimator ($\eta_{\text{sp},c}^{n,k,i} + \eta_{\text{tm},c}^{n,k,i} + \eta_{\text{lin},c}^{n,k,i} + \eta_{\text{alg},c}^{n,k,i} + \eta_{\text{rem},c}^{n,k,i}$) with $c = \text{G}$ the
<figure><img src="image_40.png" /><figcaption>Figure 37: [Section 7.11.3] Standard resolution vs. adaptive resolution: total estimator and its algebraic component (left) and number of BiCGStab iterations per time step (right)</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>Figure 38: [Section 7.11.3] Standard resolution vs. adaptive resolution: cumulated oil rate (left) and number of cells (right)</figcaption></figure>

gas component and the estimators given in (7.37)–(7.38), the algebraic estimator $\eta_{alg,c}^{n,k,i}$ of (7.37), (7.38f) (with $j = 20$ additional steps), and the relative algebraic residual given by

$$ \text{err}_{\text{alg}}^{n,k,i} := \frac{\| \mathbb{A}^{n,k-1} \mathbf{X}^{n,k,i} - \mathbb{B}^{n,k-1} \|}{\| \mathbb{B}^{n,k-1} \|} $$

with $\mathbb{A}^{n,k-1} \mathbf{X}^{n,k} = \mathbb{B}^{n,k-1}$ being the linear system resulting from the $k$-th iteration of the blended Newton/fixed-point linearization at time step $t^n$, see Section 7.6. For the standard resolution, we stop the algebraic iteration using a fixed threshold $\text{err}_{\text{alg}}^{n,k,i} \le 10^{-6}$, following our usual practice. In the adaptive resolution based on Algorithm 7.1, we fix $\gamma_{\text{alg}} = 10^{-2}$.

We remark, in the left part of Figure 37, that the algebraic estimator steadily decreases, while the total estimator almost stagnates after about a third of total number of iterations necessary to converge using the standard stopping criterion. In the right part of Figure 37, we depict the cumulated number of BiCGStab iterations at each time step (the sum of the necessary number of BiCGStab iterations at each linearization iteration of the time step). We observe a significant gain with the adaptive stopping criterion.

We compare in Figure 38 the number of cells and the cumulated rate of oil production resulting from both the standard and the adaptive resolution. We can observe in the left part of Figure 38 that the adaptive algorithm does not have any significant influence on the accuracy of production. The right part of Figure 38 shows an important reduction in the number of cells via the adaptive resolution compared with the standard one.

The necessary linear solver steps and the different CPU times of the adaptive and standard resolutions are collected in Table 2. It shows that the total number of linear solver steps is reduced by 70%, which is an important gain. An important gain is also observed in the CPU time. We in particular remark that the time spent for the adaptation (AMR and estimators evaluation) only represents a small part of the actual resolution time. Globally, a reduction factor of 3.8 of the overall CPU time is obtained by comparing the adaptive resolution with the standard one. Important increase of this reduction factor is
<table><caption>Table 2: [Section 7.11.3] Comparison between standard and adaptive resolutions</caption><thead><tr><th></th><th>Linear solver<br>steps</th><th>Resolution<br>time</th><th>AMR<br>time</th><th>Estimators<br>evaluation</th><th>Gain<br>factor</th></tr></thead><tbody><tr><td>Standard resolution</td><td>66386</td><td>1023s</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Adaptive resolution</td><td>20184</td><td>201s</td><td>42s</td><td>26s</td><td>3.8</td></tr></tbody></table>

still to be expected when also the linearization solver stopping criteria and time step adaptation are used.

## 7.12 Bibliographic resources

Problem (7.1)–(7.12) is highly involved. For this reason, the existence and uniqueness of the weak solution is supposed in Assumption 7.2. For simplified problems, existence and uniqueness of the weak solution together with existence, uniqueness, and a priori error estimates for numerical discretizations have been derived in Alt and Luckhaus [16], Alt et al. [17], Kröner and Luckhaus [169], Otto [207, 208], Knabner and Otto [166], Chen [81]–[82] and Chen and Ewing [84], Smaï [231], Caro et al. [72], Cao and Pop [71], Ruiz-Baier and Lunati [226], Feireisl et al. [133], Kou et al. [167], Murphy and Walkington [194], Beneš [46], and Brenner et al. [59, 60], see also the references therein. Some classical books on the subject are those of Bear and Bachmat [35], Chavent and Jaffré [78], and Chen et al. [86].

For model (unsteady nonlinear) advection–diffusion–reaction problems, a posteriori analyses have been initiated in Eriksson and Johnson [117], Süli [234, 235], Angermann et al. [21], Verfürth [241, 240, 244], Kröner and Ohlberger [170], Ohlberger [203, 202], Ohlberger and Rohde [204], Kröner et al. [168], Nicaise [197], and Ern et al. [122], see also the references therein. A detailed attention has been paid to time-evolution problems of parabolic type in order to derive a posteriori error estimates robust with respect to the final time. This has been achieved in Verfürth [242, 243, 238], see also Bergam et al. [47], Ern and Vohralík [123], Hilhorst and Vohralík [151], Georgoulis et al. [137, 138, 139]. Lately, locally space–time efficient a posteriori error estimates robust with respect to the final time were obtained in Ern et al. [120], see also Ern et al. [121] and Mitra and Vohralík [190].

Finally, for results on a posteriori error estimates and adaptivity for numerical discretizations of porous media flows and namely multiphase compositional flows, we refer to Saad and Zhang [228], Chen and Ewing [85], Jenny et al. [153], Chen and Liu [80], Klieber and Rivière [165], Chueh et al. [91, 90], Vohralík and Wheeler [253], El Hassouni and Mghazli [115], Bernardi et al. [51], Cancès et al. [70], Di Pietro et al. [103, 104], Faigle et al. [131, 130], Henning et al. [148], Buhr et al. [67], Ganis et al. [135], Beaude et al. [36], Ahmed et al. [10, 9], and Li et al. [179], Li and Wheeler [181], Ben Gharbia et al. [44, 45], Deucher et al. [98], see also the references therein.

## References

[1] Aavatsmark, I., Barkve, T., Bøe, Ø., and Mannseth, T. Discretization on unstructured grids for inhomogeneous, anisotropic media. I. Derivation of the methods. *SIAM J. Sci. Comput.* **19** (1998), 1700–1716. https://doi.org/10.1137/S1064827595293582.

[2] Aavatsmark, I., Eigestad, G. T., Klausen, R. A., Wheeler, M. F., and Yotov, I. Convergence of a symmetric MPFA method on quadrilateral grids. *Comput. Geosci.* **11** (2007), 333–345. http://dx.doi.org/10.1007/s10596-007-9056-8.

[3] Achdou, Y., Bernardi, C., and Coquel, F. A priori and a posteriori analysis of finite volume discretizations of Darcy’s equations. *Numer. Math.* **96** (2003), 17–42. http://dx.doi.org/10.1007/s00211-002-0436-7.

[4] Ács, G., Doleschall, S., and Farkas, É. General purpose compositional model. *Society of Petroleum Engineers Journal* **25** (1985), 543–553. https://doi.org/10.2118/10515-PA.

[5] Adams, R. A. Pure and Applied Mathematics, Vol. 65. Sobolev spaces. Academic Press [A subsidiary of Harcourt Brace Jovanovich, Publishers], New York-London, 1975.

[6] Afif, M., Bergam, A., Mghazli, Z., and Verfürth, R. A posteriori estimators for the finite volume discretization of an elliptic problem. *Numer. Algorithms* **34** (2003), 127–136. International Conference on Numerical Algorithms, Vol. II (Marrakesh, 2001), https://doi.org/10.1023/B:NUMA.0000005400.45852.f3.
[7] Agouzal, A., Baranger, J., Maître, J.-F., and Oudin, F. Connection between finite volume and mixed finite element methods for a diffusion problem with nonconstant coefficients. Application to a convection diffusion problem. *East-West J. Numer. Math.* **3** (1995), 237–254.

[8] Agouzal, A., and Oudin, F. A posteriori error estimator for finite volume methods. *Appl. Math. Comput.* **110** (2000), 239–250. https://doi.org/10.1016/S0096-3003(99)00118-6.

[9] Ahmed, E., Nordbotten, J. M., and Radu, F. A. Adaptive asynchronous time-stepping, stopping criteria, and a posteriori error estimates for fixed-stress iterative schemes for coupled poromechanics problems. *J. Comput. Appl. Math.* **364** (2020), 112312, 25. https://doi.org/10.1016/j.cam.2019.06.028.

[10] Ahmed, E., Radu, F. A., and Nordbotten, J. M. Adaptive poromechanics computations based on a posteriori error estimates for fully mixed formulations of Biot's consolidation model. *Comput. Methods Appl. Mech. Engrg.* **347** (2019), 264–294. https://doi.org/10.1016/j.cma.2018.12.016.

[11] Ahusborde, E., Azaïez, M., Ben Belgacem, F., and Bernardi, C. Automatic simplification of Darcy's equations with pressure dependent permeability. *ESAIM Math. Model. Numer. Anal.* **47** (2013), 1797–1820. http://dx.doi.org/10.1051/m2an/2013089.

[12] Ainsworth, M. Robust a posteriori error estimation for nonconforming finite element approximation. *SIAM J. Numer. Anal.* **42** (2005), 2320–2341. https://doi.org/10.1137/S0036142903425112.

[13] Ainsworth, M. A posteriori error estimation for lowest order Raviart-Thomas mixed finite elements. *SIAM J. Sci. Comput.* **30** (2007/08), 189–204. https://doi.org/10.1137/06067331X.

[14] Ainsworth, M., and Oden, J. T. *A posteriori error estimation in finite element analysis*. Pure and Applied Mathematics (New York), Wiley-Interscience [John Wiley & Sons], New York, 2000. https://doi.org/10.1002/9781118032824.

[15] Allaire, G. *Analyse numérique et optimisation*. Editions Ecole Polytechnique, Palaiseau, France, 2005.

[16] Alt, H. W., and Luckhaus, S. Quasilinear elliptic-parabolic differential equations. *Math. Z.* **183** (1983), 311–341. http://dx.doi.org/10.1007/BF01176474.

[17] Alt, H. W., Luckhaus, S., and Visintin, A. On nonstationary flow through porous media. *Ann. Mat. Pura Appl.* **(4)** **136** (1984), 303–316. http://dx.doi.org/10.1007/BF01773387.

[18] Amaziane, B., Bergam, A., El Ossmani, M., and Mghazli, Z. A posteriori estimators for vertex centred finite volume discretization of a convection-diffusion-reaction equation arising in flow in porous media. *Internat. J. Numer. Methods Fluids* **59** (2009), 259–284. https://doi.org/10.1002/fld.1456.

[19] Amaziane, B., Jurak, M., and Žgaljić Keko, A. An existence result for a coupled system modeling a fully equivalent global pressure formulation for immiscible compressible two-phase flow in porous media. *J. Differential Equations* **250** (2011), 1685–1718. http://dx.doi.org/10.1016/j.jde.2010.09.008.

[20] Angermann, L. Balanced a posteriori error estimates for finite-volume type discretizations of convection-dominated elliptic problems. *Computing* **55** (1995), 305–323. https://doi.org/10.1007/BF02238485.

[21] Angermann, L., Knabner, P., and Thiele, K. An error estimator for a finite volume discretization of density driven flow in porous media. *Appl. Numer. Math.* **26** (1998), 179–191. Proceedings of the International Centre for Mathematical Sciences Conference on Grid Adaptation in Computational PDEs: Theory and Applications (Edinburgh, 1996), https://doi.org/10.1016/S0168-9274(97)00084-6.

[22] Antonietti, P. F., Beirão da Veiga, L., Lovadina, C., and Verani, M. Hierarchical a posteriori error estimators for the mimetic discretization of elliptic problems. *SIAM J. Numer. Anal.* **51** (2013), 654–675. http://dx.doi.org/10.1137/120873157.
[23] Antontsev, S. N., Kazhikhov, A. V., and Monakhov, V. N. Translated from the Russian. *Boundary value problems in mechanics of nonhomogeneous fluids*, vol. **22** of *Studies in Mathematics and its Applications*. North-Holland Publishing Co., Amsterdam, 1990.

[24] Arioli, M. A stopping criterion for the conjugate gradient algorithms in a finite element method framework. *Numer. Math.* **97** (2004), 1–24. http://dx.doi.org/10.1007/s00211-003-0500-y.

[25] Arioli, M., Georgoulis, E. H., and Loghin, D. Stopping criteria for adaptive finite element solvers. *SIAM J. Sci. Comput.* **35** (2013), A1537–A1559. http://dx.doi.org/10.1137/120867421.

[26] Arioli, M., Liesen, J., Międlar, A., and Strakoš, Z. Interplay between discretization and algebraic computation in adaptive numerical solution of elliptic PDE problems. *GAMM-Mitt.* **36** (2013), 102–129. http://dx.doi.org/10.1002/gamm.201310006.

[27] Axelsson, O., and Kaporin, I. Error norm estimation and stopping criteria in preconditioned conjugate gradient iterations. *Numer. Linear Algebra Appl.* **8** (2001), 265–286. https://doi.org/10.1002/nla.244.

[28] Aziz, K., and Settari, A. *Petroleum Reservoir Simulation*. Applied Science Publishers, Ltd, London, 1979.

[29] Babuška, I., and Strouboulis, T. *The finite element method and its reliability*. Numerical Mathematics and Scientific Computation. The Clarendon Press Oxford University Press, New York, 2001. https://doi.org/10.1093/oso/9780198502760.001.0001.

[30] Bakushinsky, A., and Smirnova, A. A posteriori stopping rule for regularized fixed point iterations. *Nonlinear Anal.* **64** (2006), 1255–1261. http://dx.doi.org/10.1016/j.na.2005.06.031.

[31] Bangerth, W., and Rannacher, R. *Adaptive finite element methods for differential equations*. Lectures in Mathematics ETH Zürich. Birkhäuser Verlag, Basel, 2003. https://doi.org/10.1007/978-3-0348-7605-6.

[32] Baranger, J., Maître, J.-F., and Oudin, F. Connection between finite volume and mixed finite element methods. *RAIRO Modél. Math. Anal. Numér.* **30** (1996), 445–465. https://doi.org/10.1051/m2an/1996300404451.

[33] Bartels, S., and Kaltenbach, A. Exact a posteriori error control for variational problems via convex duality and explicit flux reconstruction. In *Error control, adaptive discretizations, and applications. Part 1*, vol. **58** of *Adv. Appl. Mech*. Academic Press, San Diego, CA, 2024, pp. 295–375.

[34] Bastian, P., and Rivière, B. Superconvergence and $H(\text{div})$ projection for discontinuous Galerkin methods. *Internat. J. Numer. Methods Fluids* **42** (2003), 1043–1057. https://doi.org/10.1002/fld.562.

[35] Bear, J., and Bachmat, Y. *Introduction to Modeling of Transport Phenomena in Porous Media*, vol. **4** of *Theory and Applications of Transport in Porous Media*. Kluwer Academic Publishers, Dordrecht, Holland, 1990.

[36] Beaude, L., Brenner, K., Lopez, S., Masson, R., and Smai, F. Non-isothermal compositional liquid gas Darcy flow: formulation, soil-atmosphere boundary condition and application to high-energy geothermal simulations. *Comput. Geosci.* **23** (2019), 443–470. https://doi.org/10.1007/s10596-018-9794-9.

[37] Bebendorf, M. A note on the Poincaré inequality for convex domains. *Z. Anal. Anwendungen* **22** (2003), 751–756. http://dx.doi.org/10.4171/ZAA/1170.

[38] Becker, R., Capatina, D., and Luce, R. Stopping criteria based on locally reconstructed fluxes. In *Numerical mathematics and advanced applications—ENUMATH 2013*, vol. **103** of *Lect. Notes Comput. Sci. Eng*. Springer, Cham, 2015, pp. 243–251.

[39] Becker, R., Johnson, C., and Rannacher, R. Adaptive error control for multigrid finite element methods. *Computing* **55** (1995), 271–288. http://dx.doi.org/10.1007/BF02238483.
[40] Beirão da Veiga, L., Brezzi, F., Marini, L. D., and Russo, A. Mixed virtual element methods for general second order elliptic problems on polygonal meshes. *ESAIM Math. Model. Numer. Anal.* **50** (2016), 727-747. http://dx.doi.org/10.1051/m2an/2015067.

[41] Beirão da Veiga, L., Lipnikov, K., and Manzini, G. *The mimetic finite difference method for elliptic problems*, vol. **11** of *MS&A. Modeling, Simulation and Applications*. Springer, Cham, 2014. https://doi.org/10.1007/978-3-319-02663-3.

[42] Beirão da Veiga, L. A residual based error estimator for the mimetic finite difference method. *Numer. Math.* **108** (2008), 387-406. http://dx.doi.org/10.1007/s00211-007-0126-6.

[43] Beirão da Veiga, L., and Manzini, G. An a posteriori error estimator for the mimetic finite difference approximation of elliptic problems. *Internat. J. Numer. Methods Engrg.* **76** (2008), 1696-1723. http://dx.doi.org/10.1002/nme.2377.

[44] Ben Gharbia, I., Dabaghi, J., Martin, V., and Vohralík, M. A posteriori error estimates for a compositional two-phase flow with nonlinear complementarity constraints. *Comput. Geosci.* **24** (2020), 1031-1055. https://doi.org/10.1007/s10596-019-09909-5.

[45] Ben Gharbia, I., Ferzly, J., Vohralík, M., and Yousef, S. Semismooth and smoothing Newton methods for nonlinear systems with complementarity constraints: Adaptivity and inexact resolution. *J. Comput. Appl. Math.* **420** (2023), 114765. https://doi.org/10.1016/j.cam.2022.114765.

[46] Beneš, M. Analysis of non-isothermal multiphase flows in porous media. *Math. Methods Appl. Sci.* **45** (2022), 9653-9677. https://doi.org/10.1002/mma.8328.

[47] Bergam, A., Bernardi, C., and Mghazli, Z. A posteriori analysis of the finite element discretization of some parabolic equations. *Math. Comp.* **74** (2005), 1117-1138. https://doi.org/10.1090/S0025-5718-04-01697-7.

[48] Bergam, A., Mghazli, Z., and Verfürth, R. Estimations a posteriori d'un schéma de volumes finis pour un problème non linéaire. *Numer. Math.* **95** (2003), 599-624. https://doi.org/10.1007/s00211-003-0460-2.

[49] Bernardi, C., Chacón Rebollo, T., Chacón Vera, E., and Franco Coronil, D. A posteriori error analysis for two-overlapping domain decomposition techniques. *Appl. Numer. Math.* **59** (2009), 1214-1236. http://dx.doi.org/10.1016/j.apnum.2008.06.004.

[50] Bernardi, C., Dakroub, J., Mansour, G., and Sayah, T. A posteriori analysis of iterative algorithms for a nonlinear problem. *J. Sci. Comput.* **65** (2015), 672-697. https://doi.org/10.1007/s10915-014-9980-4.

[51] Bernardi, C., El Alaoui, L., and Mghazli, Z. A posteriori analysis of a space and time discretization of a nonlinear model for the flow in partially saturated porous media. *IMA J. Numer. Anal.* **34** (2014), 1002-1036. http://dx.doi.org/10.1093/imanum/drt014.

[52] Bochev, P. B., and Hyman, J. M. Principles of mimetic discretizations of differential operators. In *Compatible spatial discretizations*, vol. **142** of *IMA Vol. Math. Appl*. Springer, New York, 2006, pp. 89-119. http://dx.doi.org/10.1007/0-387-38034-5_5.

[53] Boffi, D., Brezzi, F., and Fortin, M. *Mixed finite element methods and applications*, vol. **44** of Springer Series in Computational Mathematics. Springer, Heidelberg, 2013. https://doi.org/10.1007/978-3-642-36519-5.

[54] Boffi, D., and Di Pietro, D. A. Unified formulation and analysis of mixed and primal discontinuous skeletal methods on polytopal meshes. *ESAIM Math. Model. Numer. Anal.* **52** (2018), 1-28. https://doi.org/10.1051/m2an/2017036.

[55] Bonelle, J., and Ern, A. Analysis of compatible discrete operator schemes for elliptic problems on polyhedral meshes. *ESAIM Math. Model. Numer. Anal.* **48** (2014), 553-581. http://dx.doi.org/10.1051/m2an/2013104.

[56] Bossavit, A. "Hypercircle" et majorations d'erreur a posteriori par "corrections locales". In *Équations aux dérivées partielles et applications*. Gauthier-Villars, Éd. Sci. Méd. Elsevier, Paris, 1998, pp. 221-237.
[57] Braess, D., and Schöberl, J. Equilibrated residual error estimator for edge elements. *Math. Comp.* **77** (2008), 651–672. http://dx.doi.org/10.1090/S0025-5718-07-02080-7.

[58] Breil, J., and Maire, P.-H. A cell-centered diffusion scheme on two-dimensional unstructured meshes. *J. Comput. Phys.* **224** (2007), 785–823. http://dx.doi.org/10.1016/j.jcp.2006.10.025.

[59] Brenner, K., Chorfi, N., and Masson, R. Sequential implicit vertex approximate gradient discretization of incompressible two-phase Darcy flows with discontinuous capillary pressure. *Comput. Geosci.* **26** (2022), 147–169. https://doi.org/10.1007/s10596-021-10113-7.

[60] Brenner, K., Masson, R., Quenjel, E. H., and Droniou, J. Total velocity-based finite volume discretization of two-phase Darcy flow in highly heterogeneous media with discontinuous capillary pressure. *IMA J. Numer. Anal.* **42** (2022), 1231–1272. https://doi.org/10.1093/imanum/drab018.

[61] Brenner, S. C., and Scott, L. R. *The mathematical theory of finite element methods*, third ed., vol. **15** of *Texts in Applied Mathematics*. Springer, New York, 2008. http://dx.doi.org/10.1007/978-0-387-75934-0.

[62] Brezzi, F., Falk, R. S., and Marini, L. D. Basic principles of mixed virtual element methods. *ESAIM Math. Model. Numer. Anal.* **48** (2014), 1227–1240. https://doi.org/10.1051/m2an/2013138.

[63] Brezzi, F., and Fortin, M. *Mixed and hybrid finite element methods*, vol. **15** of Springer Series in Computational Mathematics. Springer-Verlag, New York, 1991. http://dx.doi.org/10.1007/978-1-4612-3172-1.

[64] Brezzi, F., Lipnikov, K., and Shashkov, M. Convergence of the mimetic finite difference method for diffusion problems on polyhedral meshes. *SIAM J. Numer. Anal.* **43** (2005), 1872–1896. https://doi.org/10.1137/040613950.

[65] Brezzi, F., Lipnikov, K., and Simoncini, V. A family of mimetic finite difference methods on polygonal and polyhedral meshes. *Math. Models Methods Appl. Sci.* **15** (2005), 1533–1551. http://dx.doi.org/10.1142/S0218202505000832.

[66] Bringmann, P., Feischl, M., Miraçi, A., Praetorius, D., and Streitberger, J. On full linear convergence and optimal complexity of adaptive FEM with inexact solver. *Comput. Math. Appl.* **180** (2025), 102–129. https://doi.org/10.1016/j.camwa.2024.12.013.

[67] Buhr, A., Engwer, C., Ohlberger, M., and Rave, S. ArbiLoMod, a simulation technique designed for arbitrary local modifications. *SIAM J. Sci. Comput.* **39** (2017), A1435–A1465. https://doi.org/10.1137/15M1054213.

[68] Burman, E., and Ern, A. Continuous interior penalty hp-finite element methods for advection and advection-diffusion equations. *Math. Comp.* **76** (2007), 1119–1140. https://doi.org/10.1090/S0025-5718-07-01951-5.

[69] Cancès, C., Gallouët, T., and Porretta, A. Two-phase flows involving capillary barriers in heterogeneous porous media. *Interfaces Free Bound.* **11** (2009), 239–258. http://dx.doi.org/10.4171/IFB/210.

[70] Cancès, C., Pop, I. S., and Vohralík, M. An a posteriori error estimate for vertex-centered finite volume discretizations of immiscible incompressible two-phase flow. *Math. Comp.* **83** (2014), 153–188. http://dx.doi.org/10.1090/S0025-5718-2013-02723-8.

[71] Cao, X., and Pop, I. S. Two-phase porous media flows with dynamic capillary effects and hysteresis: uniqueness of weak solutions. *Comput. Math. Appl.* **69** (2015), 688–695. https://doi.org/10.1016/j.camwa.2015.02.009.

[72] Caro, F., Saad, B., and Saad, M. Study of degenerate parabolic system modeling the hydrogen displacement in a nuclear waste repository. *Discrete Contin. Dyn. Syst. Ser. S* **7** (2014), 191–205. https://doi.org/10.3934/dcdss.2014.7.191.

[73] Carstensen, C., Lazarov, R., and Tomov, S. Explicit and averaging a posteriori error estimates for adaptive finite volume methods. *SIAM J. Numer. Anal.* **42** (2005), 2496–2521. https://doi.org/10.1137/S0036142903425422.
[74] Chaillou, A., and Suri, M. A posteriori estimation of the linearization error for strongly monotone nonlinear operators. *J. Comput. Appl. Math.* **205** (2007), 72–87. http://dx.doi.org/10.1016/j.cam.2006.04.041.

[75] Chaillou, A. L., and Suri, M. Computable error estimators for the approximation of nonlinear problems by linearized models. *Comput. Methods Appl. Mech. Engrg.* **196** (2006), 210–224. http://dx.doi.org/10.1016/j.cma.2006.03.008.

[76] Chamoin, L., and Legoll, F. An introductory review on a posteriori error estimation in finite element computations. *SIAM Rev.* **65** (2023), 963–1028. https://doi.org/10.1137/21M1464841.

[77] Chaudhry, J. H., Estep, D., and Tavener, S. J. A posteriori error analysis for Schwarz overlapping domain decomposition methods. *BIT* **61** (2021), 1153–1191. https://doi.org/10.1007/s10543-021-00864-1.

[78] Chavent, G., and Jaffré, J. Studies in Mathematics and Its Applications, Vol. 17. *Mathematical models and finite elements for reservoir simulation*. North-Holland, Amsterdam, 1986.

[79] Chavent, G., Younès, A., and Ackerer, P. On the finite volume reformulation of the mixed finite element method for elliptic and parabolic PDE on triangles. *Comput. Methods Appl. Mech. Engrg.* **192** (2003), 655–682. https://doi.org/10.1016/S0045-7825(02)00578-9.

[80] Chen, Y., and Liu, W. A posteriori error estimates of mixed methods for miscible displacement problems. *Internat. J. Numer. Methods Engrg.* **73** (2008), 331–343. http://dx.doi.org/10.1002/nme.2075.

[81] Chen, Z. Degenerate two-phase incompressible flow. I. Existence, uniqueness and regularity of a weak solution. *J. Differential Equations* **171** (2001), 203–232. http://dx.doi.org/10.1006/jdeq.2000.3848.

[82] Chen, Z. Degenerate two-phase incompressible flow. II. Regularity, stability and stabilization. *J. Differential Equations* **186** (2002), 345–376. http://dx.doi.org/10.1016/S0022-0396(02)00027-X.

[83] Chen, Z. Mathematical techniques in oil recovery. *Reservoir simulation*, vol. **77** of CBMS-NSF Regional Conference Series in Applied Mathematics. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2007. https://doi.org/10.1137/1.9780898717075.

[84] Chen, Z., and Ewing, R. E. Degenerate two-phase incompressible flow. III. Sharp error estimates. *Numer. Math.* **90** (2001), 215–240. http://dx.doi.org/10.1007/s002110100291.

[85] Chen, Z., and Ewing, R. E. Degenerate two-phase incompressible flow. IV. Local refinement and domain decomposition. *J. Sci. Comput.* **18** (2003), 329–360. http://dx.doi.org/10.1023/A:1022673427893.

[86] Chen, Z., Huan, G., and Ma, Y. *Computational methods for multiphase flows in porous media*. Computational Science & Engineering. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2006. http://dx.doi.org/10.1137/1.9780898718942.

[87] Christiansen, S. H. A construction of spaces of compatible differential forms on cellular complexes. *Math. Models Methods Appl. Sci.* **18** (2008), 739–757. http://dx.doi.org/10.1142/S021820250800284X.

[88] Christiansen, S. H. Éléments finis mixtes minimaux sur les polyèdres. *C. R. Math. Acad. Sci. Paris* **348** (2010), 217–221. http://dx.doi.org/10.1016/j.crma.2010.01.017.

[89] Christie, M. A., and Blunt, M. J. Tenth SPE comparative solution project: a comparison of upscaling techniques. *SPE Reservoir Eval. Eng.* **4** (2001), 308–317. https://doi.org/10.2118/72469-PA.

[90] Chueh, C.-C., Djilali, N., and Bangerth, W. An *h*-adaptive operator splitting method for two-phase flow in 3D heterogeneous porous media. *SIAM J. Sci. Comput.* **35** (2013), B149–B175. https://doi.org/10.1137/120866208.
[91] Chueh, C. C., Secanell, M., Bangerth, W., and Djilali, N. Multi-level adaptive simulation of transient two-phase flow in heterogeneous porous media. *Comput. & Fluids* **39** (2010), 1585–1596. https://doi.org/10.1016/j.compfluid.2010.05.011.

[92] Ciarlet, P. G. *The Finite Element Method for Elliptic Problems*, vol. **4** of *Studies in Mathematics and its Applications*. North-Holland, Amsterdam, 1978.

[93] Coats, K. H. An equation of state compositional model. *Society of Petroleum Engineers Journal* **20** (1980), 363–376. https://doi.org/10.2118/8284-PA.

[94] Coats, K. H. Implicit compositional simulation of single-porosity and dual-porosity reservoirs. In *SPE Symposium on Reservoir Simulation* (1989). Paper SPE-18427-MS, https://doi.org/10.2118/18427-MS.

[95] Cockburn, B., Di Pietro, D. A., and Ern, A. Bridging the hybrid high-order and hybridizable discontinuous Galerkin methods. *ESAIM Math. Model. Numer. Anal.* **50** (2016), 635–650. http://dx.doi.org/10.1051/m2an/2015051.

[96] Cottereau, R., Díez, P., and Huerta, A. Strict error bounds for linear solid mechanics problems using a subdomain-based flux-free method. *Comput. Mech.* **44** (2009), 533–547. https://doi.org/10.1007/s00466-009-0388-1.

[97] Destuynder, P., and Métivet, B. Explicit error bounds in a conforming finite element method. *Math. Comp.* **68** (1999), 1379–1396. http://dx.doi.org/10.1090/S0025-5718-99-01093-5.

[98] Deucher, R. H., Franc, J., Mø yner, O., and Tchelepi, H. A. Compositional reservoir simulation with a high-resolution compact stencil adaptive implicit method. *J. Comput. Phys.* **521** (2025), Paper No. 113558, 17. https://doi.org/10.1016/j.jcp.2024.113558.

[99] Deuflhard, P., and Weiser, M. *Adaptive numerical solution of PDEs*. de Gruyter Textbook. Walter de Gruyter & Co., Berlin, 2012. http://dx.doi.org/10.1515/9783110283112.

[100] Di Pietro, D. A., and Ern, A. *Mathematical aspects of discontinuous Galerkin methods*, vol. **69** of *Mathématiques & Applications (Berlin) [Mathematics & Applications]*. Springer, Heidelberg, 2012. http://dx.doi.org/10.1007/978-3-642-22980-0.

[101] Di Pietro, D. A., and Ern, A. A hybrid high-order locking-free method for linear elasticity on general meshes. *Comput. Methods Appl. Mech. Engrg.* **283** (2015), 1–21. http://dx.doi.org/10.1016/j.cma.2014.09.009.

[102] Di Pietro, D. A., Ern, A., and Lemaire, S. An arbitrary-order and compact-stencil discretization of diffusion on general meshes based on local reconstruction operators. *Comput. Methods Appl. Math.* **14** (2014), 461–472. https://doi.org/10.1515/cmam-2014-0018.

[103] Di Pietro, D. A., Flauraud, E., Vohralík, M., and Yousef, S. A posteriori error estimates, stopping criteria, and adaptivity for multiphase compositional Darcy flows in porous media. *J. Comput. Phys.* **276** (2014), 163–187. http://dx.doi.org/10.1016/j.jcp.2014.06.061.

[104] Di Pietro, D. A., Vohralík, M., and Yousef, S. An a posteriori-based, fully adaptive algorithm with adaptive stopping criteria and mesh refinement for thermal multiphase compositional flows in porous media. *Comput. Math. Appl.* **68** (2014), 2331–2347. https://doi.org/10.1016/j.camwa.2014.08.008.

[105] Ding, Y., and Lemonnier, P. Use of corner point geometry in reservoir simulation. In *International Meeting on Petroleum Engineering* (1995). Paper SPE-29933-MS, https://doi.org/10.2118/29933-MS.

[106] Dolgov, S., and Vejchodský, T. Guaranteed a posteriori error bounds for low-rank tensor approximate solutions. *IMA J. Numer. Anal.* **41** (2021), 1240–1266. https://doi.org/10.1093/imanum/draa010.

[107] Droniou, J. Finite volume schemes for diffusion equations: introduction to and review of modern methods. *Math. Models Methods Appl. Sci.* **24** (2014), 1575–1619. http://dx.doi.org/10.1142/S0218202514400041.
[108] Droniou, J., and Eymard, R. A mixed finite volume scheme for anisotropic diffusion problems on any grid. Numer. Math. 105 (2006), 35–71. https://doi.org/10.1007/s00211-006-0034-1.

[109] Droniou, J., Eymard, R., Gallouët, T., Guichard, C., and Herbin, R. The gradient discretisation method, vol. 82 of Mathématiques & Applications (Berlin) [Mathematics & Applications]. Springer, Cham, 2018. https://doi.org/10.1007/978-3-319-79042-8.

[110] Droniou, J., Eymard, R., Gallouët, T., and Herbin, R. A unified approach to mimetic finite difference, hybrid finite volume and mixed finite volume methods. Math. Models Methods Appl. Sci. 20 (2010), 265–295. http://dx.doi.org/10.1142/S0218202510004222.

[111] Droniou, J., Eymard, R., and Herbin, R. Gradient schemes: generic tools for the numerical analysis of diffusion equations. ESAIM Math. Model. Numer. Anal. 50 (2016), 749–781. http://dx.doi.org/10.1051/m2an/2015079.

[112] Durán, R., and Padra, C. An error estimator for nonconforming approximations of a nonlinear problem. In Finite element methods (Jyväskylä, 1993), vol. 164 of Lecture Notes in Pure and Appl. Math. Dekker, New York, 1994, pp. 201–205. https://doi.org/10.1201/b16924-21.

[113] Edwards, M. G. Unstructured, control-volume distributed, full-tensor finite-volume schemes with flow based grids. Comput. Geosci. 6 (2002), 433–452. Locally conservative numerical methods for flow in porous media, http://dx.doi.org/10.1023/A:1021243231313.

[114] El Alaoui, L., Ern, A., and Vohralík, M. Guaranteed and robust a posteriori error estimates and balancing discretization and linearization errors for monotone nonlinear problems. Comput. Methods Appl. Mech. Engrg. 200 (2011), 2782–2795. http://dx.doi.org/10.1016/j.cma.2010.03.024.

[115] El Hassouni, S., and Mghazli, Z. Mesh adaptation for two-phase flow. Int. J. Math. Stat. 13 (2013), 56–68.

[116] Erath, C., and Praetorius, D. A posteriori error estimate and adaptive mesh refinement for the cell-centered finite volume method for elliptic boundary value problems. SIAM J. Numer. Anal. 47 (2008/09), 109–135. http://dx.doi.org/10.1137/070702126.

[117] Eriksson, K., and Johnson, C. Adaptive streamline diffusion finite element methods for stationary convection-diffusion problems. Math. Comp. 60 (1993), 167–188, S1–S2. https://doi.org/10.2307/2153160.

[118] Ern, A., and Guermond, J.-L. Finite Elements I. Approximation and Interpolation, vol. 72 of Texts in Applied Mathematics. Springer International Publishing, Springer Nature Switzerland AG, 2021. https://doi-org/10.1007/978-3-030-56341-7.

[119] Ern, A., Nicaise, S., and Vohralík, M. An accurate $H(\text{div})$ flux reconstruction for discontinuous Galerkin approximations of elliptic problems. C. R. Math. Acad. Sci. Paris 345 (2007), 709–712. http://dx.doi.org/10.1016/j.crma.2007.10.036.

[120] Ern, A., Smears, I., and Vohralík, M. Guaranteed, locally space-time efficient, and polynomial-degree robust a posteriori error estimates for high-order discretizations of parabolic problems. SIAM J. Numer. Anal. 55 (2017), 2811–2834. https://doi.org/10.1137/16M1097626.

[121] Ern, A., Smears, I., and Vohralík, M. Equilibrated flux a posteriori error estimates in $L^2(H^1)$-norms for high-order discretizations of parabolic problems. IMA J. Numer. Anal. 39 (2019), 1158–1179. https://doi.org/10.1093/imanum/dry035.

[122] Ern, A., Stephansen, A. F., and Vohralík, M. Guaranteed and robust discontinuous Galerkin a posteriori error estimates for convection-diffusion-reaction problems. J. Comput. Appl. Math. 234 (2010), 114–130. http://dx.doi.org/10.1016/j.cam.2009.12.009.

[123] Ern, A., and Vohralík, M. A posteriori error estimation based on potential and flux reconstruction for the heat equation. SIAM J. Numer. Anal. 48 (2010), 198–223. http://dx.doi.org/10.1137/090759008.

[124] Ern, A., and Vohralík, M. Adaptive inexact Newton methods with a posteriori stopping criteria for nonlinear diffusion PDEs. SIAM J. Sci. Comput. 35 (2013), A1761–A1791. http://dx.doi.org/10.1137/120896918.
[125] Ern, A., and Vohralík, M. Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous Galerkin, and mixed discretizations. *SIAM J. Numer. Anal.* **53** (2015), 1058–1081. http://dx.doi.org/10.1137/130950100.

[126] Eymard, R., Gallouët, T., and Herbin, R. Finite volume methods. In *Handbook of Numerical Analysis, Vol. VII*. North-Holland, Amsterdam, 2000, pp. 713–1020. https://doi.org/10.1016/S1570-8659(00)07005-8.

[127] Eymard, R., Gallouët, T., and Herbin, R. Finite volume approximation of elliptic problems and convergence of an approximate gradient. *Appl. Numer. Math.* **37** (2001), 31–53. http://dx.doi.org/10.1016/S0168-9274(00)00024-6.

[128] Eymard, R., Gallouët, T., and Herbin, R. Discretization of heterogeneous and anisotropic diffusion problems on general nonconforming meshes SUSHI: a scheme using stabilization and hybrid interfaces. *IMA J. Numer. Anal.* **30** (2010), 1009–1043. http://dx.doi.org/10.1093/imanum/drn084.

[129] Eymard, R., Guichard, C., Herbin, R., and Masson, R. Vertex-centred discretization of multiphase compositional Darcy flows on general meshes. *Comput. Geosci.* **16** (2012), 987–1005. https://doi.org/10.1007/s10596-012-9299-x.

[130] Faigle, B., Elfeel, M. A., Helmig, R., Becker, B., Flemisch, B., and Geiger, S. Multi-physics modeling of non-isothermal compositional flow on adaptive grids. *Comput. Methods Appl. Mech. Engrg.* **292** (2015), 16–34. https://doi.org/10.1016/j.cma.2014.11.030.

[131] Faigle, B., Helmig, R., Aavatsmark, I., and Flemisch, B. Efficient multiphysics modelling with adaptive grid refinement using a MPFA method. *Comput. Geosci.* **18** (2014), 625–636. https://doi.org/10.1007/s10596-014-9407-1.

[132] Faille, I. A control volume method to solve an elliptic equation on a two-dimensional irregular mesh. *Comput. Methods Appl. Mech. Engrg.* **100** (1992), 275–290. https://doi.org/10.1016/0045-7825(92)90186-N.

[133] Feireisl, E., Hilhorst, D., Petzeltová, H., and Takáč, P. Mathematical analysis of variable density flows in porous media. *J. Evol. Equ.* **16** (2016), 1–19. https://doi.org/10.1007/s00028-015-0290-6.

[134] Févotte, F., Rappaport, A., and Vohralík, M. Adaptive regularization, discretization, and linearization for nonsmooth problems based on primal-dual gap estimators. *Comput. Methods Appl. Mech. Engrg.* **418** (2024), 116558. https://doi.org/10.1016/j.cma.2023.116558.

[135] Ganis, B., Pencheva, G., and Wheeler, M. F. Adaptive mesh refinement with an enhanced velocity mixed finite element method on semi-structured grids using a fully coupled solver. *Comput. Geosci.* **23** (2019), 149–168. https://doi.org/10.1007/s10596-018-9789-6.

[136] Gantner, G., Haberl, A., Praetorius, D., and Schimanko, S. Rate optimality of adaptive finite element methods with respect to overall computational costs. *Math. Comp.* **90** (2021), 2011–2040. https://doi.org/10.1090/mcom/3654.

[137] Georgoulis, E. H., Lakkis, O., and Virtanen, J. M. A posteriori error control for discontinuous Galerkin methods for parabolic problems. *SIAM J. Numer. Anal.* **49** (2011), 427–458. http://dx.doi.org/10.1137/080722461.

[138] Georgoulis, E. H., Lakkis, O., and Wihler, T. P. A posteriori error bounds for fully-discrete hp-discontinuous Galerkin timestepping methods for parabolic problems. *Numer. Math.* **148** (2021), 363–386. https://doi.org/10.1007/s00211-021-01187-7.

[139] Georgoulis, E. H., and Makridakis, C. G. Lower bounds, elliptic reconstruction and a posteriori error control of parabolic problems. *IMA J. Numer. Anal.* **43** (2023), 3212–3242. https://doi.org/10.1093/imanum/drac080.

[140] Golub, G. H., and Strakoš, Z. Estimates in quadratic formulas. *Numer. Algorithms* **8** (1994), 241–268. http://dx.doi.org/10.1007/BF02142693.
[141] Grospellier, G., and Lelandais, B. The Arcane development framework. In *Proceedings of the 8th workshop on Parallel/High-Performance Object-Oriented Scientific Computing* (New York, NY, USA, 2009), POOSC '09, ACM, pp. 4:1–4:11. http://doi.acm.org/10.1145/1595655.1595659.

[142] Haberl, A., Praetorius, D., Schimanko, S., and Vohralík, M. Convergence and quasi-optimal cost of adaptive algorithms for nonlinear operators including iterative linearization and algebraic solver. *Numer. Math.* **147** (2021), 679–725. https://doi.org/10.1007/s00211-021-01176-w.

[143] Han, W. A posteriori error analysis for linearization of nonlinear elliptic problems and their discretizations. *Math. Methods Appl. Sci.* **17** (1994), 487–508. http://dx.doi.org/10.1002/mma.1670170702.

[144] Han, W. With applications in modeling and numerical approximations. *A posteriori error analysis via duality theory*, vol. **8** of *Advances in Mechanics and Mathematics*. Springer-Verlag, New York, 2005. https://doi.org/10.1007/b101775.

[145] Harnist, A., Mitra, K., Rappaport, A., and Vohralík, M. Robust augmented energy a posteriori estimates for Lipschitz and strongly monotone elliptic problems. HAL Preprint 04033438, https://hal.inria.fr/hal-04033438, 2024.

[146] Heid, P., Praetorius, D., and Wihler, T. P. Energy contraction and optimal convergence of adaptive iterative linearized finite element methods. *Comput. Methods Appl. Math.* **21** (2021), 407–422. https://doi.org/10.1515/cmam-2021-0025.

[147] Heid, P., and Wihler, T. P. Adaptive iterative linearization Galerkin methods for nonlinear problems. *Math. Comp.* **89** (2020), 2707–2734. https://doi.org/10.1090/mcom/3545.

[148] Henning, P., Ohlberger, M., and Schweizer, B. Adaptive heterogeneous multiscale methods for immiscible two-phase flow in porous media. *Comput. Geosci.* **19** (2015), 99–114. https://doi.org/10.1007/s10596-014-9455-6.

[149] Herbin, R. An error estimate for a finite volume scheme for a diffusion-convection problem on a triangular mesh. *Numer. Methods Partial Differential Equations* **11** (1995), 165–173. https://doi.org/10.1002/num.1690110205.

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

[151] Hilhorst, D., and Vohralík, M. A posteriori error estimates for combined finite volume–finite element discretizations of reactive transport equations on nonmatching grids. *Comput. Methods Appl. Mech. Engrg.* **200** (2011), 597–613. http://dx.doi.org/10.1016/j.cma.2010.08.017.

[152] Hlaváček, I., Haslinger, J., Nečas, J., and Lovíšek, J. Translated from the Slovak by J. Jarník. *Solution of variational inequalities in mechanics*, vol. **66** of *Applied Mathematical Sciences*. Springer-Verlag, New York, 1988.

[153] Jenny, P., Lee, S. H., and Tchelep, H. A. Adaptive multiscale finite-volume method for multiphase flow and transport in porous media. *Multiscale Model. Simul.* **3** (2004/05), 50–64. http://dx.doi.org/10.1137/030600795.

[154] Jiránek, P., Strakoš, Z., and Vohralík, M. A posteriori error estimates including algebraic error and stopping criteria for iterative solvers. *SIAM J. Sci. Comput.* **32** (2010), 1567–1590. http://dx.doi.org/10.1137/08073706X.

[155] Jovanović, B. S., and Süli, E. For linear partial differential equations with generalized solutions. *Analysis of finite difference schemes*, vol. **46** of Springer Series in Computational Mathematics. Springer, London, 2014. https://doi.org/10.1007/978-1-4471-5460-0.

[156] Kaltenbacher, B. A posteriori parameter choice strategies for some Newton type methods for the regularization of nonlinear ill-posed problems. *Numer. Math.* **79** (1998), 501–528. http://dx.doi.org/10.1007/s002110050349.

[157] Karakashian, O. A., and Pascal, F. A posteriori error estimates for a discontinuous Galerkin approximation of second-order elliptic problems. *SIAM J. Numer. Anal.* **41** (2003), 2374–2399. http://dx.doi.org/10.1137/S0036142902405217.
[158] Keilegavlen, E., and Nordbotten, J. M. Conservative inexact solvers for porous media flow. In *Domain decomposition methods in science and engineering XXI*, vol. 98 of Lect. Notes Comput. Sci. Eng. Springer, Cham, 2014, pp. 333–340.

[159] Keilegavlen, E., and Nordbotten, J. M. Inexact linear solvers for control volume discretizations in porous media. *Comput. Geosci.* **19** (2015), 159–176. https://doi.org/10.1007/s10596-014-9453-8.

[160] Kelley, C. T. With separately available software. *Iterative methods for linear and nonlinear equations*, vol. **16** of *Frontiers in Applied Mathematics*. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1995. https://doi.org/10.1137/1.9781611970944.

[161] Khalil, Z., and Saad, M. Degenerate two-phase compressible immiscible flow in porous media: the case where the density of each phase depends on its own pressure. *Math. Comput. Simulation* **81** (2011), 2225–2233. http://dx.doi.org/10.1016/j.matcom.2010.12.012.

[162] Kim, K.-Y. A posteriori error analysis for locally conservative mixed methods. *Math. Comp.* **76** (2007), 43–66. http://dx.doi.org/10.1090/S0025-5718-06-01903-X.

[163] Kim, K.-Y. A posteriori error estimators for locally conservative methods of nonlinear elliptic problems. *Appl. Numer. Math.* **57** (2007), 1065–1080. http://dx.doi.org/10.1016/j.apnum.2006.09.010.

[164] Klausen, R. A., and Russell, T. F. Relationships among some locally conservative discretization methods which handle discontinuous coefficients. *Comput. Geosci.* **8** (2004), 341–377 (2005). http://dx.doi.org/10.1007/s10596-005-1815-9.

[165] Klieber, W., and Rivière, B. Adaptive simulations of two-phase flow by discontinuous Galerkin methods. *Comput. Methods Appl. Mech. Engrg.* **196** (2006), 404–419. http://dx.doi.org/10.1016/j.cma.2006.05.007.

[166] Knabner, P., and Otto, F. Solute transport in porous media with equilibrium and nonequilibrium multiple-site adsorption: uniqueness of weak solutions. *Nonlinear Anal.* **42** (2000), 381–403. https://doi.org/10.1016/S0362-546X(98)00352-6.

[167] Kou, J., Sun, S., and Wang, X. Linearly decoupled energy-stable numerical methods for multi-component two-phase compressible flow. *SIAM J. Numer. Anal.* **56** (2018), 3219–3248. https://doi.org/10.1137/17M1162287.

[168] Kröner, D., Küther, M., Ohlberger, M., and Rohde, C. A posteriori error estimates and adaptive methods for hyperbolic and convection dominated parabolic conservation laws. In *Trends in nonlinear analysis*. Springer, Berlin, 2003, pp. 289–306.

[169] Kröner, D., and Luckhaus, S. Flow of oil and water in a porous medium. *J. Differential Equations* **55** (1984), 276–288. http://dx.doi.org/10.1016/0022-0396(84)90084-6.

[170] Kröner, D., and Ohlberger, M. A posteriori error estimates for upwind finite volume schemes for nonlinear conservation laws in multidimensions. *Math. Comp.* **69** (2000), 25–39. http://dx.doi.org/10.1090/S0025-5718-99-01158-8.

[171] Kumar, K., Kyas, S., Nordbotten, J. M., and Repin, S. Guaranteed and computable error bounds for approximations constructed by an iterative decoupling of the Biot problem. *Comput. Math. Appl.* **91** (2021), 122–149. https://doi.org/10.1016/j.camwa.2020.05.005.

[172] Kuznetsov, Y., and Repin, S. New mixed finite element method on polygonal and polyhedral meshes. *Russian J. Numer. Anal. Math. Modelling* **18** (2003), 261–278. http://dx.doi.org/10.1163/156939803322380846.

[173] Kuznetsov, Y. A. Mixed finite element methods on polyhedral meshes for diffusion equations. In *Partial differential equations*, vol. **16** of *Comput. Methods Appl. Sci*. Springer, Dordrecht, 2008, pp. 27–41.

[174] Ladevèze, P. *Comparaison de modèles de milieux continus*. Ph.D. thesis, Université Pierre et Marie Curie (Paris 6), 1975.
[175] Ladevèze, P., and Leguillon, D. Error estimate procedure in the finite element method and applications. *SIAM J. Numer. Anal.* **20** (1983), 485-509.

[176] Ladevèze, P., and Pelle, J.-P. Translated from the 2001 French original by Theofanis Strouboulis. *Mastering calculations in linear and nonlinear mechanics*. Mechanical Engineering Series. Springer-Verlag, New York, 2005. https://doi.org/10.1007/b138705.

[177] Larson, M. G., and Niklasson, A. J. A conservative flux for the continuous Galerkin method based on discontinuous enrichment. *Calcolo* **41** (2004), 65-76. https://doi.org/10.1007/s10092-004-0084-7.

[178] Lazarov, R., and Tomov, S. A posteriori error estimates for finite volume element approximations of convection-diffusion-reaction equations. *Comput. Geosci.* **6** (2002), 483-503. Locally conservative numerical methods for flow in porous media, https://doi.org/10.1023/A:1021247300362.

[179] Li, H., Leung, W. T., and Wheeler, M. F. Sequential local mesh refinement solver with separate temporal and spatial adaptivity for non-linear two-phase flow problems. *J. Comput. Phys.* **403** (2020), 109074, 23. https://doi.org/10.1016/j.jcp.2019.109074.

[180] Li, H., Mu, L., and Ye, X. A posteriori error estimates for the weak Galerkin finite element methods on polytopal meshes. *Commun. Comput. Phys.* **26** (2019), 558-578. https://doi.org/10.4208/cicp.oa-2018-0058.

[181] Li, H., and Wheeler, M. F. Dynamic local coupling for multiphase flow: a compromise between efficiency and stability. *J. Comput. Phys.* **469** (2022), Paper No. 111535, 14. https://doi.org/10.1016/j.jcp.2022.111535.

[182] Liu, L., and Keyes, D. E. Approximate error bounds on solutions of nonlinearly preconditioned PDEs. *SIAM J. Sci. Comput.* **43** (2021), A2526-A2554. https://doi.org/10.1137/19M1285044.

[183] Luce, R., and Wohlmuth, B. I. A local a posteriori error estimator based on equilibrated fluxes. *SIAM J. Numer. Anal.* **42** (2004), 1394-1414. http://dx.doi.org/10.1137/S0036142903433790.

[184] Manteuffel, T. A., and White, Jr., A. B. The numerical solution of second-order boundary value problems on nonuniform meshes. *Math. Comp.* **47** (1986), 511-535, S53-S55. https://doi.org/10.2307/2008170.

[185] Mesri, Y., and Ricois, O. Construction process for an improved meshing for the simulation of a reservoir in an underground formation. Patent CA 28886110, http://brevets-patents.ic.gc.ca/ topic-cipo/cpd/fra/brevet/2886110/, 23 03, 2015.

[186] Mghazli, Z. A posteriori error analysis for finite element approximation of some groundwater models. Part I: Linear case. In *Error control, adaptive discretizations, and applications. Part 3*, vol. **60** of *Adv. Appl. Mech*. Academic Press, San Diego, CA, 2025. In press.

[187] Mghazli, Z. A posteriori error analysis for finite element approximation of some groundwater models. Part II: Richards equation and mixed finite element. In *Error control, adaptive discretizations, and applications. Part 4*, vol. **61** of *Adv. Appl. Mech*. Academic Press, San Diego, CA, 2025. In press.

[188] Mitchell, W. F. A collection of 2D elliptic problems for testing adaptive grid refinement algorithms. *Appl. Math. Comput.* **220** (2013), 350-364. https://doi.org/10.1016/j.amc.2013.05.068.

[189] Mitra, K., and Vohralík, M. Guaranteed, locally efficient, and robust a posteriori estimates for nonlinear elliptic problems in iteration-dependent norms. An orthogonal decomposition result based on iterative linearization. HAL Preprint 04156711, submitted for publication, https://hal.inria.fr/hal-04156711, 2023.

[190] Mitra, K., and Vohralík, M. A posteriori error estimates for the Richards equation. *Math. Comp.* **93** (2024), 1053-1096. https://doi.org/10.1090/mcom/3932.

[191] Morton, K. W., and Süli, E. Finite volume methods and their analysis. *IMA J. Numer. Anal.* **11** (1991), 241-260. https://doi.org/10.1093/imanum/11.2.241.
[192] Mozolevski, I., and Prudhomme, S. Goal-oriented error estimation based on equilibrated-flux re-construction for finite element approximations of elliptic problems. *Comput. Methods Appl. Mech. Engrg.* **288** (2015), 127-145. https://doi.org/10.1016/j.cma.2014.09.025.

[193] Munar, M., Cangiani, A., and Velásquez, I. Residual-based a posteriori error estimation for mixed virtual element methods. *Comput. Math. Appl.* **166** (2024), 182-197. https://doi.org/10.1016/j.camwa.2024.05.011.

[194] Murphy, T. J., and Walkington, N. J. Control volume approximation of degenerate two-phase porous flows. *SIAM J. Numer. Anal.* **57** (2019), 527-546. https://doi.org/10.1137/17M1160744.

[195] Nédélec, J.-C. Mixed finite elements in $\mathbb{R}^3$. *Numer. Math.* **35** (1980), 315-341. https://doi.org/10.1007/BF01396415.

[196] Nicaise, S. A posteriori error estimations of some cell-centered finite volume methods. *SIAM J. Numer. Anal.* **43** (2005), 1481-1503. https://doi.org/10.1137/S0036142903437787.

[197] Nicaise, S. A posteriori error estimations of some cell centered finite volume methods for diffusion-convection-reaction problems. *SIAM J. Numer. Anal.* **44** (2006), 949-978. https://doi.org/10.1137/040611483.

[198] Nochetto, R. H., Siebert, K. G., and Veeser, A. Theory of adaptive finite element methods: an introduction. In *Multiscale, nonlinear and adaptive approximation*. Springer, Berlin, 2009, pp. 409-542. http://dx.doi.org/10.1007/978-3-642-03413-8_12.

[199] Nordbotten, J. M., and Keilegavlen, E. An introduction to multi-point flux (MPFA) and stress (MPSA) finite volume methods for thermo-poroelasticity. In *Polyhedral methods in geosciences*, vol. **27** of *SEMA SIMAI Springer Ser*. Springer, Cham, [2021] ©2021, pp. 119-158. https://doi.org/10.1007/978-3-030-69363-3_4.

[200] Nordbotten, J. M., Keilegavlen, E., and Sandvin, A. Mass conservative domain decomposition for porous media flow. In *Finite Volume Method - Powerful Means of Engineering Design*, R. Petrova, Ed. IntechOpen, Rijeka, 2012, ch. 11. https://doi.org/10.5772/38700.

[201] Odeh, A. S. Comparison of solutions to a three-dimensional black-oil reservoir simulation problem (includes associated paper 9741). *J. Pet. Technol.* **33** (1981), 13-25. https://doi.org/10.2118/9723-PA.

[202] Ohlberger, M. A posteriori error estimate for finite volume approximations to singularly perturbed nonlinear convection-diffusion equations. *Numer. Math.* **87** (2001), 737-761. https://doi.org/10.1007/PL00005431.

[203] Ohlberger, M. A posteriori error estimates for vertex centered finite volume approximations of convection-diffusion-reaction equations. *M2AN Math. Model. Numer. Anal.* **35** (2001), 355-387. https://doi.org/10.1051/m2an:2001119.

[204] Ohlberger, M., and Rohde, C. Adaptive finite volume approximations for weakly coupled convection dominated parabolic systems. *IMA J. Numer. Anal.* **22** (2002), 253-280. https://doi.org/10.1093/imanum/22.2.253.

[205] Olshanskii, M. A., and Tyrtyshnikov, E. E. Theory and applications. *Iterative methods for linear systems*. Society for Industrial and Applied Mathematics, Philadelphia, PA, 2014. http://dx.doi.org/10.1137/1.9781611973464.

[206] Omnes, P., Penel, Y., and Rosenbaum, Y. A posteriori error estimation for the discrete duality finite volume discretization of the Laplace equation. *SIAM J. Numer. Anal.* **47** (2009), 2782-2807. https://doi.org/10.1137/080735047.

[207] Otto, F. $L^1$-contraction and uniqueness for quasilinear elliptic-parabolic equations. *J. Differential Equations* **131** (1996), 20-38. http://dx.doi.org/10.1006/jdeq.1996.0155.

[208] Otto, F. $L^1$-contraction and uniqueness for unstationary saturated-unsaturated porous media flow. *Adv. Math. Sci. Appl.* **7** (1997), 537-553.
[209] Papež, J. Algebraic error in numerical PDEs and its estimation. In *Error control, adaptive discretizations, and applications. Part 1*, vol. 58 of Adv. Appl. Mech. Academic Press, San Diego, CA, 2024, pp. 377–427.

[210] Papež, J., Rüde, U., Vohralík, M., and Wohlmuth, B. Sharp algebraic and total a posteriori error bounds for *h* and *p* finite elements via a multilevel approach. Recovering mass balance in any situation. *Comput. Methods Appl. Mech. Engrg.* **371** (2020), 113243. https://doi.org/10.1016/j.cma.2020.113243.

[211] Patera, A. T., and Rønquist, E. M. A general output bound result: application to discretization and iteration error estimation and control. *Math. Models Methods Appl. Sci.* **11** (2001), 685–712. http://dx.doi.org/10.1142/S0218202501001057.

[212] Payne, L. E., and Weinberger, H. F. An optimal Poincaré inequality for convex domains. *Arch. Rational Mech. Anal.* **5** (1960), 286–292.

[213] Pencheva, G. V., Vohralík, M., Wheeler, M. F., and Wildey, T. Robust a posteriori error control and adaptivity for multiscale, multinumerics, and mortar coupling. *SIAM J. Numer. Anal.* **51** (2013), 526–554. http://dx.doi.org/10.1137/110839047.

[214] Picasso, M. Adaptive finite elements for a linear parabolic problem. *Comput. Methods Appl. Mech. Engrg.* **167** (1998), 223–237. http://dx.doi.org/10.1016/S0045-7825(98)00121-2.

[215] Picasso, M. A stopping criterion for the conjugate gradient algorithm in the framework of anisotropic adaptive finite elements. *Comm. Numer. Methods Engrg.* **25** (2009), 339–355. https://doi.org/10.1002/cnm.1120.

[216] Prager, W., and Synge, J. L. Approximations in elasticity based on the concept of function space. *Quart. Appl. Math.* **5** (1947), 241–269.

[217] Rannacher, R., and Vihharev, J. Adaptive finite element analysis of nonlinear problems: balancing of discretization and iteration errors. *J. Numer. Math.* **21** (2013), 23–61. https://doi.org/10.1515/jnum-2013-0002.

[218] Raviart, P.-A., and Thomas, J.-M. A mixed finite element method for 2nd order elliptic problems. In *Mathematical aspects of finite element methods (Proc. Conf., Consiglio Naz. delle Ricerche (C.N.R.), Rome, 1975)*. Springer, Berlin, 1977, pp. 292–315. Lecture Notes in Math., Vol. 606.

[219] Repin, S. *A posteriori estimates for partial differential equations*, vol. 4 of Radon Series on Computational and Applied Mathematics. Walter de Gruyter GmbH & Co. KG, Berlin, 2008. http://dx.doi.org/10.1515/9783110203042.

[220] Rey, V., Gosselet, P., and Rey, C. Strict bounding of quantities of interest in computations based on domain decomposition. *Comput. Methods Appl. Mech. Engrg.* **287** (2015), 212–228. http://dx.doi.org/10.1016/j.cma.2015.01.009.

[221] Rey, V., Gosselet, P., and Rey, C. Strict lower bounds with separation of sources of error in non-overlapping domain decomposition methods. *Internat. J. Numer. Methods Engrg.* **108** (2016), 1007–1029. http://dx.doi.org/10.1002/nme.5244.

[222] Rey, V., Rey, C., and Gosselet, P. A strict error bound with separated contributions of the discretization and of the iterative solver in non-overlapping domain decomposition methods. *Comput. Methods Appl. Mech. Engrg.* **270** (2014), 293–303. http://dx.doi.org/10.1016/j.cma.2013.12.001.

[223] Ricois, O. Vision générale du simulateur ARCEOR. Tech. rep., IFPEN, 2011. Number 62096.

[224] Roberts, J. E., and Thomas, J.-M. Mixed and hybrid methods. In *Handbook of Numerical Analysis, Vol. II*. North-Holland, Amsterdam, 1991, pp. 523–639. https://doi.org/10.1016/S1570-8659(05)80041-9.

[225] Rose, M. E. Compact finite volume methods for the diffusion equation. *J. Sci. Comput.* **4** (1989), 261–290. https://doi.org/10.1007/BF01061058.
[226] Ruiz-Baier, R., and Lunati, I. Mixed finite element-discontinuous finite volume element discretization of a general class of multicontinuum models. *J. Comput. Phys.* **322** (2016), 666–688. https://doi.org/10.1016/j.jcp.2016.06.054.

[227] Russell, T. F., and Wheeler, M. F. Finite element and finite difference methods for continuous flows in porous media. In *The Mathematics of Reservoir Simulation*. SIAM, Philadelphia, 1983, pp. 35–106. https://doi.org/10.1137/1.9781611971071.ch2.

[228] Saad, M., and Zhang, H. Front tracking for two-phase flow in reservoir simulation by adaptive mesh. *Numer. Methods Partial Differential Equations* **13** (1997), 673–697. http://dx.doi.org/10.1002/(SICI)1098-2426(199711)13:6<673::AID-NUM5>3.0.CO;2-0.

[229] Saad, Y. *Iterative methods for sparse linear systems*, second ed. Society for Industrial and Applied Mathematics, Philadelphia, PA, 2003. https://doi.org/10.1137/1.9780898718003.

[230] Sboui, A., Jaffré, J., and Roberts, J. A composite mixed finite element for hexahedral grids. *SIAM J. Sci. Comput.* **31** (2009), 2623–2645. http://dx.doi.org/10.1137/070703703.

[231] Smaï, F. Existence of solutions for a model of multiphase flow in porous media applied to gas migration in underground nuclear waste repository. *Appl. Anal.* **88** (2009), 1609–1616. https://doi.org/10.1080/00036810902942226.

[232] Smears, I. An introduction to the a posteriori error analysis of parabolic partial differential equations. In *Error control, adaptive discretizations, and applications. Part 4*, vol. **61** of *Adv. Appl. Mech*. Academic Press, San Diego, CA, 2025. In press.

[233] Stokke, J. S., Mitra, K., Storvik, E., Both, J. W., and Radu, F. A. An adaptive solution strategy for Richards’ equation. *Comput. Math. Appl.* **152** (2023), 155–167. https://doi.org/10.1016/j.camwa.2023.10.020.

[234] Süli, E. A posteriori error analysis and global error control for adaptive finite volume approximations of hyperbolic problems. In *Numerical analysis 1995 (Dundee, 1995)*, vol. **344** of *Pitman Res. Notes Math. Ser*. Longman, Harlow, 1996, pp. 169–190.

[235] Süli, E. A posteriori error analysis and adaptivity for finite element approximations of hyperbolic problems. In *An introduction to recent developments in theory and numerics for conservation laws (Freiburg/Littenweiler, 1997)*, vol. **5** of *Lect. Notes Comput. Sci. Eng*. Springer, Berlin, 1999, pp. 123–194. https://doi.org/10.1007/978-3-642-58535-7_4.

[236] Synge, J. L. *The hypercircle in mathematical physics: a method for the approximate solution of boundary value problems*. Cambridge University Press, New York, 1957.

[237] Thomas, J.-M. *Sur l'analyse numérique des méthodes d'éléments finis hybrides et mixtes*. Ph.D. dissertation, Université Pierre et Marie Curie (Paris 6), May 1977.

[238] Tobiska, L., and Verfürth, R. Robust a posteriori error estimates for stabilized finite element methods. *IMA J. Numer. Anal.* **35** (2015), 1652–1671. https://doi.org/10.1093/imanum/dru060.

[239] van der Vorst, H. A. Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems. *SIAM J. Sci. Statist. Comput.* **13** (1992), 631–644. https://doi.org/10.1137/0913035.

[240] Verfürth, R. A posteriori error estimators for convection-diffusion equations. *Numer. Math.* **80** (1998), 641–663. https://doi.org/10.1007/s002110050381.

[241] Verfürth, R. Robust a posteriori error estimators for a singularly perturbed reaction-diffusion equation. *Numer. Math.* **78** (1998), 479–493. https://doi.org/10.1007/s002110050322.

[242] Verfürth, R. A posteriori error estimates for finite element discretizations of the heat equation. *Calcolo* **40** (2003), 195–212. http://dx.doi.org/10.1007/s10092-003-0073-2.

[243] Verfürth, R. Robust a posteriori error estimates for nonstationary convection-diffusion equations. *SIAM J. Numer. Anal.* **43** (2005), 1783–1802. http://dx.doi.org/10.1137/040604273.
[244] Verfürth, R. Robust a posteriori error estimates for stationary convection-diffusion equations. *SIAM J. Numer. Anal.* **43** (2005), 1766–1782. http://dx.doi.org/10.1137/040604261.

[245] Verfürth, R. *A posteriori error estimation techniques for finite element methods*. Numerical Mathematics and Scientific Computation. Oxford University Press, Oxford, 2013. http://dx.doi.org/10.1093/acprof:oso/9780199679423.001.0001.

[246] Vohralík, M. Equivalence between lowest-order mixed finite element and multi-point finite volume methods on simplicial meshes. *M2AN Math. Model. Numer. Anal.* **40** (2006), 367–391. http://dx.doi.org/10.1051/m2an:2006013.

[247] Vohralík, M. A posteriori error estimates for lowest-order mixed finite element discretizations of convection-diffusion-reaction equations. *SIAM J. Numer. Anal.* **45** (2007), 1570–1599. http://dx.doi.org/10.1137/060653184.

[248] Vohralík, M. A posteriori error estimation in the conforming finite element method based on its local conservativity and using local minimization. *C. R. Math. Acad. Sci. Paris* **346** (2008), 687–690. http://dx.doi.org/10.1016/j.crma.2008.03.006.

[249] Vohralík, M. Residual flux-based a posteriori error estimates for finite volume and related locally conservative methods. *Numer. Math.* **111** (2008), 121–158. http://dx.doi.org/10.1007/s00211-008-0168-4.

[250] Vohralík, M. Unified primal formulation-based a priori and a posteriori error analysis of mixed finite element methods. *Math. Comp.* **79** (2010), 2001–2032. http://dx.doi.org/10.1090/S0025-5718-2010-02375-0.

[251] Vohralík, M. Lecture notes. *A posteriori error estimates for efficiency and error control in numerical simulations*. Charles University, Prague, 2024. https://who.rocq.inria.fr/Martin.Vohralik/Enseig/APost/a_posteriori.pdf.

[252] Vohralík, M. Lecture notes. *Analysis and approximation of partial differential equations by finite elements*. ENSTA, Institut Polytechnique de Paris, 2025. https://who.rocq.inria.fr/Martin.Vohralik/Enseig/FEM/FEM.pdf.

[253] Vohralík, M., and Wheeler, M. F. A posteriori error estimates, stopping criteria, and adaptivity for two-phase flows. *Comput. Geosci.* **17** (2013), 789–812. http://dx.doi.org/10.1007/s10596-013-9356-0.

[254] Vohralík, M., and Wohlmuth, B. I. Mixed finite element methods: implementation with one unknown per element, local flux expressions, positivity, polygonal meshes, and relations to other methods. *Math. Models Methods Appl. Sci.* **23** (2013), 803–838. https://doi.org/10.1142/S0218202512500613.

[255] Vohralík, M., and Yousef, S. A simple a posteriori estimate on general polytopal meshes with applications to complex porous media flows. *Comput. Methods Appl. Mech. Engrg.* **331** (2018), 728–760. https://doi.org/10.1016/j.cma.2017.11.027.

[256] Wheeler, M., Xue, G., and Yotov, I. A multipoint flux mixed finite element method on distorted quadrilaterals and hexahedra. *Numer. Math.* **121** (2012), 165–204. http://dx.doi.org/10.1007/s00211-011-0427-7.

[257] Wheeler, M. F., and Yotov, I. A multipoint flux mixed finite element method. *SIAM J. Numer. Anal.* **44** (2006), 2082–2106. http://dx.doi.org/10.1137/050638473.

[258] Whitney, H. *Geometric integration theory*. Princeton University Press, Princeton, NJ, 1957.

[259] Younès, A., Ackerer, P., and Chavent, G. From mixed finite elements to finite volumes for elliptic PDEs in two and three dimensions. *Internat. J. Numer. Methods Engrg.* **59** (2004), 365–388. https://doi.org/10.1002/nme.874.

[260] Younès, A., Mosé, R., Ackerer, P., and Chavent, G. A new formulation of the mixed finite element method for solving elliptic and parabolic PDE with triangular elements. *J. Comput. Phys.* **149** (1999), 148–167. https://doi.org/10.1006/jcph.1998.6150.
[261] Young, L. C., and Stephenson, R. E. A generalized compositional approach for reservoir simulation. *Society of Petroleum Engineers Journal* **23** (1983), 727-742. https://doi.org/10.2118/10516-PA.

[262] Yousef, S. *A posteriori error estimates and adaptivity based on stopping criteria and adaptive mesh refinement for multiphase and thermal flows*. Application to steam-assisted gravity drainage. Ph.D. thesis, Université Pierre et Marie Curie - Paris VI, 2013. https://tel.archives-ouvertes.fr/tel-00918782.

[263] Zeidler, E. Nonlinear monotone operators, Translated from the German by the author and Leo F. Boron. *Nonlinear functional analysis and its applications*. II/B. Springer-Verlag, New York, 1990. https://doi.org/10.1007/978-1-4612-0981-2.