R Eymard

C Guichard
email: cindy.guichard@unice.fr

R Herbin

R Masson

Vertex-centred Discretization of Multiphase Compositional Darcy flows on General Meshes

. It leads to an unconditionally coercive scheme for arbitrary meshes and permeability tensors. The stencil of this vertex-centred scheme typically comprises 27 points on topologically Cartesian meshes, and the number of unknowns on tetrahedral meshes is considerably reduced, compared with usual cell-centred approaches. The efficiency of our approach is exhibited on the nearwell injection of miscible CO 2 in a saline aquifer taking into account the vaporization of H 2 O in the gas phase as well as the precipitation of salt.

Introduction

Many applications require the simulation of compositional multiphase Darcy flow in heterogeneous porous media. In oil reservoir modelling, the compositional triphase Darcy flow simulator is a key tool to predict and optimize the production of a reservoir. In sedimentary basin modelling, such models are used to simulate the migration of the oil and gas phases in a basin saturated with water at geological space and time scales. The objectives are to predict the location of the potential reservoirs as well as the quality and quantity of oil trapped therein. In CO 2 geological storage, compositional multiphase Darcy flow models are used to optimize the injection of CO 2 and to assess the long term integrity of the storage. Finally, two-phase compositional Darcy flow models are used to study the gas migration in nuclear waste repositories (which consist of porous media saturated with water), and to assess the safety of the storage.

The numerical simulation of such models first relies on a proper formulation which can account for the coupling of the mass and pore volume conservations together with the thermodynamical equilibrium. A major difficulty is the management of phase appearance and disappearance induced by the change of phase reactions assumed to be at equilibrium. Many formulations have been proposed mainly in the oil industry (see [START_REF] Cao | Development of techniques for general purpose simulators[END_REF] and the numerous references therein), and more recently for the modelling of gas migration in nuclear waste disposal (see for example [START_REF] Abadpour | Method of negative saturations for modelling two phase compositional flows with oversaturated zones[END_REF], [START_REF] Angelini | Étude de schémas numériques pour les écoulements diphasiques en milieu poreux déformable pour des maillages quelconques: application au stockage de déchets radioactifs[END_REF], [START_REF] Bourgeat | Two phase partially miscible flow and transport in porous media; application to gas migration in nuclear waste repository[END_REF]).

In the following, we propose an extension to an arbitrary number of phases of Coats formulation [START_REF] Coats | Implicit compositional simulation of singleporosity and dual-porosity reservoirs[END_REF] designed for reservoir simulation. This formulation, which is expressed in the natural variables of the thermody-namical and hydrodynamical laws, has the main advantage of dealing with a large range of models. The first resulting difficulty is the complex management of a variable number of equations and unknowns at each point of the computational domain.

The second difficulty is to exhibit a space-time discretization, coping with the strong coupling of both an elliptic (or a parabolic) unknown, the pressure, and hyperbolic (or degenerate parabolic) unknowns, the volume and mole fractions. The standard industrial answer is based on cell-centred finite volume schemes, which can be efficiently combined with an Euler implicit time integration to allow for large time steps and accurately capture the change of phases fronts.

The main drawback of cell-centred finite volume schemes is the difficulty to use them in the case of complex meshes and heterogeneous anisotropic permeabilities, which are unfortunately frequently encountered in practice to represent the basin and reservoir geometries and petrophysical properties.

Many progresses have been done during the last decade, leading to the design of cell-centred schemes which remain consistent in such situations. Among these, let us mention for example the well-known O scheme introduced in [START_REF] Aavatsmark | Discretization on non-orthogonal, quadrilateral grids for inhomogeneous, anisotropic media[END_REF], [START_REF] Aavatsmark | Discretization on unstructured grids for inhomogeneous, anisotropic media. Part I: Derivation of the methods[END_REF], [START_REF] Edwards | A flux continuous scheme for the full tensor pressure equation[END_REF], [START_REF] Edwards | Finite volume discretization with imposed flux continuity for the general tensor pressure equation[END_REF], the L scheme [START_REF] Aavatsmark | A compact multipoint flux approximation method with improved robustness[END_REF], or the Sushi scheme [START_REF] Eymard | Discretisation of heterogeneous and anisotropic diffusion problems on general non-conforming meshes, SUSHI: a scheme using stabilisation and hybrid interfaces[END_REF]. We also refer to [START_REF] Agélas | The G method for heterogeneous anisotropic diffusion on general meshes[END_REF], [START_REF] Agélas | An abstract analysis framework for nonconforming approximations of diffusion on general meshes[END_REF], and [START_REF] Agélas | Convergence of finite volume MPFA O type schemes for heterogeneous anisotropic diffusion problems on general meshes[END_REF] for the convergence analysis of cell-centred schemes in a general framework. Nevertheless it is still a challenge to design a cell-centred discretization of diffusive fluxes in the case of general meshes and heterogeneous anisotropic diffusion tensors, which is linear, unconditionally coercive and "compact", in the sense that the expression of the discrete flux at a given face of the mesh may at most involve the cells sharing a vertex with this face. For example the O and L schemes are compact but their coercivity is mesh and permeability tensor dependent. On the other hand, the Sushi scheme is unconditionally coercive but its stencil is not compact, since it involves the neighbours of the neighbours of a given control volume.

Recently, a new discretization of diffusive equations, the Vertex Approximate Gradient (VAG) scheme, using both cell and vertex unknowns, has been introduced in [START_REF] Eymard | Smallstencil 3D schemes for diffusive flows in porous media[END_REF]. The cell unknowns can be eliminated locally without any fill-in, leading to a compact vertex-centred scheme, with a typical 27 points stencil for 3D topologically Cartesian meshes. The VAG scheme is consistent, unconditionally coercive, compact, and easy to implement on general meshes and for heterogeneous anisotropic diffusion tensors. In addition, it is exact on cellwise affine solutions for cellwise constant diffusion tensors. It has exhibited a very good compromise between accuracy, robustness and CPU time in the recent FVCA6 3D benchmark [START_REF] Eymard | Benchmark 3d on discretization schemes for anisotropic diffusion problem on general grids[END_REF].

This paper is devoted to the generalization of the VAG scheme to the discretization of multiphase compositional Darcy flows. The first main idea is to remark that conservative generalized fluxes, between the centre of each cell of the mesh and its vertices, may be used as Multi-Point Flux Approximations in the compositional multiphase Darcy flow framework. The second main idea is a new procedure for assigning an amount of pore volume to the vertices, without any loss of accuracy on coarse heterogeneous meshes and taking into account different rocktypes with discontinuous capillary pressures (see subsection 4.5).

The outline of the paper is the following. Section 2 introduces a general formulation for compositional multiphase Darcy flow models, accounting for an arbitrary number of phases and phase appearance and disappearance. Section 3 details the vertex-centred discretization of compositional multiphase flow models. The VAG scheme is first recalled for diffusive equations in subsection 3.1. Then, the VAG scheme fluxes are derived and used in subsection 3.2 to discretize the compositional multiphase Darcy flow model, and the pore volume assignment procedure is detailed. Subsection 3.3 briefly discusses the algorithms to solve the nonlinear and the linear systems arising from the VAG discretization of the compositional models. The last section 4 exhibits the efficiency of the VAG discretization which is compared to the solutions obtained with the MPFA O scheme. The first test cases deal with two phase flow examples including highly heterogeneous cases and discontinuous capillary pressures. Then, the last test case considers the nearwell injection of miscible CO 2 in a saline aquifer, taking into account the vaporization of H 2 O in the gas phase as well as the precipitation of salt.

Formulation of multiphase compositional Darcy flows

We consider in this section a generalisation to an arbitrary number of phases of Coats' formulation [START_REF] Coats | Implicit compositional simulation of singleporosity and dual-porosity reservoirs[END_REF] for compositional multiphase Darcy flow models. Let P denote the set of possible phases α ∈ P, and C the set of components i ∈ C. Each phase α ∈ P is described by its non empty subset of components C α ⊂ C in the sense that it contains the chemical species X α i for i ∈ C α . It is assumed that, for any i ∈ C, the set

P i = {α ∈ P | i ∈ C α }
of phases containing the component i, is non empty.

Each phase α ∈ P is characterized by its thermodynamical properties depending on its pressure P α and its molar composition

C α = (C α i ) i∈C α .
The dependence on the temperature will not be specified in the following since, for the sake of simplicity, we only consider isothermal flows. In fact, the extension to thermal Darcy flows follows the same basic ideas and typically involves an additional unknown, the temperature T , assumed to be identical for all phases and for the porous media, as well as an additional global energy conservation equation for the fluids and the porous media.

In order to decouple the local thermodynamical equilibrium from the hydrodynamics, we shall make the simplifying assumption that they only depend on a reference pressure, denoted by P , usually defined by the pressure of a given phase or the average of the phase pressures weighted by their volume fractions. Then, for each phase α ∈ P, we shall denote by ζ α (P, C α ) its molar density, by ρ α (P, C α ) its mass density, by µ α (P, C α ) its viscosity, and by f α i (P, C α ) its fugacity coefficients for the components i ∈ C α .

The phases can exchange mass, according to the change of phases reactions

X α i X β i for all (α, β) ∈ (P i ) 2 , α = β, i ∈ C.
It results that phases can appear or disappear, and we shall denote by Q the unknown representing the set of present phases, valued in the set Q of all non empty subsets of P.

These reactions are assumed to be at equilibrium, by stating that for any i ∈ C

f α i (P, C α )C α i = f β i (P, C β )C β i ,
for any couple of present phases α and β containing the component i, i.e. such that (α, β)

∈ (Q ∩ P i ) 2 , α = β.
For a given set of present phases Q, it may occur that a component i ∈ C is absent of all present phases α ∈ Q. Hence, we define the set of absent components as a function of Q by

C Q = {i ∈ C | Q ∩ P i = ∅}.
In the following, n i will denote the number of moles of the component i ∈ C per unit volume. It is considered as an independent unknown for i ∈ C Q and is equal to

n i = φ α∈Q∩Pi ζ α (P, C α ) S α C α i ,
for i ∈ C \ C Q , where φ > 0 denotes the porosity of the porous media.

Let us also define the vector of the total molar fractions of the components

Z = n i j∈C n j i∈C .
For a prescribed reference pressure P and given total component molar fractions Z, the so called flash computes the set of present phases Q, their molar fractions θ α , and their compositions C α for α ∈ Q. It solves the following local conservation of the number of moles and the equilibrium equations

               Z i = α∈Q∩Pi θ α C α i , i ∈ C, i∈C α C α i = 1, α ∈ Q, f α i (P, C α )C α i -f β i (P, C β )C β i = 0, α = β, (α, β) ∈ (Q ∩ P i ) 2 , i ∈ C, (1)
together with the stability of the present phases Q in the sense that the system achieves a global minimum of the Gibbs free energy. It is usually obtained using negative flashes [START_REF] Whitson | The negative flash[END_REF] for the possible values of Q ∈ Q or alternatively a stability analysis [START_REF] Michelsen | The isothermal flash problem. Part I: Stability[END_REF] followed by the solution of the system (1). It will be denoted in the following by

(Q, θ α , C α , α ∈ Q) = Flash(P, Z), or by Q = flash(P, Z),
if we only retain the present phases as output.

The hydrodynamical properties of the multiphase Darcy flow system are the capillary pressures and the relative permeabilities. For simplicity in the notations (but without restrictions) they are assumed to depend only on the volume fractions or saturations of the phases denoted by S = (S α ) α∈P .

The capillary pressures P c,α (S) are such that

P α = P + P c,α (S), (2)
for all α ∈ P, and the relative permeabilities are denoted by k rα (S) for all α ∈ P.

For all present phases α ∈ Q, the multiphase Darcy velocity is defined by

U α = k rα (S) µ α (P, C α ) V α , with V α = -Λ (∇P α -ρ α (P, C α )g) ,
and where Λ denotes the permeability tensor, and g the gravity vector.

Let Ω be a bounded polyhedral subdomain of R 3 of boundary ∂Ω = Ω \ Ω, and let (0, t f ) denote the time interval. The set of unknowns first includes the set of the present phases Q, the reference pressure P , the saturations and compositions of the present phases S α , C α for α ∈ Q, and the number of moles per unit volume n i for the absent components i ∈ C Q . Note that in the remaining of the paper, it is implicitly assumed that S α = 0 for all α ∈ Q.

The system of PDEs, accounting for the conservation of the number of moles per unit volume and the multiphase Darcy laws, has to be solved together with the local algebraic closure laws accounting for the conservation of the pore volume, the thermodynamical equilibrium, and the stability of the present phases. The phase pressures P α , α ∈ Q are defined as function of the reference pressure P and the saturations S using the capillary relations [START_REF] Aavatsmark | Discretization on unstructured grids for inhomogeneous, anisotropic media. Part I: Derivation of the methods[END_REF]. We end up with the following system of equations which is set on the domain Ω × (0, t f ):

                                               ∂ t n i + div α∈Q∩Pi C α i ζ α (P, C α )k rα (S) µ α (P, C α ) V α = 0, i ∈ C, V α = -Λ ∇P α -ρ α (P, C α )g , α ∈ Q, P α = P + P c,α (S), α ∈ Q, α∈Q S α = 1, i∈C α C α i = 1, α ∈ Q, f α i (P, C α )C α i -f β i (P, C β )C β i = 0, α = β, (α, β) ∈ (Q ∩ P i ) 2 , i ∈ C, Q = flash(P, Z). (3)
To avoid redundancy in the definition of the system of equations (3), we have only retained the set of phases Q as output from the flash computation, although in practice one may also use the phase molar fractions and compositions outputs in the nonlinear solver of the discrete system. Alternatively to the flash computation, one could use a stability analysis [START_REF] Michelsen | The isothermal flash problem. Part I: Stability[END_REF] of the set of phases Q to check if additional phases may appear or not.

For the sake of simplicity, we shall only assume here boundary conditions of the Dirichlet or Neumann types. The disjoint subsets ∂Ω D and ∂Ω N of the boundary ∂Ω are such that ∂Ω D ∪ ∂Ω N = ∂Ω. The normal vector at the boundary outward the domain Ω is denoted by n.

On the Dirichlet boundary ∂Ω D , the pressure P is specified as well as the input phases α ∈ Q D , their volume fractions S α D , α ∈ Q D , and their compositions C α D , α ∈ Q D , assumed to be at equilibrium:

on ∂Ω D :    P = P D , S = S D , C α = C α D for α ∈ Q D if V α • n < 0.
On the Neumann boundary ∂Ω N , input fluxes g i ≤ 0 are prescribed for each component i ∈ C:

on ∂Ω N : α∈Q∩Pi C α i ζ α (P, C α )k rα (S) µ α (P, C α ) V α • n = g i , i ∈ C.
3 Vertex-centred discretization on generalised polyhedral meshes 3.1 Vertex-centred discretization of Darcy fluxes (VAG scheme)

For a.e. (almost every) x ∈ Ω, Λ(x) is assumed to be a 3-dimensional symmetric positive definite matrix such 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 = Ω f v dx - ∂Ω N g v dσ 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 | Smallstencil 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, for κ ∈ M and s ∈ V},
and its subspace with homogeneous Dirichlet boundary conditions on

V D W D = {v ∈ W D | v s = 0 for s ∈ V D }.

Vertex Approximate Gradient (VAG) scheme

The VAG scheme introduced in [START_REF] Eymard | Smallstencil 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 | Smallstencil 3D schemes for diffusive flows in porous media[END_REF][START_REF] Eymard | Benchmark 3D: the VAG scheme[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 for 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 ((η κ ) κ∈M , (η s ) s∈V ).

Following [START_REF] Eymard | Smallstencil 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 for all v ∈ W D ,

a D (u, v) = Ω f (x) v T (x) dx - ∂Ω g(x) v T (x) dσ,
where a D is the bilinear form defined by

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 generalized 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 generalized 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 balance equations:

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

g s = ∂Ω g(x)η s (x) dσ.
One may easily check that the VAG scheme generalized fluxes can be rewritten as follows:

F κ,s (u) = κ -Λ(x)∇u T (x) • ∇η s (x) dx. (4)
This formula provides in 2D an interpretation of the VAG scheme generalized fluxes as Control Volume Finite Element (CVFE) fluxes [START_REF] Chen | On the control volume finite element methods and their applications to multiphase flow[END_REF] on a triangular submesh. In 2D, x σ is chosen to be the mid-point of the edge σ = ss , and the interpolation is simply defined by

I σ (v) = vs+v s 2
. It results that the triangular submesh T is rather defined in 2D as the set of triangles obtained for each cell κ by joining each face σ of the cell κ to the cell centre x κ . Then, the VAG scheme reduces to the P 1 finite element scheme on the submesh T . Let a be the mid-point of sx κ (see figure 2). Using (4), one easily shows that, assuming a cellwise constant tensor field Λ, one has

F κ,s (u) = xσa ∪ x σ a -Λ(x)∇u T (x) • n κ dσ,
where x σ a (respectively x σ a) is any curved segment inside the triangle σx κ (respectively σ x κ ) and n κ is the normal outward the CVFE control volume containing x κ (see figure 2). The flexibility in the definition of the curved edges of the CVFE control volumes will be exploited in the next subsection to adapt the porous volume at the vertices and at the cell centres of the mesh in heterogeneous cases. This geometrical CVFE interpretation of the VAG generalized fluxes cannot be extended to the 3D case due to the interpolation at the face centres which has been used to avoid numerous additional unknowns at the faces.

Discretization of multiphase compositional Darcy flows

The VAG scheme was shown to have good approximation properties for single phase Darcy flows on general meshes for instance 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].

In the case of the discretization of multiphase Darcy flow models on general meshes, the first idea is then to use the generalized fluxes F κ,s (u) = -F s,κ (u) between a cell κ of the mesh and its vertices s ∈ V κ as a Multi-Point Flux Approximation. Therefore, the set of control volumes K is defined as the union of the cells and of the interior and Neumann boundary vertices

K = M ∪ V int ∪ V N ,
These generalized fluxes F κ,s (u) between the control volumes κ and s are then classically used for the approximation of the transport terms, in addition to an upwind scheme. Although these fluxes are not defined in the usual way, that is as the approximation of the continuous fluxes σ -Λ∇ū • n σ dσ on a given face σ of the mesh, the mathematical analysis developed in [START_REF] Eymard | Vertex centred Discretization of Two-Phase Darcy flows on General Meshes[END_REF] shows that they lead to a convergent scheme, at least in a particular two-phase flow case. It is proved for the decoupled case where the sum of the mobilities is independent on the saturation, that the discrete saturation converges weakly in L ∞ to the weak solution of the saturation equation. The proof follows the lines of [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] dealing with CVFE schemes (see also [START_REF] Eymard | Finite volume methods[END_REF] for TPFA schemes) using a weak BV estimate for the VAG generalized fluxes together with the finite element variational formulation of the pressure equation.

The second ingredient is the assignment of a porous volume φ K to each control volume K ∈ K such that

K∈K φ K = Ω φ(x)dx and φ K > 0 for all K ∈ K. (5)
It is achieved by a conservative redistribution to the vertices of the surrounding cell porous volumes

       φ s = ω κ∈Ms α s κ Φ κ for all s ∈ V int ∪ V N , φ κ = 1 -ω s∈Vκ\V D α s κ Φ κ for all κ ∈ M, (6)
with Φ κ = κ φ(x)dx, and α s κ ≥ 0, κ∈Ms α s κ = 1, which guarantees [START_REF] Achdou | Block preconditioning and multigrid solvers for linear systems in reservoir simulations[END_REF] provided that the parameter ω > 0 is chosen small enough.

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

α s κ = a κ,s κ ∈Ms a κ ,s , (7)
for all s ∈ V int ∪ V N and κ ∈ M s with

a κ,s = s ∈Vκ\V D 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 which are used in practical situations involving highly heterogeneous media.

Let K = K ∪ V D = M ∪ V,
denote the union of the set of control volumes K and of the Dirichlet boundary vertices, and let us set

P s = P D (s), S s = S D (s), Q s = Q D (s), C α s = C α D (s), α ∈ Q s , for all s ∈ V D .
Then, on each control volume K ∈ K, the set of discrete unknowns and Dirichlet data is denoted by

X K = Q K , P K , (S α K ) α∈Q K , (C α K ) α∈Q K , (n i,K ) i∈C Q K .
We also introduce for all κ ∈ M the vector of unknowns and Dirichlet data

X Vκ = X s , s ∈ V κ .
For all K ∈ K, let us define n i (X K ), the number of moles of the component i in the control volume K, equal to

φ K α∈Q K ∩Pi ζ α (P K , C α K ) S α K C α i,K , for all i ∈ C \ C Q K and to the independent unknown n i,K for i ∈ C Q K .
The phase pressures at each control volume K ∈ K are defined by

P α K = P K + P c,α (S K ),
for all phases α ∈ P (not only for present phases α ∈ Q K ). Then, for all s ∈ V κ , κ ∈ M, and for all phases α ∈ Q κ ∪ Q s , we define the Darcy fluxes

V α κ,s (X κ , X Vκ ) = -V α s,κ (X κ , X Vκ ) = s ∈Vκ a s κ,s P κ -P s + P c,α (S κ ) -P c,α (S s ) +ρ α κ,s (X κ , X s ) g (z κ -z s ) ,
with z κ (respectively z s ) denoting the vertical coordinate of x κ (respectively s ), and

ρ α κ,s (X κ , X s ) = S α κ S α κ + S α s ρ α P κ , C α κ + S α s S α κ + S α s ρ α P s , C α s .
Let us point out that the function ρ α κ,s (X κ , X s ) is not continuous at S α κ = S α s = 0, but its product with the relative permeability k rα (S K α κ,s ) in the upwind control volume defined below is continuous since the relative permeability k rα (S) of the phase α vanishes for S α = 0.

The transport terms are approximated using a first order upwind scheme. The upwinding is done as usual for each phase with respect to the sign of its Darcy velocity (see [START_REF] Aziz | Petroleum Reservoir Simulation[END_REF], [START_REF] Peaceman | Fundamentals of Numerical Reservoir Simulations[END_REF]). For all s ∈ V κ , κ ∈ M, and

α ∈ Q κ ∪ Q s we set K α κ,s = κ if V α κ,s (X κ , X Vκ ) ≥ 0, s if V α κ,s (X κ , X Vκ ) < 0.
The time integration scheme implicitly couples the local algebraic closure laws to the system of PDEs. For the fluxes one can use different Euler integrations, either fully implicit [START_REF] Coats | An equation of state compositional model[END_REF], or implicit in pressure, explicit in saturations and compositions (ImPES) [START_REF] Aziz | Petroleum Reservoir Simulation[END_REF], [START_REF] Peaceman | Fundamentals of Numerical Reservoir Simulations[END_REF], or implicit in pressure and saturations, explicit in compositions (ImPSat) [START_REF] Quandalle | An implicit in pressure and saturations approach to fully compositional simulation[END_REF], or also adaptive implicit (AIm) [START_REF] Thomas | Reservoir simulation using and adaptive implicit method[END_REF]. Here, for the sake of simplicity, we will restrict the presentation to the fully implicit integration case.

Let us set for all K ∈ K

C α i ζ α k rα µ α X K = C α i,K ζ α (P K , C α K ) k rα (S K ) µ α (P K , C α K )
, and let us define for all

s ∈ V κ , κ ∈ M, α ∈ Q κ ∪ Q s , the fluxes
G α i,κ,s = -G α i,s,κ = C α i ζ α k rα µ α X n K α,n κ,s V α κ,s (X n κ , X n Vκ ).
The discrete system couples the following set of equations

the discrete balance equations on each cell κ ∈ M and for each component i ∈ C:

n i (X n κ ) -n i (X n-1 κ ) ∆t n + s∈Vκ α∈Q n K α,n κ,s ∩Pi G α i,κ,s = 0, (8)
the discrete balance equations on each vertex s ∈ V int ∪ V N and for each component i ∈ C:

n i (X n s ) -n i (X n-1 s ) ∆t n + κ∈Ms α∈Q n K α,n κ,s ∩Pi G α i,s,κ + g n i,s = 0, (9)
the local closure laws for all control volume K ∈ K:

               α∈Q n K S α,n K -1 = 0, i∈C α C α,n i,K -1 = 0, α ∈ Q n K , f α i (P n K , C α,n K )C α,n i,K -f β i (P n K , C β,n K )C β,n i,K = 0, α = β, (α, β) ∈ (Q n K ∩ P i ) 2 , i ∈ C (10)
the flash for all control volume K ∈ K:

Q n K = flash(P n K , Z n K ), (11)
with

Z n K = n i (X n K ) j∈C n j (X n K ) i∈C .
It is easy to check that, at each time step n, the number of unknowns Q n K , K ∈ K matches the number of equations [START_REF] Brenner | A Convergent Finite Volume Scheme for Two-Phase Flows in Porous Media with Discontinuous Capillary Pressure Field[END_REF], and that the number of unknowns

P n K , (S α,n K ) α∈Q n K , (C α,n K ) α∈Q n K , (n n i,K ) i∈C Q n K for all K ∈ K equal to K∈K   1 + #Q n K + α∈Q n K #C α + #C Q n K   ,
matches the number of equations ( 8), ( 9), [START_REF] Aziz | Petroleum Reservoir Simulation[END_REF] equal to

K∈K   #C + 1 + #Q n K + i∈C\C Q n K (#P i ∩ Q n K -1)    , since one notices that i∈C\C Q n K (#P i ∩ Q n K -1) = α∈Q n K #C α + #C Q n K -#C.

Nonlinear solver

At each time step n = 1, • • • , N , the nonlinear system coupling the conservation equations ( 8), [START_REF] Angelini | Étude de schémas numériques pour les écoulements diphasiques en milieu poreux déformable pour des maillages quelconques: application au stockage de déchets radioactifs[END_REF] with the local closure laws ( 10) is solved using a Newton type algorithm combined with a fixed point update for the set of present phases given by [START_REF] Brenner | A Convergent Finite Volume Scheme for Two-Phase Flows in Porous Media with Discontinuous Capillary Pressure Field[END_REF]. In the following the superscript n will be dropped for the sake of the clarity of the notations.

Let us denote by

Y K = P K , (S α K ) α∈Q K , (C α K ) α∈Q K , (n i,K ) i∈C Q K the vector of unknowns in the control volume K ∈ K, by Y Vκ = (Y s , s ∈ V κ \ V D )
the vector of unknowns of the cell vertices, and by Y K = (Y K , K ∈ K) the vector of cell and vertex unknowns. The vector ȲK denotes the vector of unknowns obtained from Y K by excluding the absent components

n i,K , i ∈ C Q K .
Finally, we use similar notations for the sets of sets of phases Q Vκ and Q K .

Using these notations, the system (8), ( 9), [START_REF] Aziz | Petroleum Reservoir Simulation[END_REF] can be written as

             R κ Y κ , Y Vκ , Q κ , Q Vκ = 0, κ ∈ M, R s Y κ , Y Vκ , Q κ , Q Vκ , κ ∈ M s = 0, s ∈ V int ∪ V N , L ȲK , Q K = 0, K ∈ K, (12)
with straightforward definitions of the residual functions R κ , R s and L. Let us rewrite the full system as

R K (Y K , Q K ) = 0.
The Newton type algorithm iterates until convergence on the following steps for given stopping parameters , , and k max .

-Initial guess: usually taken as

Q (0) K = Q n-1 K , and
Y (0) K = Y n-1 K . -Compute the initial residual R K Y (0) K , Q (0) K and r (0) = R K Y (0) K , Q (0) K for a given weighted norm . . -Iterate on k = 0, • • • , k max until r (k) ≤ r (0) or r (k) ≤ or k = k max :
-Compute the Jacobian matrix

J (k) = ∂R K ∂Y K Y (k) K , Q (k) K .
-Solve the linear system

J (k) dY K = -R K Y (k) K , Q (k) K . -Update the unknowns Y K Y (k+ 1 2 ) K = Y (k) K + θ (k) dY K ,
with a full Newton step θ (k) = 1 or a possible relaxation θ (k) ∈ (0, 1).

-Update the set of phases

Q K for all K ∈ K Z (k+ 1 2 ) K = Projection {Ui≥0 | i∈C Ui=1}      n i (Y (k+ 1 2 ) K , Q (k) K ) j∈C n j (Y (k+ 1 2 ) K , Q (k) K )      i∈C , Q (k+1) K = flash P (k+ 1 2 ) K , Z (k+ 1 2 ) K .
The flash computation also provides the compositions and molar fractions of the present phases. They are used together with Y (k+ 1 2 ) K and Q

(k+1) K to update the new set of unknowns Y (k+1) K at step k + 1. -Compute the residual R K Y (k+1) K , Q (k+1) K and r (k+1) = R K Y (k+1) K , Q (k+1) K .
-If the iterations terminate before k = k max , then proceed to the next time step n + 1, otherwise recompute the time step n with a reduced ∆t n .

In view of the nonlinear system [START_REF] Bourgeat | Two phase partially miscible flow and transport in porous media; application to gas migration in nuclear waste repository[END_REF], the size of the linear system for the computation of the Newton step can be considerably reduced without fill-in first by elimination of the local closure laws [START_REF] Aziz | Petroleum Reservoir Simulation[END_REF], and second by elimination of the cell unknowns using equations [START_REF] Agélas | Convergence of finite volume MPFA O type schemes for heterogeneous anisotropic diffusion problems on general meshes[END_REF].

The elimination of the local closure laws is achieved in each control volume K ∈ K by a splitting of the unknowns ȲK into

N p K = #C -#C Q K primary unknowns Ȳ p K and N s K secondary unknowns Y s K with N s K = 1 + #Q K + α∈Q K #C α + #C Q K -#C.
For each K ∈ K, the secondary unknowns must be chosen in such a way that the square matrix

∂L ∂Y s K Ȳ p K , Y s K , Q K ,
of size N s K is non singular. This choice can be done algebraically in the general case, or defined once and for all for each set of phases Q K for specific cases such as single phase flows or simple two phase flows.

The reduced linear system is solved using an iterative solver such as GMRES or BiCGStab and a preconditioner adapted to the elliptic or parabolic nature of the pressure unknown and to the coupling with the remaining hyperbolic or degenerate parabolic unknowns. One of the most efficient preconditioners is the so-called CPR-AMG preconditioner introduced in [START_REF] Lacroix | Decoupling preconditioners in the Implicit Parallel Accurate Reservoir Simulator (IPARS)[END_REF] and [START_REF] Scheichl | Decoupling and block preconditioning for sedimentary basin simulations[END_REF]. It combines multiplicatively an Algebraic Multigrid preconditioner for a pressure block of the linear system [START_REF] Henson | BoomerAMG: A parallel algebraic multigrid solver and preconditioner[END_REF] with a more local preconditioner for the full system, such as an incomplete LU factorization. The choice of the pressure block is crucial for the efficiency of the CPR-AMG preconditioner, we refer to [START_REF] Lacroix | Decoupling preconditioners in the Implicit Parallel Accurate Reservoir Simulator (IPARS)[END_REF], [START_REF] Scheichl | Decoupling and block preconditioning for sedimentary basin simulations[END_REF], and [START_REF] Achdou | Block preconditioning and multigrid solvers for linear systems in reservoir simulations[END_REF] for a discussion of the possible choices.

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. Note that on the Cartesian meshes used below and with a permeability tensor aligned with the directions of the mesh, the MPFA O scheme reduces to the Two Point Flux Approximation scheme (TPFA). The first 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 heterogeneities and different rocktypes on the transport. The last test case focuses on phase appearance and disappearance. It simulates the nearwell injection of miscible CO 2 in a saline aquifer involving the vaporization of H 2 O and the precipitation of salt close to the well.

In all cases, the fully implicit time integration scheme described in subsection 3.2 is used and the linear systems arising from the Newton algorithm are solved using a GMRES iterative solver [START_REF] Saad | GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF] preconditioned by an ILU0 incomplete factorization with no fill-in preconditioner [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF].

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 ( 7), [START_REF] Agélas | The G method for heterogeneous anisotropic diffusion on general meshes[END_REF], the VAG scheme provides solutions which are just as accurate as the solu-tions 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 . The fluid model is a simple immiscible two-component (say 1 and 2) two-phase (say gas (g) and (w)) flow with C g = {1}, C w = {2}, no capillary effect, no gravity and the sum of the mobilities of both phases equal to one. Thus the model reduces to a hyperbolic equation for the gas saturation, still denoted by S g , coupled to a fixed elliptic equation for the pressure P . The porosity φ is constant, and the reservoir is initially saturated with water. A pressure P 1 is fixed at the left side x = 0 and a pressure P 2 at right side x = 100 such that P 1 > P 2 . The input gas saturation is set to S g D = 1 at the input boundary x = 0. Homogeneous Neumann boundary conditions g 1 = g 2 = 0 are imposed at the remaining boundaries.

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] Agélas | An abstract analysis framework for nonconforming approximations of diffusion on general meshes[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.

Two phase flow with a log normal permeability field

The next test case compares the TPFA and VAG schemes on a 2D immiscible two phase flow with a log normal permeability field and two different choices of the variance. The flow takes into account gravity effect, no capillary pressure, it considers cross relative permeabilities and a ratio of 10 of the water and gas viscosities. The gas is injected at the left end x = 0 of the (x, z) 2D rectangular domain initially saturated with water. Zero flow conditions are imposed at the top and bottom boundaries and the pressure is fixed at the left and right ends.

The mesh is a uniform Cartesian grid of size 100 × 50. Figure 6 exhibits the permeability field for the largest variance 5, the lowest variance is 1.

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

We consider an immiscible two-component (1 and 2) two-phase (gas g and water w) flow with C g = {1}, C w = {2}. The gas phase is injected into a reservoir Ω = (-100, 100)×(0, 50)×(0, 45) m 3 initially saturated with water through a well localized at the centre of the reservoir. The gas and water phases are assumed to be incompressible and the gravity and capillary effects are neglected. The ratio of the water and gas phases constant viscosities is set to 10, and we consider cross relative permeabilities k rα (S) = S α , α = g, w.

The reservoir is divided into three 15 m-thick layers as illustrated in Figure 8(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. A pressure output boundary condition P D is imposed at the boundaries x = ±100 and a homogeneous Neumann boundary condition is imposed at the remaining boundaries. The well is defined by one perforated cell at the centre of the domain with a specified gas flow rate.

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 9, 10 and 11 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, finitedifference 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 stencil of the saturation equation obtained for the VAG scheme (≤ 27 points after elimination of the cell unknowns) compared with the MPFA O scheme (≤ 7 points) on topologically Cartesian meshes.

Numerical diffusion and CPU time for a decoupled two phase flow

We consider a simple immiscible two-component (say 1 and 2) two-phase (say gas (g) and water (w)) flow with C g = {1}, C w = {2}, no capillary effect, no gravity and the sum of the mobilities of both phases equal to one. In such a case, the model reduces to a linear scalar hyperbolic equation for the gas saturation denoted by S g coupled to an elliptic equation for the pressure P . The simulation is done on the domain Ω = (0, 1) 3 with the permeability tensor Λ = I, a porosity φ = 1, and the initial gas saturation S g (x, 0) = 0. 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 1 = g 2 = 0 are imposed at the remaining boundaries. The input gas saturation is set to S g D = 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

S g (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 12 shows, for each grid, the projection of the gas saturation S g on the x-coordinate axis at the simulation time for which the gas has filled half of the reservoir. We have plotted the analytical solution S g and the discrete solutions (x κ , S g κ ) for all cells κ ∈ M obtained with the VAG and the MPFA O schemes. For the VAG scheme we use the post-processed values,

Sg κ = (1 -ω s∈Vκ\V D α s κ )S g κ + ω s∈Vκ\V D α s κ S g
s deduced from the redistribution of the volumes defined by ( 6)-( 7) . The results presented in Figure 12 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 12(a), the solutions of the VAG scheme on the Cartesian grid is plotted for both ω = 0.01 and ω = 0.3 [START_REF] Agélas | The G method for heterogeneous anisotropic diffusion on general meshes[END_REF]. The value ω = 0.3 roughly corresponds to match the pore volume of each vertex φ s with the pore volume of each surrounding cell φ κ , κ ∈ 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 TPFA scheme). On the other hand, for the tetrahedral mesh, we can notice on Figure 12(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 both ω = 0.01 and ω = 0.3.

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 with discontinuous capillary pressures

The objective of this 2D immiscible incompressible two phase flow test case is to assess the ability of the VAG scheme to deal with different rocktypes. In such cases, the main issue is how to define a capillary pressure (and similarly relative permeabilities) for each nodal control volume in order to allow for discontinuous saturations at the interface between two rocktypes.

The solution of this problem relies again on the flexibility of the VAG scheme in the definition of the porous volumes at the vertices. The rocktypes are denoted by rt κ for each cell κ ∈ M. For each node s ∈ V int ∪ V N we choose a rocktype rt s = rt κs with κ s ∈ argmax{a κ,s , κ ∈ M s } and we set

N s = {κ ∈ M s | rt κ = rt κs }.
The porous volumes ( 6) are defined with the following new definition of the weights [START_REF] Agélas | An abstract analysis framework for nonconforming approximations of diffusion on general meshes[END_REF] in such a way that the porous volume at the vertex s is taken only from the cells with the same rocktype rt κs .

α s κ =      aκ,s κ ∈Ns a κ ,s if s ∈ N s , 0 else . (13)
In the following test case, the domain Ω = (0, 100)× (0, 100) in the (x, z) plane is split into two layers Ω 1 = (0, 100) × (0, 50) and Ω 2 = (0, 100) × (50, 100), each corresponding to a rocktype rt = 1 and 2 respectively. The porosity φ = 0.2, and the permeability Λ = 10 -12 I are homogeneous and the same for both rocktypes. The relative permeabilities are set to k r,α (S α ) = (S α ) 2 , α = o, w, for both rocktypes while the capillary pressure of each rocktype is defined by

S o (x) = 0.3 if x ∈ Ω 1 , 0 if x ∈ Ω 2 .
Figures 13 and 14 compare at two different times the solutions obtained for the oil saturation along the z axis on the coarse grid 2 × 10 and on the fine grid 16 × 80 with the TPFA scheme and with three different VAG schemes. Since the permeability is homogeneous and the mesh uniform, there are two possible choices for the rocktypes at the vertices located at the interfaces between the two subdomains. The VAGs1 scheme (respectively VAGs2) is the one obtained with the choice of the rocktype 1 (respectively 2) at the interface between the two subdomains. The VAGa scheme is obtained with the weights defined by [START_REF] Agélas | An abstract analysis framework for nonconforming approximations of diffusion on general meshes[END_REF] and with a capillary pressure at the interfaces given by an harmonic average of both rocktype capillary pressures. In all cases, the parameter ω is computed such that the minimum porous volumes at vertices and at cells match. Note that the discrete solutions are in all cases independent on x although the VAG scheme does not degenerate to a 1D scheme while the TPFA scheme does for this test case. For the VAG schemes the following post-processed values of the oil saturation are plotted at the cell centre z κ along the z axis:

So κ = (1 -ω s∈Vκ\V D α s κ )S o κ + ω s∈Vκ\V D α s κ S o s .
It is known that oil can only flow by gravity to the top subdomain provided that the capillary pressures can achieve continuity at the interface, meaning here that the jump of the oil saturation at the interface must reach the value 0.5. It can be checked that it is the case for all schemes on the fine grid solution. On the coarse grid, the VAGs1 and TPFA scheme are very close, the VAGs2 scheme is slightly better, while the VAGa scheme, as could be expected, is much worse in capturing the saturation jump. 0, and the water and gas phase relative permeabilities k r,α are non decreasing functions of the reduced saturation S α = S α S w +S g , α = w, g.

k r,α ( S α ) =      S α -S r,α 1 -S r,α e α if S α ∈ [S r,α , 1] 0 if S α < S r,α ,
with e w = 5, e g = 2, S r,w = 0.3 and S r,g = 0. The reference pressure is chosen to be the gas pressure P = P g , and we set P w = P g + P c,w ( S w ) with

P c,w ( S w ) = P c,1 log S w -S r,w 1 -S r,w -P c,2 ,
for 1 ≥ S w > S r,w , where P c,1 = 2.10 4 Pa, P c,2 = 10 4 Pa. The thermodynamical equilibrium is modelled by equilibrium constants The equilibrium constants will be considered fixed in the range of pressure and temperature with the following values K H2O = 0.025, K CO2 = 0.03 and K salt = 0.39 in kg/kg. With these assumptions, the flash Q = flash P, Z admits an analytical solution, independent of the pressure P , and with entries Z in the two dimensional simplex

K i , i ∈ C such that                        C g H2O = K H2O C w
{(Z CO2 , Z salt ) | Z CO2 ≥ 0, Z salt ≥ 0, 1 -Z CO2 -Z salt ≥ 0}.
The solution is exhibited Figure 15 where we have set

E 1 = (E CO2 K CO2 , 0), E 2 = (E CO2 , 0), E 3 = (D CO2 , 0), E 4 = (0, K salt ), E 5 = (K CO2 D CO2 , K salt ), with D CO2 = 1 -K H2O (1 -K salt ) 1 -K H2O K CO2 ,
and

E CO2 = 1 -K H2O 1 -K H2O K CO2 .
The 3D nearwell grids used for the simulation are exhibited in Figure 16. The first step of the discretization is to create a radial mesh, Figure 16(a), that is exponentially refined down to the well boundary. This nearwell radial local refinement is matched with the reservoir Ω = (-15, 15) × (-15, 15) × (-7.5, 7.5) m 3 using either hexahedra (see Figure 16(b)) or both tetrahedra and pyramids, (see Figure 16(c)). The radius of the well is 10 cm and the radius of the radial zone is 5 m. The well is deviated by an angle of 20 degrees away from the vertical axis z in the x, z plane. The hexahe- The remaining of the data set and the boundary and initial conditions are the following. The porosity is set to φ = 0.2, and the permeability tensor Λ is homogeneous and isotropic equal to 1. 10 -12 m 2 . The density and the viscosity of the water phase are computed by correlations function of P and C w , and those of the gas phase by linear interpolation in the pressure P , the density of the mineral phase is fixed to ρ m = 2173 kg/l.

Homogeneous Neumann boundary conditions are imposed on the top and bottom boundaries. Along the well boundary we impose the pressure P (x, y, z) = P wellρ g g z, with P well = 300 10 +5 Pa, the input phase S g = 1 and its composition C g CO2 = 1, C g H2O = 0. On the lateral outer boundaries (resp. at initial time) the following hydrostatic pressure is imposed P (x, y, z) = P 1 -ρ l g z, with P 1 = P well -10 +5 Pa, as well as the following input (resp. initial) phase and its composition S w = 1, C w H2O = 0.84 and C w salt = 0.16.

The simulation time is fixed to 7 days in order to obtain a precipitation of the mineral up to around half of the radial zone. To avoid too small control volumes in the nearwell region, the parameter ω has been chosen in such a way that min κ∈M φ k = min s∈V φ s , leading in our case to ω ∼ 0.4.

We report in the table 1 below for both schemes and for both meshes the number of unknowns nu (#M for the MPFA O scheme and #V int ∪ V N for the VAG scheme), the number nnz of nonzeros blocks in the Jacobian after elimination of the cell unknowns for the VAG scheme, the average number of Newton iterations per time step nnl, and the average number of GMRES iterations per Newton step nl. The stopping criteriae in terms of relative residual for both the Newton convergence and the linear convergence are set to 10 -6 . Note that the simulation with the MPFA O scheme on the hybrid mesh could not be obtained due to too high memory requirement for the matrix storage. The slight increase of nonlinear iterations between the MPFA and the VAG scheme and between both meshes, as seen in table 1, is probably due to the decrease of the control volume sizes. Table 1 For both schemes and both meshes: number of unknowns nu, number nnz of nonzeros blocks in the Jacobian, average number nnl of Newton iterations per time step, and average number nl of GMRES iterations per Newton step.

Figure 17 exhibits, for the two schemes and on both grids, the rate of variation of the masses of CO 2 and of the mineral in the reservoir function of time. We can notice that the VAG scheme solutions on both meshes are almost the same and only slightly differ from the O scheme solution. The oscillations observed in Figure 17(b) are a well-known phenomenon due to the appearance of the mineral phase on each successive cells when the salt reaches its maximum solubility. where the mineral saturation higher than 0.1 percent is plotted at the end of the simulation. We can observe on Figure 19(c) the peak of mineral saturation around the well which is due to the reflux of the water phase close to the well by capillary effect. These results obtained with the VAG scheme on the hexahedral grid are similar to those obtained with the VAG scheme on the hybrid grid and to those obtained with the O scheme on the hexahedral grid.

Conclusion

A general formulation for the discretization of multiphase compositional Darcy flows is introduced allowing for an arbitrary number of components and phases and taking into account phase appearance and disappearance. The space discretization is a combination of finite element and finite volume. It applies to general meshes, which may be possibly non conforming and with non planar faces and leads to a vertex-centred compact scheme after elimination of the cell unknowns in the linear system arising from the Newton linearization. Compared with usual finite element approaches, the VAG scheme has the ability to deal with highly heterogeneous media and different rocktypes on coarse meshes due to its flexibility in the definition of the porous volumes at the vertices. The efficiency of our approach on complex meshes and for complex compositional models is exhibited on three phases three components models which simulate the nearwell injection of miscible CO 2 in a saline aquifer taking into account the vaporization of H 2 O in the gas phase as well a the deposition of the salt.

In order to better take into account discontinuous capillary pressures, a more advanced solution, following the recent work [START_REF] Brenner | A Convergent Finite Volume Scheme for Two-Phase Flows in Porous Media with Discontinuous Capillary Pressure Field[END_REF], would be to consider as vertex unknowns the capillary pressures rather than the satura-tions, in such a way that the saturations at the vertices would be allowed to be discontinuous. This approach will be developed in a future work.

Fig. 1

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

Fig. 2

 2 Fig. 2 CVFE interpretation of the fluxes F κ,s (u) in the 2D case.

Fig. 3

 3 Fig. 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 (7) 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[START_REF] Agélas | An abstract analysis framework for nonconforming approximations of diffusion on general meshes[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

Fig. 4

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

Fig. 5

 5 Fig. 5 Log normal permeability field of variance 5.

Figures 6 and 7

 7 Figures 6 and 7 exhibit the gas saturation fronts at fixed simulation times for both the TPFA and the VAG schemes. It can be checked that the TPFA and VAG solutions are very closed for both values of the variance. The VAG solution is slightly more developed which could be due to the richer flux stencil of the VAG saturation equation.

Fig. 6

 6 Fig. 6 Gas saturation, permeability field of variance 1.

Fig. 7

 7 Fig. 7 Gas saturation, permeability field of variance 5.

  (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.

Fig. 8

 8 Fig. 8 Configuration of the test case of injection of gas in an heterogeneous reservoir.

Fig. 9

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

Fig. 11

 11 Fig. 11 Cut at y = 25 m -Gas saturation at the end of the simulation.

Fig. 12

 12 Fig. 12 Propagation in the horizontal direction of the gas saturation.

P

  c,rt (S o ) = a rt S o + b rt with a 1 = a 2 = 10 +5 , b 1 = 0 and b 2 = 0.5 10 +5 . Homogeneous Neumann boundary conditions are imposed at the boundary of Ω and the flow is buoyancy driven starting from the initial oil saturation defined by

4. 6 Fig. 13

 613 Fig. 13 Oil saturation along the z axis after 15 days of simulation on the coarse and fine grids, and for the TPFA, VAGs1, VGAs2, and VAGa schemes.

H2O

  in presence of both phases w and g,C w CO2 = K CO2 C g CO2in presence of both phases w and g, C w salt = K salt in presence of both phases w and m.

(a) grid 2 × 10 (b) grid 16 × 80 Fig. 14

 108014 Fig. 14 Oil saturation along the z axis after 30 days of simulation on the coarse and fine grids, and for the TPFA, VAGs1, VGAs2, and VAGa schemes.

Fig. 15

 15 Fig. 15 Diagram of present phases in the space (Z CO 2 , Z salt ).

  Fig. 16 Nearwell meshes.

Fig. 17

 17 Fig. 17 Rate of variation of the mass of CO 2 and of mineral in the reservoir function of time.

  Fig. 17 Rate of variation of the mass of CO 2 and of mineral in the reservoir function of time.

Fig. 18

 18 Fig. 18 Trajectory of Z κ for four cells κ i , i = 1, • • • , 4 in the space (Z CO 2 , Z salt ). The four cells at all at z = -7 m and ordered according to their increasing distance to the well axis.

Figure 18

 18 Figure18exhibits the trajectory of the total mass fraction Z in the simplex (Z CO2 , Z salt ) function of time for four cells κ i , i = 1, • • • , 4. Starting from an initial state given by a single water phase and Z CO2 = 0, the mass of CO 2 increases due to the injection. The CO 2 is initially fully dissolved in the water phase, then the gas phase appears and its mass fraction θ g increases. As long as C w Salt is roughly constant, the trajectory Z is close to the line defined byZ CO2 = (1 -θ g )C w CO2 + θ g C g CO2 , Z salt = (1 -θ g )C w salt since C wCO2 and C g CO2 are both fixed by C w salt and the thermodynamical equilibrium constants. Once the water saturation S w is close to the irreducible water saturation S r,w , the composition C w salt increases rapidly due to the vaporization of H 2 O and the mineral phase appears. At the end of the simulation, it only remains the CO 2 component in the gas phase and the salt component in the mineral phase and the trajectory ends on the segment Z H2O = 0.

Figure 19 (

 19 Figure 19(a) shows the water saturation at the end of the simulation, illustrating the nearwell drying. The precipitation of the salt is exhibited in Figure 19(b)where the mineral saturation higher than 0.1 percent is plotted at the end of the simulation. We can observe on Figure19(c) the peak of mineral saturation around the well which is due to the reflux of the water phase close to the well by capillary effect. These results obtained with the VAG scheme on the hexahedral grid are similar to those obtained with the VAG scheme on the hybrid grid and to those obtained with the O scheme on the hexahedral grid.

  (a) reduced water saturation S w (b) mineral saturation such that S m >0.1% (c) Cut of the mineral saturation S m

Fig. 19

 19 Fig. 19 Saturations at the end of the simulation.

Acknowledgements This work was partially supported by GNR MoMAS.