Regular Article
Energy stable numerical methods for porous media flow type problems
INRIA, Univ. Lille, CNRS, UMR 8524 – Laboratoire Paul Painlevé, 59000
Lille, France
^{*} Corresponding author: clement.cances@inria.fr
Received:
28
February
2018
Accepted:
19
September
2018
Many problems arising in the context of multiphase porous media flows that take the form of degenerate parabolic equations have a dissipative structure, so that the energy of an isolated system is decreasing along time. In this paper, we discuss two approaches to tune a rather large family of numerical method in order to ensure a control on the energy at the discrete level as well. The first methodology is based on upwinding of the mobilities and leads to schemes that are unconditionally positivity preserving but only first order accurate in space. We present a second methodology which is based on the construction of local positive dissipation tensors. This allows to recover a second order accuracy w.r.t. space, but the preservation of the positivity is conditioned to some additional assumption on the nonlinearities. Both methods are based on an underlying numerical method for a linear anisotropic diffusion equation. We do not suppose that this building block is monotone.
© C. Cancès, published by IFP Energies nouvelles, 2018
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
1.1 Twophase porous media flows
Incompressible twophase porous media flows are often modeled by the following set of equations. In the absence of source terms, the volume of the phase α ∈ {n, w} (n and w stand for nonwetting and wetting respectively) is locally conserved along time as a consequence of(1)where ϕ denotes the porosity of the rock (supposed to be constant w.r.t. time), where s _{ α } denotes the saturation of the phase α and v _{ α } denotes the filtration speed of the phase α. It is classically assumed that the phase filtration speeds obey the generalized Darcy law(2)
The intrinsic permeability tensor of the porous medium is a definite positive and symmetric tensor field whereas the relative permeability k _{r,α } of the phase α is an increasing function of the saturation satisfying ,k _{r,α } (0) = 0 (we neglect the residual saturations). The variations of the viscosities μ _{ α } and of the densities ρ _{ α } are neglected, and g = ∇(g· x ) denotes the gravity vector. The phase pressures p = (p _{n}, p _{w}) are the main unknowns of the system together with the saturations s = (s _{n}, s _{w}). Two algebraic relations are imposed to close the problem. The first one is a constraint coming from the fact that the whole pore volume is saturated by the fluid:(3)
The second one links the capillary pressure to the nonwetting phase saturation in a monotone way:(4)where π is a maximal monotone graph from [0,1] to that may also depend on the space variable in the case where the geological environment is made of several different rocks (see for instance [1–5]).
Once complemented by noflux conditions on the boundary of the porous domain Ω, the model (1)–(4) has a very particular variational structure. As depicted in [6] (see also [7–9]), this problem can be reinterpreted as the generalized gradient flow [10] of the energy(5)with
This energy is made of the capillary energy, of the gravitational potential energy, and of a contribution related to the constraint (3):
The capillary potential is a convex function that is finite in [0,1], infinite outside of [0,1], and satisfies(6)whereis the subdifferential of Π at s. The functional is convex and its subdifferential ∂E( s ) is made of the couples h = (h _{n}, h _{w}) = (p _{n} − ρ _{n} g⋅ x , p _{w} − ρ _{w} g⋅ x ) such that the relation (4) holds.
We will not go deep into details on the description of the gradient flow structure since it is not the purpose of this paper. We refer to [7, 11] for a more complete discussion on the Wasserstein gradient flow interpretation of porous media flows and to [7, 12] for a general presentation on gradient flows in metric spaces. Let us only stress important points that will motivate our discussion on numerical schemes. First, the gradient flow structure implies that ℰ( s ) is a decreasing function of time and that the dynamics aims at maximizing this decay. This yields an energy/dissipation of the form(7)for all t ≥ 0. The energy/dissipation relation (7) can be obtained by multiplying formally the equation (1) by h _{ α } and by summing over α ∈ {n,w}. As a consequence, stable steady states s ^{∞} are local minimizers of ℰ for which each phase is hydrostatic on its support, i.e.,(8)
In order to compute in a accurate way the longtime behavior of the system (this is very important in the context of basin modeling), the numerical scheme has to be designed to make the discrete counterpart of the energy ℰ( s ) decay along time and, as much as possible, to be exact on the equilibrium (8) (see [13]). In the hydrostatic zones, the equilibrium consists in a balance between the diffusion and the convection. Thus we will avoid operator splitting and discretize the convection and the diffusion simultaneously to recover this equilibrium.
1.2 A simplified model problem
Our purpose can already be illustrated on the single scalar equation(9)
Π being a convex and coercive internal energy functional as previously and Ψ being a smooth external (possibly gravitational) potential. The mobility function η is nondecreasing on and satisfies η(0) = 0 and η(s) > 0 if s > 0. Here, we used the notation for the domain of Π and we assume that 0 ∈ D(Π). If Π is defined on the whole , i.e., , then we assume moreover that Π is superlinear:
The solutions to (9) corresponding to nonnegative initial data remain nonnegative along time. This property might be destroyed by the numerical approximation, so that we need to extend the definitions of the functions Π, π, and η for negative values of s when this is possible. In the case where ∂Π(0) contains a finite value, i.e., if(10)then Π can be artificially extended on (−∞,0) ∪ D(Π) into a convex function (still denoted by Π) with a single valued subdifferential at , i.e., ∂Π(0) = π°(0) = π(0). This can be done for instance by prescribing(11)henceif s < 0. The function η can for instance be extended by setting η(s) = −s if s < 0.
The framework of (9) is already interesting since it includes Richards equation for which s is the water saturation and Ψ = −ρ g· x . The capillary energy function Π is then defined on as the antiderivative of the capillary pressure function. In usual models (see [14], pp. 343–345), D(Π) = [0,1] and condition (10) does not hold. This also includes a linear FokkerPlanck equation (but written in a nonlinear form) when η(s) = s and π(s) = log(s) or models for crowd motions [15]. We are interested in the computation of nonnegative solutions corresponding to initial data with(12)
There is still a generalized gradient flow structure corresponding to (9), the energy being given by(13)
The counterpart to (7) is obtained by multiplying formally (9) by p + Ψ and writes(14)
Similarly to what was discussed in Section 1.1, one aims to build numerical schemes for (9) that are accurate in the longtime limit. The steady states for (9) are now given by(15)
One of the difficulties arising in the resolution of the problem (9) is that the variable p and s are related by the maximal monotone graph π that may have vertical or horizontal parts. Similar difficulties occur for instance in the context of reactive flows in porous media, cf. [16]. For Richards equation, the graph π has vertical parts as depicted on Figure 1. The solution (s, p) to the problem (9) can be deduced from the knowledge of p, but not from s. But the choice of p as a primary unknown in a naive numerical method can yield severe difficulties. Typically mass conservation can be lost, and severe troubles can be encountered in the convergence of the iterative procedure to compute the solution to the nonlinear system arising from the numerical scheme. These difficulties motivated the development of several strategies to optimize the convergence properties of the iterative procedures. An popular approach consist in making used of robust fixed point procedures with linear convergence speed rather than with Newton’s method (see for instance [17–22]). There were also important efforts carried out to fix the difficulties of Newton’s method [23–25]. Comparisons between the fixed point and the Newton’s strategies are presented for instance in [26, 27] (see also [28]). In [29], the authors combine a Picardtype fixed point strategy with Newton’s method (i.e., they perform a few fixed points iterations before running Newton’s algorithm). An alternative approach would consist in keeping both s and p as unknowns together with the additional equation p ∈ π(s) that is often rephrased as a complementary constraint and then solving the problem with a nonsmooth Newton method (see for instance [30–32]). Another classical solution consists in partitioning Ω at each time t into a part Ω_{ s } (t) where s is chosen as a primary variable and a part Ω_{ p } (t) where p is chosen as a primary variable. Switching from one variable to another is then compulsory [33] in this case. This can be seen as a particular case of the approach proposed recently in [34, 35], which consists in parametrizing the graph π in a suitable way. More precisely, the graph can be described by two functions and such that(16) cf. Figure 1. We assume moreover that the parametrization is nondegenerate in the sense that is a strictly increasing function.
Fig. 1 The maximal monotone graph π depicted on top is parametrized as with some monotone functions (red) and (blue) depicted in the bottom picture. Here, and satisfy . 
Then the knowledge of the unphysical variable w allows to reconstruct both s and p, so that on can use w as a primary variable in the computations. There are an infinity of possible choices for the parametrization of π. One example is the resolvents of π and π ^{−1}, i.e., and . The idea in [34] is to take advantage of this flexibility and to choose the parametrization in order to optimize the convergence of the Newton’s method.
Before discretizing the nonlinear transient problems (1)–(4) or (9), we need to introduce some material concerning the discretization of elliptic equations.
1.3 A diffusion building block for numerical approximation
The goal of this section is to introduce a rather general framework that can fit to many numerical methods to solve the elliptic equation(17)subject to noflux boundary conditions. There is a huge literature about the numerical resolution of the above problem. Besides classical conforming and nonconforming finite elements (see for instance [36, 37]), let us mention some approaches based on mixed finite elements [38, 39], MultiPoint Flux Approximation (MPFA) finite volumes [40–43], or Discontinuous Galerkin method [44–46]. The list of aforementioned methods and related publications is of course far from being exhaustive, and some other numerical methods will be discussed in what follows.
Let be the set of geometrical entities to which correspond the unknowns (i.e., the cells for cell centered methods, the edges for hybrid methods, or the vertices for vertex centered methods), then we are interested in numerical methods for solving (17) that reduce into a linear system of the form(18)where are such that ∑_{ K } m _{ K } = Ω. The compatibility condition on the righthand side reduces to ∑_{ K } m _{ K } f _{ K } = 0. We denote by the vector unknown and by the matrix corresponding to the system (18), i.e., where b _{ K } = m _{ K } f _{ K }. We require the transmissivity coefficients a _{ KL } (and thus the matrix ) to be symmetric(19)
In what follows, denotes the set of the couples such that the transmissivity a _{ KL } is different from 0. Thanks to the symmetry property (19), the scheme (18) is equivalent to: ,(20)
Many numerical schemes can write as (18) and (19). This is for instance the case of the classical TwoPoint Flux Approximation (TPFA) finite volume scheme [47–49], where(21)as soon as the cells K and L share an edge σ _{ KL }. In (21), x _{ K } (resp. x _{ L }) denotes the center of the cell K, whereas the coefficient m _{ K } in (18) is the Lebesgue measure of the cell K. The TPFA scheme requires strong assumptions for its consistency. The permeability tensor must be isotropic – in (21) is the harmonic mean of the permeabilities and associated to the cells K and L –, and the edge σ _{ KL } must be orthogonal to the straight line [ x _{ K }, x _{ L }] joining the cell centers K and L.
Another classical example of scheme satisfying (18) is the classical conformal finite elements with mass lumping, which can be seen as a particular box scheme [50]. In this case, the unknowns are located at the vertices of a simplicial mesh of Ω. Denoting by ϕ _{ K } the Lagrange basis function associated to the vertex , the transmissivities a _{ KL } are defined by(22)and the lumped mass associated to the vertex K is Note that in this context, the transmissivities a _{ KL } may be negative (for instance in the case where does not satisfy the Delaunay condition if ).
The HybridMixedMimetic (HMM) methods [51] containing Hybrid Finite Volumes (HFV, see [52]), Mixed Finite Volumes [53], and Mimetic Finite Differences [54, 55] also enter the framework (18)–(19), as well as Discrete Duality Finite Volumes (DDFV) [56]. We refer to Droniou’s review paper [57] and to the book [58] for a more complete overview of the methods entering our framework.
In the case where the transmissivities are nonnegative(23)the scheme (18) is monotone and it fulfills the maximum principle. The property (23) is lost as soon as the mesh and the anisotropy tensor do not fulfill restrictive conditions. But for many methods, the transmissivities remain blockwise positive. This means that there exist geometrical entities and positive definite symmetric matrices , , such that the scheme (20) rewrites(24)
The vector δ ^{ M } u represents the inner variations of u inside . Let us illustrate formula (24) on some classical methods. First, for the TPFA scheme, since a _{ KL } > 0, then one can choose , as the set of the edges, , and where M is the diamond cell corresponding to the edge (K, L). The case conformal finite elements is more interesting. In this case, denotes the set of the simplices, whereas is equal to the space dimension d. The vertices of the simplex are denoted by , andwhereas the matrix is defined by(25)
A crucial property of the local diffusion matrices is that the condition number of can be bounded by a quantity depending only on the condition number of and on the regularity (in Ciarlet’s sense [36], see also [37]) of the simplex M, i.e.,(26)
The HFV (or SUSHI [52]) also enter the framework (24). In this case, is the set of the control volumes that can be quite general (nonconvex cells having various number of edges). To each cell M, there is one cell unknown and edge unknowns where denotes the number of edges of the element M, and the inner variation vector δ ^{ M } v is given by
The matrix is built thanks to a formula similar to (25) but with an approximate gradient that is piecewise constant on the half diamonds. Here again, the conditioning of the matrix only depends on the one of and on the regularity of the mesh. The Vertex Approximate Gradient (VAG) scheme [59, 60] has a very similar structure and also enters our framework (see [61]). The main difference with SUSHI is that the unknowns are located at the vertices of the cells instead of the edges. A last example concerns the DDFV method that, at least in 2D, also enters the framework (24) as highlighted in [62, 63]. In this context, is the set of the diamond cells and once again, the local matrices have bounded conditioning numbers if has a bounded conditioning number and if the mesh fulfills some weak regularity assumption.
Before coming back to our nonlinear parabolic problem, there is a last point that we need to point out here. The coercivity of the numerical method (20) amounts to claiming that the matrix is a symmetric definite positive matrix, the lowest eigenvalue of being bounded away from 0 uniformly w.r.t. the mesh size:where
Combining this information with (26), it was proven in [64] that(27)for some C depending only on the mesh regularity and on the anisotropy of the permeability tensor. A similar property was derived for the SUSHI scheme in [52]. The abstract framework of polytopal toolboxes of Droniou et al. [58] also allows to derive inequalities of this type.
1.4 About the paper content
The (physical) energy stability of numerical methods appears to be a secondary point for a large part of the mathematical literature concerning porous media flows with capillary diffusion. This originates probably from the fact that the mathematical analysis for such problems often relies on the use of the socalled Kirchhoff transform and global pressure (see for instance [65–68]) which may mask the considerations regarding the physical energy. The question of the convergence of numerical methods is more often addressed. For degenerate parabolic scalar equations, we refer for instance to [69–79], whereas for multiphase porous media flows, we refer to [80–85]. Here again, the list is of course far from being exhaustive.
In the remaining of this paper, we will focus on the construction of schemes for (1)–(4) or more simply for (9) that make the discrete counterpart of the energy decay with time. The goal is to design schemes that capture in an accurate way the longtime behavior of the continuous model, and in particular the equilibriums (8) and (15). Two strategies will be discussed. The first one is detailed in Section 2 and consists in working with the formulation (20) and to use an appropriate upwinding of the mobilities. The second one, which is presented in Section 3, rather uses the formulation (24) and aims at preserving a similar structure even in the presence of a degenerate mobility.
For the sake of simplicity, we do not discuss here about possible source terms or other boundary conditions. Adding these terms to the problem is possible but the energy of the system would no longer decrease in general. However, under reasonable assumptions, the problem remains energy stable. We refer to [86] for a detailed treatment of source terms and inhomogeneous boundary conditions in the context masslumped finite elements.
2 Upstream mobility schemes
Before addressing the twophase flow problem (1)–(4), let us first focus on the simpler scalar problem (9).
2.1 The scheme for the simplified problem (9)
The goal is to tune the numerical method (18) to approximate the solution to (9) while preserving a good energetic behavior at the discrete level. Concerning the time discretization, we stick to the backward Euler scheme. Let be parametrizations of the graph π as in (16), then the numerical scheme writes: (28)where we have used the notation for the neighbors of K. The time step τ _{ n } is not necessarily uniform and can depend on n ≥ 1. The initial saturation is assumed to be given and once has been computed for n ≥ 1, one sets
Following [64] and[87], the mobility is chosen thank to an upwinding that takes the sign of a _{ KL } into account:(29)
Note that , so that one can write the scheme (28) under the equivalent form(30)
In the above formula and as in (20), v is an arbitrary element of .
2.2 Some properties of the scheme
As a consequence of the choice (29) for the mobility , the method preserves the nonnegativity. More precisely, this means that if is a solution of the scheme (28), then so does its projection on , where is such that and (cf. Fig. 1). So without loss of generality, we can assume that(31)
This property still holds for transmissivities a _{ KL } violating the monotonicity condition (23). Note that the inequality (31) is of interest only if .
Let us now turn to the energy estimate. Choosing in the weak formulation (30) and using the convexity inequalityfor all and for all p ∈ ∂Π(s) = π(s), we obtain the following discrete counterpart of (14):(32)where(33)
If the monotonicity condition (23) holds (like for instance for TPFA finite volume schemes), the inequality (32) is enough to ensure that the second term is nonnegative, thus the discrete energy is diminishing along time, i.e.,
Moreover, summing over n provides a control on the numerical dissipation (34)
This estimate together with (31) is enough to prove the existence of one solution ( w ^{ n }) to the scheme for all n ≥ 1. Thanks to the monotonicity of the scheme and the nondegeneracy of the parametrization , the solution is unique (see [34]). Finally, the convergence of the scheme can be proved if the approximation of the diffusion operator is consistent. This requires some classical conditions on the mesh and on the anisotropy tensor for the TPFA scheme [49].
The case where the condition (23) is not satisfied is more intricate. Indeed, the estimate (32) still hold but the fact that 0 is a lower bound in (34) is no longer true in general. deducing a control on the energy and on its dissipation from (32) is not straightforward. Indeed, the second term does not have an obvious sign and we are not able to claim that the dissipation is positive along time. However, it is possible to bound from below the dissipation rate. To do so, let us first remark that the definition (29) of the upstream mobilities implies that(35)whatever the nondecreasing function . In the above relations, denotes the interval with extremal points and and we have used the monotonicity of η. In the case where is bounded (this is natural in the porous media flow context since is), we can take advantage of this formulation with , of the inequality (27), on the Lipschitz regularity of the external potential Ψ to show (cf. [88], Sect. 3.2) that(36)
A similar inequality can be obtained even in the case where is not bounded under suitable assumptions on the nonlinearities, see for instance ([61], Sect. 3.1). Inequality (36) allows to show that the discrete energy grows at most linearly with time, as well as the existence of (at least) one solution to the scheme thanks to a topological degree argument. Moreover, it allows one to show the convergence of the scheme towards the unique solution to the continuous problem if the mobility function η is bounded and if π ^{−1} is a function (i.e., when there are no hyperbolic degeneracy in the problem (9)) when the mesh size and the time step τ _{ n } tend to 0. We refer to [88] for the details concerning the convergence analysis of the scheme (28).
The upstream mobility scheme can be tested by other quantities of the form thanks to (35). For s ∈ [0, 1], we denote by the largest interval of [−∞,+∞] such that for all . Then given a nondecreasing onto function , there corresponds a convex and coercive function such that . In particular, there holds(37)for all .
We can then take benefits of the characterization (35) of the upstream mobility to get enhanced regularity estimates. For instance, choosing(38)one gets that(39)provided Ψ is regular enough and the initial entropy is finite:
Details for estimate (39) are provided in Appendix A.1. Such an estimate is the cornerstone in the study [89] where the convergence of an upstream mobility scheme was established for the Dupuit approximation of multiphase porous media flows. It was also used in [8], or in [90] where a degenerate CahnHilliard system with a very similar mathematical structure to (1)–(4) (see [91, 92]) was considered.
2.3 Convergence of the scheme
Let us illustrate the convergence of the scheme (28) when the mesh size and the time step tend to 0. As a model problem, we consider the very simple convection diffusion equation(40)where Ω = (0, 1)^{2}, x = (x _{1}, x _{2})^{ t }, and t _{ f } = 0.05. The permeability tensor is assumed to be diagonal of the form
The equation boils down into the very simple linear equation(41)where e _{2} = (0, 1)^{ t }. We choose the initial condition and the noflux boundary condition in accordance with the exact solution(42)with α = π ^{2} + 1/4. Note that s ^{ini} ( x ) = s _{ex } ( x ,0) vanishes when x _{2} = 1.
The solution is computed thanks to a Control Volume Finite Element scheme [87, 88]. This scheme requires conformal triangulations of Ω. We use successively refined Delaunay grids from the FVCA5 benchmark [93]. In the isotropic case (corresponding to κ = 1), the scheme is monotone, i.e., condition (23) is fulfilled. This is no longer true in the anisotropic case κ = 20.
We also compute the solution to our scheme with an anisotropy ratio κ = 20. The numerical solutions are compared with those computed with the following linear scheme with centered convection: (43)
Here, Ψ_{ K } = x _{K}· e _{2}. The scheme (43) is second order accurate in space. In the isotropic case and for fine enough grids (leading to small enough Peclet numbers), the scheme turns out to be monotone. This is no longer the case in the anisotropic case κ = 20 even though the scheme remains second order accurate in space.
In order to make sure that the error corresponding to the time discretization remains small when compared to the error corresponding to the space discretization, we divide the time step by 4 when the mesh size is divided by 2. We present on Figure 2 the L _{2} (Ω × (0, t _{f})) error corresponding to the numerical solutions produced by the schemes (28) and (43) for Delaunay triangular grids from [93].
Fig. 2 L ^{2} (Ω × (0, t _{f})) error as a function of the mesh size for the scheme (28) in the isotropic case κ = 1 (solid red) and anisotropic case κ = 20 (solid blue) corresponding to upstream mobilities (29). Comparison with the solution to the linear scheme (43) in the isotropic case (dashed red) and anisotropic case (dashed blue). 
The method (28) preserves positivity whatever the anisotropy ratio. This is not the case of the linear scheme (43) that produces undershoots in the anisotropic case κ = 20. As expected because of the upwinding procedure, the scheme (28) is first order accurate in space, i.e.,(44)where u _{ h } denotes the piecewise constant reconstruction on the dual mesh and where h denotes the mesh size. But it appears on Figure 2 that the constant C appearing in (44) strongly depends on the anisotropy ratio.
To sum up, the method (28) enjoys a very strong stability when used in the nonmonotone context, but this is mainly due to the excessive numerical diffusion. Even though the method is still first order accurate w.r.t. space, the constant strongly depends on the anisotropy ratio. This makes the method inefficient in the case of large anisotropy ratios or of poor quality meshes. However, the upstream mobility finite volume scheme remains a very robust method for solving complex but isotropic problems, like for instance some degenerate CahnHilliard systems [90] or multiphase porous media flows [8, 89, 94].
2.4 Longtime behavior of the scheme
It is natural to wonder if the numerical scheme is able to reproduce in an accurate way the longtime behavior (15) of the continuous problem. This question is indeed of broad interest for porous media flows, in particular in the context where the time scales are long like for instance in basin modeling or for nuclear waste repository management. Therefore, the study of the longtime behavior of Finite Volume schemes for convectiondiffusion problems has been the purpose of several contributions (see, e.g., [13, 62, 95–97]).
We go back to the twophase flow model presented in Section 1.1 that we discretize on a Delaunay mesh thanks to a classical fully implicit upstream mobility finite volume scheme (see for instance [94] or [8]). In particular, the monotonicity relation (23) is fulfilled and the discrete saturations remain bounded between 0 and 1.
At the continuous level, the relation (8) prescribing the longtime limit boils down to the following alternative:(45)for some . The parameter γ is determined by the conservation of mass
The steady state (45) can be discretized directly into(46)with γ' fixed so that(47)
In the above relation, is the discrete porosity. In particular, in the classical case where is a (singlevalued) function, the following alternative holds:for all
Let be the discrete counterpart of the energy (5), i.e.,(48)then the energy is decreasing along time, i.e.,
One can show (cf. Appendix A.2) that is a minimizer of under the constraint (47). Therefore the relative energy is nonnegative. If the capillary pressure π is an increasing function on (0, 1), then the relative energy vanishes if and only if . This quantity can be used to illustrate the convergence of towards as n tends to ∞. In the case depicted on Figure 3, we set Ω = (−1/2,1/2)^{2}, , k _{r,α }(s) = s, μ _{n} = 10, μ _{w} = 1, ρ _{n} = 0.87, ρ _{w} = 1, g = (0, −9.81)^{ T }, ϕ ≡ 1, , and
Fig. 3 Evolution of the relative energy as a function of time. We observe the exponential convergence towards until the machine precision is reached. 
We see on Figure 3 that the relative energy converges exponentially fast towards 0 until it reaches the machine precision. This shows both that the energy is effectively dissipated along time and that given by (45) and (46) is a steady solution to the scheme.
3 Schemes with local positive dissipation tensors
The upstream mobility numerical schemes presented in the previous section enjoy very nice properties but they are merely first order accurate in space and lack robustness w.r.t. the anisotropy ratio (or to the mesh regularity) as illustrated in Figure 2. This motivates the development of alternative numerical methods with the following specifications.

No Kirchhoff transform: the scheme should only involve quantities with a clear physical meaning.

Nonlinear stability: the scheme must fulfill a discrete counterpart of the energy/energy dissipation relation (14).

Convergence: when the mesh size tends to 0 and under mild regularity assumptions, the approximate solution produced by the numerical scheme must converge towards a solution to the PDE (9).

Second order accuracy: if the solution to (9) is smooth enough (say C ^{2}), the error between the exact and the approximate solution must behave as h ^{2} where h stands for the mesh size.

Robustness: the scheme should allow general grids and general anisotropy tensors. The accuracy should not be excessively impacted by the anisotropy ratio or the mesh regularity.
3.1 Presentation of the methodology
Assume that the building block (18) is second order accurate, then a natural scheme to fulfill the specifications (i) and (iv) of the previous list for problem (9) is(49)with a centered choice for the mobility:
Multiplying the scheme by and summing over leads once again to(50)
Already in the case developed in Section 2 where was chosen thanks to upwinding, the sign of the second term of the lefthand side was unclear. It was however possible to get a sufficient control to claim that the energy was growing at most linearly, cf. (36). This conclusion does not hold any longer in general for a centered choice of the mobilities and no control on the energy can be deduced from (50). Hence the specification (ii) is not satisfied by scheme (49).
In order to correct this, we propose a scheme based on the formalism (24), that is(51)
The above relation must hold for any . The matrix is called the local dissipation tensor. It must incorporate the local diffusion but also the mobilities . In order to ensure the dissipation property for the scheme, we want to be symmetric semidefinite positive. In [61], we proposed to choose(52)with(53)for the particular choice(54)
The fact that the energy is diminishing along time is obtained by choosing and by applying a simple convexity inequality, leading to(55)the second term being nonnegative since is semidefinite positive.
Additionally, the method is globally mass conservative, i.e.,(56)where . This estimate is obtained by choosing v = 1 in (50). Let us stress that in the particular cases of the schemes studied in [61, 62, 86], the scheme is also locally conservative since fluxes can be constructed.
In general, the quantity appearing in (52) is chosen as a convex combination of . In order to assess that the scheme (50)–(52) converges, the numerical mobilities have to satisfy an additional coercivity condition, that is(57)for some uniform α > 0. The condition (56) is clearly satisfied by the choice (53) with α = 1/2, but it prohibits the choice that would have been quite natural in the context of the SUSHI [52] or VAG [59] schemes.
The implementation of the scheme (50)–(53) can appear to be too involved. An easy way to simplify it is to choose independent on i, i.e.,
The matrix then reduces to and commutes with , leading to the simpler formula . In view of the constraint (56), a natural choice for η ^{ M } ( v ) is(58)
This choice was successfully used in [86] for a method based on conformal finite elements with mass lumping (here ), and in [62, 63] for a nonlinear DDFV method (here d = 2 and ).
3.2 Conditional positivity preservation and convergence w.r.t. the grid
The scheme (50) amounts at each time step to a system of nonlinear equations of the form(59)
The functions η, and are uniformly continuous on , then is also uniformly continuous. In order to ensure the existence of a solution w ^{ n } to the system, we need some bounds on w ^{ n } (that might depend on the mesh and the time step).
In the case where (10) holds, the estimate (55) (in particular the dissipation term) provides a sufficient bound on w ^{ n } in order to apply a topological degree argument and to claim that there exists (at least) one solution to the scheme (51). In the more intricate situation where (10) is no longer satisfied, corresponding to the situation where(60)(with a slight abuse of notation), then one can show that(61)for some ζ depending on the time step τ _{ n } and on the mesh. This is done for instance in [61] of Lemma 3.7 in the context of the VAG scheme, or in [63] of Lemma 3.5 for a DDFV scheme. In both case, the proof strongly relies on the coercivity assumption (57).
As a consequence of (61), the scheme (51) preserves the positivity as soon as (60) holds. This property is unfortunately lost in general when (10) is satisfied. To illustrate this fact, we show results of [61] where the solution of the porous medium equation(62)is rewritten under the form(63)for three different choices of nonlinearities, that are

η(s) = 1 and π(s) = ss (recall that we need to extend the nonlinearities for negative saturations if (10) holds);

η(s) = 2s and π(s) = s;

η(s) = 2s ^{2} and π(s) = log(s).
The case (a) does not enter our framework since η(0) = 1 ≠ 0, but it is interesting since it corresponds to the most natural approach to solve (62). The condition (10) holds in cases (a) and (b), whereas (60) holds in case (c). Therefore, the positivity of the solutions should be guaranteed only in this last case. To illustrate this fact, let us choose Ω = {(x, y) ∈ (0, 1)^{2}}, and let us approximate the exact solution to (62) defined by(64)thanks to the VAG scheme [61] on successively refined triangular meshes from the FVCA5 benchmark [93]. The problem is here complemented with Dirichlet boundary conditions.
We observe in Tables 1–3 that second order convergence is destroyed for all the three schemes because of the lack of regularity of the exact solution. As expected, the discrete solution corresponding to the choice (c) remains positive because condition (60) is verified. This is no longer the case for the choices (a) and (b) and undershoots are observed. But the choice (b) appears to be both cheaper and more accurate than the choice (a), whereas the amplitude of the undershoots is reduced.
Let us illustrate again the ability of the approach. We consider the linear and isotropic convection diffusion equation
that we rewrite under the nonlinear form(65)
We aim to approximate the exact solution (42) thanks to the nonlinear DDFV method proposed in [62, 63]. We compute the approximate solution corresponding to the sequence of Kershaw meshes from [93], see Figure 4. The numerical results are presented in Table 4.
Fig. 4 The Kershaw meshes are highly deformed topologically cartesian grids. The coarsest mesh of the family is depicted here. 
Approximation of (42) with a nonlinear DDFV scheme [62, 63] on the Kershaw mesh family from [93], final time T = 0.25. M is the mesh index, τ is the time step, errgs and errs respectively stand for the L ^{2}(Ω × (0, T)) error on ∇ s and the L ^{∞}((0, T); L ^{2}(Ø)) error on s, whereas ordgs and ords are the corresponding convergence orders. N _{max} and N _{mean} are the maximal and mean numbers of Newton iterations, and s _{min} minimal value of the approximation of s during the whole simulation.
As expected, since π(s) = log(s) fulfills (60), the solution remains positive although the grid is very irregular. The convergence at order 2 on u is not affected by the poor mesh regularity. We observe a superconvergence for the gradient.
Finally, let us show that the methodology presented in this section is often much more robust w.r.t. anisotropy than the methodology presented in Section 2. To this end, we stick to the testcase of the linear FokkerPlanck equation written in a nonlinear form described in Section 2.3. We discretize it with the energy stable finite element scheme with mass lumping proposed in [86]. In this scheme, the mobilities η ^{ M } ( v ) are chosen according to formula (58). The final time is set to t _{f} = 0.25. The L ^{2} ((0, t _{f}) × Ω) error as a function of the mesh size is plotted on Figure 5. As expected, the method is of order 2 whatever the anisotropy ratio. But it is worth noticing that the accuracy is almost not affected by the anisotropy ratio .
Fig. 5 Convergence towards the analytical solution (42) of the FokkerPlanck equation (41) for different values of the anisotropy ratio . 
3.3 About the longtime behavior of the scheme
We aim now to illustrate the longtime behavior of the scheme. One discretizes the equation (65) complemented with noflux boundary conditions following the methodology of [61]. In particular, we take , hence we denote by s ^{ n } = w ^{ n } the vector of the unknown saturations. The scheme (51) then rewrites(66)for any . The mobilities are discretized thanks to formula (54). The mass is conserved along time, i.e.,(67)
For any ρ > 0, the vector of defined by(68)is a steady solution to the scheme (66) since for all . Then the expected longtime limit as n → ∞ is s ^{∞} where ρ has been tuned so that
To illustrate this fact, we plot on Figure 6 the evolution of the relative energy as a function of the discrete time nτ for the sequence of successively refined triangular meshes used in Tables 1–3 and for the Kershaw mesh depicted on Figure 4. The relative energy is proved to be decreasing along time and vanishes if and only if s ^{ n } = s ^{∞}.
Fig. 6 Plot of the log of the relative energy (in log scale) as a function of time. 
We show in Figure 6 that the relative energy converges exponentially fast towards 0, providing a discrete counterpart to the relationthat is deduced from (42). We observe on Figure 6 that the scheme preserves exactly (up to machine precision) the longtime behavior of the equation, i.e., the longtime limit s ^{∞} is a discretization through (67) of the exact longtime behavior of (64). Figure 6 suggests that the convergence speed is sensitive to mesh regularity (the convergence is slightly too fast on the Kershaw mesh) and to mesh size (the convergence is too slow on coarse triangular meshes). Note that the relative energy is defined only for nonnegative s ^{ n }. For a behavior like the one depicted on Figure 6, the scheme must (at least) preserve positivity and admit s ^{∞} defined by (67) as a steady state.
3.4 Application to twophase porous media flows
We show now preliminary results obtained thanks to the methodology of this section to simulate twophase porous media flows. We consider an anisotropic quarter five spot problem with injection in a topleft corner and extraction in the bottomright corner. We neglect gravity, so that the equations boil down into(69) (70)
In the above system, we eliminated the unknowns p _{n} and s _{w}. They can be deduced from p _{w} and s _{n} by the relations s _{w} = 1−s _{n} and p _{n} = p _{w} + π(s _{n}). We fix the nonlinearities η _{n}, η _{w}, π, f _{n}, and f _{w} as in [98]:and
The pure wetting phase is injected in ω _{inj} = (0, 0.05) × (0.95, 1) near the topleft corner, while the fluid occupying the area ω _{prod} = (0.98, 1) × (0, 0.02) near the bottomright corner is extracted (see Fig. 7). Defining the fractional flow functions bythen the source terms q _{n} and q _{w} are defined by
Fig. 7 Schematic representation of the injection and production wells. 
The permeability tensor is anisotropic and we take a constant saturation profile ϕ ≡ 1.
The system (69)–(70) is complemented with noflux boundary conditions and with the initial condition s _{n} (·,0) ≡ 0.9.
For the discretization, we use conforming finite elements with mass lumping on a structured triangulation made of 5000 triangles and a constant time step τ = 10^{−3}. For the phase mobilities, we make a choice of type (58). As explained in [86], the method is locally conservative: one can build equilibrated fluxes in the RaviartThomasNedelec space such thatwhere denotes the approximation of the saturation s _{ α } at time .
We plot snapshots of the saturation s_{n} on Figure 8.
Fig. 8 Approximate saturation s _{n,h } at time t = 0.2 (top left), t = 0.5 (top right), t = 1.5 (bottom left) and t = 2.5 (bottom right). 
Acknowledgments
The author warmly thanks the anonymous referees for their precious feedback. He overall wants to thank his numerous collaborators on this topics, namely Ahmed Ait Hammou Oulhaj, Konstantin Brenner, Claire ChainaisHillairet, Thomas O. Gallouët, Cindy Guichard, Stella Krell, Maxime Laborde, Léonard Monsaingeon, Flore Nabet, and Martin Vohralík. They actively contributed to the realisation of this work. This research was supported by the French National Research Agency through project GEOPOR (ANR13JS01000701) and Labex CEMPI (ANR11LABX000701).
References
 van Duijn C.J., Molenaar J., de Neef M.J. (1995) The effect of capillary forces on immiscible twophase flows in heterogeneous porous media, Transp. Porous Media 21, 71–93. [Google Scholar]
 Bertsch M., Dal Passo R., van Duijn C.J. (2003) Analysis of oil trapping in porous media flow, SIAM J. Math. Anal. 35, 1, 245–267. ISSN 00361410. [CrossRef] [MathSciNet] [Google Scholar]
 Buzzi F., Lenzinger M., Schweizer B. (2009) Interface conditions for degenerate twophase flow equations in one space dimension, Analysis 29, 299–316. [Google Scholar]
 Cancès C., Gallouët T., Porretta A. (2009) Twophase flows involving capillary barriers in heterogeneous porous media, Interfaces Free Bound. 11, 2, 239–258. [CrossRef] [Google Scholar]
 Cancès C., Pierre M. (2012) An existence result for multidimensional immiscible twophase flows with discontinuous capillary pressure field, SIAM J. Math. Anal. 44, 2, 966–992. doi: 10.1137/11082943X. URL http://hal.archivesouvertes.fr/hal00518219. [CrossRef] [MathSciNet] [Google Scholar]
 Cancès C., Gallouët T.O., Monsaingeon L. (2015) The gradient flow structure of immiscible incompressible twophase flows in porous media, C. R. Acad. Sci. Paris Ser. I Math. 353, 985–989. [CrossRef] [Google Scholar]
 Cancès C., Gallouët T.O., Monsaingeon L. (2017) Incompressible immiscible multiphase flows in porous media: a variational approach, Anal. PDE 10, 8, 1845–1876. [CrossRef] [Google Scholar]
 Cancès C., Gallouët T.O., Laborde M., MonsainGeon L. (2018) Simulation of multiphase porous media flows with minimizing movement and finite volume schemes, HAL, hal01700952. URL https://hal.archivesouvertes.fr/hal01700952/document . [Google Scholar]
 Murphy T.J., Walkington N.J. Control volume approximation of degenerate twophase porous media flows, submitted for publication. [Google Scholar]
 Mielke A. (2011) A gradient structure for reactiondiffusion systems and for energydriftdiffusion systems, Nonlinearity 24, 4, 1329–1346, ISSN 09517715. doi: 10.1088/09517715/24/4/016. URL http://dx.doi.org/10.1088/09517715/24/4/016. [Google Scholar]
 Otto F. (2001) The geometry of dissipative evolution equations: the porous medium equation, Comm. PDE 26, 1–2, 101–174, ISSN 03605302. [Google Scholar]
 Ambrosio L., Gigli N., Savaré G. (2008) Gradient flows in metric spaces and in the space of probability measures, 2nd edn, Lectures in Mathematics ETH Zürich. Birkhäuser Verlag, Basel, ISBN 9783764387211. [Google Scholar]
 BessemoulinChatard M. (2012) Développement et analyse de schémas volumes finis motivés par la préservation de comportements asymptotiques. Application à des modèles issus de la physique et de la biologie, PhD Thesis, Université Blaise Pascal – ClermontFerrand II, 2012. URL http://tel.archivesouvertes.fr/tel00836514 [Google Scholar]
 Bear J., Bachmat Y. (1990) Introduction to modeling of transport phenomena in porous media, Kluwer Academic Publishers, Dordrecht, The Netherlands. [Google Scholar]
 Maury B., RoudneffChupin A., Santambrogio F. (2010) A macroscopic crowd motion model of gradient flow type, Math. Models Methods Appl. Sci. 20, 10, 1787–1821. ISSN 02182025. doi: 10.1142/S0218202510004799. URL http://dx.doi.org/10.1142/S0218202510004799. [Google Scholar]
 Kumar K., Pop I.S., Radu F.A. (2013) Convergence analysis of mixed numerical schemes for reactive flow in a porous medium, SIAM J. Numer. Anal. 51, 4, 2283–2308. [Google Scholar]
 Zarba R.L., Bouloutas E.T., Celia M. (1990) General massconservative numerical solution for the unsaturated flow equation, Water Resour. Res. 26, 7, 1483–1496. [Google Scholar]
 Jäger W., Kacur J. (1991) Solution of porous medium type systems by linear approximation schemes, Numer. Math. 60, 3, 407–427. [Google Scholar]
 Jäger W., Kacur J. (1995) Solution of doubly nonlinear and degenerate parabolic problems by relaxation schemes, RAIRO Modél. Math. Anal. Numér 29, 5, 605–627. [Google Scholar]
 Pop I.S., Radu F.A., Knabner P. (2004) Mixed finite elements for the Richards equation: linearization procedure, J. Comput. Appl. Math. 168, 1, 365–373. [Google Scholar]
 Radu F.A., Nordbotten J.M., Pop I.S., Kumar K. (2015) A robust linearization scheme for finite volume based discretizations for simulation of twophase flow in porous media, J. Comput. Appl. Math. 289, 134–141, ISSN 03770427. URL https://doi.org/10.1016/j.cam.2015.02.051. [Google Scholar]
 Radu F.A., Kumar K., Nordbotten J.M., Pop I.S. (2018) A robust, mass conservative scheme for twophase flow in porous media including hlder continuous nonlinearities, IMA J. Numer. Anal. 38, 2, 88420. doi: 10.1093/imanum/drx032. URL http://dx.doi.org/10.1093/imanum/drx032 [Google Scholar]
 Casulli V., Zanolli P. (2010) A nested Newtontype algorithm for finite volume methods solving Richards’ equation in mixed form, SIAM J. Sci. Comp. 32, 4, 2255–2273. doi: 10.1137/100786320. URL https://doi.org/10.1137/100786320 . [CrossRef] [Google Scholar]
 Younis R., Tchelepi H.A., Aziz K. (2010) Adaptively localized continuationNewton methodnonlinear solvers that converge all the time, SPE J. 15, 2, 526–544. [CrossRef] [Google Scholar]
 Wang X., Tchelepi H.A. (2013) Trustregion based solver for nonlinear transport in heterogeneous porous media, J. Comput. Phys. 253, 114–137. [Google Scholar]
 Lehmann F., Ackerer P.H. (1998) Comparison of iterative methods for improved solutions of the fluid flow equation in partially saturated porous media, Transp. Porous Media. 31, 3, 275–292. [Google Scholar]
 Bergamaschi L., Putti M. (1999) Mixed finite elements and Newtontype linearizations for the solution of Richards’ equation, Int. J. Numer. Meth. Eng. 45, 8, 1025–1046. [CrossRef] [Google Scholar]
 Radu F.A., Pop I.S., Knabner P. (2006) Newtontype methods for the mixed finite element discretization of some degenerate parabolic equations. Numerical mathematics and advanced applications, Springer. [Google Scholar]
 List F., Radu F.A. (2016) A study on iterative methods for solving Richards’ equation, Comput. Geosci. 1–13. [Google Scholar]
 Marchand E., Müller T., Knabner P. (2012) Fully coupled generalised hybridmixed finite element approximation of twophase twocomponent flow in porous media. Part II: numerical scheme and numerical results, Comput. Geosci. 16, 3, 691–708. doi: 10.1007/s1059601292791. URL https://doi.org/10.1007/s1059601292791. [Google Scholar]
 Marchand E., Müller T., Knabner P. (2013) Fully coupled generalized hybridmixed finite element approximation of twophase twocomponent flow in porous media. Part I: Formulation and properties of the mathematical model, Comput. Geosci. 17, 2, 431–442, ISSN 15731499. doi: 10.1007/s1059601393417. URL https://doi.org/10.1007/s1059601393417. [Google Scholar]
 Ben Gharbia I. (2012) Résolution de problèmes de complémentarité : application à un écoulement diphasique dans un milieu poreux, Thesis, Université Paris Dauphine  Paris IX, December 2012. URL https://tel.archivesouvertes.fr/tel00776617 [Google Scholar]
 Diersch H.J.G., Perrochet P. (1999) On the primary variable switching technique for simulating unsaturatedsaturated flows, Adv. Water Resour. 23, 3, 271–301. [Google Scholar]
 Brenner K., Cancès C. (2017) Improving Newton’s method performance by parametrization: The case of the Richards equation, SIAM J. Numer. Anal. 55, 4, 1760–1785. doi: 10.1137/16M1083414. URL https://doi.org/10.1137/16M1083414 [Google Scholar]
 Brenner K., Groza M., Jeannin L., Masson R., Pellerin J. (2017) Immiscible twophase Darcy flow model accounting for vanishing and discontinuous capillary pressures: application to the flow in fractured porous media, Comput. Geosci. 21, 5–6, 1075–1094. [Google Scholar]
 Ciarlet P.G. (1978) The finite element method for elliptic problems, NorthHolland Publishing Co., AmsterdamNew YorkOxford, ISBN 0444850287. Studies in Mathematics and its Applications, Vol. 4. [Google Scholar]
 Ern A., Guermond J.L. (2004) Theory and Practice of Finite Elements, volume 159 of Applied Mathematical Series, Springer, New York. [CrossRef] [Google Scholar]
 Franco Brezzi and Michel Fortin (1991) Mixed and hybrid finite element methods, volume 15 of Springer Series in Computational Mathematics, SpringerVerlag, New York. ISBN 0387975829 [Google Scholar]
 Arbogast T., Wheeler M.F., Yotov I. (1997) Mixed finite elements for elliptic problems with tensor coefficients as cellcentered finite differences, SIAM J. Numer. Anal. 34, 2, 828–852. doi: 10.1137/S0036142994262585. URL https://doi.org/10.1137/S0036142994262585. [Google Scholar]
 Aavatsmark I., Barkve T., Bøe Ø., Mannseth T. (1998) Discretization on unstructured grids for inhomogeneous, anisotropic media. I. Derivation of the methods, SIAM J. Sci. Comput. 19, 5, 1700–1716. doi: 10.1137/S1064827595293582. URL http://dx.doi.org/10.1137/S1064827595293582 . [Google Scholar]
 Edwards M.G., Rogers C.F. (1998) Finite volume discretization with imposed flux continuity for the general tensor pressure equation, Comput. Geosci. 2, 4, 259–290. doi: 10.1023/A:1011510505406. URL http://dx.doi.org/10.1023/A:1011510505406. [Google Scholar]
 Edwards M.G. (2002) Unstructured, controlvolume distributed, full tensor finitevolume schemes with flow based grids, Comput. Geosci. 6, 3–4, 433–452, ISSN 14200597. doi: 10.1023/A:1021243231313. URL https://doi.org/10.1023/A:1021243231313 . [Google Scholar]
 Agelas L., Guichard C., Masson R. (2010) Convergence of finite volume MPFA O type schemes for heterogeneous anisotropic diffusion problems on general meshes, Int. J. Finite 7, 2, 33. [Google Scholar]
 Arnold D., Brezzi F., Cockburn B., Marini L. (2002) Unified analysis of discontinuous Galerkin methods for elliptic problems, SIAM J. Numer. Anal. 39, 5, 1749–1779. doi: 10.1137/S0036142901384162. URL https://doi.org/10.1137/S0036142901384162. [Google Scholar]
 Rivière B. (2008) Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations, SIAM. doi: 10.1137/1.9780898717440. URL https://epubs.siam.org/doi/abs/10.1137/1.9780898717440 . [Google Scholar]
 Di Pietro D.A., Ern A. (2012) Mathematical aspects of discontinuous Galerkin methods, volume 69 of Mathématiques & Applications (Berlin) [Mathematics & Applications], Springer, Heidelberg, ISBN 9783642229794. doi: 10.1007/9783642229800. URL http://dx.doi.org/10.1007/9783642229800. [Google Scholar]
 Herbin R. (1995) An error estimate for a finite volume scheme for a diffusiononvection problem on a triangular mesh, Numer. Methods Partial Differ. Equ. 11, 2, 165–173. doi: 10.1002/num.1690110205. URL https://doi.org/10.1002/num.1690110205. [Google Scholar]
 Eymard R., Gallouët T., Herbin R. (2000) Finite volume methods, in: Ciarlet P.G., et al. (eds), Handbook of numerical analysis, NorthHolland: Amsterdam, p. 713 1020. [Google Scholar]
 Eymard R., Gallouët T., Guichard C., Herbin R., Masson R. (2014) TP or not TP, that is the question, Comput. Geosci. 18, 285–296. [Google Scholar]
 Hackbusch W. (1989) On first and second order box schemes, Computing 41, 4, 277–296, ISSN 0010485X. doi: 10.1007/BF02241218. URL https://doi.org/10.1007/BF02241218 . [CrossRef] [MathSciNet] [Google Scholar]
 Droniou J., Eymard R., Gallouët T., Herbin R. (2010) A unified approach to mimetic finite difference, hybrid finite volume and mixed finite volume methods, Math. Models Methods Appl. Sci. 20, 2, 265–295. [Google Scholar]
 Eymard R., Gallouët T., Herbin R. (2010) Discretization of heterogeneous and anisotropic diffusion problems on general nonconforming meshes sushi: a scheme using stabilization and hybrid interfaces, IMA J. Numer. Anal. 30, 4, 1009–1043. [CrossRef] [MathSciNet] [Google Scholar]
 Droniou J., Eymard R. (2006) A mixed finite volume scheme for anisotropic diffusion problems on any grid, Numer. Math. 105, 35–71. [Google Scholar]
 Brezzi F., Lipnikov K., Simoncini V. (2005) A family of mimetic finite difference methods on polygonal and polyhedral meshes, Math. Models Methods Appl. Sci. 15, 10, 1533–1551. [Google Scholar]
 Brezzi F., Lipnikov K., Shashkov M. (2005) Convergence of the mimetic finite difference method for diffusion problems on polyhedral meshes, SIAM J. Numer. Anal. 43, 5, 1872–1896. [Google Scholar]
 Domelevo K., Omnes P. (2005) A finite volume method for the Laplace equation on almost arbitrary twodimensional grids, M2AN: Math. Model. Numer. Anal. 39, 6, 1203–1249. [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]
 Droniou J. (2014) Finite volume schemes for diffusion equations: introduction to and review of modern methods, Math. Models Methods Appl. Sci. 24, 8, 1575–1620. [Google Scholar]
 Droniou J., Eymard R., Gallouët T., Guichard C., Herbin R. (2018) The gradient discretisation method, Vol. 42, Mathématiques et Applications, Springer International Publishing, https://doi.org/10.1007/9783319790428. [Google Scholar]
 Eymard R., Guichard C., Herbin R. (2012) Smallstencil 3D schemes for diffusive flows in porous media, ESAIM: Math. Model. Numer. Anal. 46, 2, 265–290. doi: 10.1051/m2an/2011040. URL http://dx.doi.org/10.1051/m2an/2011040 . [Google Scholar]
 Eymard R., Guichard C., Herbin R. (2011) Benchmark 3D: the VAG scheme, in: Fort J., Fürst J., Halama J., Herbin R., Hubert F. (eds), Finite Volumes for Complex Applications VI Problems & Perspectives, volume 4 of Springer Proceedings in Mathematics, Springer, Berlin Heidelberg, pp. 1013–1022. ISBN 9783642206702. doi: 10.1007/9783642206719_99. URL http://dx.doi.org/10.1007/9783642206719_99 . [Google Scholar]
 Cancès C., Guichard C. (2017) Numerical analysis of a robust free energy diminishing finite volume scheme for parabolic equations with gradient structure, Found. Comput. Math. 17, 6, 1525–1584. doi: 10.1007/s1020801693286. [CrossRef] [Google Scholar]
 Cancès C., ChainaisHillairet C., Krell S. (2017) A nonlinear Discrete Duality Finite Volume Scheme for convection diffusion equations, in: Cancès C., Omnes P. (eds), FVCA8 2017 – International Conference on Finite Volumes for Complex Applications VIII, volume 199 of Springer Proceedings in Mathematics & Statistics, Lille, France, Springer International Publishing, pp. 439–447. URL https://hal.archivesouvertes.fr/hal01468811. [Google Scholar]
 Cancès C., ChainaisHillairet C., Krell S. (2017) Numerical analysis of a nonlinear freeenergy diminishing Discrete Duality Finite Volume scheme for convection diffusion equations, Comput Methods Appl. Math. doi: 10.1515/cmam20170043. URL https://hal.archivesouvertes.fr/hal01529143. Special issue on “Advanced numerical methods: recent developments, analysis and application”. [Google Scholar]
 Cancès C., Guichard C. (2016) Convergence of a nonlinear entropy diminishing Control Volume Finite Element scheme for solving anisotropic degenerate parabolic equations, Math. Comp. 85, 298, 549–580. [CrossRef] [Google Scholar]
 Chavent G., Jaffré J. (1986), Mathematical Models and Finite Elements for Reservoir Simulation, Vol. 17, Stud. Math. Appl. edition, NorthHolland, Amsterdam. [Google Scholar]
 Antontsev S.N., Kazhikhov A.V., Monakhov V.N. (1990) Boundary value problems in mechanics of nonhomogeneous fluids, vol. 22 of Studies in Mathematics and its Applications, NorthHolland Publishing Co., Amsterdam, ISBN 0444883827. Translated from the Russian. [Google Scholar]
 Gagneux G., MadauneTort M. (1996) Analyse mathématique de modèles non linéaires de l’ingénierie pétrolière, vol. 22 of Mathématiques & Applications (Berlin) [Mathematics & Applications], SpringerVerlag, Berlin, ISBN 3540605886. [Google Scholar]
 Chen Z. (2001) Degenerate twophase incompressible flow. I. Existence, uniqueness and regularity of a weak solution, J. Diff. Equ. 171, 2, 203–232. [CrossRef] [MathSciNet] [Google Scholar]
 Nochetto R.H., Verdi C. (1988) Approximation of degenerate parabolic problems using numerical integration, SIAM J. Numer. Anal. 25, 4, 784–814. doi: 10.1137/0725046. URL https://doi.org/10.1137/0725046. [Google Scholar]
 Arbogast T., Wheeler M.F., Zhang N.Y. (1996) A nonlinear mixed finite element method for a degenerate parabolic equation arising in flow in porous media, SIAM J. Numer. Anal. 33, 4, 1669–1687. doi: 10.1137/S0036142994266728. URL http://dx.doi.org/10.1137/S0036142994266728. [Google Scholar]
 Eymard R., Gallouët T., Hilhorst D., Naït Slimane Y. (1998) Finite volumes and nonlinear diffusion equations, RAIRO Modél. Math. Anal. Numér. 32, 6, 747–761. [Google Scholar]
 Eymard R., Gutnic M., Hilhorst D. (1999) The finite volume method for Richards equation, Comput. Geosci. 3, 3–4, 259–294. doi: 10.1023/A:1011547513583. URL http://dx.doi.org/10.1023/A%3A1011547513583. [Google Scholar]
 Woodward C.S., Dawson C.N. (2000) Analysis of expanded mixed finite element methods for a nonlinear parabolic equation modeling flow into variably saturated porous media, SIAM J. Numer. Anal. 37, 3, 701–724. doi: 10.1137/S0036142996311040. URL https://doi.org/10.1137/S0036142996311040. [Google Scholar]
 Eymard R., Gallouët T., Herbin R., Michel A. (2002) Convergence of finite volume schemes for parabolic degenerate equations, Numer. Math. 92, 41–82. [Google Scholar]
 Pop I.S. (2002) Error estimates for a time discretization method for the Richards’ equation, Comput. Geosci. 6, 141–160. [Google Scholar]
 Radu F.A., Pop I.S., Knabner P. (2004) Order of convergence estimates for an Euler implicit, mixed finite element discretization of Richards’ equation, SIAM J. Numer. Anal. 42, 4, 1452–1478. [Google Scholar]
 Eymard R., Hilhorst D., Vohralík M. (2006) A combined finite volumenonconforming/mixedhybrid finite element scheme for degenerate parabolic problems, Numer. Math. 105, 1, 73–131. doi: 10.1007/s002110060036z. URL http://dx.doi.org/10.1007/s002110060036z [Google Scholar]
 Radu F.A., Pop I.S., Knabner P. (2008) Error estimates for a mixed finite element discretization of some degenerate parabolic equations, Numer. Math. 109, 2, 285–311. doi: 10.1007/s0021100801399. URL http://dx.doi.org/10.1007/s0021100801399. [Google Scholar]
 Angelini O., Brenner K., Hilhorst D. (2013) A finite volume method on general meshes for a degenerate parabolic convectionreactiondiffusion equation, Numer. Math. 123, 219–257, ISSN 0029599X. doi: 10.1007/s0021101204855. URL http://dx.doi.org/10.1007/s0021101204855. [Google Scholar]
 Chen Z., Ewing R.E. (1997) Fully discrete finite element analysis of multiphase flow in groundwater hydrology, SIAM J. Numer. Anal. 34, 6, 2228–2253. [Google Scholar]
 Chen Z., Ewing R.E. (2001) Degenerate twophase incompressible flow. III. Sharp error estimates, Numer. Math. 90, 2, 215–240, ISSN 0029599X. doi: 10.1007/s002110100291. URL http://dx.doi.org/10.1007/s002110100291. [Google Scholar]
 Michel A. (2003) A finite volume scheme for twophase immiscible flow in porous media, SIAM J. Numer. Anal. 41, 4, 1301–1317. [Google Scholar]
 Epshteyn Y., Rivière B. (2009) Analysis of hp discontinuous Galerkin methods for incompressible twophase flow, J. Comput. Appl. Math. 225, 2, 487–509. doi: 10.1016/j.cam.2008.08.026. URL https://doi.org/10.1016/j.cam.2008.08.026. [Google Scholar]
 Brenner K., Masson R. (2013) Convergence of a vertex centered discretization of twophase Darcy flows on general meshes, Int. J. Finite 10, 1–37. [Google Scholar]
 Cancès C., Pop I.S., Vohralík M. (2014) An a posteriori error estimate for vertexcentered finite volume discretizations of immiscible incompressible twophase flow, Math. Comp. 83, 285, 153–188. doi: 10.1090/S002557182013027238. URL http://dx.doi.org/10.1090/S002557182013027238. [CrossRef] [MathSciNet] [Google Scholar]
 Cancès C., Nabet F., Vohralik M. (2018) Convergence and a posteriori error analysis for energystable finite element approximations of degenerate parabolic equations, in preparation. [Google Scholar]
 Forsyth P.A. (1991) A control volume finite element approach to NAPL groundwater contamination, SIAM J. Sci. Statist. Comput. 12, 5, 1029–1057. [CrossRef] [MathSciNet] [Google Scholar]
 Ait Hammou Oulhaj A., Cancès C., ChainaisHillairet C. (2018) Numerical analysis of a nonlinearly stable and positive Control Volume Finite Element scheme for Richards equation with anisotropy, ESAIM Math. Model. Numer. Anal. 52, 1532—1567. [Google Scholar]
 Ait Hammou Oulhaj A. (2018) Numerical analysis of a finite volume scheme for a seawater intrusion model with crossdiffusion in an unconfined aquifer, Numer. Methods Partial Differ. Equ. 34, 3, 857–880. doi: 10.1002/num.22234. URL https://doi.org/10.1002/num.22234. [Google Scholar]
 Cancès C., Nabet F. (2017) Finite volume approximation of a degenerate immiscible twophase flow model of CahnHilliard type, in: Cancès C., Omnes P. (eds), Finite Volumes for Complex Applications VIII  Methods and Theoretical Aspects : FVCA 8, Lille, France, June 2017, number 199 in Proceedings in Mathematics and Statistics, Cham, Springer International Publishing, pp. 431–438, ISBN 9783319573977. doi: 10.1007/9783319573977_36. http://dx.doi.org/10.1007/9783319573977_36. [CrossRef] [Google Scholar]
 Otto F., Weinan E. (1997) Thermodynamically driven incompressible fluid mixtures, J. Chem. Phys. 107, 23, 10177–10184. [Google Scholar]
 Cancès C., Matthes D., Nabet F. (2017) A twophase twofluxes degenerate CahnHilliard model as constrained Wasserstein gradient flow, HAL, hal01665338, December 2017. URL https://hal.archivesouvertes.fr/hal01665338. [Google Scholar]
 Herbin R., Hubert F. (2008) Benchmark on discretization schemes for anisotropic diffusion problems on general grids, in: Eymard R., Herard J.M. (eds), Finite Volumes for Complex Applications V, Wiley, pp. 659–692. https://www.latp.univmrs.fr/fvca5/benchmark/ [Google Scholar]
 Eymard R., Herbin R., Michel A. (2003) Mathematical study of a petroleumengineering scheme, M2AN: Math. Model. Numer. Anal. 37, 6, 937–972. [CrossRef] [Google Scholar]
 ChainaisHillairet C., Filbet F. (2007) Asymptotic behaviour of a finitevolume scheme for the transient driftdiffusion model, IMA J. Numer. Anal. 27, 4, 689–716. doi: 10.1093/imanum/drl045. URL http://dx.doi.org/10.1093/imanum/drl045. [CrossRef] [Google Scholar]
 BessemoulinChatard M., ChainaisHillairet C. (2017) Exponential decay of a finite volume scheme to the thermal equilibrium for driftdiffusion systems, J. Numer. Math. 25, 3, 147—168. [CrossRef] [Google Scholar]
 Filbet F., Herda M. (2017) A finite volume scheme for boundarydriven convectiondiffusion equations with relative entropy structure, Numer. Math. URL https://hal.archivesouvertes.fr/hal01326029. [Google Scholar]
 Ganis B., Kumar K., Pencheva G., Wheeler M., Yotov I. (2014) A global Jacobian method for mortar discretizations of a fully implicit twophase flow model, Multiscale Model. Simul. 12, 4, 1401–1423. doi: 10.1137/140952922. URL https://doi.org/10.1137/140952922. [Google Scholar]
A Some technical details
A.1 About Estimate (39)
We now give some details on the derivation of estimate (39). Note that the derivation of estimate (36) relies on similar arguments.
Let be defined by (38), then multiplying (28) by and summing over leads to(71)where
The term A _{ n } can be underestimated thanks to the convexity inequality (37):(72)
Thanks to the characterization (35) of the upstream mobility, the term B _{ n } can be underestimated bywhatever the mean value between and . Choosingwhich lies between and because of the definition (38) of , one gets thatwhereis the quantity that we want to bound from above, and
The quantitywhich is an approximation of at x _{ K }, is supposed to be bounded from below by some quantity M depending only on the regularity of the mesh. Since the method preserves the positivity of the saturations , one obtains that(73)thanks to the global conservativity of the scheme (28). As a consequence, one gets that
One concludes by noticing that Φ is bounded from below, hence
A.2 Relation (45) as an optimality condition
Before justifying why (45) can be seen as an optimality condition for the discrete energy under the mass contraint (47), let us first notice that sincethe discrete counterpartof the energy (5) can be rewritten as
The second term in the above equality does not depend on s ^{ n }, hence it can be omitted in (48) and we can write the discrete energy as a function of only, i.e., .
A second preliminary remark is the following: if ϕ _{ K } m _{ K } = 0 for some , then has no influence on the energy . Its value can be fixed arbitrarily, for instance by (45). We assume for simplicity that ϕ _{ K } m _{ K } > 0 for all even though this assumption can be easily bypassed.
Let us now go to the constrained optimization problemwhere, setting , we denoted
Note that s ^{ n } necessarily belongs to otherwise would be infinite. The problem is equivalent to the saddlepoint problem
We can swap the and the and optimize w.r.t. s _{ n }, so that we get the optimality conditionwhich is equivalent to (45).
All Tables
Approximation of (42) with a nonlinear DDFV scheme [62, 63] on the Kershaw mesh family from [93], final time T = 0.25. M is the mesh index, τ is the time step, errgs and errs respectively stand for the L ^{2}(Ω × (0, T)) error on ∇ s and the L ^{∞}((0, T); L ^{2}(Ø)) error on s, whereas ordgs and ords are the corresponding convergence orders. N _{max} and N _{mean} are the maximal and mean numbers of Newton iterations, and s _{min} minimal value of the approximation of s during the whole simulation.
All Figures
Fig. 1 The maximal monotone graph π depicted on top is parametrized as with some monotone functions (red) and (blue) depicted in the bottom picture. Here, and satisfy . 

In the text 
Fig. 2 L ^{2} (Ω × (0, t _{f})) error as a function of the mesh size for the scheme (28) in the isotropic case κ = 1 (solid red) and anisotropic case κ = 20 (solid blue) corresponding to upstream mobilities (29). Comparison with the solution to the linear scheme (43) in the isotropic case (dashed red) and anisotropic case (dashed blue). 

In the text 
Fig. 3 Evolution of the relative energy as a function of time. We observe the exponential convergence towards until the machine precision is reached. 

In the text 
Fig. 4 The Kershaw meshes are highly deformed topologically cartesian grids. The coarsest mesh of the family is depicted here. 

In the text 
Fig. 5 Convergence towards the analytical solution (42) of the FokkerPlanck equation (41) for different values of the anisotropy ratio . 

In the text 
Fig. 6 Plot of the log of the relative energy (in log scale) as a function of time. 

In the text 
Fig. 7 Schematic representation of the injection and production wells. 

In the text 
Fig. 8 Approximate saturation s _{n,h } at time t = 0.2 (top left), t = 0.5 (top right), t = 1.5 (bottom left) and t = 2.5 (bottom right). 

In the text 