Combined face based and nodal based discretizations on hybrid meshes for non-isothermal two-phase Darcy flow problems

In the last 20 years many discretization schemes have been developed to approximate the Darcy fluxes on polyhedral cells in heterogeneous anisotropic porous media. Among them, we can distinguished cell based approaches like the Two Point Flux Approximation (TPFA) or the Multi Point Flux Approximation (MPFA) schemes, face based approaches like the Hybrid Finite Volume (HFV) scheme belonging to the family of Hybrid Mimetic Mixed methods and nodal based discretizations like the Vertex Approximate Gradient (VAG) scheme. They all have their own drawbacks and advantages which typically depend on the type of cells and on the anisotropy of the medium. In this work, we propose a new methodology to combine the VAG and HFV discretizations on arbitrary subsets of cells or faces in order to choose the best suited scheme in different parts of the mesh. In our approach the TPFA discretization is considered as an HFV discretization for which the face unknowns can be eliminated. The coupling strategy is based on a node to face interpolation operator at the interfaces which must be chosen to ensure the consistency, the coercivity and the limit conformity properties of the combined discretization. ∗we would like to thank the BRGM and the Provence-Alpes-Côte d’Azur Region for co-funding the PhD of Laurence Beaude as well as the support of the CHARMS ANR project (ANR-16-CE06-0009) †Université Côte d’Azur, CNRS, Inria team Coffee, LJAD, France ‡Université Côte d’Azur, CNRS, Inria team Coffee, LJAD, France §BRGM, scientific and Technical Center, 3 avenue Claude Guillemin, BP 36009, 45060 Orléans Cedex 2 France ¶CSTJF, TOTAL S.A. Avenue Larribau, 64018 Pau, France


Introduction
The simulation of compositional multi-phase Darcy flow in heterogeneous porous media plays a major role in many applications. In the oil and gas sector, compositional multi-phase Darcy flow simulations are paramount to predict and optimize reservoir production. In sedimentary basin modelling, such models are used to simulate the migration of hydrocarbon phases, over geological space and time scales, from source rock to traps in geological formations. In CO 2 geological storage, compositional multi-phase Darcy flow models allow to optimize the injection of CO 2 and to evaluate the integrity of the storage. Two-phase compositional Darcy flow models are used to study the gas migration and to assess the long term safety of nuclear waste storages. Finally, coupling such models with the energy conservation equation lead to practical applications for both exploration and production phases of geothermal operations.
The standard industrial answer to cope with the strong coupling of both an elliptic (or parabolic) unknown, the pressure, and hyperbolic (or degenerate parabolic) unknowns, the volume and mole fractions, is based on finite volume spatial discretization, which is efficiently combined with an Euler implicit time integration to allow for sufficiently large time steps [4,22]. A major difficulty is linked to the finite volume discretization of the Darcy fluxes on the complex geometry and geology encountered in practical problems typically including fault networks, pinch-out, heterogeneities and spatially varying anisotropies of the medium. The classical cellcentred Two-Point Flux Approximation (TPFA) widely used in industrial simulators is cheap and robust but its consistency requires strong orthogonality conditions on the mesh which cannot be achieved for complex geological models. In the last 20 years, these restrictions have motivated an active research on the development of new discretization schemes to approximate the Darcy fluxes on polyhedral cells and in heterogeneous anisotropic porous media [13,18]. Still relying on the cell-centred approximation, Multi-Point Flux Approximations (MPFA) extend TPFA to consistent discretizations on general meshes with anisotropic heterogeneous media [1,14]. Yet, MPFA schemes are mesh and anisotropy conditionally stable and exhibit a very large stencil on simplectic meshes. Alternatively, nodal discretizations such as the Control Volume Finite Element (CVFE) method [19] and the Vertex Approximate Gradient (VAG) scheme [15,18,23] are unconditionally coercive and very efficient on simplectic meshes thanks to their nodal based approximation. Finally, face based discretizations such as the Hybrid Finite Volume (HFV) scheme [16] belonging to the family of Hybrid Mixed Mimetic (HMM) methods [12], or the Mixed Hybrid Finite Element method, have been developed and adapted to multi-phase Darcy flows in [2,3]. They provide accurate and unconditionally stable discretizations of the Darcy fluxes but, due to the large number of faces, remain rather expensive compared with nodal or cell-centred approaches.
Roughly speaking, all these discretizations of the Darcy fluxes have their own drawbacks and advantages which typically depend on the mesh characteristics and on the anisotropy of the medium. In this work, we propose a new methodology which combines nodal based and face based discretizations on arbitrary subsets of cells or faces in order to adapt the choice of the scheme in different parts of the mesh. In our approach, the TPFA is considered as a face based HFV scheme for which the face unknowns can be eliminated assuming that the mesh satisfies the superadmissibility property with anisotropy aligned with the mesh directions (see Lem. 2.1 of [16]). The VAG scheme is chosen as our nodal approximation since it shares a common data structure with the HFV scheme based on local to each cell transmissibility matrices. It has also the advantage, compared with more classical nodal discretizations such as CVFE, to avoid the mixing of rocktypes at nodal control volumes [18,23].
We consider two types of strategies to couple the VAG and HFV (TPFA) discretizations. The first one is based on a partition of the cells, each cell using either nodal or face unknowns. The second approach can be extended to more general partitions of the mesh based on faces, each face using either face or nodal unknowns. In both cases, the coupling is performed using a node to face interpolation operator at interfaces which must be chosen to ensure the consistency, the coercivity and the limit conformity properties of the combined VAG-HFV discretizations. The convergence analysis is performed in the gradient discretization framework [10,12,15] and convergence is proved for arbitrary cell or face partitions of the mesh. For face partitions, an additional stabilisation is required to ensure the coercivity while for cell partitions no additional stabilisation is needed and the stability is obtained at the interface thanks to the neighboring VAG cells.
At the interface, the framework preserves the discrete conservation properties of the VAG and HFV schemes with fluxes based on local to each cell transmissibility matrices which size is the number of selected nodes or/and faces on the shared boundary. This discrete conservative form leads to a natural extension of the VAG and HFV discretizations of multi-phase Darcy flow models to the combined VAG-HFV schemes.
The remainder of this paper is organized as follows. Section 2 focuses on the discretization of second order diffusion problems. It introduces our framework based, to fix ideas, on a partition of the cells into VAG, HFV and interface cells. Then, two gradient schemes are built combining the VAG and HFV schemes in their respective subset of cells coupled with two possible choices of matching discretizations in the interface cells. The stability and convergence of both discretizations are proved using the gradient discretization framework and the convergence is assessed numerically on various types of hybrid 3D meshes and compared with the standalone VAG and HFV discretizations. Using the discrete fluxes connecting each cell to its node or/and face boundary degrees of freedom, the VAG-HFV discretizations are extended to immiscible two-phase Darcy flows in Section 3. Then, numerical tests investigate, on a 1D two-phase flow reference solution, the convergence and efficiency of the VAG-HFV schemes compared with the standalone VAG and HFV discretizations. Section 4 considers the extension of the VAG-HFV discretizations to non-isothermal compositional liquid gas Darcy flows based on the formulation introduced in [5]. Finally, the model and its VAG-TPFA discretization are tested on a simplified two dimensional cross-section of the Bouillante high temperature geothermal reservoir with an hybrid cartesian triangular mesh. A reference solution, computed on a refined mesh, is compared in terms of accuracy and CPU time with the solutions obtained with the VAG scheme on a triangular mesh and the TPFA scheme on a Voronoi mesh.

Two Gradient discretizations combining the VAG and HFV schemes
Let us consider a domain Ω ⊂ R d , with d = 2, 3 the space dimension, and the following second order diffusion equation looking for the potential u ∈ H 1 0 (Ω) and the velocity q ∈ H div (Ω) such that div(q) = f on Ω, q = −Λ Λ Λ∇u on Ω. (2.1) In (2.1), f ∈ L 2 (Ω) is the source term and Λ Λ Λ ∈ L ∞ (Ω) d×d is the diffusion tensor such that there exist k ≥ k > 0 with The primal weak formulation of (2.1) amounts to find u ∈ H 1 0 (Ω) satisfying the following variational equality for all v ∈ H 1 0 (Ω) It admits a unique solution from the Lax Milgram theorem.

Polyhedral mesh and partition of the mesh
Following [15], we consider generalized polyhedral meshes of Ω. Let M be the set of cells that are disjoint open subsets of Ω such that K∈M K = Ω. For all K ∈ M, x K denotes the so-called "centre" of the cell K under the assumption that K is star-shaped with respect to x K . Let F denote the set of faces of the mesh. The faces are not assumed to be planar for the VAG discretization, hence the term "generalized polyhedral cells", but they need to be planar for the HFV discretization. We denote by V the set of vertices of the mesh. Let V K , F K , V σ respectively denote the set of the vertices of K ∈ M, faces of K and vertices of σ ∈ F. The set of edges of the mesh is denoted by E and E σ denotes the set of edges of the face σ ∈ F. Let M σ denote the set of cells sharing the face σ ∈ F. We denote by F ext the subset of faces σ ∈ F such that M σ has only one element and we set E ext = σ∈Fext E σ and V ext = σ∈Fext V σ . The mesh is assumed to be conforming in the sense that for all σ ∈ F \ F ext , the set M σ contains exactly two cells. It is assumed that, for each face σ ∈ F, there exists a so-called "centre" x σ of the face such that 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 σ .
A tetrahedral submesh of M is defined by where T K,σ,e is the tetrahedron joining the cell centre x K to the triangle T σ,e . Let ρ T denote the insphere diameter of a given tetrahedron T , h T its diameter and h T = max T ∈T h T . We will assume in the convergence analysis that the family of tetrahedral submeshes T is shape regular. Hence let us define the following shape regularity parameter of the mesh by The following combination of the VAG and HFV discretizations is based on the choice of a subset of cells M v ⊂ M on which the VAG discretization is used. Then, we define the subset of interfacial faces F hv ⊂ F \F ext by and the set of interface cells M hv by The subset of HFV cells on which the HFV discretization is used is finally defined by such that M v , M h , M hv defines a partition of the set of cells M (see Fig. 1). We also define the following subsets of nodes and faces it is assumed in the following that the face σ is planar and that x σ is the centre of gravity of σ.

Combining the VAG and HFV discretizations using the gradient discretization framework
The junction between the VAG and HFV discretizations is obtained using the gradient discretization framework introduced in [10,12,15]. This framework is based on the definition of a vector space of discrete unknowns X D , of a function reconstruction operator and of a gradient reconstruction operator The subspace of X D incorporating homogeneous Dirichlet boundary conditions is denoted by X 0 D . Then, the discretization of our model problem (2.2) is obtained by the following variational formulation: find u D ∈ X 0 D such that It admits a unique solution as soon as . D = ∇ D . (L 2 (Ω)) d defines a norm on X 0 D . As exhibited in Figure 1, our construction is based on the following set of degrees of freedom (d.o.f.) the associated vector space X D of discrete unknowns and its subspace Let us also define the subsets of d.o.f. located at the boundary of a given cell K ∈ M as The function reconstruction operator is based on an arbitrary partition {D K , D K,ν , ν ∈ dof K } of each cell K ∈ M and is defined by The gradient reconstruction operator is defined cellwise by It is based on the VAG gradient reconstruction operator ∇ v K for all VAG cells K ∈ M v and on the HFV gradient reconstruction operator ∇ h K for all HFV cells K ∈ M h . On the interface cells K ∈ M hv , the gradient reconstruction operator ∇ hv K must be built to guarantee that the gradient discretization (X D , ∇ D , Π D ) satisfies the coercivity, consistency and limit conformity properties of the gradient discretization framework, ensuring the well-posedness and convergence of the scheme (see [10,12,15] and Sect. 2.3 below).

VAG gradient reconstruction operator
Following [17], a P 1 finite element discretization is built using the tetrahedral submesh T of M and a second order interpolation at the face centres x σ , σ ∈ F \ F h defined for u D ∈ X D by For a given u D ∈ X D , we define the function Π T u D on K∈M v K as the continuous piecewise affine function on The VAG gradient reconstruction operator is obtained from this finite element discretization by setting (2.10)

HFV gradient reconstruction operator
We follow the construction presented in [16]. As shown in [11] it can be generalized as the family of Hybrid Mimetic Methods covering in the same framework Mimetic Finite Difference, Hybrid Finite Volume and Mixed Finite Volume Methods. For K ∈ M, let us set U K = u K , u σ , σ ∈ F K ∈ R #F K +1 and define where |K| is the volume of the cell K, |σ| is the surface of the face σ, and n K,σ is the unit normal vector of the face σ ∈ F K oriented outward of the cell K. Let us remark that ∇ K U K does not depend on u K since σ∈F K |σ|n K,σ = 0. Hence a stabilised gradient reconstruction is defined as follows: where K σ is the cone joining the face σ to the cell centre x K and |K σ | its d-dimensional measure. It leads to the definition of the HFV gradient reconstruction operator for u D ∈ X D as Note that the weight 1 √ d is chosen in order to recover the gradient reconstruction corresponding to the two point flux approximation in the case of a superadmissible mesh (see Lem. 2.1 of [16]).

First gradient reconstruction operator in the interface cells
Our first construction is based on a second order interpolation of the face unknown u σ at the centre of gravity x σ for each face σ ∈ F hv defined by Since x σ is the centre of gravity of the face σ, it results that Then we set for u D ∈ X D ,

Second gradient reconstruction operator in the interface cells
The second construction combines the previous interpolation at the faces σ ∈ F hv with a stabilisation of the cell gradient. As previously, for u D ∈ X D , let us set Then we rewrite the constant gradient ∇ K U K as This gradient does not actually depends on u K and must be stabilised using the residual and, for each ν ∈ dof K , the new gradient It leads to define the stabilised gradient where the weights (γ K,ν ) ν∈dof K and the partition (ω K,ν ) ν∈dof K of the cell K ∈ M hv are such that Note that, as soon as the diffusion tensor Λ Λ Λ(x) is cellwise constant, only the d-dimensional measures of the sets ω K,ν , ν ∈ dof K are used.
Remark 2.1. This second gradient reconstruction (2.13) based on interpolation and stabilisation can be applied as a standalone discretization in all cells provided that a partition of the faces σ ∈ F between those with a face unknown u σ and those with node unknowns u s , s ∈ V σ is given. The situation is different for the first gradient reconstruction (2.12) based only on interpolation which leads to a stable discretization thanks to the neighboring VAG cells (see Sect. 2.3). For example, if all faces are with node unknowns, it is clear that the first construction, if applied to all cells, will lead to an unstable discretization while the second construction reduces to the VAG discretization presented in [15] which differs from the VAG gradient reconstruction defined by (2.10).

Conservative formulation
From the cellwise definition of the gradient reconstruction, one can define the cell transmissibility symmetric positive matrices Let us define the following fluxes connecting each cell K ∈ M to its boundary d. (2.14) Then, the gradient scheme (2.4) can be formulated as a set of discrete conservation equations as follows: find (2.15) Each cell unknown u K can be eliminated from the first equation in (2.15) which depends only on u K and u ν , ν ∈ dof K . It leads to a Schur complement linear system without any fill-in depending only on the face and node unknowns u ν for ν ∈ F h ∪ V v .
Remark 2.2. For the first construction, given T h K ∈ R F K ×F K , the HFV transmissibility matrix of the cell K ∈ M hv , then the transmissibility matrix T K can be easily computed by for the second construction only. From (2.16), it is clear that T K is symmetric positive but not definite for the first construction for K ∈ M hv . Remark 2.4. In the special case for which a given cell K ∈ M hv satisfies the superadmissibility property σ ⊥ x K x σ for all σ ∈ F K , and say for Λ Λ Λ isotropic and cellwise constant, the HFV discretization transmissibility matrix T h K is diagonal leading to two point fluxes F K,σ (see Lem. 2.1 of [16]). From (2.16), it can be checked that this two point flux property is preserved by the first construction for all faces σ ∈ F K ∩ F h while it is not a priori the case for the second construction. This is one of the major advantage of the first approach when coupling the VAG and TPFA discretizations.

Gradient discretization framework
Let us recall the coercivity, consistency, and limit conformity properties for sequences of gradient discretizations introduced in [10,12,15].
Coercivity: Let C D > 0 be defined by Then, a sequence of gradient discretizations (D l ) l∈N is said to be coercive if there exist C P > 0 such that Consistency: For all u ∈ H 1 0 (Ω) and v D ∈ X 0 D let us define Then, a sequence of gradient discretizations (D l ) l∈N is said to be consistent if for all u ∈ H 1 0 (Ω) one has lim l→+∞ S D l (u) = 0.
Limit Conformity: For all q ∈ H div (Ω) and v D ∈ X 0 D , let us define and Then, a sequence of gradient discretizations (D l ) l∈N is said to be limit conforming if for all q ∈ H div (Ω) one has lim l→+∞ W D l (q) = 0.
be a gradient discretization such that . D is a norm on X 0 D , then the gradient scheme (2.4) has a unique solution u D ∈ X 0 D which satisfies the a priori estimate Let u ∈ H 1 0 (Ω) be the solution of (2.2) and let us set q = −Λ Λ Λ∇u ∈ H div (Ω). Then, one has the following error estimates:

2.3.2.
Proof of the coercivity, consistency and limit conformity properties for both constructions Proposition 2.2. Let us consider the gradient discretization D = X 0 D , ∇ D , Π D defined by (2.7), (2.8), (2.9) with the gradient reconstructions given either by (2.10), (2.11), (2.12) or by (2.10), (2.11), (2.13). Then, there exists C D depending only on θ T such that 22) and the following consistency estimate holds with C ϕ depending only on θ T and ϕ. Furthermore, the following limit conformity estimate holds with C ϕ depending only on θ T and ϕ.
Proof for the first gradient reconstruction: the consistency estimate (2.23) is a classical result already established in the case of the VAG discretization (see [7], Lem. 3.7 and 3.4) and of the HFV discretization (see [16], Lem. 4.3). The extension to our case results from the exactness of the cell gradients on affine functions as well as from the definition of θ T (2.3). Let us now prove the coercivity (2.22). Let us set for all u D ∈ X D , Π M u D (x) = u K for all x ∈ K and K ∈ M. It results from the discrete Sobolev embeddings proved in [16] Lemma 5.3, that there exists C 1 depending only on θ T such that for all u D ∈ X 0 with u σ = s∈Vσ β σ,s u s for all σ ∈ F \ F h and d K,σ = n K,σ · (x σ − x K ). For all K ∈ M v , it results from the convex combination assumption on the weights β σ,s , s ∈ V σ , the definition of θ T (2.3), and from Lemma 3.2 of [7] that there exists C 2 depending only on θ T such that For all K ∈ M h ∪ M hv , it also results from the definition of θ T (2.3), and from Lemma 4.1 of [16] that there exists C 3 depending only on θ T such that To conclude the proof of the coercivity property, let us now prove that there exists C 5 depending only on θ T such that for all It results from (2.27) and Lemma 3.4 of [7] that there exists a constant C 6 depending only on θ T such that for all u D ∈ X D and for all K ∈ M v ∪ M h , one has where h K is the diameter of the cell K. On the interface cells K ∈ M hv , from (2.27), there exists a constant C 7 depending only on θ T such that for all u D ∈ X D It is clear from (2.32) that the control of the contribution of the node s ∈ is obtained thanks to the neighboring VAG cell L.
Let us first prove the limit conformity estimate (2.24) for the gradient discretization X 0 and, from [16], it exists C depending only on θ T and ϕ such that Combining the previous identities and using that |σ|u We deduce that there exists C depending only on θ T and ϕ such that Since for σ ∈ F hv with M σ = {K, L}, K ∈ M v , there exists a constant C depending only on θ T such that we deduce that there exists a constant C depending only on θ T and ϕ such that whatever the subset M v of VAG cells one has the estimate Note that a better estimate of order h 3 2 T is obtained for |T hv | if the subset of VAG cells is such that the surface of the interface σ∈F hv |σ| remains bounded independently of the mesh.
This concludes the proof of the limit conformity estimate (2.24) for the gradient discretization X 0 D , ∇ D , Π D . The extension of this estimate to D = X 0 D , ∇ D , Π D results from the estimate for all u D ∈ X 0 D which is obtained in a similar way than the estimate (2.29) on Π M u D − Π D u D L 2 (Ω) established in the above proof of the coercivity.
Proof for the second gradient reconstruction: as stated in Remark 2.1, the second gradient reconstruction (2.13) can be used in combination with X 0 D and Π D as a standalone gradient discretization. The proof of the coercivity, consistency and limit conformity for this gradient discretization is similar to the one presented in Lemma 3.1 of [15] The limit conformity when combining this gradient discretization with the VAG gradient reconstruction (2.10) must be checked but this analysis is similar to the one performed above using that (2.9) with the gradient reconstructions given either by (2.10), (2.11), (2.12) or by (2.10), (2.11), (2.13) and such that there exists θ with θ T l ≤ θ for all l ∈ N and such that lim l→+∞ h T l = 0. Then, the sequence (D l ) l∈N is coercive, consistent and limit conforming. Therefore the gradient scheme is convergent. Furthermore, it satisfies a first order error estimate on smooth solutions.
Proof: the coercivity of the sequence of gradient discretizations results from Proposition 2.2 and from the shape regularity assumption. The consistency of the sequence of gradient discretizations results from Proposition 2.2, from lim l→+∞ h T l = 0 and from the density of C 2 (Ω) ∩ H 1 0 (Ω) in H 1 0 (Ω). The limit conformity property of the sequence of gradient discretizations results from Proposition 2.2, the density of (C 1 (Ω)) d in H div (Ω) and from the coercivity property.

Numerical tests for second order diffusion problems
In the following subsections, the VAG scheme on the full domain (vag), the HFV scheme on the full domain (hfv) and both combined VAG-HFV schemes using stabilisation (vag-hfv stab) or not (vag-hfv) are compared on various families of meshes. The objective of these test cases is to compare the accuracy of the four schemes and in particular of both combined VAG-HFV schemes. All test cases consider the exact solution u(x, y, z) = e cos(x+y+z) , with Dirichlet boundary conditions at ∂Ω. If not specified differently, the diffusion tensor Λ Λ Λ is the identity matrix.

Hexahedral meshes
Let us consider the family of uniform Cartesian grids of the domain Ω = (0, 1) 3 of size N × N × N with N = 8, 16, 32, 64. The family of hexahedral meshes is obtained by random perturbation of the Cartesian grids inside the subdomain Ω v = (0.25, 0.75) 3 as exhibited in Figure 2 for N = 8.
A convergence of order two is observed as expected in Figure 3 on the potential and of order one on the gradient for the VAG and both VAG-HFV schemes. This is not the case for the HFV scheme for which the gradient clearly does not converge due to the non planar faces in the VAG region. The combined VAG-HFV schemes solve this issue by using the consistent VAG scheme in the non planar face region. Both VAG-HFV schemes are remarkably more accurate than the VAG scheme on the potential.

Hybrid meshes with hexahedra and pyramids
Let us consider the family of uniform Cartesian grids of the domain Ω = (0, 1) 3 of size N × N × N with N = 8, 16, 32, 64. Then our family of hybrid meshes is obtained by cutting in 6 pyramids each cubic cell contained in the VAG subdomain Ω v = (0.25, 0.75) 3 as exhibited in Figure 4 for N = 8.
A convergence of order two is observed as expected in Figure 5 on the potential and of order one on the gradient for all schemes. The VAG scheme is more accurate than the HFV scheme on the gradient while it is the     contrary on the potential. The convergence of both combined VAG-HFV schemes matches with the convergence of the HFV scheme on the potential and the VAG-HFV schemes provide the best convergence of the gradient in L 2 norm.

Anisotropic test case
We consider in this subsection the previous test case with the following homogeneous anisotropic diffusion tensor As exhibited in Figure 6, the numerical results are similar as in the previous test case.   A super convergence of order two of the gradient is observed in Figure 8 for the VAG and HFV schemes on this family of uniform Cartesian meshes. This super convergence property is lost as expected for the combined VAG-HFV schemes. We also remark that the HFV scheme provides a better accuracy than the VAG scheme for this family of meshes and consequently that the convergence of the potential for the combined VAG-HFV schemes is close to the one of the VAG scheme.

Combined VAG-HFV discretization of two-phase Darcy flows
The extension of the VAG-HFV discretization to two-phase Darcy flows combines ideas presented in [18] for the VAG discretization of multi-phase Darcy flow models and in [3] for the HFV discretization of two-phase Darcy flows. It is based on the discrete fluxes F K,ν , K ∈ M, ν ∈ dof K defined in (2.14) and connecting each cell K to its nodes s ∈ V K ∩ V v and faces σ ∈ F K ∩ F h . Porous volumes are assigned as usual to all cells K ∈ M but also, following [18], to each node s ∈ V v \ V D (excluding the Dirichlet nodes V D ). Then, discrete conservation equations are derived for all K ∈ M and s ∈ V v \ V D using the porous volumes, the discrete fluxes F K,ν and an upwind approximation of the mobilities. The faces σ ∈ F h are considered as interfaces on which, following [3], the flux continuity equations are written for each phase assuming the continuity of the phase mobility. This is a natural generalisation to the HFV discretization of the harmonic transmissibility formula which is classicaly considered for Two-Point Flux Approximations [4,22].

Two-phase Darcy flow model
Let us consider the following two-phase Darcy flow model where φ(x) is the porous-medium porosity, P = {g, l} denotes the set of the non-wetting phase (denoted by g to fix ideas) and the wetting phase (denoted by l to fix ideas) and S α , α ∈ P is the phase saturation. The flow rates q α are defined by the following generalized Darcy laws for α ∈ P: where ρ α is the phase mass density, µ α is the phase dynamic viscosity, g is the gravitational acceleration, Λ Λ Λ(x) is the absolute permeability tensor, k α r (x, S α ) is the phase relative permeability and P α is the phase pressure. The model is closed by the following capillary pressure relation where P c (x, S g ) is the capillary function.
Let us consider a partition ∂Ω D ∪ ∂Ω N of the boundary ∂Ω. On ∂Ω D , we consider a Dirichlet boundary condition with prescribed pressures P α and saturations S α , α ∈ P. On ∂Ω N , homogeneous Neumann boundary conditions are imposed with q α · n = 0 for α ∈ P.

Combined VAG-HFV discretization
The mesh is assumed to be conforming with the boundary condition in the sense that there exists a subset F D ⊂ F ext such that σ∈F Dσ = ∂Ω D . Let us define the set of Dirichlet boundary HFV faces by and the set of Dirichlet boundary VAG nodes by The Neumann boundary HFV faces are defined by A rocktype rt K is assigned to each cell K ∈ M and each rocktype rt corresponds to given capillary and phase relative permeability functions denoted by P c,rt (S g ) and k α r,rt (S α ). The porosity and absolute permeability tensor are also assumed cellwise constant and denoted by φ K and Λ Λ Λ K for all K ∈ M.
Let us define the set of neighboring cells of a node s as M s = {K ∈ M | s ∈ V K }. For the VAG discretization, a single rocktype rt s is assigned to each node s ∈ V v \ V v D . It is chosen as the most permeable rocktype among all rocktypes (rt K ) K∈Ms∩M v . Then, the porous volume is distributed to the VAG nodes as follows: given a parameter ω ∈ (0, 1), we set for all Figure 9. Illustration on a 2D mesh of the distribution of the volumes α K,s |K| to each node s ∈ V K of the cell K ∈ M v in the case of a single rock type.
thus the porous volumes are for K ∈ M v , and ϕ K = φ K |K|, for K ∈ M h ∪ M hv . Note that ω is chosen small enough such that (1 − s∈V K α K,s ) > 0 (see Fig. 9). The complementary rock volume for ν ∈ M∪(V v \ V v D ) is denoted byφ ν . For HFV Dirichlet boundary faces σ ∈ F h D , let us also set rt σ = rt K with M σ = {K}.
The set of unknowns and dirichlet data, exhibited in Figure 10, is defined by and the subsets of phase pressure unknowns by for α ∈ P.
Using the combined VAG-HFV discretization, let us define the discrete Darcy fluxes for all K ∈ M and ν ∈ dof K by with G D ∈ X D such that for all ν ∈ dof D , G ν = x ν · g, where g is the gravity acceleration. The discrete generalized Darcy fluxes for all K ∈ M and ν ∈ dof K are deduced using an upwind approximation of the The time integration is based on a fully implicit Euler scheme to avoid severe restrictions on the time steps. For N t f ∈ N * , let us consider the time discretization t 0 = 0 < t 1 < · · · < t n−1 < t n · · · < t Nt f = t f of the time interval (0, t f ). We denote the time steps by ∆t n = t n − t n−1 for all n = 1, . . . , N t f . It leads to the following set of equations at each time step n = 1, · · · , N t f accounting for the cell conservation equations the VAG node conservation equations coupled with the flux continuity equations or Neumann boundary condition at HFV faces the Dirichlet boundary conditions at Dirichlet nodes and faces and with the closure laws S g,n ν + S l,n ν = 1, P g,n ν − P l,n ν = P c,rtν (S g,n ν ), ν ∈ M ∪ V v \ V v D . (3.9)

Numerical experiments on a one dimensional solution
This test case considers the domain Ω = (0, 1) 3 with homogeneous isotropic permeability Λ Λ Λ = 1 and porosity φ = 1. The gravitational acceleration g is set to zero, the relative permeabilities to k α r (S α ) = (S α ) 2 , α = g, l, the dynamic viscosities to µ g = 5 and µ l = 1, and the capillary pressure to P c (S g ) = −0.1 log(1 − S g ). Dirichlet boundary conditions are set at x = 0 with imposed non-wetting phase pressure P g = 2 and saturation S g = 0.9, as well as at x = 1 with P l = 0 and S g = 0. Homogeneous Neumann conditions are considered at the remaining boundaries. The saturation is set to S g = 0 at initial time t = 0 and the final simulation time is fixed to t f = 1.
In the following subsections, the VAG scheme on the full domain (vag), the HFV scheme on the full domain (hfv) and both combined VAG-HFV schemes using stabilisation (vag-hfv stab) or not (vag-hfv) are compared on two families of meshes. The error is computed both for the saturation and for the non-wetting phase pressure based on a numerical reference solution obtained on a one dimensional uniform grid with 1000 cells and time steps. A space time discrete L 1 norm computed from all cell and time step values is used for simplicity.
For all test cases, a uniform time stepping is used with 200 time steps on (0, t f ). The system of equations (3.5)-(3.8) is solved at each time step by using a Newton-Raphson algorithm w.r.t. the set of primary unknowns taking into account the elimination of the Dirichlet boundary conditions (3.8) and of the secondary unknowns together with the closure laws (3.9). The Jacobian matrix is computed analytically, then the VAG cell primary unknowns and conservation equations for K ∈ M v are eliminated without any fill-in by Schur complement. As mentioned in Remark 2.4, in the case of superadmissibility of the cells K ∈ M h ∪ M hv such as for Cartesian cells with isotropic permeability tensor, the HFV scheme reduces to a TPFA scheme and the face unknowns {(P g σ , P l σ ), σ ∈ F h \ F h D } are eliminated from the flux continuity equations (3.7). When using the stabilised version of the combined VAG-HFV scheme, note that the TPFA fluxes are not preserved at the faces of the interface cells. Consequently these face unknowns will not be eliminated for the stabilised VAG-HFV scheme. The reduced linear systems obtained at each Newton iteration are solved using the GMRes iterative solver combined with a CPR-AMG preconditioner [20,24].
The numerical behavior of the four schemes is reported for the two families of meshes on the finest mesh with N red the number of degree of freedom of the reduced linear systems with two primary unknowns per d.o.f., N Z red the number of non-zero 2 by 2 entries in the reduced linear systems, N newton the average number of Newton iterations per time step and N gmres the average number of GMRes iterations per Newton step. The CPU times are in seconds on 3.1 Ghz Intel Core i7 processor and 16Go RAM.

Hexahedral meshes
Let us consider the family of uniform Cartesian grids of the domain Ω = (0, 1) 3 of size N × N × N with N = 4, 8, 16, 32. The family of hexahedral meshes is obtained by random perturbation of the Cartesian grids inside the VAG subdomain Ω v = (0.25, 0.75) 3 as exhibited in Figure 2 for N = 8. Figure 11 exhibits, for the four schemes, the convergence of the error. All schemes exhibit the same order of convergence both for the saturation (lower but close to 1) and for the pressure (roughly 1). The VAG scheme is more accurate on this type of mesh since it uses more d.o.f. than the HFV scheme for the transport of the saturation accounting for the leading error term due to the first order upwind discretization of the mobilities. The combined VAG-HFV schemes both exhibit roughly the same convergence as the HFV scheme since the HFV domain is much larger than the VAG domain in this test case.   Table 1 exhibits the numerical behavior of the four schemes on the finest mesh N = 32 with 4096 hexahedra in the VAG subdomain and 28 672 cubes in the HFV subdomain. All schemes have roughly the same number of Newton iterations and hence their numerical behavior differs by the sparsity of the reduced linear systems and efficiency of the CPR-AMG preconditioner. The HFV scheme reduces to a TPFA scheme in the HFV region representing most of the domain which explains the better CPU time observed with the HFV scheme compared with the VAG scheme in this test case. The VAG-HFV scheme without stabilisation leads to the sparsest reduced system and to the lowest CPU time. The stabilisation increases the number of non-zero elements of the reduced linear systems since the face unknowns located at the interface cells are not eliminated. The CPR-AMG preconditioner is also less efficient leading to a CPU time that is twice larger than the unstabilised VAG-HFV scheme.

Hybrid meshes with hexahedra and pyramids
Let us consider the family of uniform Cartesian grids of the domain Ω = (0, 1) 3 of size N × N × N with N = 8, 16, 32. Then our family of hybrid meshes is obtained by cutting in 6 pyramids each cubic cell contained in the VAG subdomain Ω v = (0.5, 1) × (0, 1) 2 as exhibited in Figure 12 for N = 8.
As exhibited in Figure 13, the convergence of the error behaves like in the previous test case. The following Table 2 exhibits the numerical behavior of the four schemes on the finest mesh N = 32 with 16 384 cubes in the HFV subdomain and 98 304 pyramids in the VAG subdomain. The HFV scheme is roughly four times more costly than the three other schemes due the much larger number of faces in the VAG region compared with the number of nodes.

Combined VAG-HFV discretization of non-isothermal compositional two-phase Darcy flows
We consider in this section the extension of the previous combined VAG-HFV discretization to the case of a non-isothermal compositional two-phase Darcy flow. This extension is based on the formulation of the model introduced in [5]. Its main advantages compared with the related Coats' variable switch formulation [9] is to be based on a fixed set of unknowns using extended phase molar fractions and to express the thermodynamical equilibrium as complementary constraints for both phases α ∈ P. Previous works have considered the VAG discretization of isothermal and of non-isothermal compositional two-phase Darcy flows in respectively [18] and [26]. The HFV discretization of isothermal two component two-phase Darcy flows is also derived in [3] and the related Mixed-Hybrid Finite Element discretization of three-phase Darcy flows in [2]. Our extension to the combined VAG-HFV discretization follows the methodology presented in the previous section for immiscible isothermal two-phase flows which takes advantage of the cell based definition of the fluxes shared by the VAG, the HFV and by the modified scheme at interface cells.

Non-isothermal compositional two-phase Darcy flow model
Let us consider a non-isothermal compositional liquid gas Darcy flow model where each phase α ∈ P is a mixture of an arbitrary number of components with typically the water component (denoted by w) which can vaporize in the gas phase and the air component (denoted by a) which can dissolve in the liquid phase. The set of components is denoted by C. The thermodynamic properties of each phase α ∈ P depend on its pressure P α , the local equilibrium temperature of the system T and its molar fractions C α = (C α i ) i∈C . Our formulation of the model is based on the fixed set of unknowns defined by U = (P g , P l , T, S g , S l , C g , C l ).
Let us introduce the notations for the thermodynamic laws. The molar density is denoted by ζ α (P α , T, C α ) and the dynamic viscosity by µ α (P α , T, C α ) for each phase α ∈ P. The mass density is defined by ρ α (P α , T, C α ) = i∈C C α i m i ζ α where m i is the molar mass of the component i ∈ C. Let us denote by e α (P α , T, C α ) the molar internal energy, by h α (P α , T, C α ) the molar enthalpy of the phase α ∈ P, and by E r (T ) the rock energy per unit rock volume. For shorter notations, let us introduce the fluid energy per unit pore volume defined by and the number of moles of the component i ∈ C per unit pore volume denoted by In order to simplify the notations, each thermodynamic law can also be written in the following as a function of the full set of variables U still keeping the same notation for the function. Thermodynamic equilibrium between the gas and liquid phases is assumed for each component and governed by the phase fugacities denoted by f α (P α , T, C α ) = (f α i (P α , T, C α )) i∈C , α ∈ P.
Note that, as opposed to the Coats' variable switch formulation [8,9], the molar fractions C α of an absent phase α are extended by the ones at equilibrium with the present phase in the sense that the equality of the gas and liquid fugacities f g (P g , T, C g ) = f l (P l , T, C l ) always holds. This allows to fix the set of unknowns independently of the present phases.
The total molar flow rate q i of the component i ∈ C and the energy flow rate q e are obtained from the generalized Darcy velocities introduced in (3.2) such that where λ stands for the bulk thermal conductivity of the fluid and rock mixture. The molar diffusion is neglected for the sake of simplicity.
Let us write the system of equations accounting for the molar conservation of each component i ∈ C together with the energy conservation It is complemented by the capillary relation between the two phase pressures and the pore volume balance In the spirit of [21], the liquid gas thermodynamic equilibrium can be expressed as the following complementary constraints for each phase α ∈ P combined with the equality of the gas and liquid fugacities of each component To fix ideas, let us consider a Dirichlet boundary condition on ∂Ω D with prescribed phase pressures P α , molar fractions C α , saturation S α , α ∈ P and temperature T . On ∂Ω N , homogeneous Neumann boundary conditions are imposed with q α · n = 0 for α ∈ P and q e · n = 0.

Combined VAG-HFV discretization of non-isothermal compositional two-phase Darcy flows
The extension of the discrete set of unknowns (3.3) to the non-isothermal compositional two-phase flow model is defined by We also consider the subset of phase pressure unknowns P α D defined in (3.4), the subset of temperature unknowns and the subset of the physical unknowns at a given ν ∈ M ∪ V v ∪ F h D U ν = (P g ν , P l ν , T ν , S g ν , S l ν , C g ν , C l ν ).
Let us define the discrete Darcy fluxes for all K ∈ M and ν ∈ dof K by where the phase mass density is defined by Let us also introduce the discrete Fourier fluxes for all K ∈ M and ν ∈ dof K by where the thermal conductivity λ K is a cell average of the bulk thermal conductivity computed explicitly from the previous time step variables, and G K,ν is the flux function (2.14) obtained with the identity diffusion tensor in all cells. The discretization of the mobilities is obtained using the phase based upwind approximation. For each Darcy flux, let us define the phase dependent upwind control volume (Kν) α such that Let us introduce the upwind approximation of the phase molar fluxes The discrete generalised compositional Darcy and flowing enthalpy fluxes write respectively It leads to the following set of equations at each time step n = 1, . . . , N t f accounting for the cell conservation equations for K ∈ M, the VAG node conservation equations  Table 3. Primary unknowns of the degree of freedom ν ∈ M ∪ V v \ V v D depending on the active complementary constraints of the Newton-Raphson algorithm. for ν ∈ V v D ∪ F h D and with the closure laws f g i (P g,n ν , T n ν , C g,n ν ) = f l i (P l,n ν , T n ν , C l,n ν ), i ∈ C, 3) by a Newton-Raphson algorithm adapted to the complementary constraints (see [5] for details). The size of the linear system to be solved at each Newton iteration can be considerably reduced by elimination of the Dirichlet d.o.f. and by elimination for each d.o.f. ν ∈ M ∪ V v \ V v D of the local closure laws (4.8) w.r.t. to a set of #C + 4 secondary unknowns U S ν ⊂ U ν chosen in such a way that the differential of the closure laws w.r.t. U S ν is non singular. A classical choice of the set of primary unknowns U P ν = U ν \ U S ν is reported in Table 3. Furthermore, the VAG cell unknowns U K and equations (4.4), K ∈ M v can be eliminated from the linear system without any fill-in by Schur complement which considerably reduces the number of VAG d.o.f. in the case of simplectic meshes. Finally, in the case of superadmissibility of the cells K ∈ M h ∪ M hv , as in Remark 2.4, the HFV scheme reduces to a TPFA scheme and the face unknowns {(P g σ , P l σ , T σ ), σ ∈ F h \ F h D } are eliminated from the Darcy and Fourier flux continuity equations (3.7) leading to the classical harmonic transmissibilities.

Two dimensional Bouillante geothermal test case
Geothermal energy is a carbon-free non-intermittent energy source which can represent an alternative to fossil energy both for power production and direct use, heat representing a large amount of the world final energy consumption. In countries with a favorable geological context, high temperature geothermal energy can even make a significant contribution to power production. Consequently, the world installed capacity is expected to double over the present decade [6]. Compositional multi-phase thermal Darcy flow models are essential tools to provide a quantitative assessment of deep geothermal resources, from exploration phases of geothermal operations with assessment of the geothermal potential, validation of conceptual hypothesis, well siting... to production phases to achieve an optimal and sustainable exploitation and prevent resource exhaustion.
Hight temperature geothermal fields are often located in active geological settings (e.g. plate boundaries, volcanic areas...) with complex natural structures and geometries such as fault networks with discontinuous properties and fractures that act as drain or barriers on the deep transfer of mass and energy, thus controlling the distribution of geothermal resources. The geological modelling of such systems and their discretization into conforming unstructured meshes are challenging tasks and often result in meshes that are hardly tractable for flow simulations. In such situation, hybrid meshes composed of different types of cells best suited to discretize the geology and geometry in different parts of the geothermal system represent a clear asset. Then, the scheme is adapted locally to the type of mesh/cells using the methodology developed in the previous sections.
In this section, we consider a simplified geological setting which corresponds to a 2D vertical cross-section of the Bouillante high temperature geothermal field (Guadeloupe, French West Indies). The vertical cross-section is assumed to be in the plane of major fault zone acting as a regional permeable drain. Our objective is to compare the results of the simulations on different meshes and schemes of the non-isothermal compositional liquid gas Darcy flow model. The following test cases focus on the coupling between the VAG scheme and the TPFA scheme assuming that the cells K ∈ M h ∪ M hv satisfy the TPFA admissibility conditions. Only the unstabilised version of the combined VAG-TPFA scheme is considered in order to preserve the two-point fluxes at all faces σ ∈ F hv .
We consider the above non-isothermal compositional two-phase Darcy flow model with the set of water and air components C = {a, w}. The thermodynamic laws used in this test case are the following. The gas phase is assumed to have a perfect gas molar density ζ g = P g RT , R = 8.314 J K −1 mol −1 and a constant gas dynamic viscosity fixed to µ g = 2 × 10 −5 Pa s. The liquid molar enthalpy h l and the gas molar enthalpies of each component h g a , h g w are taken from [25], from which the gas molar enthalpy is deduced assuming a perfect mixture h g (P g , T, C g ) = i∈C C g i h g i (P g , T ).
The liquid molar density and viscosity are also fixed to the constant values ζ l = 1000 0.018 mol m −3 and µ l = 10 −3 Pa s in order to avoid thermal convection instabilities which would prevent the comparison of the schemes on the different types of meshes. The mass density is defined by ρ α = ζ α i∈C m i C α i with the molar masses m a = 0.029 and m w = 0.018 Kg mol −1 . The molar internal energy e α (T ) of each phase α ∈ P is considered to be equal to its enthalpy. The fugacities are defined by where the Henry constant of the air component is set to H a = 10 8 Pa and the vapour pressure P sat (T ) is given by the Clausius-Clapeyron equation In this test case, the porous medium is homogeneous with porosity φ(x) = 0.35 and isotropic permeability Λ Λ Λ(x) = K × I with K = 1 D. The relative permeabilities are defined by k α r (S α ) = (S α ) 2 for each phase α ∈ P and the capillary pressure function is given by the Corey law regularized for S g ∈ (S 1 , 1] to allow for the disappearance of the liquid phase, with b = 2×10 5 Pa and S 1 = 0.99. The rock energy per unit rock volume is fixed to E r (T ) = 2 × 10 6 T in J m −3 and the bulk thermal conductivity Figure 14. Illustration of the two dimensional geothermal reservoir and of the various conditions applied at its boundary.
of the fluid and rock mixture is fixed to λ = 3 W m 2 K −1 . Figure 14 shows the 2D vertical cross-section of the Bouillante geothermal reservoir and the conditions applied at the domain boundary. The initial and left side conditions are defined by a pure water liquid phase (S l = 1, C l w = 1, C l a = 0) at hydrostatic pressure and by a linear temperature between the fixed top and bottom temperatures. The bottom boundary is impervious (zero Darcy fluxes) with a fixed temperature of 400 K except in the interval 8000 m ≤ x ≤ 10 000 m where a pure water liquid input flux of −3 × 10 −2 mol m −2 s −1 at 550 K is imposed. The right side of the domain is supposed thermally isolated (zero Fourier flux) and impervious (zero Darcy fluxes).
The upper boundary is composed of three parts corresponding to the seabed (z ≤ 0 m and 0 ≤ x ≤ 5000 m), a sunny plain zone (0 < z ≤ 500 m and 5000 m < x ≤ 8450 m) and a rainy mountain zone (z > 500 m and 8450 m < x ≤ 11 000 m): -the seabed boundary condition is defined by a pure water liquid phase (S l = 1, C l w = 1) at hydrostatic pressure. The temperature is sea depth dependent. It is linear between the sea level z = 0 m at 300 K and z = −100 m at 278 K, then constant below z = −100 m, -the sunny plain zone is defined by the relative humidity fixed to 0.5, the temperature fixed to 300 K and the gas pressure fixed to P g = 1 atm, from which we deduce that only the gas phase is present with the molar fractions roughly equal to C g a 0.99, C g w 10 −2 , -the rainy mountain zone is characterized by a two-phase state at thermodynamic equilibrium which is defined by a fixed temperature, gas pressure and relative humidity corresponding to the following physical values:     Figures 16-19 show the temperature (on the left) and the gas saturation above the threshold of 10 −2 (on the right) at final time for the different combinations of meshes and schemes. All the numerical results are similar and no visible artefact can be seen at the interface z = −1500 m for the VAG-TPFA scheme. However, a closer comparison with the reference solution focusing on the bottom part of the domain in the temperature plots and on the gas saturation plots exhibits that the TPFA scheme (alone or combined with VAG) provides a more accurate solution. The solution is even more accurate when the TPFA scheme is associated with Cartesian cells at the bottom subdomain. It can be explained by the quasi horizontal flow lines close to the bottom boundary which is better captured by the Cartesian mesh coupled with the TPFA scheme (see Fig. 19      The vaporization occurs when the high temperature front reaches the low pressure zone close to the top boundary (at around t = 70 000 days). It confirms that the combined VAG-TPFA scheme applied on the Cartesian-triangular mesh better captures the high temperature front. Table 4 compares the numerical behavior of the simulation with N red the number of degree of freedom of the reduced linear systems with 3 primary unknowns per d.o.f., N Z red the number of non-zero 3 by 3 entries in the reduced linear systems, N t f the number of successful time steps, N chops the number of time step chops and N newton the average number of Newton-Raphson iterations per successful time step. The CPU times are in seconds on 2.9 GHz Intel Core i5 processor and 8Go RAM. This table exhibits large differences in CPU time between the different meshes and schemes whereas the numbers of cells are comparable for all meshes. Thanks to the elimination of the VAG cell unknowns for all K ∈ M v , the number of non-zero entries in the Jacobian matrix is the smallest with the VAG scheme combined with the triangular mesh. It results that the VAG scheme combined with the triangular mesh leads to the fastest solution. On the other hand, as explained above, it also leads to a less accurate solution compared with the solution obtained using the TPFA scheme combined with a Cartesian mesh on the lower part of the reservoir. All together, the Cartesian-triangular mesh combined with the VAG-TPFA scheme provides the best compromise in terms of CPU time and accuracy for this geothermal test case.

Conclusion
A new methodology is introduced in this work to combine face based (HFV or TPFA) and nodal based (VAG) discretizations on hybrid meshes in order to adapt the numerical scheme to the different types of cells and medium properties in different parts of the mesh. The stability and convergence of the combined VAG-HFV schemes is studied in the gradient discretization framework and is shown to hold on arbitrary partitions of the cells for the unstabilised version and on arbitrary partitions of the faces for the stabilised version. The framework preserves at the interface the discrete conservation properties of the VAG and HFV schemes with fluxes based on local to each cell transmissibility matrices of size the number of d.o.f. at the cell boundary. This discrete conservative form allows a natural extension of the VAG and HFV discretizations of two-phase Darcy flow models to the combined VAG-HFV schemes. Numerical results on different types of meshes show the accuracy and efficiency of the combined schemes which are compared to the standalone VAG and HFV (or TPFA) discretizations. The convergence of the schemes is first studied numerically for single-phase and two-phase Darcy flows using analytical solutions. Then, a non-isothermal compositional liquid gas Darcy flow test case representing a 2D vertical cross-section of the Bouillante geothermal reservoir is considered. For this test case, the combined VAG-TPFA scheme on an hybrid Cartesian-triangular mesh is shown to provide the best compromise between accuracy and CPU time compared with the VAG scheme on a triangular mesh and the TPFA scheme on a Voronoi mesh.