Energy stable numerical methods for porous media flow type problems

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.


Two-phase porous media flows
Incompressible two-phase porous media flows are often modeled by the following set of equations. In the absence of source terms, the volume of the phase a 2 fn; wg (n and w stand for non-wetting and wetting respectively) is locally conserved along time as a consequence of where / denotes the porosity of the rock (supposed to be constant w.r.t. time), where s a denotes the saturation of the phase a and v a denotes the filtration speed of the phase a. It is classically assumed that the phase filtration speeds obey the generalized Darcy law v a ¼ À k r;a ðs a Þ l a K $p a À q a g ð Þ : The intrinsic permeability tensor K of the porous medium is a definite positive and symmetric tensor field whereas the relative permeability k r;a of the phase a is an increasing function of the saturation satisfying k r;a ð0Þ ¼ 0 (we neglect the residual saturations). The variations of the viscosities l a and of the densities q a 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: The second one links the capillary pressure to the nonwetting phase saturation in a monotone way: where p is a maximal monotone graph from ½0; 1 to R 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][2][3][4][5]). Once complemented by no-flux conditions on the boundary of the porous domain X, the model (1)-(4) has a very particular variational structure. As depicted in [6] (see also [7][8][9]), this problem can be reinterpreted as the generalized gradient flow [10] of the energy with EðsÞ ¼ Pðs n Þ À X a2fn;wg s a q a g Á x þ vðsÞ: 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 P : R ! R [ fþ1g is a convex function that is finite in [0,1], infinite outside of [0,1], and satisfies pðsÞ ¼ @PðsÞ for s 2 ½0; 1; ð6Þ where @P s ð Þ ¼ p 2 R f jP s ð Þ À P s ð Þ À p s À s ð Þ!0 for all s 2 DðPÞg is the subdifferential of P at s. The functional E : R 2 ! R [ fþ1g is convex and its subdifferential @EðsÞ is made of the couples h ¼ h n ; h w ð Þ¼ ðp n À q n g Á x; p w À q 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 EðsÞ is a decreasing function of time and that the dynamics aims at maximizing this decay. This yields an energy/ dissipation of the form for all t ! 0. The energy/dissipation relation (7) can be obtained by multiplying formally the equation (1) by h a and by summing over a 2 fn; wg. As a consequence, stable steady states s 1 are local minimizers of E for which each phase is hydrostatic on its support, i.e., In order to compute in a accurate way the long-time 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 Eð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.

A simplified model problem
Our purpose can already be illustrated on the single scalar equation P being a convex and coercive internal energy functional as previously and É being a smooth external (possibly gravitational) potential. The mobility function g is nondecreasing on DðPÞ \ R þ and satisfies gð0Þ ¼ 0 and gðsÞ > 0 if s > 0. Here, we used the notation DðPÞ ¼ fs 2 R j PðsÞ < 1g for the domain of P and we assume that 0 2 DðPÞ. If P is defined on the whole R þ , i.e., DðPÞ \ R þ ¼ R þ , then we assume moreover that P is superlinear: PðsÞ s ¼ þ1: 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 P, p, and g for negative values of s when this is possible. In the case where @Pð0Þ contains a finite value, i.e., if then P can be artificially extended on ðÀ1; 0Þ [ DðPÞ into a convex function (still denoted by P) with a single valued subdifferential at 0, i.e., @Pð0Þ ¼ p ð0Þ ¼ pð0Þ. This can be done for instance by prescribing if s < 0. The function g can for instance be extended by setting gð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 W ¼ Àqg Á x. The capillary energy function P is then defined on R þ as the antiderivative of the capillary pressure function. In usual models (see [14], pp. 343-345), DðPÞ ¼ ½0; 1 and condition (10) does not hold. This also includes a linear Fokker-Planck equation (but written in a nonlinear form) when gðsÞ ¼ s and pðsÞ ¼ logðsÞ or models for crowd motions [15]. We are interested in the computation of non-negative solutions corresponding to initial data There is still a generalized gradient flow structure corresponding to (9), the energy being given by The counterpart to (7) is obtained by multiplying formally (9) by p þ W and writes EðsÞðtÞ þ Similarly to what was discussed in Section 1.1, one aims to build numerical schemes for (9) that are accurate in the long-time limit. The steady states for (9) are now given by 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 p 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 p 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][18][19][20][21][22]). There were also important efforts carried out to fix the difficulties of Newton's method [23][24][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 Picard-type 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 2 pðsÞ that is often rephrased as a complementary constraint and then solving the problem with a non-smooth Newton method (see for instance [30][31][32]). Another classical solution consists in partitioning X at each time t into a part X s ðtÞ where s is chosen as a primary variable and a part X 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 p in a suitable way. More precisely, the graph can be described by two functions w 7 ! sðwÞ and w 7 ! pðwÞ such that p 2 pðsÞ () 9 w s:t: s ¼ sðwÞ and p ¼ pðwÞ; cf. Figure 1. We assume moreover that the parametrization is non-degenerate in the sense that s þ p is a strictly increasing function. 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 ðs; pÞ of p. One example is the resolvents of p and p À1 , i.e., s ¼ ðid þ pÞ À1 and p ¼ ðid þ p À1 Þ À1 . 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.

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 subject to no-flux 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], Multi-Point Flux Approximation (MPFA) finite volumes [40][41][42][43], or Discontinuous Galerkin method [ [44][45][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 U 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 X where m K ð Þ K2U & R þ are such that P K m K ¼ jXj. The compatibility condition on the right-hand side the vector unknown and by A the matrix corresponding to the system (18), i.e., We require the transmissivity coefficients a KL (and thus the matrix A) to be symmetric In what follows, S denotes the set of the couples ðK; LÞ 2 U 2 such that the transmissivity a KL is different from 0. Thanks to the symmetry property (19), the scheme (18) is equivalent to: Many numerical schemes can write as (18) and (19). This is for instance the case of the classical Two-Point Flux Approximation (TPFA) finite volume scheme [47][48][49], where as soon as the cells K and L share an edge r 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 K ¼ jI must be isotropicj KL in (21) is the harmonic mean of the permeabilities j K and j L associated to the cells K and L -, and the edge r 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 P 1 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 U of a simplicial mesh T of X. Denoting by / K the Lagrange basis function associated to the vertex K 2 U, the transmissivities a KL are defined by and the lumped mass associated to the vertex K is m K ¼ R X / K dx: Note that in this context, the transmissivities a KL may be negative (for instance in the case where T does not satisfy the Delaunay condition if K ¼ I).
In the case where the transmissivities are nonnegative 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 M 2 M and positive definite symmet- The A crucial property of the local diffusion matrices A M is that the condition number of A M can be bounded by a quantity depending only on the condition number of Kdx and on the regularity (in Ciarlet's sense [36], see also [37]) of the simplex M, i.e., The HFV (or SUSHI [52]) also enter the framework (24). In this case, M is the set of the control volumes that can be quite general (non-convex cells having various number of edges). To each cell M, there is one cell unknown u M;0 and ' M edge unknowns u M;1 ; . . . ; u M;' M where ' M denotes the number of edges of the element M, and the inner variation vector d M v is given by The matrix A M 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 K 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 M 2 M 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, M is the set of the diamond cells and once again, the local matrices A M have bounded conditioning numbers if K 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 A is a symmetric definite positive matrix, the lowest eigenvalue of A being bounded away from 0 uniformly w.r.t. the mesh size: Combining this information with (26), it was proven in [64] that 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.

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 so-called Kirchhoff transform and global pressure (see for instance [65][66][67][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][70][71][72][73][74][75][76][77][78][79], whereas for multiphase porous media flows, we refer to [80][81][82][83][84][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 long-time 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.

Upstream mobility schemes
Before addressing the two-phase flow problem (1)-(4), let us first focus on the simpler scalar problem (9).

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 s; p be parametrizations of the graph p as in (16), then the numerical scheme writes: 8K 2 U; 8n ! 1; where we have used the notation N K ¼ fL j ðK ; LÞ 2 Sg for the neighbors of K . The time step s n is not necessarily uniform and can depend on n ! 1. The initial saturation s 0 K À Á K 2U is assumed to be given and once w n K has been computed for n ! 1, one sets s n K ¼ sðw n K Þ: Following [64] and [87], the mobility g n KL is chosen thank to an upwinding that takes the sign of a KL into account: gðsðw n L ÞÞ otherwise: ( ð29Þ Note that g n KL ¼ g n LK , so that one can write the scheme (28) under the equivalent form X In the above formula and as in (20), v is an arbitrary element of R U .

Some properties of the scheme
As a consequence of the choice (29) for the mobility g n KL , the method preserves the non-negativity. More precisely, this means that if w n ¼ w n K À Á K is a solution of the scheme (28), then so does its projectionw n ¼ maxðw n Fig. 1). So without loss of generality, we can assume that This property still holds for transmissivities a KL violating the monotonicity condition (23). Note that the inequality (31) is of interest only if w H > À1.
Let us now turn to the energy estimate. Choosing v ¼ pðw n K Þ þ É K À Á K in the weak formulation (30) and using the convexity inequality s À s p ! PðsÞ À Pð sÞ; for all ðs; sÞ 2 DðPÞ 2 and for all p 2 @PðsÞ ¼ pðsÞ, we obtain the following discrete counterpart of (14): where 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 non-negative, thus the discrete energy is diminishing along time, i.e., E U ðs n Þ E U ðs nÀ1 Þ: Moreover, summing over n provides a control on the numerical dissipation 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 ðs; pÞ, 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 whatever the non-decreasing function u. In the above relations, I n KL denotes the interval with extremal points w n K and w n L and we have used the monotonicity of g. In the case where g s is bounded (this is natural in the porous media flow context since s is), we can take advantage of this formulation with u ¼ p, of the inequality (27), on the Lipschitz regularity of the external potential É to show (cf. [88] A similar inequality can be obtained even in the case where g s 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 g is bounded and if p À1 is a function (i.e., when there are no hyperbolic degeneracy in the problem (9)) when the mesh size and the time step s 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 uðw n Þ thanks to (35). For s 2 ½0; 1, we denote by s À1 ðsÞ ¼ ½wðsÞ; wðsÞ the largest interval of ½À1; þ1 such that sðwÞ ¼ s for all w 2 s À1 ðsÞ. Then given a non-decreasing onto function u : R ! R, there corresponds a convex and coercive function È : ½0; 1 ! R [ fþ1g such that @ÈðsÞ ¼ u s À1 ðsÞ. a KL pðw n K Þ À pðw n L Þ À Á sðw n K Þ À sðw n L Þ À Á provided É is regular enough and the initial entropy is finite: X K2U m K Èðs 0 K Þ < 1: Details for estimate (39) are provided in Appendix A.1. Such an estimate is the corner-stone 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 Cahn-Hilliard system with a very similar mathematical structure to (1)-(4) (see [91,92]) was considered.

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 where X ¼ ð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
where e 2 ¼ ð0; 1Þ t . We choose the initial condition and the no-flux boundary condition in accordance with the exact solution with a ¼ p 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 X. We use successively refined Delaunay grids from the FVCA5 benchmark [93]. In the isotropic case (corresponding to j ¼ 1), the scheme is monotone, i.e., condition (23) is fulfilled. This is no longer true in the anisotropic case j ¼ 20.
We also compute the solution to our scheme with an anisotropy ratio j ¼ 20. The numerical solutions are compared with those computed with the following linear scheme with centered convection: 8K 2 U; 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 j ¼ 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 ðX Â ð0; t f ÞÞ error corresponding to the numerical solutions produced by the schemes (28) and (43) for Delaunay triangular grids from [93].
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 j ¼ 20. As expected because of the upwinding procedure, the scheme (28) is first order accurate in space, i.e., jju À u h jj L 2 ðXÂð0;t f ÞÞ Ch; ð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 non-monotone 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 Cahn-Hilliard systems [90] or multiphase porous media flows [8,89,94].

Long-time behavior of the scheme
It is natural to wonder if the numerical scheme is able to reproduce in an accurate way the long-time 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 long-time behavior of Finite Volume schemes for convection-diffusion problems has been the purpose of several contributions (see, e.g., [13,62,[95][96][97]).
We go back to the two-phase 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 long-time limit boils down to the following alternative: either s 1 n 2 f0; 1g or pðs 1 for some c 2 R: The parameter c is determined by the conservation of mass Z X /s 1 The steady state (45) can be discretized directly into with c 0 fixed so that X K2U m K / K s 1 n;K ¼ In the above relation, / K ¼ /ðx K Þ is the discrete porosity. In particular, in the classical case where p : ð0; 1Þ ! R is a (single-valued) function, the following alternative holds: either s 1 n;K 2 f0; 1g or pðs 1 n;K Þ ¼ ðq n À q w Þg Á x K À c 0 for all K 2 U. Let E U be the discrete counterpart of the energy (5), i.e., Pðs n n;K Þ þ s n n;K ðq w À q n Þg Á x K m K ; ð48Þ then the energy is decreasing along time, i.e., E U ðs n n Þ E U ðs nÀ1 n Þ; n ! 1: One can show (cf. Appendix A.2) that s 1 is a minimizer of E U under the constraint (47). Therefore the relative energy E U ðs n n Þ À E U ðs 1 n Þ is non-negative. If the capillary pressure p is an increasing function on ð0; 1Þ, then the relative energy vanishes if and only if s n n ¼ s 1 n . This quantity can be used to illustrate the convergence of s n n towards s 1 n as n tends to 1. In the case depicted on Figure 3, we set X ¼ ðÀ1=2; 1=2Þ 2 , K ¼ I, k r;a ðsÞ ¼ s, l n ¼ 10, l w ¼ 1, q n ¼ 0:87, q w ¼ 1, g ¼ ð0; À9:81Þ T , / 1, s ini ðxÞ ¼ e À4jxj 2 , and 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 s 1 n given by (45) and (46) is a steady solution to the scheme.

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.
(i) No Kirchhoff transform: the scheme should only involve quantities with a clear physical meaning. (ii) Nonlinear stability: the scheme must fulfill a discrete counterpart of the energy/energy dissipation relation (14). (iii) 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). (iv) Second order accuracy: if the solution to (9) is smooth enough (say C 2 ), the error between the exact

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 with a centered choice for the mobility: Multiplying the scheme by s n pðw n K Þ þ É K À Á and summing over K 2 U leads once again to Already in the case developed in Section 2 where g n KL was chosen thanks to upwinding, the sign of the second term of the left-hand 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 The above relation must hold for any v 2 U. The matrix B M ðw n Þ 2 R ' M Â' M is called the local dissipation tensor. It must incorporate the local diffusion A M but also the mobilities gðsðw n ÞÞ. In order to ensure the dissipation property for the scheme, we want B M ðw n Þ to be symmetric semi-definite positive. In [61], we proposed to choose with for the particular choice The fact that the energy is diminishing along time is obtained by choosing v ¼ pðw n Þ þ W and by applying a simple convexity inequality, leading to the second term being non-negative since B M ðw n Þ is semidefinite positive. Additionally, the method is globally mass conservative, i.e., X where s n K ¼ sðw n K Þ. 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 g M i ðvÞ appearing in (52) is chosen as a convex combination of gðsðv K M j ÞÞ In order to assess that the scheme (50)-(52) converges, the numerical mobilities have to satisfy an additional coercivity condition, that is for some uniform a > 0. The condition (56) is clearly satisfied by the choice (53) with a ¼ 1=2, but it prohibits the choice g M i ðvÞ ¼ gðsðv K M 0 ÞÞ 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 g M i ðvÞ independent on i, i.e.,  Fig. 3. Evolution of the relative energy E U ðs n Þ À E U ðs 1 Þ as a function of time. We observe the exponential convergence towards 0 until the machine precision is reached.
The matrix H M ðvÞ then reduces to ffiffiffiffiffiffiffiffiffiffiffiffi g M ðvÞ p I ' M and commutes with A M , leading to the simpler formula B M ðvÞ ¼ g M ðvÞA M . In view of the constraint (56), a natural choice for g M ðvÞ is This choice was successfully used in [86] for a method based on conformal P 1 finite elements with mass lumping (here ' M ¼ d), and in [62,63] for a nonlinear DDFV method (here d ¼ 2 and ' M ¼ 3).

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 The functions g, s and p are uniformly continuous on R, then F 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 lim s!0 þ pðsÞ ¼ À1; ð60Þ (with a slight abuse of notation), then one can show that for some f depending on the time step s 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 is rewritten under the form for three different choices of nonlinearities, that are (a) gðsÞ ¼ 1 and pðsÞ ¼ jsjs (recall that we need to extend the nonlinearities for negative saturations if (10) holds); (b) gðsÞ ¼ 2jsj and pðsÞ ¼ s; (c) gðsÞ ¼ 2s 2 and pðsÞ ¼ logðsÞ.
The case (a) does not enter our framework since gð0Þ ¼ 1 6 ¼ 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 ¼ fðx; yÞ 2 ð0; 1Þ 2 g, K ¼ 1 0 0 100 and let us approximate the exact solution to (62) defined by sðx; y; tÞ ¼ maxð0; 2t À xÞ ð 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 @ t s þ $ Á ðse 2 À $sÞ ¼ 0; 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.
As expected, since pð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 super-convergence 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 test-case of the linear Fokker-Planck equation written in a nonlinear form described in Section 2.3. We discretize it with the energy stable P 1 finite element scheme with mass lumping proposed in [86]. In this scheme, the mobilities g M ðvÞ are chosen according to formula (58). The final time is set to t f ¼ 0:25. The L 2 ðð0; t f Þ Â XÞ 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 j.

About the long-time behavior of the scheme
We aim now to illustrate the long-time behavior of the scheme. One discretizes the equation (65) complemented with no-flux boundary conditions following the methodology of [61]. In particular, we take s ¼ Id, hence we denote by s n ¼ w n the vector of the unknown saturations. The scheme (51) then rewrites for any v 2 U. The mobilities are discretized thanks to formula (54). The mass is conserved along time, i.e., X For any q > 0, the vector s 1 ¼ s 1 is a steady solution to the scheme (66) since d M logðs n Þ þ W ð Þ¼0 for all M 2 M. Then the expected long-time limit as n ! 1 is s 1 where q has been tuned so that X To illustrate this fact, we plot on Figure 6 the evolution of the relative energy E U ðs n Þ À E U ðs 1 Þ as a function of the discrete time ns 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 1 .
We show in Figure 6 that the relative energy converges exponentially fast towards 0, providing a discrete counterpart to the relation s ex ðx; tÞ À e ÀWðxÞ Ce Àat that is deduced from (42). We observe on Figure 6 that the scheme preserves exactly (up to machine precision) the long-time behavior of the equation, i.e., the long-time limit s 1 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 non-negative s n . For a behavior like the one depicted on Figure 6, the scheme must (at least) preserve positivity and admit s 1 defined by (67) as a steady state.

Application to two-phase porous media flows
We show now preliminary results obtained thanks to the methodology of this section to simulate two-phase porous media flows. We consider an anisotropic quarter five spot problem with injection in a top-left corner and extraction in the bottom-right corner. We neglect gravity, so that the equations boil down into /@ t s n þ $ Á g n ðs n ÞK$ p w þ pðs n Þ ð Þ ð Þ ¼ q n ðs n ; xÞ; ð69Þ À/@ t s n þ $ Á g w ðs n ÞK$p w ð Þ¼q w ðs n ; xÞ: 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   Table 4. 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, s is the time step, errgs and errs respectively stand for the L 2 ðX Â ð0; T ÞÞ error on rs and the L 1 ðð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.  relations s w ¼ 1 À s n and p n ¼ p w þ pðs n Þ. We fix the nonlinearities g n , g w , p, f n , and f w as in [98]:  Fig. 7). Defining the fractional flow functions by f n ðsÞ ¼ g n ðsÞ g n ðsÞ þ g w ðsÞ and f w ðsÞ ¼ g w ðsÞ g n ðsÞ þ g w ðsÞ ; then the source terms q n and q w are defined by q a ðs; xÞ ¼ f a ð0Þ1 x inj ðxÞ À f a ðsÞ1 x prod ðxÞ; a 2 fn; wg: The permeability tensor K ¼ 1 0 0 5 is anisotropic and we take a constant saturation profile / 1. The system (69)- (70) is complemented with no-flux boundary conditions and with the initial condition s n ðÁ; 0Þ 0:9.
For the discretization, we use conforming P 1 finite elements with mass lumping on a structured triangulation T made of 5000 triangles and a constant time step s ¼ 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 / n a;h in the Raviart-Thomas-Nedelec space RTN 1 ðT Þ such that s n a;h À s nÀ1 a;h s þ $ Á / n a ¼ q n a;h ; a 2 fn; wg; where s n a;h denotes the P 1 approximation of the saturation s a at time t n ¼ P n k¼1 s k . We plot snapshots of the saturation s n on Figure 8.

Injection zone
Extraction zone Fig. 7. Schematic representation of the injection and production wells. the energy E U ðs n n Þ. Its value can be fixed arbitrarily, for instance by (45). We assume for simplicity that / K m K > 0 for all K 2 U even though this assumption can be easily bypassed.
Let us now go to the constrained optimization problem s 1 n 2 argmin s n 2X E U ðs n n Þ; where, setting m ¼ / K m K ð Þ2R U , we denoted Note that s n necessarily belongs to ½0; 1 U otherwise E U would be infinite. The problem is equivalent to the saddle-point problem min sn2R U max c 0 2R E U ðs n Þ þ c 0 ðs n À s 0 n Þ Á m È É : We can swap the min and the max and optimize w.r.t. s n , so that we get the optimality condition which is equivalent to (45).