Robert Eymard
email: robert.eymard@univ-mlv.fr

Cindy Guichard
email: cindy.guichard@unice.fr

Raphaele Herbin
email: raphaele.herbin@latp.univ-mrs.fr

Roland Masson
email: roland.masson@unice.fr

Vertex centred Discretization of Two-Phase Darcy flows on General Meshes

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

The simulation of multiphase flow in porous media is a complex task, which has been the object of several works over a long period of time, see the reference books [START_REF] Aziz | Petroleum Reservoir Simulation[END_REF] and [START_REF]Fundamentals of Numerical Reservoir Simulations[END_REF]. Several types of numerical schemes have been proposed in the past decades. Those which are implemented in industrial codes are mainly built upon cell centred approximations and discrete fluxes, in a framework which is also that of the method we propose here. Let us briefly sketch this framework. The 3D simulation domain Ω is meshed by disjoint control volumes K ∈ M. Let us denote by Λ the diffusion matrix (which is a possibly full matrix depending on the point of the domain).

For each control volume K ∈ M, the set of neighbours N K ⊂ M is the set of all control volumes L involved in the mass balance in K, which means that the following approximation formula is used: -K div Λ∇ P dx the flux F K,L (P ), between control volumes K and L, is a linear function of the components of P which ensures the following conservativity property: F K,L (P ) = -F L,K (P ).

(

) 1
Such a linear function may be defined by

F K,L (P ) = J∈M K,L a J K,L p J , (2)
where the family (a J K,L ) J∈M K,L and M K,L ⊂ M are such that J∈M K,L a J K,L = 0, in order that F K,L (P ) vanishes on constant families.

The conservation of N c components within N α fluid phases typically couples the conservation equations of the components

Φ ∂ t A i + Nα α=1 div M α i F α = 0, i = 1, . . . , N c , (3)
together with the Darcy laws

F α = -Λ ∇ P α -ρ α g , α = 1, . . . , N α , (4)
and additional closure laws such as capillary pressure relations, thermodynamical equilibrium and pore volume conservation. In equations ( 3) and ( 4), A i represents the accumulation of component i per unit pore volume (assumed to take into account the dependence of the porosity with respect to the pressure), M α i is the amount of component i transported by the phase α, P α is the phase pressure, ρ α is the bulk density of the phase α, Φ is the porous volume, and g is the gravity acceleration.

The cell centred finite volume discretization of (3) and (4) reads

Φ K ∆t n (A n K,i -A n-1 K,i ) + Nα α=1 L∈N K
M n,α K,L,i F n,α K,L = 0, for all i = 1, . . . , N c , F n,α K,L = F K,L (P n,α ) -ρ n,α K,L g • (x L -x K ), for all α = 1, . . . , N α ,

where n is the time index, ∆t n = t n -t n-1 is the time step between time t n and time t n-1 , Φ K is the porous volume of the control volume K ∈ M, A K,i represents the accumulation of component i in the control volume K per unit pore volume, M α K,L,i is the amount of component i transported by the phase α from the control volume K to the control volume L (generally computed by taking the upstream value with respect to the sign of F α K,L ), P α is the family of the pressure unknowns of the phase α in all the control volumes, ρ α K,L is the bulk density of the phase α between the control volumes K and L, and x K is the centre of control volume K (which is not necessarily the barycenter of K).

When applying scheme [START_REF] Edwards | Finite volume discretization with imposed flux continuity for the general tensor pressure equation[END_REF], one should be very wary of the use of conformal finite elements in the case of highly heterogeneous media. Indeed, assuming that the control volumes are vertex centred with vertices located at the interfaces between different media, then the porous volume concerned by the flow of very permeable medium includes that of non permeable medium. This may lead to surprisingly wrong results on the component velocities. A possible interpretation of these poor results is that, when seen as a set of discrete balance laws, the finite element method provides the same amount of impermeable and permeable porous volume for the accumulation term for a node located at a heterogeneous interface.

We present in this paper the extension of a new scheme, called Vertex Approximate Gradient (VAG) scheme [START_REF] Eymard | Small-stencil 3d schemes for diffusive flows in porous media[END_REF], to multiphase flows in porous media. This extension provides discrete balances which may also be expressed by [START_REF] Edwards | Finite volume discretization with imposed flux continuity for the general tensor pressure equation[END_REF], thanks to a special choice of the control volumes and of the discrete fluxes, which can be written by [START_REF] Aavatsmark | A compact multipoint flux approximation method with improved robustness[END_REF]. As a result, the component velocities are correctly approximated. The purpose of respecting the form (2), [START_REF] Edwards | Finite volume discretization with imposed flux continuity for the general tensor pressure equation[END_REF] is to be able to plug it easily into an existing reservoir code, commonly using Multi-Point Flux Approximation (MPFA) [START_REF] Aavatsmark | Discretization on non-orthogonal, quadrilateral grids for inhomogeneous, anisotropic media[END_REF][START_REF] Aavatsmark | A compact multipoint flux approximation method with improved robustness[END_REF][START_REF] Edwards | Finite volume discretization with imposed flux continuity for the general tensor pressure equation[END_REF] data structure, by simply redefining the control volumes and the coefficients a J K,L of the discrete flux.

Although part of this scheme is vertex centred, we show that the solution, obtained in the case of a very heterogeneous medium, remains accurate on coarse meshes. This is a great advantage of this scheme, which is also always coercive, symmetric, and leads to a 27-stencil on hexahedral structured meshes. Moreover, the VAG scheme is very efficient on tetrahedral meshes, since the scheme can then be written with the nodal unknowns only, thus inducing a reduction of the number of degrees of freedom by a factor 5 compared with cell centred finite volume schemes such as MPFA schemes.

This paper is organized as follows. In Section 2, we present the extension of the one-phase VAG scheme to two-phase flows. We then provide the convergence analysis of the scheme in a simplified case in Section 3. Then numerical examples show the efficiency and the accuracy of the scheme in Section 4.

Vertex centred Discretization on generalised polyhedral meshes

Vertex centred discretization of Darcy fluxes (VAG scheme)

Let Ω be a bounded polyhedral subdomain of R 3 of boundary ∂Ω = Ω\Ω. The normal vector at the boundary outward the domain Ω is denoted by n. The disjoint subsets ∂Ω D and ∂Ω N of the boundary ∂Ω (each boundary face is assumed to be included in one of them) are such that ∂Ω D ∪ ∂Ω N = ∂Ω.

For a.e. (almost every) x ∈ Ω, Λ(x) denotes a 3-dimensional symmetric positive definite matrix such that that there exist

β 0 ≥ α 0 > 0 with α 0 ξ 2 ≤ ξ t Λ(x)ξ ≤ β 0 ξ 2 ,
for all ξ ∈ R 3 and for a.e. x ∈ Ω.

We consider the following diffusion equation

   div (-Λ∇ū) = f in Ω, ū = u d on ∂Ω D , -Λ∇ū • n = g on ∂Ω N .
Its variational formulation: find ū ∈ H 1 (Ω) such that ū = u d on ∂Ω D , and

Ω Λ∇ū • ∇v dx + ∂Ω N g vdσ = Ω f dx for all v ∈ H 1 D (Ω) = {w ∈ H 1 (Ω) | w = 0 on ∂Ω D }, admits a unique solution ū provided that the measure of ∂Ω D is nonzero, f ∈ L 2 (Ω), u d ∈ H 1/2 (∂Ω D ) and g ∈ L 2 (∂Ω N )
, which is assumed in the following.

Following [START_REF] Eymard | Small-stencil 3d schemes for diffusive flows in porous media[END_REF], we consider generalised polyhedral meshes of Ω. Let M be the set of cells that are disjoint open subsets of Ω such that κ∈M κ = Ω. For all κ ∈ M, x κ denotes the so called "centre" of the cell κ under the assumption that κ is star-shaped with respect to x κ . Let F denote the set of faces of the mesh which are not assumed to be planar, hence the term "generalised polyhedral cells". We denote by V the set of vertices of the mesh. Let V κ , F κ , V σ respectively denote the set of the vertices of κ ∈ M, faces of κ, and vertices of σ ∈ F. For any face σ ∈ F κ , we have V σ ⊂ V κ . Let M s denote the set of the cells sharing the vertex s. The set of edges of the mesh is denoted by E and E σ denotes the set of edges of the face σ ∈ F. It is assumed that for each face σ ∈ F, there exists a so called "centre" of the face x σ such that

x σ = s∈Vσ β σ,s s, with s∈Vσ β σ,s = 1,
where β σ,s ≥ 0 for all s ∈ V σ . The face σ is assumed to match with the union of the triangles T σ,e defined by the face centre x σ and each of its edge e ∈ E σ .

It is assumed that ∂Ω D = σ∈F D σ and that ∂Ω N = σ∈F N σ for a partition F = F D ∪ F N of F. Let V int = V \ ∂Ω denote the set of interior vertices, and V ext = V ∩ ∂Ω the set of boundary vertices. Let us then define the partition

V ext = V D ∪ V N of V ext with V D = σ∈F D V σ and V N = V ext \ V D .
The previous discretization is denoted by D and we define the discrete space

W D = {v κ ∈ R, v s ∈ R, κ ∈ M, s ∈ V},
and its subspace with homogeneous Dirichlet boundary conditions on

V D W D D = {v κ ∈ R, v s ∈ R, κ ∈ M, s ∈ V | v s = 0 for s ∈ V D }.

Vertex Approximate Gradient (VAG) scheme

The VAG scheme introduced in [START_REF] Eymard | Small-stencil 3d schemes for diffusive flows in porous media[END_REF] is based on a piecewise constant discrete gradient reconstruction for functions in the space W D . Several constructions are proposed based on different decompositions of the cell. Let us recall the simplest one based on a conforming finite element discretization on a tetrahedral sub-mesh, and we refer to [START_REF] Eymard | Benchmark 3d: the vag scheme[END_REF][START_REF] Eymard | Small-stencil 3d schemes for diffusive flows in porous media[END_REF] for two other constructions sharing the same basic features.

For all σ ∈ F, the operator

I σ : W D → R such that I σ (v) = s∈Vσ β σ,s v s ,
is by definition of x σ a second order interpolation operator at point x σ .

Let us introduce the tetrahedral sub-mesh T = {T κ,σ,e , e ∈ E σ , σ ∈ F κ , κ ∈ M} of the mesh M, where T κ,σ,e is the tetrahedron defined by the cell centre x κ and the triangle T σ,e as shown by Figure 1. For a given v ∈ W D , we define the function v T ∈ H 1 (Ω) as the continuous piecewise affine function on each tetrahedron T of T such that v T (x κ ) = v κ , v T (s) = v s , and v T (x σ ) = I σ (v) for all κ ∈ M, s ∈ V, σ ∈ F. The nodal basis of this finite element discretization will be denoted by η κ , η s , κ ∈ M, s ∈ V.

Following [START_REF] Eymard | Small-stencil 3d schemes for diffusive flows in porous media[END_REF], the Vertex Approximate Gradient (VAG) scheme is defined by the discrete variational formulation: find u ∈ W D such that u s = u d s for all s ∈ V D and

a D (u, v) + ∂Ω g(x) v T (x) dσ = Ω f (x) v T (x) dx for all v ∈ W D D ,
with a D the bilinear form such that

a D (u, v) = Ω ∇u T (x) • Λ(x) ∇v T (x) dx for all (u, v) ∈ W D × W D , and
u d s = 1 ∂Ω D η s (x) dσ ∂Ω D u d (x)η s (x) dσ for all s ∈ V D .

Conservative fluxes

Let us define for all κ ∈ M and s, s

∈ V κ a s κ,s = κ ∇η s (x) • Λ(x)∇η s (x) dx. One has a D (u, v) = κ∈M s∈Vκ s ∈Vκ a s κ,s (u s -u κ )(v s -v κ ),
leading to the definition of the following conservative fluxes between a given cell κ ∈ M and its vertices s

∈ V κ F κ,s (u) = s ∈Vκ a s κ,s (u κ -u s ),
and

F s,κ (u) = -F κ,s (u).
The VAG scheme is equivalent to the following discrete system of conservation laws:

                     s∈Vκ F κ,s (u) = κ f (x) η κ (x) dx for all κ ∈ M, κ∈Ms F s,κ (u) + g s = Ω f (x) η s (x) dx for all s ∈ V int ∪ V N , u s = u d s for all s ∈ V D ,
where

g s = ∂Ω g(x)η s (x) dσ.

Extension of the VAG scheme to multiphase flow in porous media

The VAG scheme has shown its good approximation properties for single phase Darcy flows on general meshes as exhibited in the benchmark results [START_REF] Eymard | Benchmark 3d: the vag scheme[END_REF][START_REF] Eymard | Benchmark 3d on discretization schemes for anisotropic diffusion problem on general grids[END_REF]. This suggests the idea to use the fluxes F κ,s (u) = -F s,κ (u) between a cell κ of the mesh and its vertices s ∈ V κ to discretize multiphase Darcy flow models on general meshes.

To achieve this, we first define the set of control volumes K ∈ M as the union of the cells and of the vertices

M = M ∪ V,
and we associate a volume m K to each control volume

K ∈ M such that K∈M m K = m Ω and m K > 0 for all K ∈ M, (6)
with

m Ω = Ω dx.
It is achieved by a conservative redistribution to the vertices of the surrounding cell volumes

       m s = ω κ∈Ms α s κ m κ for all s ∈ V, m κ = 1 -ω s∈Vκ α s κ m κ for all κ ∈ M, (7)
with m κ = κ dx, and α s κ ≥ 0, κ∈Ms α s κ = 1, which guarantees [START_REF] Eymard | convergence d'un schéma de type éléments finis-volume finis pour un système formé d'une équation elliptique et d'une équation hyperbolique[END_REF] provided that the parameter ω > 0 is chosen small enough.

In practice, the weights α s κ are chosen in such a way that the volumes m s at the vertices are mainly taken from the surrounding cells with the highest permeabilities, using the formula:

α s κ = a κ,s κ ∈Ms a κ ,s , for all s ∈ V and κ ∈ M s , (8)
with a κ,s = s ∈Vκ a s κ,s > 0. This choice of the weights is the key ingredient to obtain an accurate approximation of the saturations and compositions, on the coarse meshes used in practical situations involving highly heterogeneous media.

Then, the fluxes F κ,s (u) = -F s,κ (u) are used to connect the control volumes κ and s for the discretization of multiphase Darcy flows in the framework defined by ( 1), ( 2), [START_REF] Edwards | Finite volume discretization with imposed flux continuity for the general tensor pressure equation[END_REF] with

N K = V κ for K = κ, N K = M s for K = s and M κ,s = M s,κ = V κ .
Although the fluxes F κ,s (u) and F s,κ (u) are not as usual defined as the approximation of the continuous fluxes σ -Λ∇ P • n σ on a given face σ of the mesh, the mathematical analysis given in the following section shows that they lead to a convergent scheme, at least in a particular two-phase flow case.

Convergence analysis

We develop the convergence analysis in the case of a simplified two-phase Darcy flow model coupling an elliptic equation for the pressure P with an advection equation for the saturation denoted by ū. This model can be obtained from a general two-phase Darcy flow model, assuming that both phases are incompressible, that the sum of the mobilities is constant, and neglecting gravity and capillary effects:

                     div -Λ∇ P = 0 on Ω × (0, t f ), -Λ∇ P • n = g on ∂Ω × (0, t f ), Ω P (x)dx = 0, ∂ t ū + div -ū Λ∇ P = 0 on Ω × (0, t f ), ū = u -on ∂Ω -× (0, t f ), ū| t=0 = u 0 on Ω, (9)
where

∂Ω -= {x ∈ ∂Ω | g(x) < 0}, g ∈ L 2 (∂Ω) is such that ∂Ω g(x)dσ = 0, u 0 ∈ L ∞ (Ω), and u -∈ L ∞ ∂Ω -× (0, t f ) . Let ∂Ω + = {x ∈ ∂Ω | g(x)
≥ 0}, and let us define the function space

D = {ϕ ∈ C ∞ c (R 3 × R, R) | ϕ(x, t) = 0 on R 3 × (t f , +∞) and on ∂Ω + × R}.
The system (9) admits the following weak formulation: find ū ∈ L ∞ Ω × (0, t f ) , P ∈ H 1 (Ω) such that Ω P (x)dx = 0, and for all Q ∈ H 1 (Ω) and ϕ ∈ D one has

               Ω ∇ P (x) • ∇Q(x) dx + ∂Ω g(x) Q(x) dσ = 0, - t f 0 Ω ū(x, t) ∂ t ϕ(x, t) dx dt - Ω u 0 (x) ϕ(x, 0) dx + t f 0 Ω ū(x, t) Λ(x)∇ P (x) • ∇ϕ(x, t) dx dt + t f 0 ∂Ω - g(x, t) u -(x, t)ϕ(x, t) dσ dt = 0. ( 10
)
The existence of a weak solution to [START_REF] Eymard | Benchmark 3d on discretization schemes for anisotropic diffusion problem on general grids[END_REF] is shown in [START_REF] Eymard | convergence d'un schéma de type éléments finis-volume finis pour un système formé d'une équation elliptique et d'une équation hyperbolique[END_REF] through the convergence of a finite volume scheme. The uniqueness is shown in [START_REF] Caginalp | Nonlinear equations with coefficients of bounded variations in two space variables[END_REF] assuming enough regularity on the data and on the domain Ω such that the pressure solution P is in C 2 (Ω).

The time discretization of the interval (0, t f ) is denoted by S = t 0 , t 1 , • • • , t N with t 0 = 0 and t N = t f and we set ∆t n = t n -t n-1 for all n = 1, • • • , N . All the unknowns at time t n will be denoted in the following with the superscript n.

Applying to (9) the space time discretization issued from Section 2, we look for

P ∈ W D and u n ∈ W D , n = 1, • • • , N such that                                  s∈Vκ F κ,s (P ) = 0 for all κ ∈ M, κ∈Ms F s,κ (P ) + g + s + g - s = 0 for all s ∈ V, Ω P T (x) dx = 0 , m κ u n κ -u n-1 κ ∆t n + s∈Vκ F κ,s (P ) u n κ,s = 0 for all κ ∈ M, n = 1, • • • , N, m s u n s -u n-1 s ∆t n + κ∈Ms F s,κ (P ) u n κ,s + g + s u n s + g - s u -,n s = 0 for all s ∈ V, n = 1, • • • , N, (11)
where

g + s = ∂Ω g + (x)η s (x) dσ, g - s = ∂Ω g -(x)η s (x) dσ, setting y + = 1 2 (y + |y|), y -= 1 2 (y -|y|) for y ∈ R, u -,n s =    1 ∆t n ∂Ω η s (x)dx t n t n-1 ∂Ω u -(x, t)η s (x)dσdt if s ∈ V ext , 0 if s ∈ V int , u 0 κ = 1 Ω η κ (x)dx Ω u 0 (x)η κ (x)dx, u 0 s = 1 Ω η s (x)dx Ω u 0 (x)η s (x)dx
and using the upwind approximation

u k,s = u κ if F κ,s (P ) ≥ 0, u s if F κ,s (P ) < 0.
It is easy to show (see e.g. [START_REF] Eymard | Finite volume methods[END_REF]) that there exists a unique solution to the system [START_REF] Hubert | Basic error estimates for elliptic problems: Handbook of Numerical Analysis II[END_REF] and that the solution u n satisfies a discrete maximum principle. The discrete solution ũD,S is reconstructed from u n , n = 1, • • • , N on the whole space time domain using again the weights defined by ( 7)-( 8) in order to ensure the global conservation and to be accurate for highly heterogeneous media on coarse meshes:

ũD,S (x, t) = (1 - s∈Vκ ωα s κ )u n κ + s∈Vκ ωα s κ u n s on κ × (t n-1 , t n ) for all κ ∈ M, n = 1, • • • , N.
From the discrete maximum principle, ũD,S satisfies the estimate

ũD,S L ∞ (Ω×(0,t f )) ≤ max u 0 L ∞ (Ω) , u - L ∞ (∂Ω -×(0,t f )) . ( 12
)
For a given discretization D, let us define for each tetrahedron T of the sub-mesh T , its diameter h T and the diameter ρ T of the insphere of T . Then, we define

h D = max T ∈T h T , θ D = max T ∈T h T ρ T , γ D = max T ∈T h D h T , M D = max κ∈M #V κ ,
where #A denotes the cardinal of a given set A. The maximum time step of the time interval discretization S of (0, t f ) is denoted by ∆t S .

Let us consider a family of space and time discretizations D l and S l , l ∈ N of respectively Ω and (0, t f ). It is assumed that the family of tetrahedral sub-meshes T l is quasi-uniform and shape regular in the sense that there exist reals θ, γ, M , such that θ D l ≤ θ, γ D l ≤ γ, and M D l ≤ M for all l ∈ N. Then we can prove the following theorem.

Theorem 3.1. Assuming that the discretization parameters h D l , ∆t S l , and ω l in [START_REF] Eymard | Finite volume methods[END_REF] tend to zero when l → +∞, then the pressure solution P T l converges to P in H 1 (Ω), and there exists a subsequence of ũD l ,S l , l ∈ N which converges for the L ∞ (Ω × (0, t f )) weak-topology to a weak solution ū of [START_REF] Eymard | Benchmark 3d on discretization schemes for anisotropic diffusion problem on general grids[END_REF]. If the weak solution of (10) is unique, then the full sequence converges.

The convergence of the pressure P T l is a classical result of finite element approximation using the assumption of shape regularity θ D l ≤ θ for all l ∈ N (see e.g. [START_REF] Hubert | Basic error estimates for elliptic problems: Handbook of Numerical Analysis II[END_REF]).

The weak convergence proof of ũD l ,S l up to a subsequence follows similar arguments than in [START_REF] Eymard | convergence d'un schéma de type éléments finis-volume finis pour un système formé d'une équation elliptique et d'une équation hyperbolique[END_REF] and [START_REF] Eymard | Finite volume methods[END_REF] with some adaptations. We shall first prove the following weak BV estimate. Lemma 3.1. Let D and S denote space and time discretizations of Ω × (0, t f ), then the discrete solution u n ∈ W D , n = 1, . . . , N satisfies the following weak BV estimate

N n=1 ∆t n κ∈M s∈Vκ |F κ,s (P )||u n κ -u n s | + s∈V -g - s |u -,n s -u n s | ≤ C 1 √ h D ,
with a constant C depending only on the discretization parameters θ D , γ D , M D , and on the data u -, u 0 , g, Λ, Ω, and t f .

Proof of Lemma 3.1: Let us denote the equations in ( 11) by (11.1) to (11.5) and consider the following expression

2 κ∈M (11.4) u n κ - 1 2 (11.1) (u n κ ) 2 + 2 s∈V (11.5) u n s - 1 2 (11.2) (u n s ) 2 .
Multiplying all terms by ∆t n and summing over n = 1, • • • , N , we obtain after some computations that

κ∈M m κ (u N κ ) 2 -(u 0 κ ) 2 + s∈V m s (u N s ) 2 -(u 0 s ) 2 + N n=1 κ∈M m κ (u n κ -u n-1 κ ) 2 + s∈V m s (u n s -u n-1 s ) 2 + N n=1 ∆t n κ∈M s∈Vκ |F κ,s (P )|(u n κ -u n s ) 2 + s∈V -g - s (u -,n s -u n s ) 2 + N n=1 ∆t n s∈V g + s (u n s ) 2 + g - s (u -,n s ) 2 = 0.
Hence one has

N n=1 ∆t n κ∈M s∈Vκ |F κ,s (P )|(u n κ -u n s ) 2 + s∈V -g - s (u -,n s -u n s ) 2 ≤ m Ω u 0 2 L ∞ (Ω) + t f u -2 L ∞ (∂Ω -×(0,t f )) ∂Ω -g -(x)dσ,
and by the Cauchy Schwarz inequality and λ M (A κ ) its maximum eigenvalue, then one has

N n=1 ∆t n κ∈M s∈Vκ |F κ,s (P )||u n κ -u n s | + s∈V -g - s |u -,n s -u n s | ≤ m Ω u 0 2 L ∞ (Ω) + t f u -2 L ∞ (∂Ω -×(0,t f )) ∂Ω -g -(x)dσ
E ≤ κ∈M (#V κ ) 1 2 s∈Vκ |F κ,s (P )| 2 1 2 ≤ M 1 2 D κ∈M λ M (A κ ) 1 2
s,s ∈Vκ a s κ,s (P κ -P s )(P κ -P s )

1 2 ≤ M 1 2 D κ∈M λ M (A κ ) 1 2
κ∈M s,s ∈Vκ a s κ,s (P κ -P s )(P κ -P s )

1 2 ≤ M 1 2 D κ∈M λ M (A κ ) 1 2 Ω Λ∇P T • ∇P T dx 1 2 .
We conclude the proof from the estimates

κ∈M λ M (A κ ) ≤ C m Ω β 0 θ 2 D γ 2 D h 2 D ,
with C depending only on M D , and

Ω Λ∇P T • ∇P T dx ≤ Ω Λ∇ P • ∇ P dx.
Proof of Theorem 3.1: From Estimate (12), there exists a subsequence of ũD l ,S l , l ∈ N (still denoted by l ∈ N for simplicity) which converges to a function w ∈ L ∞ (Ω × (0, t f )) for the weak-topology. The aim of the proof is now to show that w is a weak solution of [START_REF] Eymard | Benchmark 3d on discretization schemes for anisotropic diffusion problem on general grids[END_REF].

Let ϕ ∈ D and let us set ϕ n κ = ϕ(x κ , t n-1 ) and ϕ n s = ϕ(s, t n-1 ). We consider the following expression

E 1 l + E 2 l = 0, with E 1 l = N n=1 κ∈M m κ (u n κ -u n-1 κ ) ϕ n κ + s∈V m s (u n s -u n-1 s ) ϕ n s , and
E 2 l = - N n=1 ∆t n κ∈M s∈Vκ (u n κ -u n s )(F κ,s (P ) -ϕ n κ + F κ,s (P ) + ϕ n s ) + s∈V -g - s (u -,n s -u n s ) ϕ n s .
Let us then prove that E 1 l converges to T 1 and E 2 l to T 2 with

T 1 = - t f 0 Ω w(x, t) ∂ t ϕ(x, t) dx dt - Ω u 0 (x) ϕ(x, 0) dx, and
T 2 = t f 0 Ω w(x, t) Λ(x)∇ P (x) • ∇ϕ(x, t) dx dt + t f 0 ∂Ω - g(x, t) u -(x, t)ϕ(x, t) dσ dt. Let ũn κ = (1 -ω l s∈Vκ α s κ )u n κ + ω l s∈Vκ α s κ u n s . One has E 1 l = F 1 l + ω l F 2 l with F 1 l = N n=1 κ∈M m κ (ũ n κ -ũn-1 κ ) ϕ n κ = - κ∈M m κ N n=1 ũn κ t n t n-1 ∂ t ϕ(x κ , t)dt + ũ0 κ ϕ(x κ , 0) , and
F 2 l = N n=1 κ∈M m κ s∈Vκ α s κ (u n s -u n-1 s )(ϕ n s -ϕ n κ ),
which is such that

F 2 l = - κ∈M m κ s∈Vκ α s κ N n=1 u n s t n t n-1 (∂ t ϕ(s, t) -∂ t ϕ(x κ , t))dt + u 0 s (ϕ(s, 0) -ϕ(x κ , 0)) . Hence |E 1 l -F 1 l | → 0 with F 1 l = N n=1 κ∈M m κ (ũ n κ -ũn-1 κ ) ϕ n κ = - t f 0 Ω ũD l ,S l (x, t)∂ t ϕ(x, t)dxdt - Ω ũ0 D l (x) ϕ(x, 0)dx,
where ϕ n κ = 1 mκ κ ϕ(x, t n-1 )dx and ũ0 D l is the piecewise constant function equal to ũ0 κ on each cell κ ∈ M. Since F 1 l → T 1 as l → ∞, this concludes the proof that E 1 l converges to T 1 .

Let us set

E 3 l = - N n=1 ∆t n κ∈M s∈Vκ F κ,s (P ) (u n κ -u n s ) ϕ n s + s∈V -g - s (u -,n s -u n s ) ϕ n s .
One has

|E 3 l -E 2 l | = - N n=1 ∆t n κ∈M s∈Vκ F κ,s (P ) -(u n κ -u n s ) (ϕ n s -ϕ n κ ),
which tends to zero thanks to Lemma 3.1.

Let us add to E 3 l the expression

N n=1 ∆t n κ∈M (11.1) u n κ ϕ n κ + s∈V (11.2) u n s ϕ n s .
We thus obtain that

E 3 l = N n=1 ∆t n κ∈M s∈Vκ F κ,s (P ) (ϕ n κ -ϕ n s ) u n κ + s∈V (g - s u -,n s + g + s u n s ) ϕ n s = N n=1 ∆t n Ω u n D l (x)Λ(x)∇P T l (x) • ∇ϕ n T l (x)dx + s∈V (g - s u -,n s + g + s u n s ) ϕ n s ,
where u n D l is the piecewise constant function equal to u n κ on each cell κ ∈ M, and ϕ n T l = κ∈M ϕ n κ η κ + s∈V ϕ n s η s . We conclude that E 3 l converges to T 2 since P T l converges to P in H 1 (Ω), ϕ vanishes on the support of g + and since the sequence u D l ,S l , l ∈ N defined by

u D l ,S l (x, t) = u n D l (x) for a.e. (x, t) ∈ Ω × (t n-1 , t n ) for all n = 1, • • • , N l ,
also converges to w for the weak-topology of L ∞ (Ω × (0, t f )), thanks to the assumption that ω l > 0 converges to zero.

Numerical examples

The numerical solutions computed by the VAG scheme applied to multiphase flow are compared with the solutions resulting from the cell centred MPFA O scheme [START_REF] Aavatsmark | Discretization on non-orthogonal, quadrilateral grids for inhomogeneous, anisotropic media[END_REF][START_REF] Edwards | Finite volume discretization with imposed flux continuity for the general tensor pressure equation[END_REF]. The first three test cases are designed to better understand the properties of the VAG scheme for two-phase flows regarding the sensitivity of the solution to the parameter ω used in the redistribution of the porosity, and the effect of large heterogeneities on the transport. The last test case is a more realistic example which simulates the near-well injection of miscible CO 2 in an aquifer.

Decoupled two phase flow

We consider the simplified two phase flow model of subsection 3 on the domain Ω = (0, 1) 3 with Λ = I, u 0 = 0, and modified boundary conditions. Let (x, y, z) denote the Cartesian coordinates of x. We specify a pressure P 1 at the left side x = 0 and a pressure P 2 at right side x = 1 such that P 1 > P 2 . Homogeneous Neumann boundary conditions g = 0 are imposed at the remaining boundaries. The input gas saturation is set to u -= 1 at the input boundary x = 0. The system admits an analytical solution given by P (x) = (P 2 -P 1 )x + P 1 , and ū(x, t) = 1 if x (P 1 -P 2 ) t, 0 else. We consider two different grids for this test. The first one is a uniform Cartesian grid of size 32×32×32. The second grid is composed of 15266 tetrahedra. Both meshes are extracted from the FVCA6 3D Benchmark [START_REF] Eymard | Benchmark 3d on discretization schemes for anisotropic diffusion problem on general grids[END_REF].

Figure 2 shows, for each grid, the gas saturation u in function of the x-coordinate at the time of the simulation for which the gas has filled half of the reservoir. We have plotted the analytical solution ū and the discrete solutions (x κ , u κ ) for all cells κ ∈ M obtained with the VAG scheme and the MPFA O scheme. For the VAG scheme, we use the post-processed values

ũκ = (1 -ω s∈Vκ α s κ )u κ + ω s∈Vκ α s κ u s ,
deduced from the redistribution of the volumes defined by ( 7)-( 8) . The results presented in Figure 2 clearly show that, for each grid, the discrete solutions of both schemes intersect the analytical solution at the point ( 12 , 1 2 ), which exhibits that the velocity of the flow is well approximated.

On Figure 2(a), the solutions of the VAG scheme on the Cartesian grid is plotted for both ω = 0.01 and ω = 0.3 [START_REF] Eymard | Finite volume methods[END_REF]. The value ω = 0.3 roughly corresponds to match the volume of each vertex m s with the volume of each surrounding cell m κ , κ ∈ M s . As expected, the choice ω = 0.3 leads to a slightly less diffusive scheme than the VAG scheme with ω = 0.01. We note also that the VAG scheme is slightly less diffusive on such meshes than the MPFA O scheme (which degenerates for Cartesian grids to the Two Points Flux Approximation scheme denoted by TPFA). On the other hand, for the tetrahedral mesh, we can notice on Figure 2(b) that the VAG scheme is slightly more diffusive than the MPFA O scheme. This has been observed for both values of ω = 0.01 and 0.3.

Note also that for both type of meshes the convergence of the VAG scheme has been obtained numerically for a fixed value of ω and for both ω = 0.01 and ω = 0.3. This may suggest that the assumption that ω tends to zero, made in the convergence analysis, is not necessary in practical cases.

In terms of CPU time a ratio of 15 is observed between the simulation time obtained with the MPFA O scheme on the tetrahedral mesh and the VAG scheme on the same mesh. This huge factor is due to the reduced size of the linear system obtained with the VAG scheme after elimination of the cell unknowns compared with the MPFA O scheme, both in terms of number of unknowns (around five time less for the VAG scheme than for the MPFA O scheme) and in terms of number of non zero elements per line.

Two phase flow for a strongly heterogeneous test case on a coarse mesh

The aim of the following test case is to show that, thanks to the redistribution of the porous volume at the vertices defined by ( 8), [START_REF] Eymard | Finite volume methods[END_REF], the VAG scheme provides solutions which are just as accurate as the solutions given by cell centred schemes in the case of large jumps of the permeability tensor on coarse meshes.

Let us consider a stratified reservoir Ω = (0, 100) × (0, 50) × (0, 100) m 3 with five horizontal layers l = 1, • • • , 5 of thickness 20 m, and numbered by their increasing vertical position. The even layers are drains of constant high isotropic permeability K d and odd layers are barriers of constant isotropic low permeability K b with K d K b = 10 4 . We consider again the simplified two phase flow model of subsection 3. We specify a pressure P 1 at the left side x = 0 and a pressure P 2 at right side x = 100 such that P 1 > P 2 . Homogeneous Neumann boundary conditions g = 0 are imposed at the remaining boundaries. The input gas saturation is set to u -= 1 at the input boundary x = 0. The reservoir is initially saturated with water.

The mesh is a coarse uniform Cartesian grid of size 100 × 1 × 5 with only one cell in the width of each layer as shown in Figure 3. It clearly shows that the solution provided by the VAG scheme is independent on the parameter ω and matches the solution of the TPFA scheme for the choice of the weights [START_REF] Eymard | Benchmark 3d: the vag scheme[END_REF]. On the contrary, the gas breakthrough obtained by the VAG scheme with the uniform weights is clearly delayed when the parameter ω, i.e. the pore volume at the vertices, increases. This is due to the fact that the total pore volume defined by the cells of the drains plus the vertices at the interface between the drains and the barriers is roughly independent of ω in the first case but increases with the parameter ω in the second case.

Grid orientation effect for a two phase flow on a Cartesian grid

The reservoir Ω = (-100, 100) × (0, 50) × (0, 45) m 3 is divided into three 15 m-thick layers as illustrated in Figure 5(a). The top and bottom layers are assumed to be geological barriers with very low permeability, whereas the medium layer is considered as a highly permeable drain. The ratio of the permeabilities between the drain and the barriers is set to 10 4 and the media is assumed to be isotropic. We consider an immiscible two-phase (gas g and water w) flow described by the following conservation equations

φ∂ t (ρ w (1 -S g )) -∇ • ρ w k rw (1 -S g ) µ w Λ∇P = 0, φ∂ t (ρ g S g ) -∇ • ρ g k rg (S g ) µ g Λ∇P = 0
, where • φ is the porosity, • S g is the saturation of the gas phase,

• ρ α and µ α are respectively the constant density and viscosity of the phase α = w, g • k rα is the relative permeability of the phase α = w, g.

The gravity and capillary effects are neglected. The ratio of the water and gas phases viscosities is set to 10, and we consider cross relative permeabilities k rα (S) = S α , α = g, w. To close the problem, a pressure boundary condition P ext is imposed at the boundaries x = ±100 and an homogeneous Neumann boundary condition is imposed at the remaining boundaries. The gas phase is injected into a reservoir initially saturated with water through a well localized at the centre of the reservoir. The well is defined by one perforated cell at the centre of the domain and is discretized using the Peaceman well index [START_REF]Fundamentals of Numerical Reservoir Simulations[END_REF] with a fixed well pressure P well such that P well > P ext .

A uniform Cartesian grid of size 100 × 10 × 15 has been used for the simulation. The results obtained with the VAG scheme are compared in Figures 6, 7 and 8 with those obtained with the O scheme, which degenerates to the TPFA scheme in the case of a Cartesian grid and an isotropic permeability. We clearly notice that the gas saturation front computed by the TPFA scheme is spread in the horizontal and vertical directions x, y of the grid, while the front exhibited by the VAG scheme is clearly radial as could be expected. This numerical diffusion along the axes of the mesh is a well-known phenomena called the Grid Orientation Effect (GOE) (see [START_REF] Vinsome | One approach to the grid orientation problem in reservoir simulation[END_REF], [START_REF] Yanosik | A nine-point, finite-difference reservoir simulator for realistic prediction of adverse mobility ratio displacements[END_REF]). The GOE appears due to viscous instabilities when a fluid with low viscosity is injected into a viscous fluid as it is the case here. The better behaviour of the VAG scheme compared with the TPFA scheme with respect to the GOE is due to the richer flux stencil provided by the VAG scheme.

Near-well injection of CO 2 in an aquifer

In this last example, we consider the numerical simulation of a two-component (H 2 O and CO 2 ) two-phase (water w and gas g) Darcy flow with solubility of the CO 2 component in the water phase. The media is anisotropic and homogeneous with diagonal permeability tensor in the (x, y, z) coordinates, and we consider an injection of CO 2 from a deviated well into the reservoir saturated with water. The capillary effects are neglected.

We use two types of 3D near-well grids. Each of them is based on a hexahedral radial mesh that is exponentially refined down to the well boundary in order to take into account the singular pressure distribution close to the well. This radial local refinement implies to build a matching mesh between the radial grid and the reservoir CPG grid using either hexahedra or both tetrahedra and pyramids as seen in Figure 9. This transition zone enables to couple the near-well simulation and a global reservoir simulation. The number of cells is roughly the same for both grids as shown in Table 1. cell unknowns nodal unknowns hexahedral mesh 74 679 74 800 hybrid mesh 77 599 -with 28 704 tetrahedra 37 883 Table 1. Numbers of cell and nodal unknowns for both near-well meshes.

There is to our knowledge no analytical solution to this problem. Thus, this work focuses on the comparison of the solution obtained by the VAG scheme and the O scheme on each grid, both in terms of accuracy and efficiency.

The model accounts for the following conservation equations of the two components

H 2 O : φ∂ t ρ w S w C w H 2 O -∇ • ρ w C w H 2 O k rw (S w ) µ w Λ(∇P -ρ w g) = 0, CO 2 : φ∂ t ρ w S w C w CO 2 + ρ g S g -∇ • ρ w C w CO 2 k rw (S w ) µ w Λ(∇P -ρ w g) + ρ g k rg (S g ) µ g Λ(∇P -ρ g g) = 0
, where • φ is the porosity with φ = 0.1,

• g is the gravity with g = (0, 0, 10) T m.s -2 ,

• C α i is the mass fraction of the component i in the phase α, • S α , ρ α and µ α are respectively the saturation, the density and the viscosity of the phase α, with ρ w = 1000 kg.m -3 , ρ g = 500 kg.m -3 , µ w = 3. 10 -4 Pa.s and µ g = 3. 10 -5 Pa.s, • k rα (S α ) is the relative permeability of the phase α defined here by k rα (S α ) = S α , • Λ is the diagonal permeability tensor with diagonal entries λ x = λ y = 1. 10 -13 m 2 and λ z = λx 10 m 2 .

The problem is closed by the volume balance and the sum to one of the mass fractions

S w + S g = 1, C w H 2 O + C w CO 2 = 1,
as well as the thermodynamic equilibrium defined by the solubility C and the following complementary conditions (C w CO 2 -C) ≤ 0, S g ≥ 0, (C w CO 2 -C)S g = 0.

In the following experiments the solubility C is set to the constant value 0.05. The CO 2 is injected from the deviated well in gaseous state inside the reservoir saturated with water of composition C w The final simulation time is fixed to 4 hours so that the gas phase has reached all sides of the external domain. In order to obtain a good convergence behavior in time, the time step is fixed to 150 seconds for both families of meshes and both schemes.

On Figure 10 has been plotted, in function of time, the rate of variation of the total mass of CO 2 in the reservoir. About the post-processing, the mass of CO 2 is computed by considering all the control volumes. Thus it means all the cells for the O scheme, and all the cells plus all the vertices for the VAG scheme. It can be observed that the VAG scheme solution is less sensitive to the type of the grid than the MPFA O scheme solution. Also note that several values of ω (8) have been tested without significant influence on the solution of the VAG scheme.

Conclusion

We have extended in this paper the Vertex Approximate Gradient (VAG) scheme to multiphase Darcy flows and proved its convergence for a simplified two phase Darcy flow model coupling an elliptic equation for the pressure and a linear hyperbolic equation for the saturation.

The VAG scheme has several advantages compared with usual cell centred approaches for multiphase Darcy flows. The scheme is coercive for general meshes and general permeability tensors, it leads to a compact 27 points stencil on topologically Cartesian meshes, and it reduces considerably the number of unknowns in the case of tetrahedral meshes. Compared with usual finite element approaches, the VAG scheme has the ability to deal with highly heterogeneous media on coarse meshes due to its flexibility in the definition of the porous volumes, and it can be easily implemented for more complex models in existing reservoir simulators based on a graph of transmissibilities data structure.

Figure 1 .

 1 Figure 1. Tetrahedron T κ,σ,e of the sub-mesh T .

2 .

 2 s∈Vκ |F κ,s (P )| + s∈V -g -s 1 It remains to estimate E = κ∈M s∈Vκ |F κ,s (P )|. For κ ∈ M, let A κ denote the symmetric positive definite matrix a s κ,s s,s ∈Vκ

κ∈M ( 11 . 4 )

 114 ϕ n κ -(11.1) u n κ ϕ n κ + s∈V (11.5) ϕ n s -(11.2) u n s ϕ n s Multiplying all terms by ∆t n and summing over n = 1, • • • , N , we obtain that

Figure 2 .

 2 Figure 2. Propagation in the horizontal direction of the gas saturation.

Figure 3 .

 3 Figure 3. Mesh and layers. In red the drains, and in green the barriers.

Figure 4

 4 Figure 4 exhibits the evolution of the cumulative gas flow rate at the right boundary, using either the weights α s κ defined by (8) in subfigure 4(a) or the uniform weights α s κ = 1 #Ms in subfigure 4(b). It is compared with the solution obtained with the TPFA scheme on both subfigures.It clearly shows that the solution provided by the VAG scheme is independent on the parameter ω and matches the solution of the TPFA scheme for the choice of the weights (8). On the contrary, the gas breakthrough obtained by the VAG scheme with the uniform weights is clearly delayed when the parameter ω, i.e. the pore volume at the vertices, increases. This is due to the fact that the total pore volume defined by the cells of the drains plus the vertices at the interface between the drains and the barriers is roughly independent of ω in the first case but increases with the parameter ω in the second case.

Figure 4 .

 4 Figure 4. Cumulative gas flow rate at the right boundary function of time.

  (a) Mesh and layers. In red the drain, and in green the barriers. (b) xz 2D view at y = 25. The well perforated cell is depicted by the blue square.

Figure 5 .

 5 Figure 5. Configuration of the test case of injection of gas in an heterogeneous reservoir.

Figure 6 .

 6 Figure 6. Gas saturation such that S g > 0.1% after a short time of injection.

Figure 7 .

 7 Figure 7. Cut at z = 22.5 m -Gas saturation at the end of the simulation.

Figure 8 .

 8 Figure 8. Cut at y = 25 m -Gas saturation at the end of the simulation.

  (a) unstructured mesh with only hexahedra (b) hybrid mesh with hexahedra, tetrahedra and pyramids

Figure 9 .

 9 Figure 9. Near-well meshes.

H 2 O

 2 = 1. Homogeneous Neumann boundary conditions are set on the north and south sides for both phases and Dirichlet boundary conditions are set for all other sides. Hydrostatic pressure condition is imposed on the outer Dirichlet boundary and in the reservoir at the initial state P (x, y, z) = P 1 + x 30 (P 2 -P 1 ) -ρ w g • x, with P 1 = 240 10 +5 Pa, P 2 = P 1 -1 10 +4 Pa and (x, y, z) ∈ [-15, 15] × [-15, 15] × [-7.5, 7.5]. The initial saturations and compositions of the reservoir are those of pure water : S w = 1 and C w H 2 O = 1. The input fluid at the outer Dirichlet boundary is also pure water. To complete the description of the test case, hydrostatic pressure condition is also imposed at the well boundary P (x, y, z) = P well -ρ g g • x, with P well = 250 10 +5 Pa. As mentioned above the injected fluid at the well boundary is pure gas : S g = 1 and C w CO 2 = 1.

Figure 10 .

 10 Figure 10. Total mass flow rate of CO 2 in the reservoir in function of time.