Simulation of multiphase porous media flows with minimizing movement and finite volume schemes

The Wasserstein gradient flow structure of the PDE system governing multiphase flows in porous media was recently highlighted in [C. Canc\`es, T. O. Gallou\"et, and L. Monsaingeon, {\it Anal. PDE} 10(8):1845--1876, 2017]. The model can thus be approximated by means of the minimizing movement (or JKO) scheme, that we solve thanks to the ALG2-JKO scheme proposed in [J.-D. Benamou, G. Carlier, and M. Laborde, {\it ESAIM Proc. Surveys}, 57:1--17, 2016]. The numerical results are compared to a classical upstream mobility Finite Volume scheme, for which strong stability properties can be established.


Multiphase porous media flows as Wasserstein gradient flow
Because of their wide range of interest in the applications, multiphase flows in porous media have been the object of countless scientific studies. In particular, there has been an extensive effort in order to develop reliable and efficient tools for the simulation of such flows. In many practical situations, the characteristic size of the pores (typically of the order the micrometer for regular sandstones) is much smaller than the characteristic size of the domain of interest. The direct numerical simulation of fluid flows at the pore scale is therefore not tractable. The use of homogenised models of Darcy type is therefore commonly used to simulate porous media flows. The derivation of such models is the purpose of a very extended literature. We refer for instance to [3] for an extended introduction to the modelling of porous media flows. But let us stress that, as far as we know, there is no rigorous mathematical derivation of homogenised models for multiphase porous media flows.
Because of the very large friction of the fluid with the porous matrix, the energy is dissipated and inertia is often naturally neglected in the Darcy-type models. The resulting models therefore have a formal gradient flow structure, as highlighted in [10] for immiscible incompressible multiphase porous media flows. This was then rigorously established in [11] that the equations governing such flows can be reinterpreted as a gradient flow in some appropriate Wasserstein space. The goal of this paper is to explore how this new point of view can be used to simulate multiphase flows in porous media.

Incompressible immiscible multiphase flows
As a first step, let us recall the equations governing multiphase porous media flows. We remain synthetic here and refer to the monograph [3] for a rather complete presentation of the models. The porous medium is represented by a convex-bounded open subset of R d (d ≤ 3). Within this porous medium, N + 1 phases are flowing. Denoting by s = (s 0 , . . . , s N ) the saturations, that is, the volume ratios of the various phases in the fluid, the following total saturation relations have to be fulfilled: In what follows, we denote by = s ∈ R N + s 0 + s 1 + · · · + s N = 1 , and by X = s : → R N s(x) ∈ for a.e. x ∈ .
As a consequence of (1.1a), the composition of the fluid is fully characterised by the knowledge of s * = (s 1 , . . . , s N ) ∈ * = (s 1 , . . . , Concerning the evolution, each phase is convected with its own speed: where ω stands for the porosity of the medium and is assumed to be constant in the sequel for simplicity. Then a straightforward rescaling in time allows to choose ω = 1. We further assume a no flux condition across the boundary ∂ for each phase, hence the mass is conserved along time. This motivates the introduction of the set where s 0 = s 0 i : → is a prescribed initial data. The phase speeds v i are prescribed by the Darcy law [14] v i = − κ μ i (∇p i − ρ i g), i ∈ {0, . . . , N}.
( 1 . 1 c) In (1.1c), κ denotes the permeability of the porous medium. For simplicity, it is assumed to be constant and positive. We refer to [11] for the case of space-dependent anisotropic permeability tensors. The fluid viscosity and density are denoted by μ i > 0 and ρ i ≥ 0, respectively, whereas g = −ge z denotes the gravity. The unknown phase pressures p = ( p i ) 0≤i≤N are related to the saturations by N capillary pressure relations: The capillary pressure functions π = (π i ) 1≤i≤N are assumed to derive from a strictly convex and -concave potential : * → R + for some > 0, that is, ∀s * , s * ∈ * , with s * = s * . (1.2) This implies that π : * → R N is strictly monotone (thus one-to-one) and Lipschitz continuous: 0 < π( s * ) − π (s * ) · s * − s * ≤ | s * − s * | 2 , ∀s * , s * ∈ * , with s * = s * , and thus The last inequalities have to be understood in the sense of the symmetric matrices. The function is extended by +∞ outside of * . As established in [11], the problem (1.1) can be interpreted as the Wasserstein gradient flow of the energy: In formula (1.3), the exterior gravitational potential = ( i ) 0≤i≤N is given by 3) one can consider a large class of arbitrary potential ; see [11] for details.
The constraint (1.1a) is incorporated in the energy rather than in the geometry thanks to the term We refer to [5,8,28] for a presentation of the multiphase optimal transportation problem for which the constraint (1.1a) is directly incorporated in the geometry. In order to be more precise in our statements, we need to introduce some extra material concerning the Wasserstein distance to be used to equip A. This is the purpose of the next section.

Remark 1.2
In the previous work [11], the uniform convexity of the capillary potential was required. In (1.2), we relax this assumption into a mere strict convexity requirement. This can be done by slightly adapting the proofs of [11].

Wasserstein distance
For i ∈ {0, . . . , N} we define Given s i , s i ∈ A i , the set of admissible transport plans between s i and s i is given by where M + ( × ) stands for the set of Borel measures on × , and γ (k) i is the kth marginal of the measure γ i . The quadratic Wasserstein distance W i on A i is then defined as Equivalently, the continuity equation (1.1b) allows to give the following dynamical characterisation.
where the minimum runs over curves of measures t → s i,t ∈ A i with endpoints s i,0 , s i,1 and velocity fields t in the sense of distributions. [4], the right variables to be used in the Benamou-Brenier formula (1.6) are not the velocity v, but in fact the momentum m = sv, since the action A(s, m) = |m| 2 s = s|v| 2 becomes then jointly convex in both arguments.

Remark 1.4 As originally developed in
A third equivalent formulation is the Kantorovich dual problem: Proposition 1.5 There holds where the maximum runs over all pairs Any maximiser is called a (pair of optimal) Kantorovich potential.
The viscosity μ i and permeability κ appear in (1.5)-(1.7) as scaling factors in the cost function μ i |x − y| 2 /κ, and this is required for consistency with Darcy's law (1.1c). For more general heterogeneous permeability tensors K(x), one could use instead the intrinsic distance d 2 i (x, y) induced on by the Riemannian tensor μ i K −1 (x); see [27] for a general approach of Wasserstein distances with variable coefficients and [11] in the particular context of multiphase flows in porous media.
With the phase Wasserstein distances (W i ) 0≤i≤N at hand, we can define the global Wasserstein distance W on A := A 0 × · · · × A N by setting

Approximation by minimisation scheme
As already mentioned, the problem (1.1) is the Wasserstein gradient flow of our singular energy (1.3); see our earlier works [10,11]. Rather than discussing the meaning of gradient flows in the Wasserstein setting, we refer to the monograph [2] for an exposition of gradient flows in abstract metric spaces [2] and to [35,36] for a detailed overview. As is now well understood from the work of Jordan et al. [25], one possible way to formalise this gradient flow structure is to implement the JKO scheme (named after Jordan, Kinderlehrer and Otto [25], also referred to as DeGiorgi's minimising movement, see [15]). Given an initial datum s 0 ∈ A with energy E(s 0 ) < ∞ and a time step τ > 0, the strategy consists in This is a variant in the Wasserstein space of the implicit variational Euler scheme: indeed, in Euclidean spaces x ∈ R d and for smooth functions E : R d → R, the Euler-Lagrange equation corresponding to minimising x → 1 2τ |x − x n | 2 + E(x) is nothing but the finite difference approximation x n+1 −x n τ = −∇E(x n+1 ). We refrain from giving more details at this stage and refer again to [2,35,37].
Due to lower semi-continuity and convexity, it is easy to prove that the minimisation problem (1.8) is well posed; hence the discrete solution s τ is uniquely and unambiguously defined. But we still need to construct approximate phase pressures p τ = ( p 1,τ , p 2,τ ). Their construction makes use of the backward Kantorovich potentials (see [11,Section 3] From classical optimal transport theory [35], v n+1 should be interpreted as the discrete velocity driving the ith phase, which will automatically give ∂ t s i + ∇· (s i v i ) = 0 in the limit τ → 0. Hence (1.9) is a discrete counterpart of Darcy law (1.1c). The capillary relation (1.1d) hold as well at the discrete level thanks to relations (1.10), whereas the total saturation constraint (1.1a) is automatically enforced in (1.8) thanks to E(s n+1 ) < ∞. For the sake of brevity we omit the details and refer again to [11].

Main properties of the approximation
Since our system (1.1) of partial differential equations (PDEs) is highly nonlinear, taking the limit s(t) = lim τ →0 s τ (t) will require sufficient compactness both in time and space. In this section we sketch the main arguments leading to such compactness.
Compactness in time is derived from the classical total square distance estimate below, which is a characteristic feature of any JKO variational discretisation. Testing s = s n as a competitor in (1.8) gives first 1 2τ This implies of course the energy monotonicity E(s n+1 ) ≤ E(s n ), but summing over n, we also get the total square distance estimate in the form By definition of the piecewise-constant interpolation, an easy application of the Cauchy-Schwarz inequality gives then the approximate equicontinuity: uniformly in τ , which yields the desired compactness in time (see [2,Proposition 3.10] or [18,Theorem C.10]). Compactness in space will be obtained exploiting the flow interchange technique from [29]. Roughly speaking, this amounts to estimating the dissipation of the driving functional E along a well-behaved auxiliary gradient flow, driven by an auxiliary functional and starting from the minimiser s n+1 . More explicitly, we define the ε-perturbations = (s 0, , . . . ,s N, ) as solutions to the independent heat equations: The key observation is that, for each i = 0, . . . , N, the above heat equation is a gradient flow in the Wasserstein space (A i , W i ) with driving functional μ i H, where the Boltzmann entropy In addition to the usual regularising effects, this heat equation is particularly well behaved here in the sense that it preserves the total saturation constraint denotes the JKO functional, then by optimality of the minimiser s n+1 in (1.8) we must have lim sup The energy term E(s ) = (s * ) can easily be differentiated under the integral sign (with respect to ), while the variation of the first W 2 (s , s n ) term can be estimated using the evolution variational inequality [2] for the well-behaved H-flows (this metric characterisation precisely requires some displacement convexity of the auxiliary flow, see [19,Theorem 2.23]). Omitting again the details, one gets in the end the dissipation estimate see [11,Section 2.2] for the details. Exploiting the previous total square distance estimate and summing over n = 0, . . . , T/τ (or equivalently integrating in time), we control next: for arbitrary T > 0 and fixed initial datum s 0 . It is worth recalling at this stage that, due to our assumption (1.2), π (s * ) = ∇ s * (s * ) is a strictly monotone thus invertible map of s* due to the strict convexity of . The compactness with respect to the space variable of (s τ ) τ >0 then follows from (1.13).

Remark 1.7 A formal but more PDE-oriented explanation of the above flow interchange simply
consists in taking log(s i ) as a test function in the weak formulation of system (1.1). The delicate technical part is to justify this computation and mimic this formal chain rule in the discrete time setting in order to retrieve enhanced regularity of the JKO minimisers.
Exploiting the above compactness, one can argue as in [11] and finally prove the following convergence results. The existence of a weak solution to the problem (1.1) is a direct byproduct.
and the limit ((s, p)) is a weak solution of (1.1).

Numerical approximation of the flow
We present here the ALG2-JKO scheme and the upstream mobility finite volume (FV) scheme. The first method is based on the variational JKO scheme (1.8) described in Section 1.3, whereas the second method is based on the PDE formulation of the problem (1.1) given by (1.1a)-(1.1d). Both methods are well adapted for gradient flows equations, and more precisely we will check the following key properties for the numerical solutions: • preservation of the positivity; • conservation of the mass and saturation constraints; and • energy dissipation along solutions.

The ALG2-JKO scheme
This algorithm relies on the seminal work of Benamou and Brenier [4], where an augmented Lagrangian approach was used to compute Wasserstein distances. In [6], this approach was extended to the computation of Wasserstein gradient flows. The method is very well suited for computing solutions to constrained gradient flows, as it will appear in the numerical simulations presented in Section 3.

The augmented Lagrangian formulation
Roughly speaking, the ALG2-JKO scheme consists in rewriting the single JKO step (1.8) as a more fashionable (and effectively implementable) convex minimisation problem. In order to do so, let us first introduce the convex lower-semicontinuous 1-homogeneous action function given, for all (s, m) ∈ R × R d , by if s > 0, 0 ifs = 0 and m = 0, +∞ otherwise.
We recall that m = sv is the momentum variable in the continuity equation ∂ t s + ∇· (sv) = 0 and |m| 2 /s = s|v| 2 is a kinetic energy, see Remark 1.4. As originally observed in [4], the function A can be seen as the support function Taking advantage of the Benamou-Brenier formula (1.6), and given the previous JKO step s n , (1.8) can be recast as where the infimum runs over curves of measures t → s t = (s 0,t , . . . , s N,t ) ∈ A and momenta t → Note that only the initial endpoint s t=0 = s n is prescribed for the curve (s t ) t∈ [0,1] . The terminal endpoint is free and contributes to the objective functional (2.4) through the E(s t=1 ) term, and the JKO minimiser will be retrieved as s n+1 = s t=1 . Note also that the minimising curve (s t ) t∈ [0,1] in (2.4) and (2.5) will automatically be a Wasserstein geodesic between the successive JKO minimisers s t=0 = s n and s t=1 = s n+1 . As a first step towards a Lagrangian formulation, we rewrite the constraint (2.5) as a sup problem with multipliers φ i (t, x) : the problems (2.4) and (2.5) finally become after a few elementary manipulations Here E * τ denotes the Legendre transform of E τ := τ E. This dual problem can be reformulated as and χ K 2μ i /κ stands again for the characteristic function of K 2μ i /κ . Introducing a Lagrange multiplier σ = (s, m,s 1 ) for the constraint φ = q, finding a minimiser s n+1 in the JKO scheme (1.8) is thus equivalent to finding a saddle-point of the Lagrangian Here we slightly abuse the notations: 1] are time-depending curves, whiles 1 ∈ A is independent of time. The scalar product in (2.7) is We stress that the free variables 1 is a priori independent of the curve (s t ) t∈[0, 1] , but that the saddle-point will ultimately satisfy s t=1 =s 1 . In the Lagrangian (2.7), the original unknowns (s, m,s 1 ) become the Lagrange multipliers for the constraint q = φ,, that is, respectively: For some fixed regularisation parameter r > 0, we introduce now the augmented Lagrangian where the extra regularising term is given by the L 2 norm: Observe that being a saddle-point of (2.7) is equivalent to being a saddle-point of (2.8); see for instance [22]. Thus in order to solve one step of the JKO scheme (1.8), it suffices to find a saddle-point of the augmented Lagrangian L r .

Algorithm and discretisation
The augmented Lagrangian algorithm ALG2 aims at finding a saddle-point of L r and consists in a splitting scheme. Starting from (φ 0 , q 0 , σ 0 ), we generate a sequence (φ k , q k , σ k ) k≥0 by induction as follows: Step 1: minimise with respect to φ: Step 2: minimise with respect to q: Step 3: maximise with respect to σ , which amounts here to updating the multiplier by the gradient ascent formula: Since step 3 is a mere point-wise update we only describe in details the first two steps. In order to keep the notations light we sometimes write s i (t, x) = s i,t (x), and likewise for any other variable depending on time.
• The first step corresponds to solving N + 1 independent linear elliptic problems in time and space, namely • The second step splits into two convex point-wise subproblems. The first one corresponds to projections onto the parabolas K 2μ i /κ : This projection P K 2μ i /κ onto K 2μ i /κ is explicitly given by (see [33]) where λ is the largest real root of the cubic equation The second subproblem should update c. To this end, we need to solve the point-wise proximal problem: for each x ∈ , where E * τ (x, ·) is the Legendre transform of the energy density E τ (x, ·) = τ E(x, ·) in its second argument (E being implicitly defined as E(s) = E(x, s(x)) dx).
Notice that the energy functional E only plays a role in the minimisation with respect to the internal c variable, namely the second subproblem (2.9) in step 2. In Section 3, we will try to make this step explicit for our two particular applications.
In order to implement this algorithm in a computational setting we use P2 finite elements in time and space for φ, and P1 finite elements for σ and q. The variables ∇ t, ) are understood as the projection onto P1 finite elements and the algorithm was implemented using FreeFem++ [24]. The convergence of this algorithm is known in finite dimension [22], that is, the iterates (φ k , q k , σ k ) are guaranteed to converge to a saddle-point (φ, q, σ ) as k → ∞. Once the saddle-point is reached, the output σ = (s, m,s 1 ) is a minimiser for the problems (2.4) and (2.5) and the solution of the JKO scheme (1.8) is simply recovered as s n+1 =s 1 = s| t=1 .
Numerically, the Benamou-Brenier formula involves an additional time dimension to be effectively discretised in each elementary JKO step, and this can be seen as a drawback. However, the successive JKO densities are close due to the small time step τ → 0 (indeed W (s n+1 , s n ) = O( √ τ ) from the total square distance estimate (1.11)) and, in practice, only a very few inner time steps are needed.

Some properties of the approximate solution
As previously mentioned, the above Lagrangian framework can be practically implemented by simply projecting the (infinite dimensional) problem onto P1/P2 finite elements. Provided that the iteration procedure (steps 1-3 in Section 2.1.2) converges as k → ∞, as guaranteed from [22], the saddle-point σ = (s, m,s 1 ) satisfies by construction: As a consequence of (i), the scheme preserves the positivity, that is, s n+1 i ≥ 0, whereas (ii) ensures the mass conservation s n+1 i = s n i . Moreover, the fully discrete ALG2-JKO scheme preserves by construction the gradient flow structure; hence the scheme is automatically energy diminishing. Since the energy functional (1.3) includes the χ term accounting for the saturation constraint s i = 1, one can and should include this convex indicator term in the discretised energy. This constraint is then passed on to the proximal operator to be used in the implementation; see Section 3 for details. As a result, the saturation constraint is satisfied.

Upstream mobility finite volume scheme
The ALG2-JKO scheme described in the previous section will be compared to the widely used upstream mobility finite volume scheme [7,21,34]. As a first step, let us detail how is discretised.

The finite volume mesh
The domain is assumed to be polygonal. Then following [20], an admissible mesh consists in a triplet T, E, (x K ) K∈T . The elements K of T are open polygonal convex subsets of called control volumes. Their boundaries are made of elements σ ∈ E of codimension 1 (edges if d = 2 or faces if d = 3). Let K, L be two distinct elements of T, then K ∩ L is either empty, or reduced to a point (a vertex), or there exists σ ∈ E denoted by σ = K|L such that K ∩ L = σ . In particular, two control volumes share at most one edge. We denote by E K = σ ∈ E | σ ∈E k σ = ∂K the set of the edges associated to an element K ∈ T, and by N K = {L ∈ T there exists σ = K|L ∈ E}, the set of the neighboring control volumes to K. We also denote by The last element (x K ) K∈T of the triplet corresponds to the so-called cell centers. To each control volume K ∈ T, we associate an element x K ∈ such that for all L ∈ N K , the straight line (x K , x L ) is orthogonal to the edge K|L. This implicitly requires that x K and x L are distinct, and we denote by d σ = |x K − x L | for σ = K|L the distance between the cell centers of the neighboring control volumes K and L. For σ ∈ E K ∩ E ext , we denote by x σ the projection of x K on the hyperplane containing σ, and by d σ = |x K − x σ |. We also require that the vector x L − x K is oriented in the same sense as the normal n K,σ to σ ∈ E K outward with respect to K. We refer to Figure 1 for an illustration of the notations used hereafter. Beyond cartesian grids, there are two classical ways to construct admissible meshes in the above sense when d = 2. The first one consists in the classical Delaunay triangulation, the cell center x K of K ∈ T being the center of the circumcircle of K. The second classical construction consists in choosing the cell centers (x K ) at first and then to construct T as the associated Voronoï diagram.
In what follows, we denote by m K the d-dimensional Lebesgue measure of the control volume K ∈ T, while m σ denotes the (d−1)-dimensional Lebesgue measure of the edge σ ∈ E. We also denote by a σ = m σ d σ the transmissivity of the edge σ . In order to simplify the presentation, we restrict our presentation to the case of uniform time discretisations with time step τ > 0. The extension to the case of time discretisations with varying time steps does lead to any particular difficulty.

Definition of the finite volume scheme
The finite volume scheme relies on the discretisation of the Euler-Lagrange equations (1.1) rather than on the minimising movement scheme (1.8). The main unknowns to the problems are located at the cell centers (x K ) K∈T . They consist in discrete saturations s n i,K s i (x K , nτ ) and discrete pressures p n i,K p i (x K , nτ ). In what follows, we denote by s n K = s n i,K 0≤i≤N (resp. p n K = p n i,K 0≤i≤N ) and s n T = s n K K∈T .
The first equation of the scheme is a straightforward consequence of (1.1a), that is, This motivates the introduction of the discrete counterpart X T of X defined by X T = {s T |s K ∈ for all K ∈ T}; so that (2.10a) amounts to requiring that s n T belongs to X T for all n (the non-negativity of the saturations will be established later on). The capillary pressure relations (1.1d) are discretised into Integrating (1.1b) over the control volume K ∈ T (recall here that the porosity ω was artificially set to 1) and using Stokes' formula, one gets the natural approximation: where v i is related to p i through to Darcy law (1.1c). Thanks to the orthogonality condition on the mesh, the choice is consistent -we use the shortened notation i,K = i (x K ). In accordance with the no-flux boundary conditions, we impose that v n i,K,σ = 0, ∀σ ∈ E K ∩ E ext , ∀n ≥ 1.
It remains to define the approximate saturations s n i,σ for σ ∈ E int . We use here the very classical upwind choice [7,21,34], that is, Note that even though the mapping (s n T , p n T ) → s n E = ((s n i,σ ) 0≤i≤n ) σ ∈E int is discontinuous, the quantity s n i,σ v n i,K,σ depends on a continuous way of the main unknowns. The scheme (2.10) amounts to a nonlinear system of equations to be solved a each time step. This will be practically done thanks to Newton-Raphson method. But before, we establish some properties of the FV scheme, namely the energy decay, the entropy control, the non-negativity of the saturations, or the existence of a solution (s n T , p n T ) to the scheme.

Some properties of the approximate solution
The first key property of the FV scheme that we point out is the non-negativity of the saturations: . . , N}, ∀K ∈ T, ∀n ≥ 1.
In order to establish this estimate, it suffices to rewrite (2.10c) as thanks to (2.10e). In the previous expression, we used the convention a − = max(0, −a) ≥ 0. Assume for contradiction that s n i,K is negative, then so does the left-hand side, while the right-hand side is non-negative by induction. Together with (2.10a), this shows that The scheme is mass conservative for the N + 1 phases since Together with the no-flux boundary conditions, this shows that the mass is conserved along time: Then the discrete solution s n T remains in the discrete counterpart A T of A defined as the elements The second term in the above expression is clearly non-negative. Concerning the first term, one can use the constraint (2.10a) to rewrite as the last inequality being a consequence of the convexity of . This establishes that the scheme is energy diminishing: denoting by one has The last a priori estimate we want to point out is the discrete counterpart of the flow interchange estimate. It is obtained by multiplying (2.10c) by τ μ i log(s n i,K ) and by summing over i ∈ {0, . . . , N} and K ∈ T, leading to As already discussed in Remark 1.7 this corresponds to taking log s i as a test function in the weak formulation of the continuous PDEs. The first term of (2.14) can be estimated thanks to an elementary convexity inequality: Note that the entropy functional H is bounded on X T . The second term of (2.14) can be estimated as follows. First, the concavity of s → log(s) yields so that the upwind choice (2.10e) for s n i,σ ensures that Using the expression (2.10d) of v n i,K,σ and the relation (2.10a) on the saturations, one gets that Recalling the definition (1.4) of the external potential and denoting by i,σ = i (x σ ), one has Since 0 ≤ s n i,K ≤ 1, this implies that On the other hand, the assumption (1.2) on the capillary pressure potential ensures that Hence collecting the previous inequalities in (2.14) provides the following discrete L 2 loc (H 1 )estimate one the capillary pressures Clearly, (2.15) is the discrete counterpart of the estimate (1.13) obtained thanks to the flow interchange technique. The derivation of a discrete L 2 loc (H 1 ) estimate on the phase pressures from (2.15) and (2.13) requires one additional assumption on the capillary pressure functions (π i ) 1≤i≤n . More precisely, we assume that π i only depends on s i : Since is convex, the functions π i are increasing. Assumption (2.16) is needed to establish that, at least for fine enough grids, there holds for some uniform α. Thanks to this estimate, one can follow the lines of [11, Proposition 3.4 and Corollary 3.5] (see also [13]) to derive the estimate The phase pressures being defined up to an additive constant (recall that they are related to Kantorovich potentials), one has to fix this degree of freedom. This can be done by enforcing Based on the a priori estimates (2.11) and (2.17), we can make use of a topological degree argument (see for instance [16]) to claim that there exists (at least) one solution to the scheme. Moreover, assuming some classical regularity on the mesh T (see for instance [1]), one can prove the piecewise-constant approximate solutions converge towards a weak solution when the size of the mesh T and the time step τ tend to 0. This convergence results together with the properties (2.11)-(2.17) as well as the wide popularity of this scheme in the engineering community makes this scheme a reference for solving (1.1). In the next section, we show that the ALG2-JKO scheme presented in Section 2.1 produces very similar results: same qualitative results, conservation of the mass of each phase and preservation of the positivity.

Numerical experiments
In this section, we compare the numerical results produced by the ALG2-JKO scheme presented in Section 2.1 with the upstream mobility finite volume scheme of Section 2.2. In the sequel the regularisation parameter r introduced in the augmented Lagrangian formulation (2.8) is fixed to r = 1 for simplicity, which gives satisfactory numerical results. The case of a three-phase flow (typically water, oil and gas) is presented in Section 3.2, whereas a two-phase flow is simulated in Section 3.1. In both cases, we do not have analytical solutions at hand and the results are compared thanks to snapshots. Note the both time discretisations are of order 1. The extension to order two methods is a challenging task. Concerning the ALG2-JKO scheme, one possibility could be to use the order 2 approximation based on the midpoint rule proposed in [26], but there is no rigorous foundation to this work up to now as far as we know. An alternative approach would be to use the variational BDF2 approach proposed in [30]. But the variational problem to be solved at each time step is no longer convex-concave, so that its practical resolution becomes more involving. Concerning the finite volume scheme, there is (up to our knowledge) no time integrator of order 2 that ensures the decay of a general energy. Going to higher order time discretisations yield also difficulties concerning the preservation of the positivity. This explains why the backward Euler scheme is very popular in the context of the simulation of multiphase porous media flows.

Two-phase flow with Brooks-Corey capillarity
As a first example we consider a two-phase flow, where water (s 0 ) and oil (s 1 ) are competing within the background porous medium. For the capillary pressure, we choose the very classical Brooks-Corey (or Leverett) model: The capillary pressure is depicted on Figure 2. We refer to [3] for an overview of the classical capillary pressure relation for two-phase flows. As in Section 1.1, the corresponding energy reads explicitly E(s 0 , s 1 ) = 0 s 0 + 1 s 1 − 2α (1 − s 1 ) 1/2 + χ (s 0 , s 1 ).
As already mentioned, only the second subproblem (2.9) in step 2 of the ALG2-JKO algorithm depends on the choice of the energy functional. For the above particular case, this reads: for each x ∈ and setting c : This minimisation problem is equivalent to computing where the proximal operator Prox f of a given convex, lower semicontinuous function f : Thanks to Moreau's identity it suffices to compute Prox E τ in order to determine Prox E * τ , and we never actually compute the Legendre transform E * τ (x, ·). Computing the proximal operator c k+1 (x) = Prox E * τ (x,·) (c) thus amounts to evaluating Finally, (c 0 ,c 1 ) := Prox E τ (x,·) (c 0 , c 1 ) is computed by solving and then settingc 0 = 1 −c 1 . More explicitly,c 1 is the positive part of the root on (−∞, 1) of To conclude, we set (c n+1 0 (x), c n+1 1 (x)) = (c 0 −c 0 , c 1 −c 1 ). We start from the initial profile depicted in Figure 2. In Figure 3, we compare the numerical solutions of problem (1.1) with Brooks-Corey capillarity (3.1) obtained thanks to the ALG2-JKO scheme and to the upstream mobility finite volume scheme. Simulations with the ALG2-JKO scheme are carried using a structured grid with 5000 triangles and 2601 vertices in space and a single inner time step, and with 200 JKO steps (τ = 0.05). Simulations with the upstream mobility finite volume scheme are performed on the corresponding Cartesian grid with 2500 squares. The time step τ appearing in 28 can be also set to 0.05 here since Newton's method converges rather easily in this test case.  As expected, the results produced by the two schemes are very similar. The dense phase (the water) is instantaneously diffused in the whole domain because of the singularity of π 1 near 1. When time goes, oil slowly moves to the top because of buoyancy.

Three-phase flow with quadratic capillary potential
In the second test case, we consider the case of a three-phase flow where water (s 0 ), oil (s 1 ) and gas (s 2 ) are in competition within the porous medium. Here we assume that the capillary pressure functions π 1 and π 2 are linear: The corresponding capillary potential is then given by (s * ) = α 1 2 (s 2 1 ) + α 2 2 (s 2 2 ).
Assumptions (1.2) and (2.16) are fulfilled, so that we are in the theoretical framework of our statements, that is, convergence of the minimising movement scheme and of the finite volume scheme. However, the problem is difficult to simulate because of the rather large ratios on the viscosities. Indeed, the phase 0 represents water, the phase 1 corresponds to oil and the phase 2 corresponds to gas, and we set μ 0 = 1, μ 1 = 50, μ 2 = 0.1, and ρ 0 = 1, ρ 1 = 0.87, ρ 2 = 0.1.
The resulting energy in the JKO scheme (1.8) is given by and we denote accordingly, for x ∈ and c = (c 0 , c 1 , Setting again c = −φ k+1 (1, x) +s k 1 (x) and taking advantage of Moreau's identity (3.2), the second subproblem (2.9) of step 2 is equivalent to, for all x ∈ R d ,  Evaluating the proximal operatorc := Prox E τ (x,·) (c) is equivalent to solving (c 1 ,c 2 ) = argmin The solution (u 1 , u 2 ) of the unconstrained version of (3.3) is explicitly given by  Figures 5-7 show the evolution of the three phases with quadratic capillarity potential. Again, the simulation with the ALG2-JKO scheme is carried out using a 50 × 50 discretisation in space, with a single inner time step. There are 200 JKO steps (τ = 0.05). The convergence of the augmented Lagrangian iterative method is rather slow: it took around 10 h on a laptop to produce the results with FreeFem++. But because of the large viscosity ratio, Newton's method had severe difficulties to converge for the upstream mobility scheme. A very small time step (τ = 10 −4 ) was needed, so that approximately 1 days of computation on a server were needed to produce the results with Matlab. Again, simulations with the upstream mobility finite volume scheme are performed on a Cartesian grid with 2500 squares. Once again, both methods produce similar results, as highlighted in Figures 5-7.
Due to the large viscosity ratios, two distinct time scale appear in the numerical results. Since water and gas have smaller mobilities, they move much faster than oil. This quick phenomenon is not well captured by the ALG2-JKO scheme. The interface between oil and gas is already almost horizontal at t = 0.1. This horizontal interface is captured by the finite volume scheme but not by the ALG2-JKO scheme that encounters difficulties to converge for the early time steps. The finite volume scheme also has difficulties to converge, enforcing us to consider very small time steps. Oil is much less mobile and its interface with the two other phases remains almost vertical at that time. Then oil evolves slowly towards its equilibrium state that consists in a horizontal layer trapped between gas above and water below. This longtime equilibrium is not yet reached for t = 10.

Energy dissipation
As already highlighted, both schemes dissipate the energy along time. The goal of this test case is to compare the energy dissipation. To this end, we consider a test case proposed in [9]. We consider a two-phase flow with oil (i = 1) and water (i = 0) with ρ 1 = 0.87, ρ 0 = 1, μ 1 = 10, and μ 0 = 1, while κ = 1 and ω = 1. The capillary pressure law is given by p 1 − p 0 = π 1 (s 1 ) = s 1 2 ; so that the energy is defined by We consider the initial data s 0 1 (x) = e −4|x| 2 . At equilibrium, the saturation s ∞ 1 minimises E under the constraints s ∞ 1 ∈ [0, 1] and It is, therefore, given by either s ∞ 1 ∈ {0, 1} or π 1 (s 1 ) = (ρ 1 − ρ 0 )g · x + γ , (3.5) the constant γ being fixed thanks to (3.4). Similar calculations can be performed in the discrete settings, both for the ALG2-JKO scheme and the finite volume scheme. Then one computes for both scheme the relative energy E(s 1 ) − E(s ∞ 1 ) ≥ 0 that we plot as a function of time in Figure 8. The convergence towards the equilibrium appears to be exponential in both cases.

Conclusion
We proposed to apply the ALG2-JKO scheme of [6] to simulate multiphase porous media flows. The results have been compared to the widely used upstream mobility finite volume scheme. The ALG2-JKO scheme appears to be robust with respect to the capillary pressure function and overall with respect to the viscosity ratios. The method is parameter free (the only parameter r has a rather low influence and is chosen equal to 1 in the computations) and is unconditionally converging whatever the time step. This is a great advantage when compared to the Newton method that may require very small time steps in the presence of large viscosity ratios. Moreover, the ALG2-JKO scheme preserves the positivity of the saturations, the constraint on the sum of the saturations, and it is locally conservative. Its main drawback concerns the restriction to linear mobility function so that formulae (2.2) and (2.3) hold (this can probably be extended to the non-physical case of concave mobilities [17] but we did not push into this direction). Finally, let us stress that the code depends only at stage (2.9) of the energy. Therefore, the extension of the ALG2-JKO approach to multiphase models with different energies (like, for instance, degenerate Cahn-Hilliard models [12,32]) is not demanding once the code is written. A natural extension to this work would be to add source terms corresponding, for instance, to production wells. This would, for instance, require to adapt the material of [23] to our context.

Conflicts of interest
None.