Martin Vohralík
email: martin.vohralik@inria.fr

Mary F Wheeler

‡ M F Wheeler

A posteriori error estimates, stopping criteria, and adaptivity for two-phase flows *

Keywords: two-phase flow, a posteriori error estimate, general framework, discretization error, linearization error, iterative coupling error, algebraic solver error

This paper develops a general abstract framework for a posteriori estimates for immiscible incompressible two-phase flows in porous media. We measure the error by the dual norm of the residual and, for mathematical correctness, employ the concept of global and complementary pressures in the analysis. Our estimators allow to estimate separately the different error components, namely the spatial discretization error, the temporal discretization error, the linearization error, the iterative coupling error, and the algebraic solver error. We propose an adaptive algorithm wherein the different iterative procedures (iterative linearization, iterative coupling, iterative solution of linear systems) are stopped when the corresponding errors do not affect significantly the overall error, and wherein the spatial and temporal errors are equilibrated. Consequently, important computational savings can be achieved while guaranteeing a user-given precision. The developed framework covers fully implicit, implicit pressureexplicit saturation, or iterative coupling formulations; conforming spatial discretization schemes such as the vertex-centered finite volume method or the finite element method, and nonconforming spatial discretization schemes such as the cell-centered finite volume method, the mixed finite element method, or the discontinuous Galerkin method; linearizations such as the Newton or the fixed-point one; and general linear solvers. Numerical experiments for a model problem are presented to illustrate the theoretical results. Only by stopping timely the linear and nonlinear solvers, speed-ups by a factor between 10 and 20 in terms of the number of total linear solver iterations are achieved.

Introduction

Let an open bounded polygonal (polyhedral) domain Ω ⊂ R d , d = 2, 3, and a time interval (0, T ), T > 0, be given and set Q := Ω × (0, T ). We consider the immiscible incompressible two-phase flow in porous media in the space-time cylinder Q in the form

∂ t (φs α ) -∇• k r,α (s w ) µ α K(∇p α + ρ α g∇z) = q α , (1.1a)
s n + s w = 1, (1.1b) p np w = p c (s w ).

(1.1c)

Here the unknowns are s α , the phase saturations, and p α , the phase pressures, α ∈ {n, w}. The subscripts n, w stand for nonwetting and wetting, respectively. Typically, the nonwetting phase is oil and the wetting one is water. For the sake of simplicity, we suppose that the porosity φ, as well as the phase viscosities µ α and the phase densities ρ α are all constant. The permeability tensor K and the phase sources q α , α ∈ {n, w}, are only supposed to depend on the space coordinate x and on the time t. For the sake of simplicity, we suppose q α piecewise constant in time on time mesh defined below. In (1.1a)-(1.1c), z stands for the vertical coordinate and g for the gravitation acceleration constant. The system (1.1a)-(1.1c) is nonlinear and coupled because of the presence of p c , the capillary pressure, and of k r,α , the phase relative permeabilities, which are all given functions of the wetting phase saturation s w . For example, in the Brooks-Corey [START_REF] Brooks | Hydraulic properties of porous media[END_REF] model, Here p d is the entry pressure and s rw and s rn are respectively the wetting and nonwetting residual saturations. System (1.1a)-(1.1c) is degenerate as the phase relative permeabilities k r,α can become zero. Define the phase Darcy velocities u α , α ∈ {n, w}, by

u α := - k r,α (s w ) µ α K(∇p α + ρ α g∇z). (1.2)
The initial condition is imposed through

s w (•, 0) = s 0 w in Ω. (1.3)
For the sake of simplicity, we suppose homogeneous Neumann and Dirichlet boundary conditions

u α •n Ω = 0 on Γ N × (0, T ), α ∈ {n, w}, (1.4a)
s w = p w = 0 on Γ D × (0, T ).

(1.4b)

Here ∂Ω = Γ N ∪ Γ D , Γ N ∩ Γ D = ∅, both Γ N and Γ D are simply connected, and Γ D is of nonzero measure.

Conditions (1.4) can be replaced by more realistic ones. Problem (1.1a)-(1.4b) is of fundamental importance in petroleum engineering. Many results on this problem and on its numerical approximation have been derived in the past. The analysis of (1.1a)-(1.4b) including the existence, uniqueness, and well-posedness results has been carried out in [START_REF] Kröner | Flow of oil and water in a porous medium[END_REF][START_REF] Chavent | Mathematical models and finite elements for reservoir simulation[END_REF][START_REF] Antontsev | Boundary value problems in mechanics of nonhomogeneous fluids[END_REF][START_REF] Arbogast | The existence of weak solutions to single porosity and simple dual-porosity models of two-phase incompressible flow[END_REF][START_REF] Chen | Degenerate two-phase incompressible flow. I. Existence, uniqueness and regularity of a weak solution[END_REF][START_REF] Chen | Degenerate two-phase incompressible flow. II. Regularity, stability and stabilization[END_REF][START_REF] Chen | Degenerate two-phase incompressible flow. III. Sharp error estimates[END_REF][START_REF] Cancès | Two-phase flows involving capillary barriers in heterogeneous porous media[END_REF][START_REF] Khalil | Degenerate two-phase compressible immiscible flow in porous media: the case where the density of each phase depends on its own pressure[END_REF][START_REF] Amaziane | An existence result for a coupled system modeling a fully equivalent global pressure formulation for immiscible compressible two-phase flow in porous media[END_REF], see also [START_REF] Alt | Quasilinear elliptic-parabolic differential equations[END_REF][START_REF] Otto | L 1 -contraction and uniqueness for quasilinear elliptic-parabolic equations[END_REF][START_REF] Van Duijn | Effective equations for two-phase flow with trapping on the micro scale[END_REF] for degenerate problems, and the references therein. For the use and analysis of mixed finite element methods for the numerical approximation of (1.1a)-(1.4b) we refer to, e.g., [START_REF] Ewing | The approximation of the pressure by a mixed method in the simulation of miscible displacement[END_REF][START_REF] Arbogast | A nonlinear mixed finite element method for a degenerate parabolic equation arising in flow in porous media[END_REF][START_REF] Wheeler | Efficient algorithms for multiscale modeling in porous media[END_REF] and the references therein, for discontinuous Galerkin methods, to, e.g., [START_REF] Epshteyn | Analysis of hp discontinuous Galerkin methods for incompressible twophase flow[END_REF][START_REF] Ern | Accurate velocity reconstruction for discontinuous Galerkin approximations of two-phase porous media flows[END_REF][START_REF] Ern | Discontinuous Galerkin approximation of two-phase flows in heterogeneous porous media with discontinuous capillary pressures[END_REF][START_REF] Arbogast | A discontinuous Galerkin method for two-phase flow in a porous medium enforcing H(div) velocity and continuous capillary pressure[END_REF] and the references therein, for cell-centered finite volume methods to, e.g., [START_REF] Eymard | Mathematical study of a petroleum-engineering scheme[END_REF] and the references therein, and for vertexcentered finite volume methods to, e.g., [START_REF] Huber | Node-centered finite volume discretizations for the numerical simulation of multiphase flow in heterogeneous porous media[END_REF] and the references therein. Multiscale and mortar techniques, efficient parallelization, and multinumerics and multiphysics formulations have been investigated in [START_REF] Peszyńska | Mortar upscaling for multiphase flow in porous media[END_REF]. First adaptive mesh refinement strategies were prosed and tested in, e.g., [START_REF] Saad | Front tracking for two-phase flow in reservoir simulation by adaptive mesh[END_REF][START_REF] Chen | Degenerate two-phase incompressible flow. IV. Local refinement and domain decomposition[END_REF][START_REF] Klieber | Adaptive simulations of two-phase flow by discontinuous Galerkin methods[END_REF][START_REF] Chen | A posteriori error estimates of mixed methods for miscible displacement problems[END_REF][START_REF] Chueh | Multi-level adaptive simulation of transient two-phase flow in heterogeneous porous media[END_REF]. Linearization, linear solver techniques, and stopping criteria for multiphase flows are discussed in, e.g., [START_REF] Wallis | Constrained residual acceleration of conjugate residual methods[END_REF][START_REF] Lacroix | Decoupling preconditioners in the implicit parallel accurate reservoir simulator (IPARS)[END_REF][START_REF] Lacroix | Iterative solution methods for modeling multiphase flow in porous media fully implicitly[END_REF].

The purpose of the present paper is to derive a posteriori estimates for numerical approximations of the problem (1.1a)-(1.4b). Our estimates give a guaranteed and fully and easily computable upper bound on the selected error measure, the dual norm of the residual augmented by the distance of the approximate global and complementary pressures to proper function spaces. Recall that such error measure leads to the energy error for linear problems (cf., e.g., [START_REF] Ern | A posteriori error estimation based on potential and flux reconstruction for the heat equation[END_REF]), and it is shown in [START_REF] Cancès | An a posteriori error estimate for vertex-centered finite volume discretizations of immiscible incompressible two-phase flow[END_REF] that this is an upper bound on the error between the exact and approximate saturations, global pressures, and complementary pressures for conforming discretizations. Our estimates also allow to distinguish, estimate separately, and compare different error components. The principal error component is the discretization error, due to the numerical scheme chosen, the local space mesh size, and the local time mesh size. This can be decomposed into space discretization error and time discretization error. The subsidiary error component is the error due to various iterative procedures involved in the calculation. This includes linearization error, iterative coupling error, or linear solver error. We next devise adaptive algorithms where all the iterative procedures on a given time level are stopped whenever the individual errors drop to the level at which they do not affect significantly the overall error. Simultaneously, the space and time discretization errors are adjusted so that they are of similar size.

The benefits of such a procedure are twofold. Firstly, the overall error is controlled and strategies for obtaining a user-given final precision at the end of the simulation can be devised. Secondly, it is likely to lead to important computational savings, as performing an excessive number of unnecessary linearization/iterative coupling/linear solver iterations and using too fine (with respect to the other components of the error) space or time meshes can be avoided. These concepts have been known for long time in the engineering practice but only recently, rigorous mathematical analysis has been started in model cases. In particular, linear solver error estimation and linear solver stopping criteria have been developed in, e.g., [START_REF] Becker | Adaptive error control for multigrid finite element methods[END_REF][START_REF] Patera | A general output bound result: application to discretization and iteration error estimation and control[END_REF][START_REF] Arioli | Stopping criteria for iterations in finite element methods[END_REF], nonlinear solver error estimation and nonlinear solver stopping criteria are treated in, e.g., [START_REF] Han | A posteriori error analysis for linearization of nonlinear elliptic problems and their discretizations[END_REF][START_REF] Chaillou | Computable error estimators for the approximation of nonlinear problems by linearized models[END_REF], and spatial and temporal errors are estimated and balanced in, e.g., [START_REF] Picasso | Adaptive finite elements for a linear parabolic problem[END_REF][START_REF] Repin | Estimates of deviations from exact solutions of initial-boundary value problem for the heat equation[END_REF][START_REF] Verfürth | A posteriori error estimates for finite element discretizations of the heat equation[END_REF][START_REF] Makridakis | Elliptic reconstruction and a posteriori error estimates for parabolic problems[END_REF]. Inexact Newton methods have been studied in, e.g., [START_REF] Bank | Analysis of a multilevel iterative method for nonlinear finite element equations[END_REF][START_REF] Moret | A Kantorovich-type theorem for inexact Newton methods[END_REF][START_REF] Eisenstat | Globally convergent inexact Newton methods[END_REF][START_REF] Brown | Convergence theory of nonlinear Newton-Krylov algorithms[END_REF][START_REF] Eisenstat | Choosing the forcing terms in an inexact Newton method[END_REF][START_REF] Deuflhard | Global inexact Newton multilevel FEM for nonlinear elliptic problems[END_REF][START_REF] Deuflhard | Affine invariance and adaptive algorithms[END_REF]. Herein, we build upon the results of [START_REF] Jiránek | A posteriori error estimates including algebraic error and stopping criteria for iterative solvers[END_REF][START_REF] Ern | A posteriori error estimation based on potential and flux reconstruction for the heat equation[END_REF][START_REF] El Alaoui | Guaranteed and robust a posteriori error estimates and balancing discretization and linearization errors for monotone nonlinear problems[END_REF][START_REF] Hilhorst | A posteriori error estimates for combined finite volume-finite element discretizations of reactive transport equations on nonmatching grids[END_REF][START_REF] Ern | Adaptive inexact Newton methods with a posteriori stopping criteria for nonlinear diffusion PDEs[END_REF] which give guaranteed and robust a posteriori estimates with error components distinction.

The present paper gives a posteriori error estimates in a general setting without a specification of the underlying numerical treatment. Examples of the application of this abstract framework to different discrete formulations, spatial and temporal discretizations, linearizations, and linear solvers are given in Section 6, with some further examples in [START_REF] Di Pietro | An a posteriori error estimator for a finite volume discretization of the two-phase flow[END_REF][START_REF] Cancès | An a posteriori error estimate for vertex-centered finite volume discretizations of immiscible incompressible two-phase flow[END_REF]. In order to unify the presentation, we have chosen once and for all as the primary unknowns the saturation and pressure of the wetting phase. Adjustments to all other choices are easily possible.

Preliminaries

We specify here the notation and function spaces used, characterize the weak solution, give our assumptions on the approximate solutions, and define the error measure.

Function spaces and space and time meshes

We denote by H 1 (Ω) the Sobolev space of functions from L 2 (Ω) which admit a weak gradient in [L 2 (Ω)] d . The subspace H 1 D (Ω) of H 1 (Ω) contains functions with zero trace on Γ D . H 1 (Ω) and H 1 D (Ω) functions are continuous in the sense of traces. Let H(div, Ω) be the space of vector-valued functions from [L 2 (Ω)] d which admit a weak divergence in L 2 (Ω). H(div, Ω) functions have a continuous normal trace, so that H(div, Ω) is appropriate for representing mass-conservative vectors. Denote

X := L 2 ((0, T ); H 1 D (Ω)); (2.1)
for v ∈ X, we set

v X := T 0 ∇v 2 dt 1 2
and observe that X is the usual energy space for parabolic problems and that ϕ X is the associated energy norm. Below, (•, •) stands for the L 2 (Ω) scalar product and •, • for the duality pairing of (

H 1 D (Ω)) ′ and (H 1 D (Ω)).
We consider a strictly increasing sequence of discrete times {t n } 0≤n≤N such that t 0 = 0 and t N = T , together with a set of meshes {T n h } 0≤n≤N . For all 1 ≤ n ≤ N , we define the time interval I n := (t n-1 , t n ] and the time step τ n := t nt n-1 . For all 0 ≤ n ≤ N , we assume that T n h covers exactly Ω. The meshes T n h can be composed of general polygonal (polyhedral) elements. For all K ∈ T n h , h K denotes the diameter of the mesh element K. The discrete times and meshes can be constructed by a space-time adaptive time-marching algorithm such as those of Sections 4.1, 4.3 and 5.1, 5.2, 5.4 below.

Let 0 ≤ n ≤ N . We first define the broken Sobolev space H 1 (T n h ) as the space of all functions v ∈ L 2 (Ω) such that v| K ∈ H 1 (K) for all K ∈ T n h . The symbol ∇ denotes the corresponding broken gradient, i.e., a gradient of the function restricted to a mesh element K. We define P 1 τ (H 1 (T )) as the space of functions v continuous and piecewise affine in time, given by v n ∈ H 1 (T n h ) for every discrete time t n , 0 ≤ n ≤ N , i.e.,

{v n = v(•, t n )} 0≤n≤N in H 1 (T n h
). We will use a similar notation P 1 τ (H 1 D (Ω)) for functions continuous and piecewise affine in time with values in H 1 D (Ω) and denote

∂ n t v := ∂ t v| In = 1 τ n (v n -v n-1 ). (2.2)
Similarly, P 0 τ (H(div, Ω)) stands for functions piecewise constant in time with values in H(div, Ω); for v ∈ P 0 τ (H(div, Ω)), we set

v n := v(•, t)| I n .

Weak solution definition via the global and complementary pressures

In order to characterize the error in an approximate solution to (1.1a)-(1.4b), we first need to define the weak solution of (1.1a)-(1.4b). Following [START_REF] Antontsev | Boundary value problems in mechanics of nonhomogeneous fluids[END_REF][START_REF] Arbogast | The existence of weak solutions to single porosity and simple dual-porosity models of two-phase incompressible flow[END_REF][START_REF] Chen | Degenerate two-phase incompressible flow. I. Existence, uniqueness and regularity of a weak solution[END_REF][START_REF] Chen | Degenerate two-phase incompressible flow. II. Regularity, stability and stabilization[END_REF][START_REF] Chen | Degenerate two-phase incompressible flow. III. Sharp error estimates[END_REF][START_REF] Amaziane | An existence result for a coupled system modeling a fully equivalent global pressure formulation for immiscible compressible two-phase flow in porous media[END_REF], we introduce the global and complementary pressures. We would like to stress that these mathematical quantities only appear in order to describe the weak solution and to measure errors but are not applied in the numerical schemes.

Let the phase mobilities be denoted by

λ α (a) := k r,α (a) µ α α ∈ {n, w}. (2.3)
We define the global pressure

p(s w , p w ) := p w + sw 0 λ n (a) λ w (a) + λ n (a) p ′ c (a)da (2.4a)
and the complementary pressure

q(s w ) := - sw 0 λ w (a)λ n (a) λ w (a) + λ n (a) p ′ c (a)da. (2.4b)
Next, in order to simplify the arguments below, we define the functions v α , α ∈ {n, w}, of wetting saturation and pressure (s w , p w ), by v w (s w , p w ) := -K λ w (s w )∇p(s w , p w ) + ∇q(s w )

+ λ w (s w )ρ w g∇z , (2.5a)

v n (s w , p w ) := -K λ n (s w )∇p(s w , p w ) -∇q(s w ) + λ n (s w )ρ n g∇z . (2.5b)
Note that v α (s w , p w ) are formally equivalent to the phase velocities u α given by (1.2) and are always well defined.

We suppose that the data are regular enough so that the weak solution (s w , p w ) to (1.1a)-(1.4b), setting s n := 1s w , can be characterized by

s w ∈ C([0, T ]; L 2 (Ω)), s w (•, 0) = s 0 w , ∂ t s w ∈ L 2 ((0, T ); (H 1 D (Ω)) ′ ), (2.6a) p(s w , p w ) ∈ X, (2.6b) q(s w ) ∈ X, (2.6c) T 0 ∂ t (φs α ), ϕ -(v α (s w , p w ), ∇ϕ) -(q α , ϕ) dt = 0 ∀ϕ ∈ X, α ∈ {n, w}. (2.6d)
We refer to [START_REF] Chen | Degenerate two-phase incompressible flow. I. Existence, uniqueness and regularity of a weak solution[END_REF][START_REF] Eymard | Mathematical study of a petroleum-engineering scheme[END_REF][START_REF] Cancès | An a posteriori error estimate for vertex-centered finite volume discretizations of immiscible incompressible two-phase flow[END_REF] for details.

Approximate saturations and pressures

Our a posteriori error estimates will be given for general approximate wetting saturations s w,hτ and general approximate wetting pressures p w,hτ , not linked to any particular numerical scheme. More precisely, recalling the definition of the space P 1 τ (H 1 (T )) from Section 2.1, we merely require s w,hτ , p w,hτ ∈ P 1 τ (H 1 (T )). Thus, cases where s w,hτ and p w,hτ are nonconforming in the sense that p(s w,hτ , p w,hτ ) ∈ X and q(s w,hτ ) ∈ X are included. In general, the notation v hτ stands for a space-time function continuous and piecewise affine in time and piecewise polynomial in space on the meshes T n h and v n h := v hτ (•, t n ), 0 ≤ n ≤ N , for the piecewise polynomial in space. We also assume for simplicity that the initial and Dirichlet boundary conditions are satisfied exactly, i.e. s 0 w,h = s 0 w and s w,hτ = p w,hτ = 0 on Γ D × (0, T ).

Error measure

The first question in a posteriori error estimates is that of the error "measure". In linear problems, one usually chooses the energy norm for a global error measure. In nonlinear problems, the situation is more difficult. One approach consists in taking the dual norm of the residual, i.e., of the difference of the nonlinear operator applied to the exact and approximate solutions, cf. [START_REF] Chaillou | Computable error estimators for the approximation of nonlinear problems by linearized models[END_REF][START_REF] El Alaoui | Guaranteed and robust a posteriori error estimates and balancing discretization and linearization errors for monotone nonlinear problems[END_REF][START_REF] Dolejší | A framework for robust a posteriori error control in unsteady nonlinear advection-diffusion problems[END_REF][START_REF] Ern | Adaptive inexact Newton methods with a posteriori stopping criteria for nonlinear diffusion PDEs[END_REF]. We also refer to [START_REF] Angermann | Balanced a posteriori error estimates for finite-volume type discretizations of convection-dominated elliptic problems[END_REF][START_REF] Verfürth | Robust a posteriori error estimates for nonstationary convection-diffusion equations[END_REF] for the use of dual norms in singularly perturbed linear problems. The advantage is that such a measure is dictated by the problem at hand; it simplifies the analysis and leads to sharper (and possibly robust, as in [START_REF] Vohralík | Guaranteed and fully robust a posteriori error estimates for conforming discretizations of diffusion problems with discontinuous coefficients[END_REF][START_REF] El Alaoui | Guaranteed and robust a posteriori error estimates and balancing discretization and linearization errors for monotone nonlinear problems[END_REF][START_REF] Dolejší | A framework for robust a posteriori error control in unsteady nonlinear advection-diffusion problems[END_REF][START_REF] Ern | Adaptive inexact Newton methods with a posteriori stopping criteria for nonlinear diffusion PDEs[END_REF]) estimates. Let s w , p w and s w,hτ , p w,hτ be respectively the exact and approximate wetting saturations and pressures as described in Sections 2.2 and 2.3. Let s n := 1s w and s n,hτ := 1s w,hτ . We define the error measure by

|||(s w -s w,hτ , p w -p w,hτ )||| := α∈{n,w} sup ϕ∈X, ϕ X =1 T 0 ∂ t (φs α ) -∂ t (φs α,hτ ), ϕ -(v α (s w , p w ) -v α (s w,hτ , p w,hτ ), ∇ϕ) dt 2 1 2 + inf p∈X T 0 K(λ w (s w,hτ ) + λ n (s w,hτ )) ∇(p(s w,hτ , p w,hτ ) -p) 2 dt 1 2 + inf q∈X T 0 K∇(q(s w,hτ ) -q) 2 dt 1 2 . (2.7)
For 1 ≤ n ≤ N , a local-in-time version on the time interval I n , consisting in replacing in (2.7) the time integrals T 0 by In and the space X by X| In , is denoted by |||(s ws w,hτ , p wp w,hτ )||| In . The first term of the error measure (2.7) represents the dual norm of the residual; for p(s w,hτ , p w,hτ ) ∈ X and q(s w,hτ ) ∈ X, it equals to zero if and only if s w,hτ coincides with s w and p w,hτ with p w , whenever the weak solution is unique. The second and third terms measure the nonconformity; the terms K(λ w (s w,hτ ) + λ n (s w,hτ )) and K in front of the broken gradients represent weights with appropriate physical units and are deduced from (2.5a)-(2.5b).

A general a posteriori error estimate

We present here a general a posteriori error estimate giving a guaranteed and fully computable upper bound on the error measure (2.7).

Pressures and velocities reconstructions

In order to proceed generally, without any specification of the numerical treatment used, we now make the following assumption: Assumption 3.1 (Reconstructions). We assume that there exist scalar functions phτ and qhτ and vector functions u α,hτ , α ∈ {n, w}, such that phτ , qhτ ∈ X, u α,hτ ∈ P 0 τ (H(div, Ω)), and u α,hτ •n Ω = 0 on Γ N × (0, T ), α ∈ {n, w}. Moreover, we suppose that u α,hτ satisfy

(q n α -∂ n t (φs α,hτ ) -∇•u n α,h , 1) K = 0 ∀K ∈ T n h , ∀1 ≤ n ≤ N, α ∈ {n, w}. (3.1)
We will call phτ the global pressure, qhτ the complementary pressure, and u α,hτ , α ∈ {n, w}, the phase velocities reconstructions, respectively.

Remark 3.2 (Pressures and velocities reconstructions). In the continuous setting, the global pressure p(s w , p w ) and the complementary pressure q(s w ) belong to the space X, see (2.6b)-(2.6c). Furthermore, the normal traces of the phase velocities v α (s w , p w ) (or u α ), α ∈ {n, w}, are continuous and the equilibrium condition

∂ t (φs α ) + ∇•u α = q α (3.2)
holds, cf. (1.1a) and (1.2), which expresses the mass balance and local conservativity for the phase fluxes. These properties are not necessarily satisfied at a discrete level. The pressure reconstructions phτ , qhτ and the velocity reconstructions u α,hτ , α ∈ {n, w}, of Assumption 3.1 restore the properties of the continuous level at the discrete one.

A posteriori error estimate

We now describe our estimators. Let a time step n, 1 ≤ n ≤ N , and a mesh element K ∈ T n h be given. Recall first the Poincaré inequality:

ϕ -ϕ K K ≤ C P,K h K ∇ϕ K ∀ϕ ∈ H 1 (K), (3.3)
where ϕ K is the mean value of the function ϕ on the element K and C P,K = 1/π whenever the element K is convex [START_REF] Payne | An optimal Poincaré inequality for convex domains[END_REF][START_REF] Bebendorf | A note on the Poincaré inequality for convex domains[END_REF]. For α ∈ {n, w}, define the residual estimators

η n R,K,α := C P,K h K q n α -∂ n t (φs α,hτ ) -∇•u n α,h K (3.4a)
and the flux estimators

η n F,K,α (t) := u n α,h -v α (s w,hτ , p w,hτ )(t) K . (3.4b)
The nonconformity estimators are given by

η n NC,K,1 (t) := K(λ w (s w,hτ ) + λ n (s w,hτ )) ∇(p(s w,hτ , p w,hτ ) -phτ ) K (t), (3.4c)
η n NC,K,2 (t) := K∇(q(s w,hτ ) -qhτ ) K (t). (3.4d)
We then have the following result:

Theorem 3.3 (A posteriori
+ N n=1 In K∈T n h η n NC,K,1 (t) 2 dt 1 2 + N n=1 In K∈T n h η n NC,K,2 (t) 2 dt 1 2
.

(3.5)

Proof. The proof is straightforward using the definition of error measure (2.7) and Assumption 3.1. The second and third terms in (3.5) follow immediately from (2.7). We thus only have to establish that the first term is an upper bound on the corresponding term in (2.7). Let α ∈ {n, w} and ϕ ∈ X with ϕ X = 1 be given. Then (2.6d) implies that

T 0 ∂ t (φs α ) -∂ t (φs α,hτ ), ϕ -(v α (s w , p w ) -v α (s w,hτ , p w,hτ ), ∇ϕ) dt = T 0 (q α -∂ t (φs α,hτ ), ϕ) + (v α (s w,hτ , p w,hτ ), ∇ϕ) dt.
Let now 1 ≤ n ≤ N be given. Adding and subtracting (u n α,h , ∇ϕ), applying the Green theorem, the assumption (3.1), the Poincaré inequality (3.3), and the Cauchy-Schwarz inequality, (q n α -∂ n t (φs α,hτ ), ϕ) + (v α (s w,hτ , p w,hτ ), ∇ϕ)

= (q n α -∂ n t (φs α,hτ ) -∇•u n α,h , ϕ) + (v α (s w,hτ , p w,hτ ) -u n α,h , ∇ϕ) = K∈T n h (q n α -∂ n t (φs α,hτ ) -∇•u n α,h , ϕ -ϕ K ) K + (v α (s w,hτ , p w,hτ ) -u n α,h , ∇ϕ) K ≤ K∈T n h (η n R,K,α + η n F,K,α (t)) ∇ϕ K . Thus N n=1 In ∂ t (φs α ) -∂ n t (φs α,hτ ), ϕ -(v α (s w , p w ) -v α (s w,hτ , p w,hτ ), ∇ϕ) dt ≤ N n=1 In K∈T n h (η n R,K,α + η n F,K,α (t)) ∇ϕ K dt.
The theorem follows by the Cauchy-Schwarz inequality and using ϕ X = 1.

Concept of application to different numerical methods

For the theoretical analysis of this paper, we only need Assumption 3.1. The practical application of the present framework to different numerical methods consists in specifying the construction of phτ , qhτ , and u α,hτ , α ∈ {n, w} that we outline now.

In conforming methods one obtains p(s w,hτ , p w,hτ ) ∈ X and q(s w,hτ ) ∈ X, so that we can immediately set phτ := p(s w,hτ , p w,hτ ) and qhτ := q(s w,hτ ). In nonconforming numerical methods, we typically choose phτ ∈ P

1 τ (H 1 D (Ω)), pn h := I av (p(s n w,h , p n w,h )), 0 ≤ n ≤ N , and qhτ ∈ P 1 τ (H 1 D (Ω)), qn h := I av (q(s n w,h )), 0 ≤ n ≤ N .
Here I av is a postprocessing averaging operator which sets the Lagrangian degrees of freedom inside Ω to the average of the values from the different elements sharing this degree of freedom, see [START_REF] Achdou | A priori and a posteriori analysis of finite volume discretizations of Darcy's equations[END_REF][START_REF] Karakashian | A posteriori error estimates for a discontinuous Galerkin approximation of second-order elliptic problems[END_REF].

The choice of u α,hτ , α ∈ {n, w}, is more involved. In mixed finite element methods, in addition to the approximate wetting saturation s w,hτ and pressure p w,hτ described in Section 2.3, one also directly obtains phase velocity approximations u α,hτ ∈ P 0 τ (H(div, Ω)), α ∈ {n, w}, satisfying (3.1). More precisely, for every time interval

I n , 1 ≤ n ≤ N , typically u n α,h ∈ RTN(T n h )
, where RTN(T n h ) is the Raviart-Thomas-Nédélec finite-dimensional subspace of H(div, Ω), cf. [START_REF] Brezzi | Mixed and hybrid finite element methods[END_REF]. In other numerical methods, obtaining u α,hτ ∈ P 0 τ (H(div, Ω)) satisfying (3.1) is possible by means of local postprocessing. In the context of linear elliptic equations, we refer the reader to [START_REF] Eymard | Finite volume approximation of elliptic problems and convergence of an approximate gradient[END_REF][START_REF] Vohralík | Residual flux-based a posteriori error estimates for finite volume and related locally conservative methods[END_REF] for cell-centered finite volume methods, to [START_REF] Kim | A posteriori error estimators for locally conservative methods of nonlinear elliptic problems[END_REF][START_REF] Ainsworth | A posteriori error estimation for discontinuous Galerkin finite element approximation[END_REF][START_REF] Ern | An accurate H(div) flux reconstruction for discontinuous Galerkin approximations of elliptic problems[END_REF] for discontinuous Galerkin methods, and to [START_REF] Destuynder | Explicit error bounds in a conforming finite element method[END_REF][START_REF] Luce | A local a posteriori error estimator based on equilibrated fluxes[END_REF][START_REF] Braess | Equilibrated residual error estimator for edge elements[END_REF][START_REF] Vohralík | Guaranteed and fully robust a posteriori error estimates for conforming discretizations of diffusion problems with discontinuous coefficients[END_REF] for vertex-centered finite volume and finite element methods. For nonlinear elliptic equations, such constructions are unified for different numerical methods in [START_REF] Ern | Adaptive inexact Newton methods with a posteriori stopping criteria for nonlinear diffusion PDEs[END_REF]. In the context of two-phase flows, the constructions of u α,hτ , α ∈ {n, w}, can be found in [START_REF] Di Pietro | An a posteriori error estimator for a finite volume discretization of the two-phase flow[END_REF] for cell-centered finite volume methods, in [START_REF] Cancès | An a posteriori error estimate for vertex-centered finite volume discretizations of immiscible incompressible two-phase flow[END_REF] for vertex-centered finite volume methods, and in [START_REF] Ern | Accurate velocity reconstruction for discontinuous Galerkin approximations of two-phase porous media flows[END_REF][START_REF] Ern | Discontinuous Galerkin approximation of two-phase flows in heterogeneous porous media with discontinuous capillary pressures[END_REF][START_REF] Arbogast | A discontinuous Galerkin method for two-phase flow in a porous medium enforcing H(div) velocity and continuous capillary pressure[END_REF] for the discontinuous Galerkin method.

Stopping criteria and adaptivity for fully implicit discretizations

In this section we show how the abstract a posteriori error estimate of Section 3 can be adopted to fully implicit discretizations of (1.1a)-(1.4b). We also describe how to take into account the additional error from iterative linearization and iterative solution of algebraic linear systems and distinguish the different error components. We finally propose stopping criteria for the various iterations and design a fully adaptive algorithm.

A fully implicit formulation

Keeping p w and s w as unknowns and expressing s n as a function of s w from (1.1b) and p n as a function of p w and s w from (1.1c), we arrive at the following equivalent form of (1.1a)-(1.1c):

∂ t (φs w ) -∇•(λ r,w (s w )K(∇p w + ρ w g∇z)) = q w , (4.1a) -∂ t (φs w ) -∇•(λ r,n (s w )K(∇(p w + p c (s w )) +ρ n g∇z)) = q n . (4.1b)
Let us now suppose some discretization of the above system in both space and time, starting from s 0 w,h ∈ H 1 (T 0 h ). We suppose implicit (backward Euler) discretization in time. This leads, on a time level n, 1 ≤ n ≤ N , to a system of nonlinear algebraic equations that can be schematically written in the form

SS n w SP n w SS n n SP n n S n w P n w = D n w D n n , (4.2)
where S n w is the algebraic vector of discrete unknowns corresponding to the wetting saturation s n w,h and P n w is the algebraic vector of discrete unknowns corresponding to the wetting pressure p n w,h . In practice, (4.2) is solved using an iterative linearization such as the Newton-Raphson method combined with an iterative algebraic solver, by means of the following algorithm:

1. Let the initial wetting saturation S 0 w (and pressure P 0 w ) be given. Set n = 1.

2. Set up the system of nonlinear algebraic equations leading to (4.2).

(a)

Choose some initial wetting saturation S n,0 w and pressure P n,0 w . Typically, these are the saturation and pressure from the last time step, S n-1

w and P n-1 w . Set k := 1.
(b) Set up the following linear system: find S n,k w and P n,k w , the solutions to ii. Perform a step of a chosen iterative algebraic solver for the linear system (4.3), starting from S n,k,i-1

SS n,k-1 w SP n,k-1 w SS n,k-1 n SP n,k-1 n S n,k w P n,k w = D n,k-1 w D n,k-1 n . (4.3) Here SS n,k-1 w , SP n,k-1 w , SS n,k-1
w and P n,k,i-1 w
. This gives approximations S n,k,i w and P n,k,i w . iii. Build piecewise polynomial representations of the wetting saturations and pressures s n,k,i

w,h ∈ H 1 (T n h ) and p n,k,i w,h ∈ H 1 (T n h ) from S n,k,i
w and P n,k,i w , according to the given numerical method. Define the space-time functions s n,k,i w,hτ and p n,k,i w,hτ ; these are affine in time on the time interval I n , given by s n-1 w,h and p n-1 w,h at time t n-1 and by s n,k,i w,h and p n,k,i w,h at time t n . We use this refined notation in place of s w,hτ , p w,hτ of Section 2.3. iv. This step only concerns nonconforming methods. From s n,k,i w,h and p n,k,i w,h , prescribe the reconstructions pn,k,i v. From the discretization, build the phase velocities reconstructions u n,k,i α,h ∈ RTN(T n h ), α ∈ {n, w} (cf. Assumption 3.1). More precisely, the goal is to obtain the decompositions, α ∈ {n, w},

h := I av (p(s n,k,i w,h , p n,k,i w,h )) and qn,k,i h := I av (q(s n,k,i w,h )).
u n,k,i α,h = d n,k,i α,h + l n,k,i α,h + a n,k,i α,h , (4.4a)
d n,k,i α,h , l n,k,i α,h , a n,k,i α,h ∈ RTN(T n h ). (4.4b)
Herein, a n,k,i α,h will be used to monitor the error in the solution of the linear algebraic system (4.3), l n,k,i α,h will be used to monitor the error in the linearization of (4.2) by (4.3), and d n,k,i α,h will be used to monitor the discretization error. Structurally, this can be achieved as follows: A. From the given numerical method, reconstruct locally d n,k,i α,h , α ∈ {n, w}. Typically, the degrees of freedom of d n,k,i α,h are directly prescribed using the available functions v α (s n,k,i w,h , p n,k,i w,h ); in any case, this construction should be independent of the linearization used to obtain (4.3) and of the iterative algebraic solver used to solve (4.3). B. From S n,k,i

w and P n,k,i w , compute the algebraic residual vectors R n,k,i w and R n,k,i n of (4.3): R n,k,i w R n,k,i n := - SS n,k-1 w SP n,k-1 w SS n,k-1 n SP n,k-1 n S n,k,i w P n,k,i w + D n,k-1 w D n,k-1 n . (4.5)
C. For simplicity, in this and in the next step, we suppose s n,k,i w,h and p n,k,i w,h piecewise constant on T n h . Generalizations follow as in [START_REF] Ern | Adaptive inexact Newton methods with a posteriori stopping criteria for nonlinear diffusion PDEs[END_REF]. From the given numerical method, define implicitly l n,k,i α,h , α ∈ {n, w}, such that

(q n α -∂ n t (φs n,k,i α,hτ ) -∇•(d n,k,i α,h + l n,k,i α,h ), 1) K = R n,k,i α | K ∀K ∈ T n h , α ∈ {n, w}. (4.6)
It is crucial to ensure that l n,k,i α,h go to zero when S n,k w , P n,k w , the solutions of (4.3), converge to S n w , P n w , the solution of (4.2). D. Construct locally a n,k,i α,h , α ∈ {n, w}, such that

(∇•a n,k,i α,h , 1) K = R n,k,i α | K ∀K ∈ T n h , α ∈ {n, w}, (4.7

An a posteriori error estimate distinguishing the space, time, linearization, and algebraic errors

We now further develop the framework of Section 3 in order to distinguish the space, time, linearization, and algebraic errors. Fix α ∈ {n, w} and consider the algorithm of Section 4.1 on the n-th time step, linearization step k, and iterative algebraic solver step i. Observe from (4.4a), (4.6), and (4.7) that u n,k,i α,h satisfies

(q n α -∂ n t (φs n,k,i α,hτ ) -∇•u n,k,i α,h , 1) K = 0 ∀K ∈ T n h , (4.8)
i.e., u n,k,i α,h satisfies assumption (3.1). Rewriting (3.4a)-(3.4d) with these notations, we pose, for K ∈ T n h ,

η n,k,i R,K,α := C P,K h K q n α -∂ n t (φs n,k,i α,hτ ) -∇•u n,k,i α,h K , (4.9a) η n,k,i F,K,α (t) := u n,k,i α,h -v α (s n,k,i w,hτ , p n,k,i w,hτ )(t) K , (4.9b)
η n,k,i NC,K,1 (t) := K(λ w (s n,k,i w,hτ ) + λ n (s n,k,i w,hτ )) ∇(p(s n,k,i w,hτ , p n,k,i w,hτ ) -pn,k,i hτ ) K (t), (4.9c) η n,k,i NC,K,2 (t) := K∇(q(s n,k,i w,hτ ) -qn,k,i hτ ) K (t). (4.9d)
We then have, as in Section 3.

η n,k,i R,K,α + η n,k,i F,K,α (t) 2 dt 1 2 + In K∈T n h η n,k,i NC,K,1 (t) 2 dt 1 2 + In K∈T n h η n,k,i NC,K,2 (t) 2 dt 1 2
.

(4.10)

We now distinguish the different error components. Define the spatial estimator

η n,k,i sp,K (t) := α∈{n,w} d n,k,i α,h -v α (s n,k,i w,h , p n,k,i w,h ) K + η n,k,i R,K,α 2 + η n,k,i NC,K,1 (t) 2 + η n,k,i NC,K,2 (t) 2 1 2 (4.11a)
and, for α ∈ {n, w}, the temporal estimators

η n,k,i tm,K,α (t) := v α (s n,k,i w,hτ , p n,k,i w,hτ )(t) -v α (s n,k,i w,h , p n,k,i w,h ) K , (4.11b)
the linearization estimators η n,k,i lin,K,α := l n,k,i α,h K , (4.11c)

and the algebraic estimators η n,k,i alg,K,α := a n,k,i α,h K . (4.11d)

Define also global versions of these estimators as Then

η n,k,i sp :=    3 In K∈T n h η n,k,i sp,K (t) 2 dt    1 2 , (4.12a) η n,k,i tm :=    α∈{n,w} In K∈T n h η n,k,i tm,K,α (t) 2 dt    1 2 , (4.12b) η n,k,i lin :=    α∈{n,w} τ n K∈T n h η n,k,i lin,K,α 2    1 2 , (4.12c) η n,k,i alg :=    α∈{n,w} τ n K∈T n h η n,k,i alg,K,α 2    1 2 . ( 4
|||(s w -s n,k,i w,hτ , p w -p n,k,i w,hτ )||| In ≤ η n,k,i sp + η n,k,i tm + η n,k,i lin + η n,k,i alg .

Stopping criteria and optimal balancing of the different error components

We now discuss how to balance the error components of Corollary 4.2.

In step 3(c)vi of the algorithm of Section 4.1, we evaluate η n,k,i sp , η n,k,i tm , η n,k,i lin , and η n,k,i alg and apply

η n,k,i alg ≤ γ alg (η n,k,i sp + η n,k,i tm + η n,k,i lin ) (4.13a)
as the stopping criterion for the iterative solution of the linear system (4.3). Here 0 < γ alg ≤ 1 is userspecified weight. Criterion (4.13a) expresses that there is no need to continue with the linear solver iterations if the overall error is dominated by the other components.

Similarly, in step 3d, we evaluate η n,k,i sp , η n,k,i tm , and η n,k,i lin and stop the iterative linearization of (4.2) whenever

η n,k,i lin ≤ γ lin (η n,k,i sp + η n,k,i tm ), (4.13b)
where 0 < γ lin ≤ 1 is user-specified weight. Criterion (4.13b) expresses that there is no need to continue with the linearization iterations if the overall error is dominated by the other components. Finally, in step 4, we evaluate η n,k,i sp , η n,k,i tm , and η n,k,i sp,K for all K ∈ T n h . The purpose is to achieve

η n,k,i sp ≈ η n,k,i tm , (4.13c)
In

η n,k,i sp,K (t) 2 dt 1 2
are comparable for all

K ∈ T n h , (4.13d)
and

η n ≤ ε n . (4.13e)
Here η n is given by (4.10) and ε n is user-specified criterion for the maximal error allowed on the time interval I n .

Remark 4.3 (Local stopping criteria). Following [START_REF] Jiránek | A posteriori error estimates including algebraic error and stopping criteria for iterative solvers[END_REF], [START_REF] El Alaoui | Guaranteed and robust a posteriori error estimates and balancing discretization and linearization errors for monotone nonlinear problems[END_REF], and [START_REF] Ern | Adaptive inexact Newton methods with a posteriori stopping criteria for nonlinear diffusion PDEs[END_REF], versions of (4.13a) and (4.13b) localized on each element of the mesh T n h should be used whenever one intends to refine adaptively the meshes T n h . Remark 4.4 (Evaluation cost). The evaluation of the different estimators of Corollaries 4.1 and 4.2 and of the stopping criteria (4.13a)-(4.13e) has linear cost in terms of the number of the elements of the meshes T n h . Moreover, it can be done completely in parallel. In practice, the estimators maybe only evaluated always after several iterations and various computational simplifications may be devised. Examples are given in Section 6 below.

Stopping criteria and adaptivity for implicit pressure-explicit saturations-type discretizations

We describe here our a posteriori error estimates and stopping criteria for implicit pressure-explicit saturations-type discretizations.

Iterative coupling for the pressure-saturation formulation

We first proceed as in Section 4 to obtain (4.1a)-(4.1b). We keep the wetting phase saturation equation (4.1a), and we replace the nonwetting phase saturation equation (4.1b) by the sum of (4.1a) and (4.1b).

We thus arrive at the following equivalent "pressure-saturation" formulation of (1.1a)-(1.1c):

-∇•((λ r,w (s w )+λ r,n (s w ))K(∇p w +ρ w g∇z) +λ r,n (s w )K(∇p c (s w )+ρ n g∇z -ρ w g∇z)) = q w + q n , (5.1a)

∂ t (φs w ) -∇•(λ r,w (s w )K(∇p w +ρ w g∇z)) = q w . (5.1b)
This formulation leads to the following solution algorithm:

1. Let the initial wetting saturation s 0 w,h ∈ H 1 (T 0 h ) (and pressure p 0 w,h ∈ H 1 (T 0 h )) be given. Set n = 1.

(a)

Choose some initial wetting saturation s n,0 w,h (S n,0 w is the corresponding algebraic vector). Typically, this is the approximate saturation from the last time step, s n-1 w,h . Similarly, for the pressure, typically define p n,0 w,h := p n-1 w,h . Set k := 1. (b) Set up the following linear elliptic problem, arising from (5.1a), with p w as the unknown:

-∇•((λ r,w (s n,k-1 w,h ) + λ r,n (s n,k-1 w,h )) ×K(∇p w + ρ w g∇z) + λ r,n (s n,k-1 w,h ) ×K(∇p c (s n,k-1 w,h ) + ρ n g∇z -ρ w g∇z)) = q w + q n .
(5.2)

After a spatial discretization, this problem corresponds to, in matrix form,

P n,k-1 wn P n,k w = D n,k-1 wn , (5.3)
where the matrix P n,k-1 wn and the right-hand side vector D n,k-1 wn depend on S n,k-1 w .

(c) i. Choose some initial pressure P n,k,0 w ; typically, P n,k,0

w = P n,k-1 w
. Set i := 1. ii. Perform a step of a chosen iterative algebraic solver for the solution of (5.3), starting from P n,k,i-1 w . At the present stage, we have approximations S n,k-1 w and P n,k,i w . iii. Build piecewise polynomial representations of the wetting saturations s n,k-1

w,h ∈ H 1 (T n h ) and pressures p n,k,i w,h ∈ H 1 (T n h ) from S n,k-1
w and P n,k,i w , according to the given numerical method.

Define the space-time function p n,k,i w,hτ ; it is affine in time on the time interval I n , given by p n-1

w,h at time t n-1 and by p n,k,i w,h at time t n . iv. Set up the following problem which arises from (5.1b): find s w such that -∇• λ r,w (s n,k-1 w,h )K(∇p n,k,i w,h + ρ w g∇z) = q w -∂ t (φs w ).

(5.4) Discretize (5.4) in space and in time. The temporal discretization is explicit. This gives the vector S n,k,i w . Build s n,k,i w,h ∈ H 1 (T n h ) and the space-time approximation s n,k,i w,hτ , given by s n-1

w,h at t n-1 and by s n,k,i w,h at t n . v. This step only concerns nonconforming methods. From s n,k,i w,h and p n,k,i w,h , prescribe the reconstructions pn,k,i h := I av (p(s n,k,i w,h , p n,k,i w,h )) and qn,k,i h := I av (q(s n,k,i w,h )). Define the global pressure reconstruction pn,k,i vi. From the given discretization, reconstruct the phase velocities u n,k,i α,h ∈ RTN(T n h ), α ∈ {n, w} (cf. Assumption 3.1). More precisely, we need to obtain the decompositions (4.4a)-(4.4b). Here, a n,k,i α,h , l n,k,i α,h , and d n,k,i α,h are used to monitor the algebraic, iterative coupling, and discretization error, respectively. Structurally, this can be achieved as follows: A. From the given numerical method, reconstruct locally the fluxes d n,k,i α,h , α ∈ {n, w}, as in the step 3(c)vA of the algorithm in Section 4.1, using the functions v α (s n,k,i w,h , p n,k,i w,h ). B. From S n,k-1 w and P n,k,i w , compute the algebraic residual vector R n,k,i wn :

R n,k,i wn := -P n,k-1

wn P n,k,i w + D n,k-1 wn . (5.5)
C. From the given method, define a vector field l n,k,i w,h + l n,k,i n,h ∈ RTN(T n h ) such that

(q n w + q n n -∇•(d n,k,i n,h + d n,k,i w,h + l n,k,i w,h + l n,k,i n,h ), 1) K = R n,k,i wn | K ∀K ∈ T n h .
(5.6)

In contrast to Section 4, where l n,k,i α,h is readily obtained for each individual phase α ∈ {n, w}, we have the sum of the total "iterative coupling error". In order to obtain l n,k,i w,h ∈ RTN(T n h ), we apply local mass balance from (5.4):

(q n w -∂ t (φs n,k,i w,hτ ) -∇•(d n,k,i w,h + l n,k,i w,h ), 1) K = 0 ∀K ∈ T n h .
(5.7)

Combined with (5.6), this also sets l n,k,i n,h . D. Construct a vector a n,k,i 3. Check whether the spatial and temporal errors are comparable (see (4.13c)), whether the spatial errors are equally distributed in the computational domain (see (4.13d)), and whether the total error is small enough (see (4.13e)); if this is the case, and t n < T , set n := n + 1 and go to step 2a. If not, refine the time step τ n and/or the space mesh T n h and go to step 2a.

n ∈ RTN(T n h ) such that (∇•a n,k,i n , 1) K = R n,k,i wn | K ∀K ∈ T n h , (5.8

Implicit pressure-explicit saturation formulation

Implicit pressure-explicit saturation discretization (IMPES) corresponds to the iterative coupling algorithm of Section 5.1 where only one step in k (k = 1) is done.

An a posteriori error estimate distinguishing the space, time, iterative coupling, and algebraic errors

We now use the framework of Section 3, or more precisely that developed in Section 4.2, in order to distinguish the space, time, iterative coupling, and algebraic errors. Fix α ∈ {n, w} and consider the algorithm of Section 5.1 on the n-th time step, iterative coupling step k, and algebraic solver step i. The approximate wetting saturation and pressure at our disposal at the present stage are thus (s n,k,i w,hτ , p n,k,i w,hτ ). Define

u n,k,i α,h := d n,k,i α,h + l n,k,i α,h + a n,k,i α,h , α ∈ {n, w}.
It follows from (5.6)-(5.8) that the reconstructions of the phases fluxes u n,k,i α,h satisfy (4.8). Replacing the terminology "k-th linearization step" by "k-th iterative coupling step", Corollary 4.1 holds true for this case. Analogously, from the definitions (4.11a)-(4.12d) (η n,k,i lin,K,α and η n,k,i lin represent the iterative coupling errors), Corollary 4.2 also holds true for this case.

Stopping criteria and optimal balancing of the different error components

Stopping criteria to be used in steps 2(c)vii, 2d, and 3 of the algorithm of Section 5.1 for optimal balancing of the different error components and overall error control are exactly the same as in Section 4.3.

Applications to finite volumes and numerical experiments

We present here the application of our theoretical results and numerical experiments for two finite volume discretizations of (1.1a)-(1.4b). We neglect the gravity terms (set z := 0), use the Brooks-Corey relations, see Section 1, and consider a test case from [START_REF] Klieber | Adaptive simulations of two-phase flow by discontinuous Galerkin methods[END_REF] with Ω = (0, 300)m × (0, 300)m, T = 4•10 6 s,

φ = 0.2, K = 10 -11 I m 2 , q w = q n = 0 s -1 , µ w = 5•10 -4 kg m -1 s -1 , µ n = 2•10 -3 kg m -1 s -1 , s rw = s rn = 0, p d = 5•10 3 kg m -1 s -2 , s 0 w = 0.2 on all K ∈ T 0 h , K ∈ K, s 0 w = 0.95 on K ∈ T 0 h , K ∈ K,
where I is the identity tensor and K the 18m × 18m block in the lower left corner. Let K be the 18m × 18m block in the upper right corner. We assume homogeneous Neumann boundary conditions (1.4a) everywhere except ∂ K ∩ ∂Ω and ∂ K ∂Ω, representing respectively the injection and production wells. Here, we impose the (inhomogeneous) Dirichlet conditions (1.4b) for the wetting phase saturation and pressure s w , p w . These are respectively equal to 0.95 and 3.45•10 6 kg m -1 s -2 in the lower left corner and to 0.2 and 2.41•10 6 kg m -1 s -2 in the upper right corner. This is a classical setting for the quarter five spot problem, where the flow is driven by the pressure gradient from the injection to the production well. Contrarily to [START_REF] Klieber | Adaptive simulations of two-phase flow by discontinuous Galerkin methods[END_REF], we use the wetting residual saturation s rw equal to 0 and not to 0.15, in order to treat a possibly degenerate problem. We consider two classical discretization approaches. First, we test fully implicit cell-centered finite volumes on uniform meshes consisting of rectangular parallelepipeds. Here, the approximations of the phase velocities u α in (1.2) are H(div, Ω)-conforming but the approximations of the global and complementary pressures p(s w , p w ) and q(s w ) in (2.4a), (2.4b), respectively, are not H 1 D (Ω)-conforming. The second scheme employs a vertex-centered finite volume method. Here the phase velocities are not H(div, Ω)-conforming but the global and complementary pressure approximations are H 1 D (Ω)-conforming. In this latter case, we employ iterative coupling and use adaptive meshing.

Fully implicit cell-centered finite volumes on regular rectangular parallelepipeds

We consider here one fixed mesh T h , consisting of rectangular parallelepipeds, and focus on adaptive stopping of the linear solver using (4.13a) and adaptive stopping of the nonlinear solver using (4.13b) in the context of Section 4.

Scheme definition

Let E int K stand for the faces e of an element K ∈ T h not lying on ∂Ω. We consider the following cell-centered finite volume discretization of (4.1a)-(4.1b): for all K ∈ T h , define s 0 w,K := s 0 w | K . Then, for all 1 ≤ n ≤ N , look for s n w,h , pn w,h , piecewise constant on T h , such that

φ s n w,K -s n-1 w,K τ n |K| + eKL∈E int K F w,eKL (s n w,h , pn w,h ) = 0, (6.1a) -φ s n w,K -s n-1 w,K τ n |K| + eKL∈E int K F n,eKL (s n w,h , pn w,h ) = 0, (6.1b)
for all elements K ∈ T h not contained in K and K; in K and K, the Dirichlet boundary values are imposed.

The face fluxes are defined as

F w,eKL (s n w,h , pn w,h ) := - λ r,w (s n w,K ) + λ r,w (s n w,L ) 2 |K| × pn w,L -pn w,K |x K -x L | |e KL |, (6.2a) F n,eKL (s n w,h , pn w,h ) := - λ r,n (s n w,K ) + λ r,n (s n w,L ) 2 |K| × pn w,L + p c (s n w,L ) -(p n w,K + p c (s n w,K )) |x K -x L | |e KL |. (6.2b)
Above, |K| stands for the measure of K ∈ T h , |e KL | for the measure of e KL , |x Kx L | for the distance of the barycenters x K and x L of K, L ∈ T h , and |K| for the absolute value of K. (6.1a)-(6.1b) represents the system of nonlinear algebraic equations (4.2).

Linearization and linear systems solution

Following Section 4, we obtain on the Newton-Raphson linearization step k and algebraic solver step i applied to (6.1a)-(6.1b)

φ s n,k,i w,K -s n-1 w,K τ n |K| + eKL∈E int K F k-1 w,eKL (s n,k,i w,h , pn,k,i w,h ) = -R n,k,i w,K , (6.3a)
-φ s n,k,i w,K -s n-1 w,K τ n |K| + eKL∈E int K F k-1 n,eKL (s n,k,i w,h , pn,k,i w,h ) = -R n,k,i n,K , (6.3b)
where the linearized face fluxes take the analytical form, for α ∈ {n, w},

F k-1 α,eKL (s n,k,i w,h , pn,k,i w,h ) (6.4) :=F α,eKL (s n,k-1 w,h , pn,k-1 w,h ) + M∈{K,L} ∂F α,eKL ∂s w,M (s n,k-1 w,h , pn,k-1 w,h )•(s n,k,i w,M -s n,k-1 w,M ) + M∈{K,L} ∂F α,eKL ∂ pw,M (s n,k-1 w,h , pn,k-1 w,h )•(p n,k,i w,M -pn,k-1 w,M ).
System (6.3a)-(6.3b) represents a cell-centered finite volume description of the system (4.5). We initialize p0 w,h = 2.41•10 6 kg m -1 s -2 everywhere except K, where we set p0 w,h := 3.45•10 6 kg m -1 s -2 .

Flux reconstructions and evaluation of the estimators

We now identify the error estimators of Section 4.2. We begin by identifying the fluxes d n,k,i α,h , l n,k,i α,h , and a n,k,i α,h of (4.4a). We construct them in the lowest-order Raviart-Thomas-Nédélec space RTN 0 (T h ), cf. [START_REF] Brezzi | Mixed and hybrid finite element methods[END_REF].

As in [START_REF] Di Pietro | An a posteriori error estimator for a finite volume discretization of the two-phase flow[END_REF][START_REF] Ern | Adaptive inexact Newton methods with a posteriori stopping criteria for nonlinear diffusion PDEs[END_REF], following [START_REF] Eymard | Finite volume approximation of elliptic problems and convergence of an approximate gradient[END_REF] in the linear case, the degrees of freedom of d n,k,i α,h and l n,k,i α,h are simply fixed by the finite volume fluxes of (6.2a)-(6.2b) and (6.4): for all K ∈ T h not contained in K, K and for all

e KL ∈ E int K , set ((d n,k,i α,h + l n,k,i α,h )•n K , 1) eKL := F k-1 α,eKL (s n,k,i w,h , pn,k,i w,h ) (6.5) and (d n,k,i α,h •n K , 1) eKL := F α,eKL (s n,k,i w,h , pn,k,i w,h ), (6.6)
α ∈ {n, w}. On the boundary faces of K, the value is set to zero in accordance with the homogeneous Neumann boundary condition (1.4a).

We now identify a n,k,i α,h . These are obtained, following [START_REF] Ern | Adaptive inexact Newton methods with a posteriori stopping criteria for nonlinear diffusion PDEs[END_REF]Section 4], by performing ν additional iterations of the algebraic solver,

a n,k,i α,h := d n,k,i+ν α,h + l n,k,i+ν α,h -(d n,k,i α,h + l n,k,i α,h ), (6.7)
where

d n,k,i+ν α,h +l n,k,i+ν α,h
are given by (6.5) with i replaced by i+ν. This definition does not give exactly (4.7) but is simple and sufficiently precise in practice for large enough ν (set to 15 in the present numerical experiments). Setting u n,k,i α,h := d n,k,i α,h + l n,k,i α,h + a n,k,i α,h , we obtain (3.1), up to the neglected misfit from the construction (6.7) of a n,k,i α,h . By this construction, the residual estimators of (4.9a) take very small values and are neglected henceforth.

In a cell-centered finite volume discretization, the piecewise constant pressure pn,k,i w,h does not have sufficient smoothness required by our framework. Thus, following [70, Section 3.2.2], we introduce a postprocessed wetting phase pressure p n,k,i w,h , quadratic on all K ∈ T h :

-λ r,w (s n,k,i w,K )K∇(p n,k,i w,h | K ) = d n,k,i w,h | K , (6.8a) p n,k,i w,h (x K ) = pn,k,i w,K . (6.8b)
This postprocessing links the pressure p n,k,i w,h to the reconstructed flux d n,k,i w,h by an equivalent of the Darcy law (1.2), whereas p n,k,i w,h takes the value of the original finite volume pressure pn,k,i w,K at the barycenters. One particular consequence of (6.8a) is that

v w (s n,k,i w,h , p n,k,i w,h )| K = d n,k,i w,h | K (6.9)
on all K ∈ T h , so that the first terms in (4.11a) disappear for the wetting phase. Repeating the same procedure for the nonwetting phase, we are lead to define the piecewise quadratic postprocessing p n,k,i

n,h -λ r,n (s n,k,i w,K )K∇(p n,k,i n,h | K ) = d n,k,i n,h | K , (6.10a)
p n,k,i n,h (x K ) = pn,k,i w,K + p c (s n,k,i w,K ). (6.10b)

In the nonwetting case, since we only have (6.10b) and not p n,k,i n,h | K = p c (s n,k,i w,K ) + p n,k,i w,h | K which would match fully with (1.1c), we do not have a complete equivalent of the phase velocities link (6.9). We, however, consider the difference of v n (s n,k,i w,h , p n,k,i w,h ) and d n,k,i n,h as negligible. Consequently, the nonwetting phase first terms in (4.11a) disappear as well. Moreover,

η n,k,i tm ≈    α∈{n,w} τ n 3 K∈T h d n,k,i α,h -d n-1 α,h 2 K    1 2 (6.11)
in a simplification of (4.12b), exploiting the piecewise affine behavior in time. Thus, reconstructions (6.8a)-(6.8b) and (6.10a)-(6.10b) of p n,k,i α,h need not be constructed in practice. We are now left with evaluating/approximating the nonconformity estimators η n,k,i NC,K,1 (t) and η n,k,i NC,K,2 (t) of (4.9c)-(4.9d). We want to avoid the difficult constructions p(s n,k,i w,hτ , p n,k,i w,hτ ) and q(s n,k,i w,hτ ) according to (2.4a)-(2.4b). We thus, following (2.5a)-(2.5b), define a piecewise quadratic pn,k,i

h by -λ w (s n,k,i w,K ) + λ n (s n,k,i w,K ) ×K∇(p n,k,i h | K ) = (d n,k,i w,h + d n,k,i n,h )| K , (6.12a) pn,k,i h (x K ) = p(s n,k,i w,K , pn,k,i w,K ), (6.12b)
and a piecewise quadratic qn,k,i

h by K∇(q n,k,i h | K ) = λ n (s n,k,i w,K )K∇(p n,k,i h | K ) + d n,k,i n,h | K , (6.13a) qn,k,i h (x K ) = q(s n,k,i w,K ). ( 6
:= I av (q n,k,i h ).
We then approximate the space-time integrals In η n,k,i NC,K,1 (t) 2 dt (6.14a)

≈ τ n 3 K(λ w (s n,k,i w,K ) + λ n (s n,k,i w,K ))∇(p n,k,i h -pn,k,i h ) 2 K + K(λ w (s n-1 w,K ) + λ n (s n-1 w,K ))∇(p n-1 h -pn-1 h ) 2 K + (K(λ w (s n,k,i w,K ) + λ n (s n,k,i w,K ))∇(p n,k,i h -pn,k,i h ), K(λ w (s n-1 w,K ) + λ n (s n-1 w,K ))∇(p n-1 h -pn-1 h )) K and In η n,k,i NC,K,2 (t) 2 dt (6.14b) ≈ τ n 3 K∇(q n,k,i h -qn,k,i h ) 2 K + K∇(q n-1 h -qn-1 h ) 2 K + (K∇(q n,k,i h -qn,k,i h ), K∇(q n-1 h -qn-1 h )) K , cf. [42, proof of Lemma 6.1].
Summarizing the above developments of estimators:

1. Identify the fluxes d n,k,i

. 5. Approximate η n,k,i sp of (4.12a) while omitting the first two terms in (4.11a) and using (6.14a)-(6.14b).

Construct the pressures pn,k

A further simplification of the estimators

The key for our a posteriori error estimates are the flux reconstructions d n,k,i α,h , l n,k,i α,h , and a n,k,i α,h . It is crucial to observe that in order to evaluate our estimators, we only need the [L 2 ] d norms involving these reconstructions. These can, however, be evaluated or approximated by quadrature formulas on each mesh element, so that it is not physically necessary to construct the reconstructions; knowledge of their values at quadrature points is sufficient. One can proceed similarly for the pressures pn,k,i h , qn,k,i h and their averagings pn,k,i h and qn,k,i h ; typically the quadrature points can be chosen equal to the localization of the Lagrange degrees of freedom, so that we can only compute the coefficients in the Lagrange bases. In Section 6.2.3 below, we give a concrete example of such a simplification in the context of vertex-centered finite volumes.

Computational performances

We evaluate here the computational performance of the above cell-centered finite volume approach. We consider a uniform 50 × 50 spatial mesh, fixed in time. The initial time step τ 1 was chosen as τ 1 = 10 5 s and was not changed by the criterion (4.13c). In order to improve the numerical treatment, we have employed a scaling of all the length units by 10 3 , i.e., we express all length units in millimeters instead of meters. The values of the wetting pressure than in particular range between 2410 and 3450, instead of 2.41•10 6 and 3.45•10 6 in the original units. We have approximated the differentiation in (6.4) numerically with a parameter ǫ = 10 -12 . The systems (6.3a)-(6.3b) arise here from an exemplary algebraic Figure 1: Wetting saturation (top), spatial a posteriori error estimates η n,k,i sp,K (t) (middle), and wetting pressure (bottom) at times 1.3•10 6 s, 2.6•10 6 s, and 4•10 6 s solver, the generalized minimum residual method (GMRes) [START_REF] Saad | GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF] without restarts, with the Jacobi (diagonal) preconditioning.

We compare two computational approaches. In the classical one, the GMRes iteration in (6.3a)-(6.3b) is continued until the relative algebraic residual becomes smaller than 10 -13 , and the Newton linearization is only stopped when

S n,k w -S n,k-1 w P n,k w -P n,k-1 w w,∞ ≤ 10 -11 . (6.15)
Here, • w,∞ stands for the weighted • ∞ norm, where the weights for the saturation unknowns are 1 and the weights for the pressure unknowns are 10 -3 . Together with the scaling of the length units, all the quantities in (6.15) are of order 1. In the adaptive approach proposed in this paper, we rely instead on the stopping criteria (4.13a) and (4.13b) with γ alg = γ lin = 0.001. In order to evaluate the algebraic error fluxes, we choose ν = 15 in (6.7), i.e., we perform 15 additional GMRes steps. If (4.13a) is not satisfied, then we directly set i := i + ν. Thus we only evaluate (4.13a) each 15 GMRes steps and the additional steps are not wasted.

In Figure 1, we illustrate the evolution of wetting saturation and pressure during the simulation. Tiny oscillations can be remarked, since no upwinding or other stabilization has been used. We also plot in Figure 1 the evolution of the spatial a posteriori error estimators η n,k,i sp,K (t) of (4.11a) approximated using (6.14a)-(6.14b) (divided by τ n ). We can see that they nicely follow the saturation front, while also detecting errors at the inflow and outflow. The results of Figure 1 come from the case where the adaptive stopping criteria (4.13a) and (4.13b) have been used, but practically undistinguishable results are obtained using the classical stopping criteria.

2: Estimators η n,k,i sp /τ n , η n,k,i tm /τ n , η n,k,i lin /τ n , η n,k,i alg /τ n ,
We next investigate the behavior of the estimators of the different error components. In order to present the results, we fix one time, 2.6•10 6 s. In the left part of Figure 2, we track the dependence of the different estimators on the GMRes iterations for the first Newton step on this time. The classical approach requires 1530 iterations until the relative algebraic residual becomes smaller than 10 -13 , whereas the adaptive criterion (4.13a) only requires 435 GMRes iterations. In the right part of Figure 2, we plot the different estimators as function of the Newton iteration. 11 iterations are necessary to reach (6.15), whereas only 3 iterations are sufficient to arrive at (4.13b).

Figures 3 and4 then assess the overall computational performances of the two approaches. In Figure 3 we plot the GMRes relative algebraic residual on each time and Newton step (left) and the number of Newton iterations on each time step (right). We can see that much higher (lower, respectively) numbers are sufficient with adaptive stopping criteria. We note that in particular much fewer Jacobian matrix assemblies are necessary in our approach. Figure 4 then gives the number of GMRes iterations for each linear system solved (left) and the cumulative number of GMRes iterations as function of time (right). From the last graph, we can conclude that in the adaptive approach the number of cumulative GMRes iterations is approximately 12-times smaller compared to that in the classical one. In addition, this ratio is growing Finally, Figure 5 gives an example of the reconstructed pressures which are at the heart of our spatial estimators: the global pressure pn,k,i h of (6.12a)-(6.12b) and the complementary one qn,k,i h of (6.13a)-(6.13b), as well as their averagings pn,k,i h = I av (p n,k,i h ) and qn,k,i h = I av (q n,k,i h ). The plots are given on an example of a rough 10 × 10 mesh at time 1.3•10 6 s.

Iteratively coupled implicit pressure-explicit saturation vertex-centered finite volumes on adaptively refined meshes with hanging nodes

Here, following Section 5.1, we consider adaptive stopping of the linear solver (4.13a), adaptive stopping of the iterative coupling (4.13b), adaptive choice of the time step (4.13c), and adaptive mesh refinement (4.13d).

Scheme definition

In addition to the simplicial meshes T n h , let D n h , 0 ≤ n ≤ N , be the dual meshes with dual volumes D formed around the vertices of T n h . In our setting, the dual meshes D n h consist of squares and the vertices of T n h are given by the square barycenters, see Figure 6, left and middle. We let V n h stand for the space of continuous, piecewise affine functions on T n h . In the wetting saturation space V n h,sw ⊂ V n h , the values in the vertices of T n h lying in K are fixed to 0.95, and in those lying in K to 0.2. Similarly, in the wetting pressure space V n h,pw ⊂ V n h , the vertex values in K are fixed to 3.45•10 6 kg m -1 s -2 and in K to 2.41•10 6 kg m -1 s -2 . Finally, D int,n h stands for all the elements of D n h not contained in K and K. The vertex-centered finite volume scheme reads as follows. Define s 0 w,h ∈ V 0 h,sw by the values of s 0 w at the vertices of T 0 h . Then, for all 1 ≤ n ≤ N and on a given iterative coupling step k ≥ 1, with given s n,k-1 w,h ∈ V n h,sw , we look for p n,k w,h ∈ V n h,pw discretizing the pressure equation (5.2) and realizing (5.3), On i-th step of an iterative algebraic solver applied to (6.16a), we then obtain the realization of (5.5) in the form

-λ r,w (s n,k-1 w,h ) + λ r,n (s n,k-1 w,h ) K∇p n,k w,h •n D +λ r,n (s n,k-1 w,h )K∇p c (s n,k-1 w,h )•n D , 1 ∂D\∂Ω = 0 ∀D ∈ D int,n h . ( 6
-λ r,w (s n,k-1 w,h ) + λ r,n (s n,k-1 w,h ) K∇p n,k,i w,h •n D +λ r,n (s n,k-1 w,h )K∇p c (s n,k-1 w,h )•n D , 1 ∂D\∂Ω = -R n,k,i wn,D ∀D ∈ D int,n h , (6.17a
) whereas (6.16b) becomes the realization of (5.4),

s n,k,i w,D := τ n φ|D| λ r,w (s n,k-1 w,h )K∇p n,k,i w,h •n D ,1 ∂D\∂Ω + s n-1 w,D , D ∈ D int,n h . ( 6
.17b)

Flux reconstructions and evaluation of the estimators

Let S n-1,n h , 1 ≤ n ≤ N , be a simplicial mesh, a conforming refinement of all D n-1 h , D n h , T n-1 h
, and T n h , see Figure 6, right. For D ∈ D n h , we denote by S D the restriction of S n-1,n h to D. All the flux reconstructions will be constructed on RTN 0 (S n-1,n h ). Following Section 5, proceeding in an equivalent but slightly different way, we first construct the fluxes d n,k,i wn,h , l n,k,i wn,h , and a n,k,i wn , the sums of the individual phase fluxes. Let 1 ≤ n ≤ N , k ≥ 1, and i ≥ 1. For ((d n,k,i wn,h + l n,k,i wn,h )

•n D , 1) e := -λ r,w (s n,k-1 w,h ) + λ r,n (s n,k-1 w,h ) K∇p n,k,i w,h •n D + λ r,n (s n,k-1 w,h )K∇p c (s n,k-1 w,h )•n D , 1 e
and, similarly but evaluating the saturation terms on the saturation s n,k,i w,h obtained from (6.17b), . For this purpose, we follow the equilibration of [START_REF] Vohralík | Guaranteed and fully robust a posteriori error estimates for conforming discretizations of diffusion problems with discontinuous coefficients[END_REF], see the construction of t 2 in Section 4.3.3 in this reference. This gives

(d n,k,i wn,h •n D , 1) e := -λ r,w (s n,k,i w,h ) + λ r,n (s n,k,i w,h ) K∇p n,k,i w,h •n D + λ r,n (s n,k,i w,h )K∇p c (s n,k,i w,h )•n D , 1
-(∇•(d n,k,i wn,h + l n,k,i wn,h ), 1) K = 1 |D| (R n,k,i wn,D , 1) K (6.18) and -(∇•d n,k,i wn,h , 1) K = 1 |D| (R n,k,i wn,D , 1) K for all subsimplices K ∈ S D , where we set R n,k,i wn,D := -(d n,k,i wn,h •n D , 1
) ∂D\∂Ω . This corresponds to the localization of the mass balance equation (6.17a) over the individual elements K ∈ S D and simultaneously ensures that l n,k,i wn,h converges to zero as the iterative coupling converges. In particular, (5.6) follows. Similar to (6.7), we finally define

a n,k,i wn,h := d n,k,i+ν wn,h + l n,k,i+ν wn,h -(d n,k,i wn,h + l n,k,i wn,h ).
We now split the total fluxes to the phase fluxes d n,k,i α,h , l n,k,i α,h , and a n,k,i α,h , α ∈ {n, w}. We employ the saturation/wetting phase balance equation (6.17b). As above, we first prescribe d n,k,i w,h + l n,k,i w,h . For D ∈ D int,n h and any face e of S n-1,n h included in ∂D but not in ∂Ω, we set

((d n,k,i w,h + l n,k,i w,h )•n D , 1) e := -λ r,w (s n,k-1 w,h )K∇p n,k,i w,h •n D , 1 e and (d n,k,i w,h •n D , 1) e := -λ r,w (s n,k,i w,h )K∇p n,k,i w,h •n D , 1 e .
We then once again following [71, construction of t 2 in Section 4.3.3] equilibrate the remaining degrees of freedom of (d n,k,i w,h + l n,k,i w,h )| D and d n,k,i w,h | D in such a way that

-φ s n,k,i w,D -s n-1 w,D τ n -∇•(d n,k,i w,h + l n,k,i w,h ), 1 K = 0 (6.19)
and that

-φ s n,k,i w,D -s n-1 w,D τ n -∇•d n,k,i w,h , 1 K = 1 |D| (R n,k,i w,D , 1) K for all K ∈ S D , where R n,k,i w,D := -φ s n,k,i w,D -s n-1 w,D τ n , 1 D -(d n,k,i w,h •n D , 1) ∂D\∂Ω .
As for a n,k,i w,h , we simply set it to 0, as there is no algebraic error associated with the computation of the wetting phase approximations. The nonwetting phase fluxes are then fixed by

u n,k,i n,h := u n,k,i wn,h -u n,k,i w,h (6.20)
and similarly for all d, l, a.

We are now ready to evaluate/approximate the estimators. We use directly (4.11d) and (4.11c). We next approximate v α (s n,k,i w,h , p n,k,i w,h ) by -λ r,α (s n,k,i w,h )K∇p n,k,i α,h , with p n,k,i n,h := p n,k,i w,h + p c (s n,k,i w,h ). Then, similar to (6.11), we approximate the time integral

η n,k,i tm ≈ α∈{n,w} τ n 3 K∈T n h λ r,α (s n,k,i w,h )K∇p n,k,i α,h -λ r,α (s n-1 w,h )K∇p n-1 α,h 2 K 1 2 . (6.21)
Finally, as there is no nonconformity, both η n,k,i NC,K,1 and η n,k,i NC,K,2 are zero, whereas it follows from the equilibrations (6.18) and (6.19) that η n,k,i R,K,α are close to zero. Thus we can approximate

η n,k,i sp ≈ 3 α∈{n,w} τ n K∈T n h d n,k,i α,h + λ r,α (s n,k,i w,h )K∇p n,k,i α,h 2 K 1 2 . (6.22)
It is to be noted that these estimators take the same form as those obtained in [START_REF] Cancès | An a posteriori error estimate for vertex-centered finite volume discretizations of immiscible incompressible two-phase flow[END_REF] (therein, no simplifications have been made).

A simplified evaluation of the estimators

The above evaluation of our estimators relies on the flux reconstructions d n,k,i α,h , l n,k,i α,h , a n,k,i α,h , α ∈ {n, w}. As in Section 6.1.4, but here still in a more direct way, the physical constructions of these fluxes can be completely avoided. All these flux reconstructions are in the RTN 0 (S n-1,n h ) space. Let us likewise approximate (λ r,α (s n,k,i w,h )K∇p n,k,i α,h )| K and (λ r,α (s n-1 w,h )K∇p n-1 α,h )| K in the space RTN 0 (K), for all K ∈ S n-1,n h . We do so by preserving the normal fluxes over all faces of every K (these are the degrees of freedom in RTN 0 (K)). Then all we need to evaluate are the L 2 (K) norms of vectors from RTN 0

(K), K ∈ S n-1,n h . Let v h ∈ RTN 0 (K), v h = e∈EK β e v e .
Here β e are the degrees of freedom, the face fluxes (v h •n K , 1) e , and v e are using that the three-point face-barycentric formula is exact for polynomials of order two. Thus, it is enough to obtain the degrees of freedom (face normal fluxes) for all the concerned fluxes, which is straightforward, and then to use (6.23) in order to evaluate all the estimators.

Computational performances, uniform meshes

We test here the computational performances of our iterative coupling approach. We start by a setting similar to Section 6.1.5. The spatial mesh is thus fixed, uniform 50 × 50, and the time steps likewise fixed and of size 10 4 s; bigger time steps seem to block the convergence of the iterative coupling. The unrestarted GMRes with Jacobi preconditioning is employed. The adaptive approach, with adaptive stopping criteria with γ alg = γ lin = 0.001 and this time ν = 5, is compared with the classical one, where the GMRes iteration in (6.17a) is continued until the relative algebraic residual becomes smaller than 10 -13 , and the iterative coupling is only stopped when (6.15) is satisfied. We use the simplified evaluation of the estimators of Section 6.2.3.

In the left part of Figure 7, we track the dependence of the different estimators on the GMRes iterations for the second iterative coupling step at time 2.6•10 6 s. 215 iterations are necessary in the classical case versus 20 in the adaptive one. In the right part of Figure 7, the different estimators are then plotted as function of the iterative coupling iteration. Only 3 iterations are needed in the adaptive case, versus 10 in the classical one.

Figures 8 and9 then assess the overall computational performances of the two approaches. Figure 8 plots the GMRes relative algebraic residual on each time and iterative coupling step (left) and the number of iterative coupling iterations on each time step (right). We can see that much higher (lower, respectively) numbers are sufficient in the adaptive case. Figure 9 then gives the number of GMRes iterations for each linear system solved (left) and the cumulative number of GMRes iterations as function of time (right). Overall, 18-times less total GMRes iterations are needed in the tested case.

Computational performances, adaptive meshes

We finally test the iterative coupling vertex-centered finite volume approach in the fully adaptive case, relying on all (4.13a)-(4.13d). In fact, the local version of the stopping criteria (4.13a)-(4.13b) of Remark 4.3, with γ alg = γ lin = 0.01, was used. As before, the unrestarted GMRes with Jacobi preconditioning was employed. The initial time step was τ 1 = 5 • 10 3 s; after 20 time steps, it has been brought to 10 4 s. The time error is not dominating; these time steps are chosen as they ensure the convergence of the iterative coupling. The initial spatial mesh was uniform 15 × 15 and was refined and derefined at each time step in order to achieve (4.13d). More precisely, it turns out that better results are achieved when the adaptivity is driven by the water phase components only, α = w in (6.22). At most two levels of refinement, where each square can be refined into 9 subsquares, are allowed for this test case. Examples of the resulting adaptive meshes are given in Figure 10. We can see that they follow very nicely the saturation front. More importantly, the final result matches in quality the one with the finest meshes employed uniformly over the whole space-time domain, for a several times smaller computational price. A detailed study of the adaptive mesh refinement case will be presented elsewhere.

Conclusions

In the present paper, we have developed a comprehensive framework for a posteriori error control and adaptivity for the immiscible incompressible two-phase flow in porous media. This framework covers various formulations, numerical methods, linearizations, and linear algebraic solvers. An important use of our results is in adequate stopping of the iterative algebraic solvers and iterative linearizations. Still relying on fixed spatial and temporal meshes, speed-ups by factors of tens can be achieved. Unlike the usual stopping criteria employing the L 2 -norm of the relative algebraic residual or the L ∞ -norm of the differences of the nonlinear iterates, our estimates and criteria are expressed in a unified way relying on the common language of the phase fluxes. With spatial and temporal meshes adaptivity, where error localization plays a crucial role, still more important computational savings can be achieved through a wise usage of the computer resources, while controlling the overall error. Applications to more complex settings, elaborate solution strategies, and further model problems such as compositional flows will be presented elsewhere.

ks w -s rw 1 -

 1 r,w (s w ) = s 4 e , k r,n (s w ) = (1s e ) 2 (1s 2 e ) and p c (s w ) = p d s s rws rn .

n , and SP n,k- 1 n 1 wn 1 w.

 111 are matrices formed from S n,k-are vectors likewise formed from S n,k-1 w and P n,k-1 w . (c) i. Choose some initial saturation S n,k,0 w and pressure P n,k,0 w . Typically, S n,k,0 w = S n,k-1 w and P n,k,0 w = P n,k-Set i := 1.

  Define the global pressure reconstruction pn,k,i hτ and the complementary pressure reconstruction qn,k,i hτ (cf. Assumption 3.1) affine in time on the time interval I n by pn-1 h and qn-1 h at time t n-1 (constructed from s 0 w,h , p 0 w,h for n = 1) and by pn,k,i h and qn,k,i h at time t n .

  .12d) Corollary 4.1 and the triangle and Cauchy-Schwarz inequalities yield: Corollary 4.2 (An a posteriori error estimate distinguishing the space, time, linearization, and algebraic errors). Let the assumptions of Corollary 4.1 be satisfied. Let the estimators be given by (4.12a)-(4.12d).

  hτ and the complementary pressure reconstruction qn,k,i hτ (cf. Assumption 3.1) affine in time on the time interval I n by pn-1 h and qn-1 h at time t n-1 and by pn,k,i h and qn,k,i h at time t n .

  -(6.13b) and set up pn,k,i h and qn,k,i h by the postprocessing averaging operator I av . 3. Evaluate η n,k,i lin from (4.11c), (4.12c) and η n,k,i alg from (4.11d) and (4.12d). 4. Approximate η n,k,i tm by (6.11).

Figure

  Figure 2: Estimators η n,k,i sp /τ n , η n,k,i tm /τ n , η n,k,i lin /τ n , η n,k,i alg /τ n ,and their sum as function of the GMRes iterations on the first Newton iteration (left) and as function of the Newton iterations (right) at time 2.6•10 6 s

Figure 3 :

 3 Figure 2: Estimators η n,k,i sp /τ n , η n,k,i tm /τ n , η n,k,i lin /τ n , η n,k,i alg /τ n ,and their sum as function of the GMRes iterations on the first Newton iteration (left) and as function of the Newton iterations (right) at time 2.6•10 6 s

Figure 4 :

 4 Figure 4: Number of GMRes iterations on each time and Newton step (left) and cumulative number of GMRes iterations as a function of time (right)

Figure 5 :

 5 Figure 5: Reconstructed global pressure pn,k,i h of (6.12a)-(6.12b) (top left), reconstructed complementary pressure qn,k,i h

Figure 6 :

 6 Figure 6: A dual square mesh D n h (left), the triangular mesh T n h (in red, middle), and the fine triangular mesh S n-1,n h

  e . The fluxes through ∂D ∩ ∂Ω are set to zero, in accordance with the Neumann boundary condition (1.4a). This fixes the face normal fluxes over the boundaries of all D ∈ D int,n h but not those in the interior of D ∈ D int,n h

Figure 7 :e∈EK β e (x f -a e ) 2 + e∈EK β e (y f -b e ) 2 ,

 722 Figure 7: Estimators η n,k,isp /τ n , η n,k,i tm /τ n , η n,k,i lin /τ n , η n,k,i alg /τ n , and their sum as function of the GMRes iterations on the second iterative coupling iteration (left) and as function of the iterative coupling iterations (right) at time 2.6•10 6 s

Figure 8 :

 8 Figure 8: GMRes relative algebraic residual on each time and iterative coupling step (left, part of the time interval) and number of iterative coupling iterations on each time step (right)

Figure 9 :

 9 Figure 9: Number of GMRes iterations on each time and iterative coupling step (left, part of the time interval) and cumulative number of GMRes iterations as a function of time (right)

Figure 10 :

 10 Figure 10: Wetting saturation on adaptively refined meshes at times 1.3•10 6 s, 2.6•10 6 s, and 4•10 6 s

  Let the pressures and velocities reconstructions phτ , qhτ , and u α,hτ , α ∈ {n, w}, satisfy Assumption 3.1. Let the estimators be given by (3.4a)-(3.4d). Then |||(s ws w,hτ , p wp w,hτ )|||

		1
	N	2
	estimate of the overall error). Let (s w , p w ) be the weak wetting saturation and pressure characterized by (2.6a)-(2.6d). Let (s w,hτ , p w,hτ ) ∈ [P 1 τ (H 1 (T ))] 2 be the approximate wetting α∈{n,w} n=1 In K∈T n h saturation and pressure. ≤ η n R,K,α + η n F,K,α (t) 2 dt

  If not, k := k + 1 and go back to step 3b.[START_REF] Amaziane | An existence result for a coupled system modeling a fully equivalent global pressure formulation for immiscible compressible two-phase flow in porous media[END_REF]. Check whether the spatial and temporal errors are comparable (see (4.13c) below), whether the spatial errors are equally distributed in the computational domain (see(4.

					)
	using, for instance, the algorithm of [49, Section 7.3] or the simplification of [43, Sec-
	tion 4]. It is crucial to ensure that a n,k,i α,h	go to zero when R n,k,i α	go to zero.
	vi. Check the convergence criterion for the linear solver (see (4.13a) below); if this criterion is
	reached, set S n,k w := S n,k,i w	and P n,k w	:= P n,k,i w	. If not, set i := i + 1 and go back to step 3(c)ii.
	(d) Check the convergence criterion for the nonlinear solver (see (4.13b) below); if this criterion is
	reached, set S n w := S n,k w , P n w := P n,k w , and, for nonconforming methods, pn h := pn,k,i

h

, qn h := qn,k,i h . 13d) below), and whether the total error is small enough (see (4.13e) below); if this is the case, and t n < T , set n := n + 1 and go to step 2. If not, refine the time step τ n and/or the space mesh T n h and go to step 2.

  , and u n,k,i α,h be as specified in Section 4.1. Let the estimators be given by (4.9a)-(4.9d). Then |||(s w -s n,k,i w,hτ , p wp n,k,i w,hτ )||| In

	≤ η n :=
	α∈{n,w} In K∈T n h

2, the local-in-time iterative-algorithms-running version of Theorem 3.3: Corollary 4.1 (Local-in-time estimate for linearization and algebraic iterates). Let (s w , p w ) be the weak wetting saturation and pressure characterized by (2.6a)-(2.6d). Consider the n-th time step, k-th linearization step, and i-th algebraic solver step of the algorithm of Section 4.1. Let s n,k,i w,hτ and p n,k,i w,hτ , pn,k,i hτ and qn,k,i hτ

  .13b)In practice, we approximate the integrals in (6.12b) and (6.13b) by a quadrature formula. As usual, we define pn,k,i at time t n . Following Section 4.1, we also define pn,k,i

	pn,k,i h	and qn,k,i	h	:= I av (p n,k,i h	at time t n-1 and by ) and qn,k,i h h

hτ and qn,k,i hτ affine in time on the time interval I n , given by pn-1 h and qn-1 h

The first author was supported by the ERT project "Enhanced oil recovery and geological sequestration of CO 2 : mesh adaptivity, a posteriori error control, and other advanced techniques" and by the GNR MoMaS. The second author was partially supported by the NSF-CDI under contract number DMS 0835745 and the Center for Frontiers of Subsurface Energy Security under Contract No. DE-SC0001114.